Magnetohydrodynamic Code Research Articles

RELAP5/Mod3.3 MHD module development and validation: WCLL-TBM mock-up model

Magnetohydrodynamic (MHD) phenomena significantly impact the design of magnetic-confinement fusion reactors, particularly in the context of breeding blankets (BB) utilizing liquid metals (LMs) as working fluids. These effects arise from the interaction between the electro-conductive flowing metal and the magnetic field used for plasma confinement in the reactor chamber. Induced electrical currents in the liquid generate Lorentz forces, thereby altering flow behaviour in comparison to standard hydrodynamic conditions. For example, MHD effects modify velocity distribution and mass transport within ducts, amplify pressure losses, and influence heat transfer mechanisms. Accurate estimation of these impacts is crucial for the effective design of a liquid metal breeding blanket. While computational tools are essential for fusion-related physical analyses, no comprehensive MHD code currently exists for simulating all relevant phenomena in a liquid metal blanket. In this context, models predicting both distributed and concentrated MHD pressure drops have been integrated into the thermal-hydraulic system code RELAP5/Mod3.3. The Verification and Validation (V&V) process compares code results to direct numerical simulations and experimental data. For validation, RELAP5 recreates experimental results of a Water Cooled Lithium Lead test blanket module at magnetic field intensities ranging from Ha=500−3000, confirming the reliability of the newly implemented MHD subroutines for predicting pressure drops within this parameter range.

MHD stability of spherical tokamak equilibria with non-monotonic q-profiles

We use the 3D magnetohydrodynamic (MHD) code M3D-C1 [Jardin et al., Comput. Sci. Discovery 5, 014002 (2012)] to examine the MHD stability and subsequent evolution of NSTX shot 129169. This discharge had a period with a non-monotonic safety factor profile, q (reversed shear), which was terminated by a MHD event that abruptly lowered the central safety factor, q0, and greatly reduced the peakedness of the pressure profile. We show that the equilibrium just before the MHD event occurred was linearly unstable to many pressure-driven infernal modes. Modes with toroidal mode number n≥3 all had rational surfaces very close to the minimum value of q. However, a non-resonant pressure-driven (1, 1) mode was also present, and this dominated the nonlinear evolution. The final state in the simulation, after the MHD activity subsided, had a reduced and flattened pressure profile and a nearly monotonic q-profile, in qualitative agreement with experimental results. The initial state was also unstable to the resistive interchange criteria in the reversed-shear region, but the final state was stable everywhere. The “double tearing mode” (DTM) does not appear to play a role in the MHD activity of this discharge. In Appendix A, we show that in a torus, the DTM is strongly stabilized by pressure, but it is destabilized in cylindrical geometry (which has been the most extensively analyzed in the literature).

Magnetic and thermal acceleration in extragalactic jets

Aims. Relativistic jets launched from active galactic nuclei accelerate up to highly relativistic velocities within a length scale of between a few parsecs and tens of parsecs. The precise way in which this process takes place is still unclear. While magnetic acceleration is known to be able to accelerate relativistic outflows, little attention has been paid to the role of thermal acceleration. The latter has been assumed to act only on compact regions very close to the central engine, and to become negligible on parsec scales. However, this holds under the assumption of small internal energies relative to the magnetic ones, and whether or not this assumption is valid and what happens when we drop this assumption are open questions. Methods. We used a 2D relativistic magnetohydrodynamical code to explore jet acceleration from subparsec to parsec scales. As initial conditions for our models, we used observational constraints on jet properties derived by means of very long-baseline interferometry observations for a Fanaroff Riley I radio galaxy, NGC 315. We investigated the parameter space established for this source and performed a number of simulations of magnetically, thermally, or kinetically dominated jets at injection, and compared our results with the observations. Furthermore, we employed different models to characterize our jets, involving different magnetic field configurations (i.e., force-free vs. nonforce-free) and varying shear layer thicknesses. Results. Our simulated jets show that when thermal energy is comparable to or exceeds magnetic energy, thermal acceleration becomes significant at parsec scales. This result has important consequences, potentially extending the acceleration region far beyond the collimation scales, as thermal acceleration can effectively operate within a conically expanding jet. In all the models, we find acceleration to be driven by expansion, as expected. A number of our models allow us to reproduce the acceleration and opening angles observed in NGC 315. Finally, our results indicate that disk-launched winds might play an important role in jet propagation. Namely, when the jet has an initial force-free magnetic field configuration, thicker shear layers are needed to shield the internal spine from the action of the external medium and thus delay the growth of instabilities.

