MHD by adamdempsey90 · Pull Request #123 · lanl/artemis

adamdempsey90 · 2026-03-31T16:02:53Z

Background

This adds ideal MHD. It should work with AMR and all coordinate systems. I have tested on the magnetized blast problem in cartesian/cylindrical coordinates with/without AMR (modified the existing pgen), the Orszag-Tang vortex (new pgen), and the Brio-Wu shock tube (new input file).

The outline of new features/changed things:

Added current unit scaling parameter and $\mu_0$.
LLF/HLLE modified to compute the MHD fluxes.
Added HLLD (this is horrible btw) basically from Athena++
Two new tasks in the driver to compute the edge EMFs from the magnetic face fluxes and then a task to update the face fields from the edge EMFs.
New prolongation/restriction operators for edge and face fields.
New PostInitialization hook that sets primitive B from face B. Face B is set in the pgens.
New edge geometry quantities.
Changes are threaded through dual energy.
New documentation

Overall, it's not too many new things and there isn't much in the input file that needed to change.

Some things I'm not too happy with that we might want to change:

For speed I choose to just use the Alfven speed in the timestep determination, rather that the max of all the fast magnetosonic speeds one could calculate at all the faces. But obviously, if we take fewer timesteps at the cost of a more expensive EstimateTimeStep kernel, it is probably worth it.
HLLD is a mess but I don't know of any other way to do it.
I'm sure there are BC things we need to add, but I'm not sure at the moment. It's on the pgen to do what it needs to do.

This was done in concert with GPT.

I haven't tested on GPUs yet or anything high res. Initial testing looks pretty good though, I added a diagnostic div(B) field at that remains at machine epsilon even through AMR in cylindrical/cartesian coordinates.

Description of Changes

Closes #122

Checklist

New features are documented
Tests added for bug fixes and new features
(@lanl.gov employees) Update copyright on changed files
Any contribution that was created or modified with the assistance of generative AI must have a comment disclosing this such as // This file was created in part or in whole by generative AI

…olvers.

…he LLF face fluxes. Only do mhd fluxes for n=0 species. Add a divB diagnostic variable.

…e cell field from the face field

…mpere now

…ion. This comes at increased complexity.

adamdempsey90 · 2026-07-30T22:33:22Z


        // Add the "flux source terms"
        for (int n = 0; n < nspecies; ++n) {
+          const bool mhd = (F == Fluid::gas) && do_mhd && (n == 0);


Suggested change

const bool mhd = (F == Fluid::gas) && do_mhd && (n == 0);

const bool mhd = do_mhd && (F == Fluid::gas) && (n == 0);

adamdempsey90 · 2026-08-03T17:37:21Z

This is now ready for review @pdmullen

pdmullen

I'd like to see the linear modes test adapted to MHD and the field loop advection test before this goes in.

pdmullen · 2026-08-13T00:54:28Z

@@ -0,0 +1,2050 @@
+# Athena++ data at time=8.000000e-02  cycle=1559  variables=prim 


Does this need an Athena++ copyright?

pdmullen · 2026-08-13T01:36:26Z

+          const auto xv = coords.GetCellCenter(vg, b, k, j, i);
+          const auto &bnds = coords.GetBounds();
+          const Real bx =
+              ((bnds.x1[1] - xv[0]) * vmesh(b, TE::F1, field::face::B(), k, j, i) +


Maybe we should make a helper function for this. Maybe it can be templated on coordinates so that we can just do arithmetic average for cartesian.

pdmullen · 2026-08-13T01:40:00Z

+                          vmesh(b, TE::F3, field::face::B(), k + threed, j, i)) /
+                     (bnds.x3[1] - bnds.x3[0]))
+                  : vmesh(b, TE::F3, field::face::B(), k, j, i);
+          emag = MHD::MagneticEnergyDensity(bx, by, bz, mu0_code);


Maybe a helper function for cell-centered fields would be nice too. Would clean up duplicated above with below in this file.

