Astrophysical Magnetohydrodynamics Research Articles

Magnetically induced gluing bifurcations of three-tori in spherical Couette flows

The global dynamics of three-tori associated with homoclinic/heteroclinic global (infinite period) bifurcations are investigated for the magnetized spherical Couette problem, a paradigmatic model in geo- and astrophysical magnetohydrodynamics (MHD). A novel homoclinic bifurcation, involving collision between three-tori, is described. In addition, a heteroclinic bifurcation connecting two unstable two-tori with a stable three-torus is also analyzed. The role of the flow's spatial symmetries in this bifurcation scenario is also investigated. This bifurcation scenario gives rise to MHD flows that combine small with extremely large time scales.

Parthenon—a performance portable block-structured adaptive mesh refinement framework

On the path to exascale the landscape of computer device architectures and corresponding programming models has become much more diverse. While various low-level performance portable programming models are available, support at the application level lacks behind. To address this issue, we present the performance portable block-structured adaptive mesh refinement (AMR) framework Parthenon, derived from the well-tested and widely used Athena++ astrophysical magnetohydrodynamics code, but generalized to serve as the foundation for a variety of downstream multi-physics codes. Parthenon adopts the Kokkos programming model, and provides various levels of abstractions from multidimensional variables, to packages defining and separating components, to launching of parallel compute kernels. Parthenon allocates all data in device memory to reduce data movement, supports the logical packing of variables and mesh blocks to reduce kernel launch overhead, and employs one-sided, asynchronous MPI calls to reduce communication overhead in multi-node simulations. Using a hydrodynamics miniapp, we demonstrate weak and strong scaling on various architectures including AMD and NVIDIA GPUs, Intel and AMD x86 CPUs, IBM Power9 CPUs, as well as Fujitsu A64FX CPUs. At the largest scale on Frontier (the first TOP500 exascale machine), the miniapp reaches a total of 1.7 × 1013 zone-cycles/s on 9216 nodes (73,728 logical GPUs) at [Formula: see text] weak scaling parallel efficiency (starting from a single node). In combination with being an open, collaborative project, this makes Parthenon an ideal framework to target exascale simulations in which the downstream developers can focus on their specific application rather than on the complexity of handling massively-parallel, device-accelerated AMR.

GR-Athena++: Puncture Evolutions on Vertex-centered Oct-tree Adaptive Mesh Refinement

Abstract Numerical relativity is central to the investigation of astrophysical sources in the dynamical and strong-field gravity regime, such as binary black hole and neutron star coalescences. Current challenges set by gravitational-wave and multimessenger astronomy call for highly performant and scalable codes on modern massively parallel architectures. We present GR-Athena++, a general-relativistic, high-order, vertex-centered solver that extends the oct-tree, adaptive mesh refinement capabilities of the astrophysical (radiation) magnetohydrodynamics code Athena++. To simulate dynamical spacetimes, GR-Athena++ uses the Z4c evolution scheme of numerical relativity coupled to the moving puncture gauge. We demonstrate stable and accurate binary black hole merger evolutions via extensive convergence testing, cross-code validation, and verification against state-of-the-art effective-one-body waveforms. GR-Athena++ leverages the task-based parallelism paradigm of Athena++ to achieve excellent scalability. We measure strong-scaling efficiencies above 95% for up to ∼1.2 × 104 CPUs and excellent weak scaling is shown up to ∼105 CPUs in a production binary black hole setup with adaptive mesh refinement. GR-Athena++ thus allows for the robust simulation of compact binary coalescences and offers a viable path toward numerical relativity at exascale.

Towards performance portability in the Spark astrophysical magnetohydrodynamics solver in the Flash-X simulation framework