Development of a tokamak magnetohydrodynamic code with the discontinuous Galerkin and Weighted Essentially Non-Oscillatory methods

In this study, we present the development of a new initial-value magnetohydrodynamic (MHD) code for toroidal geometry using discontinuous Galerkin (DG) and Weighted Essentially Non-Oscillatory (WENO) methods. The code utilizes a triangular mesh based on the flux of the fixed boundary equilibrium in the poloidal plane, which is uniformly divided in the toroidal direction. By solving the conservative perturbed generalized Lagrange multiplier (GLM) MHD model, the code can simulate both ideal and resistive toroidal plasmas. We validate the code by performing the linear calculations of the internal kink mode and tearing mode, and our results show good agreement with calculations conducted using an eigen-value MHD code. Furthermore, nonlinear simulations of the resistive internal kink mode and tearing mode are also carried out. We discuss the challenges and limitations of the code, as well as the future efforts to improve its performance. The development of this new initial-value MHD code opens up new possibilities for studying and understanding the dynamics of toroidal plasmas, contributing to the advancement of fusion research.

Mancha3D Code: Multipurpose Advanced Nonideal MHD Code for High-Resolution Simulations in Astrophysics

The Mancha3D code is a versatile tool for numerical simulations of magnetohydrodynamic (MHD) processes in solar/stellar atmospheres. The code includes nonideal physics derived from plasma partial ionization, a realistic equation of state and radiative transfer, which allows performing high-quality realistic simulations of magnetoconvection, as well as idealized simulations of particular processes, such as wave propagation, instabilities or energetic events. The paper summarizes the equations and methods used in the Mancha3D (Multifluid (-purpose -physics -dimensional) Advanced Non-ideal MHD Code for High resolution simulations in Astrophysics 3D) code. It also describes its numerical stability and parallel performance and efficiency. The code is based on a finite difference discretization and a memory-saving Runge–Kutta (RK) scheme. It handles nonideal effects through super-time-stepping and Hall diffusion schemes, and takes into account thermal conduction by solving an additional hyperbolic equation for the heat flux. The code is easily configurable to perform different kinds of simulations. Several examples of the code usage are given. It is demonstrated that splitting variables into equilibrium and perturbation parts is essential for simulations of wave propagation in a static background. A perfectly matched layer (PML) boundary condition built into the code greatly facilitates a nonreflective open boundary implementation. Spatial filtering is an important numerical remedy to eliminate grid-size perturbations enhancing the code stability. Parallel performance analysis reveals that the code is strongly memory bound, which is a natural consequence of the numerical techniques used, such as split variables and PML boundary conditions. Both strong and weak scalings show adequate performance up to several thousands of processors (CPUs).

Synthetic Mirnov diagnostic for the validation of experimental observations.

A synthetic Mirnov diagnostic has been developed to investigate the capabilities and limitations of an arrangement of Mirnov coils in terms of a mode analysis. Eight test cases have been developed, with different coil arrangements and magnetic field configurations. Three of those cases are experimental configurations of the stellarator Wendelstein 7-X. It is observed that, for a high triangularity of the flux surfaces, the arrangement of the coils plays a significant role in the exact determination of the poloidal mode number. For the mode analysis, torus and magnetic coordinates have been used. In most cases, the reconstruction of the poloidal mode number of a prescribed mode was found to be more accurate in magnetic coordinates. As an application, the signal of an Alfvén eigenmode, which has been calculated with a three-dimensional magnetohydrodynamics code, is compared to experimental observations at Wendelstein 7-X. For the chosen example, the calculated and measured mode spectra agree very well and additional information on the toroidal mode number and localization of the mode has been inferred.