pdmullen · 2026-08-13T01:41:56Z

            // Compute SIE via dual energy formalism and apply floor
-            Real sieg = ArtemisUtils::DualEnergySIE(vmesh, b, n, k, j, i, de_switch, hx);
+            const Real emag =
+                (do_mhd && (n == 0)) ? vmesh(b, field::cell::energy(), k, j, i) : 0.0;


why only n==0? I suspect that this is just a placeholder bc not all physics packages respect multiple gas species... but the magnetic field should inform all fields if all are coupled.

pdmullen · 2026-08-13T01:42:27Z

          // Compute SIE via dual energy formalism and apply floor
-          Real sie = ArtemisUtils::DualEnergySIE(vmesh, b, n, k, j, i, de_switch, hx);
+          const Real emag =
+              (do_mhd && (n == 0)) ? vmesh(b, field::cell::energy(), k, j, i) : 0.0;


same comments as above.... won't repeat hereafter.

pdmullen · 2026-08-13T02:16:59Z

  const bool do_moment = do_radiation && pin->DoesBlockExist("radiation/moment");
  const bool update_fluxes = pin->GetOrAddBoolean("gas", "update_fluxes", true);

  // Check configuration selection compatibility


Shouldn't mhd require do_gas and update_fluxes (with the latter so that we can get upwind direction)?

pdmullen · 2026-08-13T02:19:19Z

+              const Real t1 = SQR(vp_(b, IVX, k, j, i) + rfv[0]) -
+                              (mhd ? SQR(vp_(b, field::cell::B(0), k, j, i)) / mu0 : 0.0);
+              const Real t2 = SQR(vp_(b, IVY, k, j, i) + rfv[1]) -
+                              (mhd ? SQR(vp_(b, field::cell::B(1), k, j, i)) / mu0 : 0.0);
+              const Real t3 = SQR(vp_(b, IVZ, k, j, i) + rfv[2]) -
+                              (mhd ? SQR(vp_(b, field::cell::B(2), k, j, i)) / mu0 : 0.0);


should these divide out rho since wdt contains rho?

pdmullen · 2026-08-13T02:29:22Z

              // a = (a < 0.0) ? 0.0 : sqrt(gm1 * a);
              // Einfeldt (1988)
              const Real ngam = 0.5 * sqrtdl * sqrtdr * SQR(isdlpdr);
              Real a =


I am confused about the dimensionality of this... can you check?

qa is a fast magnetosonic speed and the second term has units of velocity squared?

pdmullen · 2026-08-13T02:30:16Z

            const Real ke = 0.5 * w_d * (SQR(vel1) + SQR(vel2) + SQR(vel3));
