MHD interface reconstruction and HLLD flux path skip auxiliary passives (`QFX/UFX`)

[WeiqunZhang]

Summary

MHD reconstruction (plm, ppm_mhd) reconstructs species (QFS) but not auxiliary passives (QFX), yet later reads QFX for EOS. The HLLD solver similarly builds/pass-through species states but does not initialize flux/star-state entries for UFX.

Location

• Source/mhd/mhd_plm.cpp:279
• Source/mhd/mhd_plm.cpp:317
• Source/mhd/mhd_ppm.cpp:330
• Source/mhd/mhd_ppm.cpp:375
• Source/mhd/hlld.cpp:112
• Source/mhd/hlld.cpp:195
• Source/mhd/hlld.cpp:319

Problem Details

• plm/ppm_mhd: only QFS interface values are populated; QFX remains unset while EOS aux is read from interface states.
• hlld: only species are explicitly initialized for advection in FL/FR, Us*, and Uss*; auxiliary conserved components can remain undefined in states
  used to form final fluxes.

Impact

• Undefined aux values can enter EOS calls for interface thermodynamics.
• Auxiliary conserved fluxes can be uninitialized or inconsistent when NumAux > 0.
• Can produce non-deterministic updates and corrupted auxiliary fields.

Suggested Patch

diff --git a/Source/mhd/mhd_plm.cpp b/Source/mhd/mhd_plm.cpp
--- a/Source/mhd/mhd_plm.cpp
+++ b/Source/mhd/mhd_plm.cpp
@@
+    for (int n = 0; n < NumAux; n++) {
+      Real un{};
+      Real A[nslp];
+      load_stencil(s, idir, i, j, k, QFX+n, A);
+      if (idir == 0) un = s(i,j,k,QU);
+      else if (idir == 1) un = s(i,j,k,QV);
+      else un = s(i,j,k,QW);
+      Real dA = uslope(A, flatn(i,j,k), false, false);
+      if (idir == 0) qleft(i+1,j,k,QFX+n) = q_zone(QFX+n) + 0.5_rt*(1.0_rt - dtdx*un) * dA;
+      else if (idir == 1) qleft(i,j+1,k,QFX+n) = q_zone(QFX+n) + 0.5_rt*(1.0_rt - dtdx*un) * dA;
+      else qleft(i,j,k+1,QFX+n) = q_zone(QFX+n) + 0.5_rt*(1.0_rt - dtdx*un) * dA;
+      qright(i,j,k,QFX+n) = q_zone(QFX+n) - 0.5_rt*(1.0_rt + dtdx*un) * dA;
+    }

diff --git a/Source/mhd/mhd_ppm.cpp b/Source/mhd/mhd_ppm.cpp
--- a/Source/mhd/mhd_ppm.cpp
+++ b/Source/mhd/mhd_ppm.cpp
@@
+    for (int n = 0; n < NumAux; n++) {
+      Real Ips, Ims;
+      Real un = (idir == 0) ? q_arr(i,j,k,QU) : ((idir == 1) ? q_arr(i,j,k,QV) : q_arr(i,j,k,QW));
+      load_stencil(q_arr, idir, i, j, k, QFX+n, s);
+      ppm_reconstruct(s, flat, sm, sp);
+      ppm_int_profile_single(sm, sp, s[i0], un, dtdx, Ips, Ims);
+      if (idir == 0) qleft(i+1,j,k,QFX+n) = Ips;
+      else if (idir == 1) qleft(i,j+1,k,QFX+n) = Ips;
+      else qleft(i,j,k+1,QFX+n) = Ips;
+      qright(i,j,k,QFX+n) = Ims;
+    }

diff --git a/Source/mhd/hlld.cpp b/Source/mhd/hlld.cpp
--- a/Source/mhd/hlld.cpp
+++ b/Source/mhd/hlld.cpp
@@
     for (int n = 0; n < NumSpec; n++) {
       FL(UFS+n) = qL(QVELN) * uL(UFS+n);
     }
+#if NAUX_NET > 0
+    for (int n = 0; n < NumAux; n++) {
+      FL(UFX+n) = qL(QVELN) * uL(UFX+n);
+    }
+#endif
@@
     for (int n = 0; n < NumSpec; n++) {
       FR(UFS+n) = qR(QVELN) * uR(UFS+n);
     }
+#if NAUX_NET > 0
+    for (int n = 0; n < NumAux; n++) {
+      FR(UFX+n) = qR(QVELN) * uR(UFX+n);
+    }
+#endif
@@
     for (int n = 0; n < NumSpec; n++) {
       UsL(UFS+n) = qL(QFS+n) * UsL(URHO);
       UsR(UFS+n) = qR(QFS+n) * UsR(URHO);
     }
+#if NAUX_NET > 0
+    for (int n = 0; n < NumAux; n++) {
+      UsL(UFX+n) = qL(QFX+n) * UsL(URHO);
+      UsR(UFX+n) = qR(QFX+n) * UsR(URHO);
+    }
+#endif
@@
     for (int n = 0; n < NumSpec; n++) {
       UssL(UFS+n) = UsL(UFS+n);
       UssR(UFS+n) = UsR(UFS+n);
     }
+#if NAUX_NET > 0
+    for (int n = 0; n < NumAux; n++) {
+      UssL(UFX+n) = UsL(UFX+n);
+      UssR(UFX+n) = UsR(UFX+n);
+    }
+#endif

Prepared by Codex