TokaMaker: An open-source time-dependent Grad-Shafranov tool for the design and modeling of axisymmetric fusion devices

In this paper, we present a new static and time-dependent MagnetoHydroDynamic (MHD) equilibrium code, TokaMaker, for axisymmetric configurations of magnetized plasmas, based on the well-known Grad-Shafranov equation. This code utilizes finite element methods on an unstructured triangular grid to enable capturing accurate machine geometry and simple mesh generation from engineering-like descriptions of present and future devices. The new code is designed for ease of use without sacrificing capability and speed through a combination of Python, Fortran, and C/C++ components. A detailed description of the numerical methods of the code, including a novel formulation of the boundary conditions for free-boundary equilibria, and validation of the implementation of those methods using both analytic test cases and cross-code validation is shown. Results show expected convergence across tested polynomial degree for analytic and cross-code test cases.

Enhancement of ECCD by the current condensation effect for stabilizing large magnetic islands caused by neoclassical tearing modes in tokamak plasmas

The radio frequency current condensation effect reported in Reiman and Fisch (2018 Phys. Rev. Lett. 121 225001) is modeled in the nonlinear resistive magnetohydrodynamic code. A series of numerical investigations have been performed to investigate the enhancement of electron cyclotron current drive (ECCD) by the current condensation effect during the control of neoclassical tearing mode (NTM) in tokamak plasmas. In the numerical model, both the parallel transport and the perpendicular transport of electron temperature are considered. The EC driven current and driven perturbed electron temperature can nonlinearly evolve within the given magnetic configuration and eventually reach saturation states. The input power threshold of ECCD and the fold bifurcation phenomenon are numerically verified via nonlinear simulations. The numerical results show good agreements with the analytical results. Moreover, spatial distributions of EC current for the two solutions at different condensed level are displayed. The control effectiveness of ECCD for large NTM islands has been evaluated while considering the current condensation effect. While taking into account current condensation effect, for a sufficiently large input power, a larger island can be more effectively stabilized than a smaller one, which suggests a reassessment of the previous idea that the ECCD should always be turned on as early as possible. The potential physics mechanism behind the ECCD control have all been discussed in detail. Furthermore, the condensation effect is found to have favorable effects on the radial misalignment of ECCD. In the consideration of the situation for extremely localized control needs, a highly peaked heating profile is adopted to verify that the fold bifurcation phenomenon still exists and the current condensation effect can still take effect in this extreme condition.

Star-forming environments in smoothed particle magnetohydrodynamics simulations II: re-simulating isolated clumps to determine equivalence of extracted clumps and parent simulations

ABSTRACT What is the numerical reproducibility of a stellar system (including its discs) when evolving only a subset of (partially-evolved) smoothed particle hydrodynamics (SPH) particles? To investigate this, we modelled the evolution of 29 star-forming clumps that were extracted from our previous simulations that investigated the formation and early evolution of low-mass star clusters. These clumps were evolved using a three-dimensional smoothed particle radiation magnetohydrodynamics code, where we included or excluded non-ideal magnetohydrodynamics to match the cluster simulation. While star formation proceeded as expected, we were unable to identically reproduce any of the systems present at the end of the cluster simulations. However, the final distributions of stellar mass, stellar system mass, disc mass, and disc radii were reproduced statistically; unfortunately, the distribution of average magnetic field strengths in the discs was not reproduced statistically, but this may be a result of our updated algorithms governing the evolution of the magnetic field. Therefore, given that our clumps yield stellar masses that are statistically similar to those in the original low-mass star clusters, we have demonstrated that we can statistically reproduce systems (aside from their magnetic field strength) by evolving a subset of SPH particles. Therefore, clumps such as these can be used as initial conditions to investigate the formation of isolated stars from less-contrived initial environments.

Are Electric-field-driven Magnetohydrodynamic Simulations of the Solar Corona Sensitive to the Initial Condition?

Magnetohydrodynamic (MHD) simulations of the solar corona are often started from a potential field initial condition, which may be far from the true state of the coronal magnetic field. To test the sensitivity of the simulations to the initial condition, an electric-field-driven MHD code was initiated from three different initial conditions: the commonly used potential field, a nonlinear force-free field (NLFFF) extrapolation, which is believed to be a better representation of the true coronal magnetic field, and a purely vertical field. After relaxation, the magnetic field in the MHD simulations was compared. It was found that there is a tendency for all of the initial conditions to converge to a similar connectivity, particularly in areas of strong fields that are not close to the boundary, suggesting that the final state of the simulation is not particularly sensitive to the initial condition. The relaxed solution more closely resembles the NLFFF initial condition than the other initial conditions, but still exhibits substantial differences from it.

Testing MURaM and MPS-ATLAS against the quiet solar spectrum

Context. Three-dimensional (3D) radiative magnetohydrodynamic (MHD) simulations are the only way to model stellar atmospheres without any ad hoc parameterisations. Several 3D radiative MHD codes have achieved good quantitative agreement with observables for our Sun. Aims. We aim to validate the most recent version of the MURaM code by comparing MURaM results to well-established quiet-Sun measurements, in particular spatially averaged measurements that are relevant for stellar studies. This extends the number of solar observables that MURaM can reproduce with high precision. Our validation is an essential condition to ensure that MURaM can be used to accurately calculate the spectra of other cool stars. Methods. We simulated the solar photosphere and upper convection zone, which harbours a small-scale-dynamo. Using time series of 3D snapshots, we calculated the spectral irradiance, limb darkening, and selected spectral lines, which we compared to observations. Results. The computed observables agree well with the observations; in particular, the limb darkening of the quiet Sun is reproduced remarkably well.

MUSES: A nonlinear magnetohydrodynamics discontinuous Galerkin code for fusion plasmas

The article describes initial steps toward the development of our nonlinear magnetohydrodynamics code, MUSES. Unstructured grids and general coordinates are employed so that MUSES code is expected to have a capability to represent realistic geometries of fusion devices. Discontinuous Galerkin method is chosen as a code framework because of its high spatial accuracy and numerical stability on unstructured grids. The governing equation of the ideal magnetohydrodynamics in this work is Powell 8-wave model which can mitigate the divergence-free issue. Roe scheme is used as a numerical flux since the Powell 8-wave model is a nonconservative system. A cross code comparison is carried out with a linear magnetohydrodynamics code, MINERVA. Owing to discontinuous Galerkin method, a linear theory of an ideal internal kink mode is reproduced quantitatively although the upwind method is employed for numerical stability.

Bridging Scales in Black Hole Accretion and Feedback: Magnetized Bondi Accretion in 3D GRMHD

Fueling and feedback couple supermassive black holes (SMBHs) to their host galaxies across many orders of magnitude in spatial and temporal scales, making this problem notoriously challenging to simulate. We use a multi-zone computational method based on the general relativistic magnetohydrodynamic (GRMHD) code KHARMA that allows us to span 7 orders of magnitude in spatial scale, to simulate accretion onto a non-spinning SMBH from an external medium with a Bondi radius of R B ≈ 2 × 105 GM •/c 2, where M • is the SMBH mass. For the classic idealized Bondi problem, spherical gas accretion without magnetic fields, our simulation results agree very well with the general relativistic analytic solution. Meanwhile, when the accreting gas is magnetized, the SMBH magnetosphere becomes saturated with a strong magnetic field. The density profile varies as ∼r −1 rather than r −3/2 and the accretion rate is consequently suppressed by over 2 orders of magnitude below the Bondi rate . We find continuous energy feedback from the accretion flow to the external medium at a level of . Energy transport across these widely disparate scales occurs via turbulent convection triggered by magnetic field reconnection near the SMBH. Thus, strong magnetic fields that accumulate on horizon scales transform the flow dynamics far from the SMBH and naturally explain observed extremely low accretion rates compared to the Bondi rate, as well as at least part of the energy feedback.

Study accretion and ejection using a new GPU-accelerated GRMHD code

AbstractWe study disks and jets in various accretion states (SANE and MAD) using novel, GPU-accelerated general-relativistic magneto-hydrodynamic (GR-MHD) code which we developed, based on HARM. This code, written in CUDA-c and uses OpenMP to parallelize multi-GPU setups, allows high resolution simulations of accretion disks and the formation and structure of jets without the need of multi-node supercomputer infrastructure. A 2563 simulation is well within the reach of an Nvidia DGX-V100 server, with the computation being a factor about 100 times faster if only the CPU was used.We use this code to examine several disk structures, wind and jet properties in the MAD and SANE states. In the MAD state, we find that the magnetic flux threading the horizon mostly depends on the spin of the BH. This implies that the jet structure and power are strong functions of the spin, with non-spinning BHs have the widest jets.

Weibel-like instability in magnetohydrodynamics

In magnetohydrodynamics (MHD), a density perturbation perpendicular to an electron temperature gradient generates a magnetic field around itself that acts to increase the perturbation, which can lead to instability. An MHD dispersion relation is obtained for perturbations perpendicular to a fixed electron temperature gradient with an initial in-plane magnetic field, including resistivity, viscosity, and the electrothermal coefficient. Instability occurs for sufficiently small electron temperature-gradient scale lengths determined by the ion collisionless skin depth. Both viscosity and resistivity are required to prevent growth at arbitrarily small spatial scales and to give a physical result for the fastest growing mode. The perpendicular electrothermal coefficient is only significant for a narrow range of low electron Hall parameters, causing a modest reduction in magnetic field growth and modifying the criteria for instability in the presence of viscosity. If the definition of the Weibel instability [E. S. Weibel, Phys. Rev. Lett. 2, 83 (1959)] is extended to include all instabilities due to anisotropy in the electron velocity distribution, then this is a Weibel-like instability because an electron temperature gradient implies an anisotropic electron velocity distribution. The implications for the formation of filaments in laser-produced plasmas and for the verification of MHD codes are considered.

Magnetohydrodynamic features and particle transport during thermal quench in HL-2A with density source

Diverse disruption mitigation strategies based on massive material injection have been developed in recent decades, aiming to uniformly deplete the thermal energy stored within plasmas during the thermal quench (TQ) while simultaneously elevating electron density to facilitate runaway electron suppression. Irrespective of the detailed dynamics of the material delivery scheme, deposition location and subsequent density mixing are pivotal in achieving highly efficient mitigation, however, which are markedly influenced by magnetohydrodynamic (MHD) activities. In order to assess the influence of MHD-induced transport on disruption mitigation, a simulation of TQ triggered by pure deuterium (D) deposition is conducted using a three-dimensional (3D) nonlinear reduced MHD code, JOREK. Steady density sources (the deposition rate of 1024 D atoms per second is greater than in real experiments) are introduced at various locations to explore the dynamics. The findings distinctly reveal two types of TQ processes, contingent on locations of deposition (LoD) of the neutral D source. Evidently, the results underscore the effectiveness of proper density mixing and moderated MHD in disruption mitigation. Nonlinear mode coupling emerges as a significant factor in shaping the final outcomes of TQ. Specifically, the 5/2 mode contributes to edge collapse, whereas the 3/2 mode is instrumental in core collapse. Additionally, the investigation indicates that the rise in core density is contingent on LoD, exhibiting threshold behavior. This threshold is observed within the q = 2/1 surface of equilibria, and a rapid increase in core density is witnessed when the density source crosses this threshold. The outcomes point toward the important role of E × B convection due to the 1/1 mode evolution in the process.

Self-consistent coupling of JOREK and CARIDDI: On the electromagnetic interaction of 3D tokamak plasmas with 3D volumetric conductors

The self-consistent electromagnetic integration of tokamak plasma magnetohydrodynamic (MHD) models with external conductors is considered in this manuscript. We integrate the extended MHD code JOREK and the eddy current code CARIDDI, via a scheme based on the virtual casing principle enabling two-way coupled fully implicit simulations. The robustness of this new integrated tool is validated against JOREK-STARWALL and CarMa0NL simulations for vertical displacement events in a simplified ASDEX Upgrade-like geometry and for n = 1 tearing modes in a circular high-aspect-ratio tokamak. The impact of 3D wall structures on the plasma evolution is highlighted introducing ports at different locations of a simplified axisymmetric wall and characterizing its effect on the plasma evolution. In addition to realistically capturing the 3D interactions between plasma and wall structures, the coupled codes will allow in the future 1:1 comparisons to experiments using virtual magnetic diagnostics that capture screening effects accurately.

Modelling solar coronal magnetic fields with physics-informed neural networks

ABSTRACT We present a novel numerical approach aiming at computing equilibria and dynamics structures of magnetized plasmas in coronal environments. A technique based on the use of neural networks that integrates the partial differential equations of the model, and called physics-informed neural networks (PINNs), is introduced. The functionality of PINNs is explored via calculation of different magnetohydrodynamic (MHD) equilibrium configurations, and also obtention of exact two-dimensional steady-state magnetic reconnection solutions. Advantages and drawbacks of PINNs compared to traditional numerical codes are discussed in order to propose future improvements. Interestingly, PINNs is a meshfree method in which the obtained solution and associated different order derivatives are quasi-instantaneously generated at any point of the spatial domain. We believe that our results can help to pave the way for future developments of time dependent MHD codes based on PINNs.

PHARE: Parallel hybrid particle-in-cell code with patch-based adaptive mesh refinement

Modeling multi-scale collisionless magnetized processes constitutes an important numerical challenge. By treating electrons as a fluid and ions kinetically, the so-called hybrid Particle-In-Cell (PIC) codes represent a promising intermediary between fully kinetic codes, limited to model small scales and short durations, and magnetohydrodynamic codes used large scale. However, simulating processes at scales significantly larger than typical ion particle dynamics while resolving sub-ion dissipative current sheets remain extremely difficult. This paper presents a new hybrid PIC code with patch-based adaptive mesh refinement. Here, hybrid PIC equations are solved on a hierarchy of an arbitrary number of Cartesian meshes of incrementally finer resolution dynamically mapping regions of interest, and with a refined time stepping. This paper presents how the hybrid PIC algorithm is adapted to evolve such mesh hierarchy and the validation of the code on a uniform mesh, fixed refined mesh and dynamically refined mesh.

GRaM-X: a new GPU-accelerated dynamical spacetime GRMHD code for Exascale computing with the Einstein Toolkit

We present GRaM-X (General Relativistic accelerated Magnetohydrodynamics on AMReX), a new GPU-accelerated dynamical-spacetime general relativistic magnetohydrodynamics (GRMHD) code which extends the GRMHD capability of Einstein Toolkit to GPU-based exascale systems. GRaM-X supports 3D adaptive mesh refinement (AMR) on GPUs via a new AMR driver for the Einstein Toolkit called CarpetX which in turn leverages AMReX, an AMR library developed for use by the United States DOE’s Exascale Computing Project. We use the Z4c formalism to evolve the Einstein equations and the Valencia formulation to evolve the equations of GRMHD. GRaM-X supports both analytic as well as tabulated equations of state. We implement TVD and WENO reconstruction methods as well as the HLLE Riemann solver. We test the accuracy of the code using a range of tests on static spacetime, e.g. 1D magnetohydrodynamics shocktubes, the 2D magnetic rotor and a cylindrical explosion, as well as on dynamical spacetimes, i.e. the oscillations of a 3D Tolman-Oppenheimer-Volkhof star. We find excellent agreement with analytic results and results of other codes reported in literature. We also perform scaling tests and find that GRaM-X shows a weak scaling efficiency of ∼40%–50% on 2304 nodes (13824 NVIDIA V100 GPUs) with respect to single-node performance on OLCF’s supercomputer Summit.