+            Real me = 0.0;
+            if (do_mhd && (n == 0)) {
+              const Real bx = 0.5 * (vmesh(b, TE::F1, field::face::B(), k, j, i) +


no fancy centroid stuff here?

pdmullen · 2026-08-13T02:32:11Z

+              coords.template GetScaleFactorsFace<X2DIR>(vg, b, k, j, i - 1);
+          const Real h3e = coords.template GetEdgeScaleFactor<X3DIR>(vg, b, k, j, i);
+          Real &emf = v.flux(b, TE::E3, field::face::B(), k, j, i);
+          emf = h3e * 0.25 *


straight averaging? no upwinding a la Gardiner & Stone?

adamdempsey90 and others added 23 commits March 27, 2026 11:38

Initial commit of mhd. Not working. Need geometry terms and riemann s…

2183c60

…olvers.

format

f5b879f

fix indexing issue

11d642a

deep copy conserved data correctly for conserved face fields

002f72d

Add new routine to assemble edge EMF from face fluxes. Correctly do t…

dcf5b73

…he LLF face fluxes. Only do mhd fluxes for n=0 species. Add a divB diagnostic variable.

Correct some errors

c5c5ccf

include magnetic energy in total energy

50a06a1

Add Brio-Wu and Orszag-Tang tes5t problems

b99ec56

Actually update the out-of-plane fields

bdf93bf

correctly handle the boundaries when mhd is on

bcf0b5d

hlle + mhd

a1016fd

add hlld. It's a mess and should be refactored

b7e8705

make changes necessary for mhd in generic coordinates

5e5682a

assemble the edge emfs before flux correction

8dcde76

set both cell and face fields

321f63e

format

8f14247

Add edge support to coordinates

5dbd6e2

non-cartesian coordinate support for mhd

4ef97d1

amr support for face and edge fields

8529742

Move mhd things from utils to mhd.hpp. Add a post init hook to set th…

84830e5

…e cell field from the face field

Proper mhd unit support. User can set the current scaling factor to a…

2f58a0c

…mpere now

mhd docs

eaa494e

remove unused mhd variables

9dd4f29

adamdempsey90 requested a review from pdmullen March 31, 2026 16:02

adamdempsey90 added 6 commits April 15, 2026 17:20

Capture

806f1b9

Merge remote-tracking branch 'origin/develop' into dempsey/mhd

b8d45cb

PAR_VAR

1b31713

FluxSource depends on update_mhd

5f63729

Only initialize fields if they are non-zero

91e35ff

Use b/sqrt(mu) everywhere

fdac7a1

adamdempsey90 added 5 commits June 16, 2026 13:32

Add div(B) monitoring. Fix issues with div(B) growth during prolongat…

78dbc13

…ion. This comes at increased complexity.

Format

2ea6262

Merge remote-tracking branch 'origin/develop' into dempsey/mhd

c1cb9f9

See if this fixes sie recovery

fb8f292

Merge remote-tracking branch 'origin/develop' into dempsey/mhd

e4146e4

adamdempsey90 commented Jul 30, 2026

View reviewed changes

adamdempsey90 added 17 commits July 30, 2026 16:47

Use isfinite for the comparisons to Null...

4514b76

format

7337901

Fix pressure issue in hlld

8c78000

Add Brio-Wu shock tube to test framework

f357014

Fix energy flux issue with mhd

d560dd0

Test all riemann solvers and format

24dfbd7

format

5b638dd

Fix energy using stale magnetic field

8825c1e

Correct format specifier

5559bd2

Parallel brio-wu

8197f05

Add orszag-tang

a592574

Add mhd tests to parallel suite

587d0ee

format

f1fe23d

Fix script and lower tolerance

c06a631

Update script to output a plot

9d217cf

Update script again, add athena reference

1e2d693

format

af4cc68

pdmullen requested changes Aug 13, 2026

View reviewed changes

Navigation Menu

Uh oh!

MHD - #123

MHD#123
adamdempsey90 wants to merge 51 commits into
developfrom
dempsey/mhd

adamdempsey90 commented Mar 31, 2026

Uh oh!

adamdempsey90 Jul 30, 2026

Uh oh!

adamdempsey90 commented Aug 3, 2026

Uh oh!

pdmullen left a comment

Uh oh!

pdmullen Aug 13, 2026

Uh oh!

pdmullen Aug 13, 2026

Uh oh!

pdmullen Aug 13, 2026

Uh oh!

pdmullen Aug 13, 2026

Uh oh!

pdmullen Aug 13, 2026

Uh oh!

pdmullen Aug 13, 2026

Uh oh!

pdmullen Aug 13, 2026

Uh oh!

pdmullen Aug 13, 2026

Uh oh!

pdmullen Aug 13, 2026

Uh oh!

pdmullen Aug 13, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

	const bool mhd = (F == Fluid::gas) && do_mhd && (n == 0);
	const bool mhd = do_mhd && (F == Fluid::gas) && (n == 0);

		@@ -0,0 +1,2050 @@
		# Athena++ data at time=8.000000e-02 cycle=1559 variables=prim

Uh oh!

Conversation

adamdempsey90 commented Mar 31, 2026

Background

Description of Changes

Checklist

Uh oh!

Choose a reason for hiding this comment

Uh oh!

adamdempsey90 commented Aug 3, 2026

Uh oh!

pdmullen left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants