mhd: prim_half() uses Step 6 intermediate fluxes instead of final Step 9 fluxes in unsplit CTU/CT update

[zingale]

Description

construct_ctu_mhd_source() appears to compute the half-time
cell-centered primitive state q2D with the wrong fluxes.

Summary

In the unsplit MHD update, q2D is built via:

prim_half(obx, q2D_arr, q_arr,
          flxx1D_arr, flxy1D_arr, flxz1D_arr, dt);

at mhd/Castro_mhd.cpp.

However, at that point flxx1D_arr, flxy1D_arr, and flxz1D_arr no
longer contain the original 1D predictor fluxes. They were overwritten
in Step 6 with the averaged 2D intermediate fluxes:

• mhd/Castro_mhd.cpp lines 481-496

After that, Steps 7-9 compute the final half-time interface states and
final Godunov fluxes into flux[0], flux[1], and flux[2]:

• mhd/Castro_mhd.cpp lines 545, 557, 569

prim_half() is a straight flux-divergence update:

• mhd/mhd_util.cpp lines 130-149

So q2D is currently being updated from the Step 6 intermediate
fluxes, not from the final Step 9 fluxes.

Why This Is A Bug

The final edge-centered electric field update uses q2D:

• mhd/Castro_mhd.cpp lines 588, 595, 602

That means the final CT electric fields are based on a half-time
cell-centered state built from an earlier flux stage than the one used
for the final conservative update. In an unsplit CTU/CT MHD scheme,
this is inconsistent and can introduce stage errors in the magnetic
update and final EMFs.

Expected Behavior

q2D should be constructed using the final half-time fluxes in
flux[0..2], not the Step 6 averaged intermediate fluxes in
flxx1D_arr, flxy1D_arr, and flxz1D_arr.

Suggested Fix

Change the prim_half() call to use the final flux arrays:

prim_half(obx, q2D_arr, q_arr,
          flxx_arr, flxy_arr, flxz_arr, dt);

instead of:

prim_half(obx, q2D_arr, q_arr,
          flxx1D_arr, flxy1D_arr, flxz1D_arr, dt);

Notes

This is a stage/differencing bug in the unsplit algorithm, not an
x/y/z permutation bug.