MHD

doc/src/physics.rst

@@ -81,6 +81,63 @@ An example that triggers refinement when the gas density is larger than ``10.0``
    deref_thr = 0.02
 
 
+MHD
+"""
+
+Ideal magnetohydrodynamics is enabled by setting ``mhd = true`` in the ``<physics>`` block.
+When MHD is active, |code| evolves a face-centered magnetic field using constrained transport (CT).
+This keeps the magnetic update tied to the divergence-preserving CT representation.
+
+At minimum, an MHD run should enable both gas and MHD:
+
+::
+
+   <physics>
+   gas = true
+   mhd = true
+
+If there are multiple gas species defined, only the first species couples to the magnetic field, i.e., its total energy contains the magnetic energy and only its velocity induces an electric fied.
+Magnetic fluxes are computed through the gas Riemann solver selected in ``<gas>``.
+The available solver choices are:
+
+* ``llf``
+* ``hlle``
+* ``hlld``
+
+The ``hlld`` solver is the most specialized MHD option currently available.
+An example input block that enables MHD is,
+
+::
+
+   <physics>
+   gas = true
+   mhd = true
+
+   <gas>
+   riemann = hlle
+   reconstruct = plm
+
+   <gas/eos/ideal>
+   gamma = 1.4
+
+In dimensional runs, magnetic quantities use the same unit-conversion system as the rest of the code.
+When ``artemis/physical_units = cgs`` and ``artemis/unit_conversion = base``, the base-unit specification now includes electrical current:
+
+::
+
+   <artemis>
+   physical_units = cgs
+   unit_conversion = base
+   length = 1.0
+   time = 1.0
+   mass = 1.0
+   temperature = 1.0
+   current = 1.0
+
+The ``current`` parameter is used together with length, time, and mass to derive the magnetic-field conversion and the code-unit value of :math:`\mu_0`.
+In ``scalefree`` runs, these conversions reduce to unity.
+
+
 Gasses support several microphysics models including `Cooling`_, `Viscosity`_, and `Conduction`_.
 These are controlled by adding additional nodes under the ``<gas>`` node, e.g., ``<gas/viscosity>``.
 We describe each of these below.

inputs/orszag_tang/orszag_tang.in

@@ -0,0 +1,84 @@
+# ========================================================================================
+#  (C) (or copyright) 2023-2026. Triad National Security, LLC. All rights reserved.
+#
+#  This program was produced under U.S. Government contract 89233218CNA000001 for Los
+#  Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+#  for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+#  in the program are reserved by Triad National Security, LLC, and the U.S. Department
+#  of Energy/National Nuclear Security Administration. The Government is granted for
+#  itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+#  license in this material to reproduce, prepare derivative works, distribute copies to
+#  the public, perform publicly and display publicly, and to permit others to do so.
+# ========================================================================================
+
+<artemis>
+problem = orszag_tang    # name of the pgen
+coordinates = cartesian  # coordinate system
+
+<parthenon/job>
+problem_id = orszag_tang  # problem ID: basename of output filenames
+
+<parthenon/output1>
+variables = gas.prim.density,  &
+            gas.prim.pressure, &
+            gas.prim.velocity, &
+            field.cell.B,      &
+            field.cell.divB
+file_type = hdf5  # HDF5 data dump
+dt        = 0.05  # time increment between outputs
+
+<parthenon/time>
+nlim       = -1    # cycle limit
+tlim       = 0.5   # time limit
+integrator = rk2   # time integration algorithm
+ncycle_out = 1    # interval for stdout summary info
+
+<parthenon/mesh>
+nghost = 2
+refinement = none
+numlevel = 1
+
+nx1    = 64       # Number of zones in X1-direction
+x1min  = 0.0       # minimum value of X1
+x1max  = 1.0       # maximum value of X1
+ix1_bc = periodic  # inner-X1 boundary flag
+ox1_bc = periodic  # outer-X1 boundary flag
+
+nx2    = 64       # Number of zones in X2-direction
+x2min  = 0.0       # minimum value of X2
+x2max  = 1.0       # maximum value of X2
+ix2_bc = periodic  # inner-X2 boundary flag
+ox2_bc = periodic  # outer-X2 boundary flag
+
+nx3    = 1         # Number of zones in X3-direction
+x3min  = -0.5      # minimum value of X3
+x3max  = 0.5       # maximum value of X3
+ix3_bc = periodic  # inner-X3 boundary flag
+ox3_bc = periodic  # outer-X3 boundary flag
+
+<parthenon/meshblock>
+nx1 = 64  # Number of cells in each MeshBlock, X1-dir
+nx2 = 64  # Number of cells in each MeshBlock, X2-dir
+nx3 = 1   # Number of cells in each MeshBlock, X3-dir
+
+<physics>
+gas = true
+mhd = true
+dust = false
+
+<gas>
+cfl = 0.4
+reconstruct = plm
+riemann = llf
+update_fluxes = true
+dfloor = 1.0e-10
+siefloor = 1.0e-10
+
+<gas/eos/ideal>
+gamma = 1.66666666667
+
+<problem>
+rho0 = 0.22104853207207686
+p0   = 0.13262911924324611
+v0   = 1.0
+b0   = 0.28209479177387814

inputs/shock/brio_wu.in

@@ -0,0 +1,92 @@
+# ========================================================================================
+#  (C) (or copyright) 2026. Triad National Security, LLC. All rights reserved.
+#
+#  This program was produced under U.S. Government contract 89233218CNA000001 for Los
+#  Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+#  for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+#  in the program are reserved by Triad National Security, LLC, and the U.S. Department
+#  of Energy/National Nuclear Security Administration. The Government is granted for
+#  itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+#  license in this material to reproduce, prepare derivative works, distribute copies to
+#  the public, perform publicly and display publicly, and to permit others to do so.
+# ========================================================================================
+
+<artemis>
+problem = shock          # name of the pgen
+coordinates = cartesian  # coordinate system
+
+<parthenon/job>
+problem_id = brio_wu     # problem ID: basename of output filenames
+
+<parthenon/output1>
+variables = gas.prim.density,  &
+            gas.prim.velocity, &
+            gas.prim.pressure, &
+            field.cell.B,      &
+            field.cell.divB
+file_type = hdf5  # HDF5 data dump
+dt        = 0.02  # time increment between outputs
+
+<parthenon/time>
+nlim       = -1    # cycle limit
+tlim       = 0.1   # time limit
+integrator = rk2   # time integration algorithm
+ncycle_out = 10
+
+<parthenon/mesh>
+nghost = 2
+refinement = none
+numlevel = 1
+
+nx1    = 512      # Number of zones in X1-direction
+x1min  = -0.5     # minimum value of X1
+x1max  = 0.5      # maximum value of X1
+ix1_bc = ic       # inner-X1 boundary flag
+ox1_bc = ic       # outer-X1 boundary flag
+
+nx2    = 1        # Number of zones in X2-direction
+x2min  = 0.0      # minimum value of X2
+x2max  = 1.0      # maximum value of X2
+ix2_bc = periodic # inner-X2 boundary flag
+ox2_bc = periodic # outer-X2 boundary flag
+
+nx3    = 1        # Number of zones in X3-direction
+x3min  = 0.0      # minimum value of X3
+x3max  = 1.0      # maximum value of X3
+ix3_bc = periodic # inner-X3 boundary flag
+ox3_bc = periodic # outer-X3 boundary flag
+
+<parthenon/meshblock>
+nx1 = 512  # Number of cells in each MeshBlock, X1-dir
+nx2 = 1    # Number of cells in each MeshBlock, X2-dir
+nx3 = 1    # Number of cells in each MeshBlock, X3-dir
+
+<physics>
+gas = true
+mhd = true
+dust = false
+
+<gas>
+cfl = 0.4
+reconstruct = plm
+riemann = llf
+update_fluxes = true
+dfloor = 1.0e-10
+siefloor = 1.0e-10
+
+<gas/eos/ideal>
+gamma = 2.0
+
+<problem>
+rhol = 1.0
+rhor = 0.125
+vxl = 0.0
+vxr = 0.0
+pl = 1.0
+pr = 0.1
+bx = 0.75
+byl = 1.0
+byr = -1.0
+bzl = 0.0
+bzr = 0.0
+xdisc = 0.0

src/CMakeLists.txt

@@ -48,6 +48,9 @@ set (SRC_LIST
   gravity/point_mass.cpp
   gravity/uniform.cpp
 
+  mhd/mhd.cpp
+  mhd/mhd.hpp
+
   nbody/nbody.cpp
   nbody/nbody.hpp
   nbody/nbody_advance.cpp

src/artemis.cpp

@@ -21,6 +21,7 @@
 #include "gas/gas.hpp"
 #include "geometry/geometry.hpp"
 #include "gravity/gravity.hpp"
+#include "mhd/mhd.hpp"
 #include "nbody/nbody.hpp"
 #include "radiation/moments/moments.hpp"
 #include "radiation/radiation.hpp"
@@ -102,6 +103,7 @@ Packages_t ProcessPackages(std::unique_ptr<ParameterInput> &pin) {
   const bool do_radiation = pin->GetOrAddBoolean("physics", "radiation", false);
   const bool do_coagulation = pin->GetOrAddBoolean("physics", "coagulation", false);
   const bool do_raytrace = pin->GetOrAddBoolean("physics", "raytrace", false);
+  const bool do_mhd = pin->GetOrAddBoolean("physics", "mhd", false);
   const bool do_orbital_advection =
       pin->GetOrAddBoolean("physics", "orbital_advection", false);
 
@@ -144,6 +146,7 @@ Packages_t ProcessPackages(std::unique_ptr<ParameterInput> &pin) {
   artemis->AddParam("do_moment", do_moment);
   artemis->AddParam("do_orbital_advection", do_orbital_advection);
   artemis->AddParam("do_raytrace", do_raytrace);
+  artemis->AddParam("do_mhd", do_mhd);
   artemis->AddParam("update_fluxes", update_fluxes);
 
   // Set coordinate system
@@ -152,6 +155,7 @@ Packages_t ProcessPackages(std::unique_ptr<ParameterInput> &pin) {
   Coordinates coords = geometry::CoordSelect(sys, ndim);
   artemis->AddParam("coords", coords);
   artemis->AddParam("coord_sys", sys);
+  artemis->AddParam("ndim", ndim);
 
   geometry::CoordParams cpars(pin.get());
   artemis->AddParam("coord_params", cpars);
@@ -202,6 +206,12 @@ Packages_t ProcessPackages(std::unique_ptr<ParameterInput> &pin) {
     }
   }
   if (do_raytrace) packages.Add(RT::Initialize(pin.get(), units, constants));
+  if (do_mhd) packages.Add(MHD::Initialize(pin.get(), units, constants, packages));
+
+  if (do_mhd && pin->GetOrAddBoolean("mhd", "monitor_divb", false)) {
+    artemis->PreStepDiagnosticsMesh = ArtemisUtils::PreStepDiagnosticsRemeshDivB;
+    artemis->PostStepDiagnosticsMesh = ArtemisUtils::PostStepDiagnosticsRemeshDivB;
+  }
 
   // Assign geometry-specific FillDerived functions
   if (do_gas || do_dust) {

src/artemis.hpp

@@ -96,6 +96,17 @@ PAR_VAR(grav, phi);
 PAR_VAR(grav, rhs);
 } // namespace grav
 
+namespace field {
+namespace face {
+PAR_VAR(field.face, B);
+} // namespace face
+namespace cell {
+PAR_VAR(field.cell, B);
+PAR_VAR(field.cell, energy);
+PAR_VAR(field.cell, divB);
+} // namespace cell
+} // namespace field
+
 namespace geom {
 PAR_VAR(geom, x1v);
 PAR_VAR(geom, x2v);
@@ -112,9 +123,15 @@ PAR_VAR(geom, hx3f2);
 PAR_VAR(geom, hx1f3);
 PAR_VAR(geom, hx2f3);
 PAR_VAR(geom, hx3f3);
+PAR_VAR(geom, hx1e1);
+PAR_VAR(geom, hx2e2);
+PAR_VAR(geom, hx3e3);
 PAR_VAR(geom, dx1);
 PAR_VAR(geom, dx2);
 PAR_VAR(geom, dx3);
+PAR_VAR(geom, dl1);
+PAR_VAR(geom, dl2);
+PAR_VAR(geom, dl3);
 PAR_VAR(geom, vol);
 PAR_VAR(geom, ax1);
 PAR_VAR(geom, ax2);
@@ -159,7 +176,7 @@ enum class Coordinates {
 };
 
 // ...Riemann solvers
-enum class RSolver { hllc_general, hlle, llf, hllc_gamma, null };
+enum class RSolver { hllc_general, hlle, llf, hllc_gamma, hlld, null };
 // ... Upwinding (left vs right state)
 enum class Upwind { l, r, null };
 // ...Reconstruction algorithms

src/artemis_driver.cpp

@@ -75,6 +75,7 @@ ArtemisDriver<GEOM>::ArtemisDriver(ParameterInput *pin, ApplicationInput *app_in
   do_moment = artemis_pkg->template Param<bool>("do_moment");
   do_coagulation = artemis_pkg->template Param<bool>("do_coagulation");
   do_raytrace = artemis_pkg->template Param<bool>("do_raytrace");
+  do_mhd = artemis_pkg->template Param<bool>("do_mhd");
 
   // Update fluxes option--gas fields are needed for radiation temperature updates but for
   // rad-only test problems turn off advection
@@ -273,9 +274,14 @@ TaskCollection ArtemisDriver<GEOM>::StepTasks() {
         diff_flx = vflx | tflx;
       }
 
+      TaskID edge_emf = none;
+      if (do_mhd) {
+        edge_emf = tl.AddTask(gas_flx, MHD::AssembleEdgeEMF, u0.get());
+      }
+
       // Communicate and set fluxes
       auto send_flx =
-          tl.AddTask(gas_flx | dust_flx | diff_flx,
+          tl.AddTask(gas_flx | dust_flx | diff_flx | edge_emf,
                      parthenon::SendBoundBufs<parthenon::BoundaryType::flxcor_send>, u0);
       auto recv_flx = tl.AddTask(start_flx_recv, parthenon::ReceiveFluxCorrections, u0);
       auto set_flx = tl.AddTask(recv_flx, parthenon::SetFluxCorrections, u0);
@@ -284,10 +290,16 @@ TaskCollection ArtemisDriver<GEOM>::StepTasks() {
       auto update =
           tl.AddTask(gas_flx | dust_flx | set_flx, ArtemisUtils::ApplyUpdate<GEOM>,
                      u0.get(), u1.get(), g0, g1, bdt);
+      auto update_mhd = gas_flx | set_flx;
+      if (do_mhd) {
+        update_mhd = tl.AddTask(edge_emf | set_flx, ArtemisUtils::ApplyFaceUpdate<GEOM>,
+                                u0.get(), u1.get(), g0, g1, bdt);
+      }
 
       // Apply "coordinate source terms"
-      TaskID gas_coord_src = update, dust_coord_src = update;
-      if (do_gas) gas_coord_src = tl.AddTask(update, Gas::FluxSource, u0.get(), bdt);
+      TaskID gas_coord_src = update | update_mhd, dust_coord_src = update;
+      if (do_gas)
+        gas_coord_src = tl.AddTask(update | update_mhd, Gas::FluxSource, u0.get(), bdt);
       if (do_dust) dust_coord_src = tl.AddTask(update, Dust::FluxSource, u0.get(), bdt);
 
       // Apply (gas) diffusion sources

src/artemis_driver.hpp

@@ -57,7 +57,7 @@ class ArtemisDriver : public EvolutionDriver {
   bool do_gravity, do_self_gravity, do_nbody, do_rotating_frame;
   bool do_orbital_advection;
   bool do_cooling, do_drag, do_viscosity, do_conduction, do_diffusion;
-  bool do_coagulation, do_raytrace;
+  bool do_coagulation, do_raytrace, do_mhd;
   bool update_fluxes;
   const bool is_restart;
 };