Simulations of core-collapse supernovae, and other astrophysical phenomena, are quintessential extreme-scale computing challenges. For core-collapse supernova simulations to be carried out by the ExaStar project under the Exascale Computing Project umbrella, a robust, efficient, and state-of-the-art magnetohydrodynamics solver is a critical requirement. In Flash-X, the primary software instrument for ExaStar, a new magnetohydrodynamics solver has been designed and implemented from the ground up to achieve accuracy and efficiency for simulations of complex astrophysical flows. This new solver, dubbed Spark, uses high-order spatial reconstruction, Runge–Kutta time integration, and an efficient cell-centered approach to satisfying the divergence-free condition for the magnetic fields. Spark was written to be optimized for data locality in cache hierarchy of CPUs. Since data locality optimizations for cache hierarchy are not directly compatible with those of accelerators, we have taken the approach of using program synthesis to avoid massive amounts of code replication that would be necessary if we were to maintain two different versions of the solver. Our program synthesis relies on a simple key-dictionary approach, implemented in python, that enables us to assemble the version of the solver suitable for the target hardware from code fragments identified by specific keys. In this paper, we describe the data locality optimizations of the solver for CPUs and accelerators and the program synthesis tools that enable this portability. We also detail the parallel performance of Spark for both CPUs and accelerators.

The Pencil Code, a modular MPI code for partial differential equations and particles: multipurpose and multiuser-maintained

The Pencil Code is a highly modular physics-oriented simulation code that can be adapted to a wide range of applications. It is primarily designed to solve partial differential equations (PDEs) of compressible hydrodynamics and has lots of add-ons ranging from astrophysical magnetohydrodynamics (MHD) to meteorological cloud microphysics and engineering applications in combustion. Nevertheless, the framework is general and can also be applied to situations not related to hydrodynamics or even PDEs, for example when just the message passing interface or input/output strategies of the code are to be used. The code can also evolve Lagrangian (inertial and noninertial) particles, their coagulation and condensation, as well as their interaction with the fluid.

An Extension of the Athena++ Framework for Fully Conservative Self-gravitating Hydrodynamics

Numerical simulations of self-gravitating flows evolve a momentum equation and an energy equation that account for accelerations and gravitational energy releases due to a time-dependent gravitational potential. In this work, we implement a fully conservative numerical algorithm for self-gravitating flows, using source terms, in the astrophysical magnetohydrodynamics framework Athena++. We demonstrate that properly evaluated source terms are conservative when they are equivalent to the divergence of a corresponding "gravity flux" (i.e., a gravitational stress tensor or a gravitational energy flux). We provide test problems that demonstrate several advantages of the source-term-based algorithm, including second order convergence and round-off error total momentum and total energy conservation. The fully conservative scheme suppresses anomalous accelerations that arise when applying a common numerical discretization of the gravitational stress tensor that does not guarantee curl-free gravity.

Fourth-order accurate finite-volume CWENO scheme for astrophysical MHD problems

In this work, a simple fourth-order accurate finite volume semi-discrete scheme is introduced to solve astrophysical magnetohydrodynamics (MHD) problems on Cartesian meshes. Hydrodynamic quantities like density, momentum and energy are discretised as volume averages. The magnetic field and electric field components are discretised as area and line averages respectively, so as to employ the constrained transport technique, which preserves the solenoidality of the magnetic field to machine precision. The present method makes use of a dimension-by-dimension approach employing a 1-D fourth-order accurate centrally weighted essentially non-oscillatory (1D-CWENO4) reconstruction polynomial. A fourth-order accurate, strong stability preserving (SSP) Runge-Kutta method is used to evolve the semi-discrete MHD equations in time. Higher-order accuracy of the scheme is confirmed in various linear and nonlinear multi-dimensional tests and the robustness of the method in avoiding unphysical numerical artifacts in the solution is demonstrated through several complex MHD problems.

Kinematic dynamo action of a precession-driven flow based on the results of water experiments and hydrodynamic simulations.

ABSTRACTThe project DRESDYN (DREsden Sodium facility for DYNamo and thermohydraulic studies) conducted at Helmholtz–Zentrum Dresden–Rossendorf (HZDR) provides a new platform for a variety of liquid sodium experiments devoted to problems of geo- and astrophysical magnetohydrodynamics. The most ambitious experiment within this project is a precession-driven dynamo experiment that currently is under construction. It consists of a cylinder filled with liquid sodium that simultaneously rotates around two axes. The experiment is motivated by the idea of a precession-driven flow as a complementary energy source for the geodynamo or the ancient lunar dynamo. In the present study, we address numerical and experimental examinations in order to identify parameter regions where the onset of magnetic field excitation will be most probable. Both approaches show that in the strongly nonlinear regime the flow is essentially composed of the directly forced primary Kelvin mode and higher modes in terms of standing inertial waves that arise from nonlinear self-interactions. A peculiarity is the resonance-like emergence of an axisymmetric mode that represents a double roll structure in the meridional plane, which, however, only occurs in a very limited range of the precession ratio. This axisymmetric mode turns out to be beneficial for dynamo action, and kinematic simulations of the magnetic field evolution induced by the time-averaged flow exhibit magnetic field excitation at critical magnetic Reynolds numbers around , which is well within the range of the planned liquid sodium experiment.

Topological data analysis and diagnostics of compressible magnetohydrodynamic turbulence

The predictions of mean-field electrodynamics can now be probed using direct numerical simulations of random flows and magnetic fields. When modelling astrophysical magnetohydrodynamics, it is important to verify that such simulations are in agreement with observations. One of the main challenges in this area is to identify robust quantitative measures to compare structures found in simulations with those inferred from astrophysical observations. A similar challenge is to compare quantitatively results from different simulations. Topological data analysis offers a range of techniques, including the Betti numbers and persistence diagrams, that can be used to facilitate such a comparison. After describing these tools, we first apply them to synthetic random fields and demonstrate that, when the data are standardized in a straightforward manner, some topological measures are insensitive to either large-scale trends or the resolution of the data. Focusing upon one particular astrophysical example, we apply topological data analysis to H iobservations of the turbulent interstellar medium (ISM) in the Milky Way and to recent magnetohydrodynamic simulations of the random, strongly compressible ISM. We stress that these topological techniques are generic and could be applied to any complex, multi-dimensional random field.

WOMBAT: A Scalable and High-performance Astrophysical Magnetohydrodynamics Code

Abstract We present a new code for astrophysical magnetohydrodynamics specifically designed and optimized for high performance and scaling on modern and future supercomputers. We describe a novel hybrid OpenMP/MPI programming model that emerged from a collaboration between Cray, Inc. and the University of Minnesota. This design utilizes MPI-RMA optimized for thread scaling, which allows the code to run extremely efficiently at very high thread counts ideal for the latest generation of multi-core and many-core architectures. Such performance characteristics are needed in the era of “exascale” computing. We describe and demonstrate our high-performance design in detail with the intent that it may be used as a model for other, future astrophysical codes intended for applications demanding exceptional performance.

Pegasus: A new hybrid-kinetic particle-in-cell code for astrophysical plasma dynamics

We describe Pegasus, a new hybrid-kinetic particle-in-cell code tailored for the study of astrophysical plasma dynamics. The code incorporates an energy-conserving particle integrator into a stable, second-order–accurate, three-stage predictor–predictor–corrector integration algorithm. The constrained transport method is used to enforce the divergence-free constraint on the magnetic field. A δf scheme is included to facilitate a reduced-noise study of systems in which only small departures from an initial distribution function are anticipated. The effects of rotation and shear are implemented through the shearing-sheet formalism with orbital advection. These algorithms are embedded within an architecture similar to that used in the popular astrophysical magnetohydrodynamics code Athena, one that is modular, well-documented, easy to use, and efficiently parallelized for use on thousands of processors. We present a series of tests in one, two, and three spatial dimensions that demonstrate the fidelity and versatility of the code.

Divergence-free interpolation of vector fields from point values — exact ∇ ⋅B = 0 in numerical simulations

Abstract In astrophysical magnetohydrodynamics (MHD) and electrodynamics simulations, numerically enforcing the ∇·B= 0 constraint on the magnetic field has been difficult. We observe that for point-based discretization, as used in finite-difference type and pseudo-spectral methods, the ∇·B= 0 constraint can be satisfied entirely by a choice of interpolation used to define the derivatives of B. As an example we demonstrate a new class of finite-difference-type derivative operators on a regular grid which has the ∇·B= 0 property. This principle clarifies the nature of ∇·B≠ 0 errors. The principles and techniques demonstrated in this Letter are particularly useful for the magnetic field, but can be applied to any vector field. This Letter serves as a brief introduction to the method and demonstrates an implementation showing convergence.

Laboratory astrophysics and non-ideal equations of state: the next challenges for astrophysical MHD simulations

Laboratory astrophysics holds great promise not only as a highly effective validation tool for astrophysical magneto-hydrodynamics (MHD) codes but it also presents a unique challenge for these codes. The high-density plasmas found in these experiments are not well modeled by the ideal equations of state (EOS) found in most astrophysical simulation codes. To solve this problem, we replaced the ideal EOS scheme in an existing MHD code, AstroBEAR, with a non-ideal EOS method and validated our implementation with van der Waals shock tube tests. The improved code is also able to model flows that contain more than one material, as required in laboratory experiments. Simulations of jet experiments performed at the OMEGA Laser reproduce the morphology of the jet much better than when the code used a single material and an ideal EOS.

Athena: A New Code for Astrophysical MHD

A new code for astrophysical magnetohydrodynamics (MHD) is described. The code has been designed to be easily extensible for use with static and adaptive mesh refinement. It combines higher order Godunov methods with the constrained transport (CT) technique to enforce the divergence-free constraint on the magnetic field. Discretization is based on cell-centered volume averages for mass, momentum, and energy, and face-centered area averages for the magnetic field. Novel features of the algorithm include (1) a consistent framework for computing the time- and edge-averaged electric fields used by CT to evolve the magnetic field from the time- and area-averaged Godunov fluxes, (2) the extension to MHD of spatial reconstruction schemes that involve a dimensionally split time advance, and (3) the extension to MHD of two different dimensionally unsplit integration methods. Implementation of the algorithm in both C and FORTRAN95 is detailed, including strategies for parallelization using domain decomposition. Results from a test suite which includes problems in one-, two-, and three-dimensions for both hydrodynamics and MHD are given, not only to demonstrate the fidelity of the algorithms, but also to enable comparisons to other methods. The source code is freely available for download on the web.

Magnetorotational-type Instability in Couette-Taylor Flow of a Viscoelastic Polymer Liquid

We describe an instability of viscoelastic Couette-Taylor flow that is directly analogous to the magnetorotational instability (MRI) in astrophysical magnetohydrodynamics, with polymer molecules playing the role of magnetic field lines. By determining the conditions required for the onset of instability and the properties of the preferred modes, we distinguish it from the centrifugal and elastic instabilities studied previously. Experimental demonstration and investigation should be much easier for the viscoelastic instability than for the MRI in a liquid metal. The analogy holds with the case of a predominantly toroidal magnetic field such as is expected in an accretion disk, and it may be possible to access a turbulent regime in which many modes are unstable.

A high order Godunov scheme with constrained transport and adaptive mesh refinement for astrophysical magnetohydrodynamics

Aims. In this paper, we present a new method to perform numerical simulations of astrophysical MHD flows using the Adaptive Mesh Refinement framework and Constrained Transport. Methods. The algorithm is based on a previous work in which the MUSCL-Hancock scheme was used to evolve the induction equation. In this paper, we detail the extension of this scheme to the full MHD equations and discuss its properties. Results. Through a series of test problems, we illustrate the performances of this new code using two different MHD Riemann solvers (Lax-Friedrich and Roe) and the need of the Adaptive Mesh Refinement capabilities in some cases. Finally, we show its versatility by applying it to two completely different astrophysical situations well studied in the past years: the growth of the magnetorotational instability in the shearing box and the collapse of magnetized cloud cores.Conclusions. We have implemented a new Godunov scheme to solve the ideal MHD equations in the AMR code RAMSES. We have shown that it results in a powerful tool that can be applied to a great variety of astrophysical problems, ranging from galaxies formation in the early universe to high resolution studies of molecular cloud collapse in our galaxy.

Constrained Transport Algorithms for Numerical Relativity. I. Development of a Finite‐Difference Scheme

A scheme is presented for accurately propagating the gravitational field constraints in finite difference implementations of numerical relativity. The method is based on similar techniques used in astrophysical magnetohydrodynamics and engineering electromagnetics, and has properties of a finite differential calculus on a four-dimensional manifold. It is motivated by the arguments that 1) an evolutionary scheme that naturally satisfies the Bianchi identities will propagate the constraints, and 2) methods in which temporal and spatial derivatives commute will satisfy the Bianchi identities implicitly. The proposed algorithm exactly propagates the constraints in a local Riemann normal coordinate system; {\it i.e.}, all terms in the Bianchi identities (which all vary as $\partial^3 g$) cancel to machine roundoff accuracy at each time step. In a general coordinate basis, these terms, and those that vary as $\partial g\partial^2 g$, also can be made to cancel, but differences of connection terms, proportional to $(\partial g)^3$, will remain, resulting in a net truncation error. Detailed and complex numerical experiments with four-dimensional staggered grids will be needed to completely examine the stability and convergence properties of this method. If such techniques are successful for finite difference implementations of numerical relativity, other implementations, such as finite element (and eventually pseudo-spectral) techniques, might benefit from schemes that use four-dimensional grids and that have temporal and spatial derivatives that commute.

New numerical tools to study waves and instabilities of flowing plasmas

Studying plasma waves and instabilities is an indispensable part of present thermonuclear fusion and astrophysical magnetohydrodynamics (MHD). Up till recently, spectral analysis was mostly restricted to static plasmas. However, the assumption of a static plasma is unrealistic not only for astrophysical but also for modern fusion research. Plasmas with flow have been shown to have spectra essentially different from those of static plasmas [Phys. Rev. Lett. 84 (2000) 2865]. We present two new numerical tools for spectral studies of plasmas with flow. The first one, a program called FINESSE (Finite Element Solver for Stationary Equilibria), computes equilibria of non-static plasmas for a variety of fusion and astrophysical configurations (tokamaks, solar loops, solar winds, etc.). Ideal and resistive spectra of the computed equilibria are studied with another tool, a program called PHOENIX. In PHOENIX, the large-scale generalized eigenproblems are solved with the recently proposed iterative Jacobi–Davidson method [SIAM J. Matrix Anal. Appl. 17 (1996) 401]. Our numerical examples show how FINESSE and PHOENIX can be used to study the effect of the poloidal flows on Toroidal Alfven Eigenmodes.

Criteria for Scaled Laboratory Simulations of Astrophysical MHD Phenomena

We demonstrate that two systems described by the equations of the ideal magnetohydrodynamics (MHD) evolve similarly, if the initial conditions are geometrically similar and certain scaling relations hold. The thermodynamic properties of the gas must be such that the internal energy density is proportional to the pressure. The presence of the shocks is allowed. We discuss the applicability conditions of the ideal MHD and demonstrate that they are satisfied with a large margin both in a number of astrophysical objects, and in properly designed simulation experiments with high-power lasers. This allows one to perform laboratory experiments whose results can be used for quantitative interpretation of various effects of astrophysical MHD.

The ZEUS code for astrophysical magnetohydrodynamics: new extensions and applications

A brief review of the numerical methods implemented in the widely used ZEUS code for astrophysical magnetohydrodynamics (MHD) is given. Extensions of the numerical methods to treat problems in nonideal MHD in a variety of regimes are discussed. In particular, methods for treating Ohmic dissipation, and methods for studying partially ionized plasmas in which the ion and neutral components are weakly coupled through a collisional drag term are considered. Recent application of the methods to the study of the dynamics of accretion disks is described.