change how we do MHD prim_half

Source/mhd/mhd_util.H

@@ -21,60 +21,6 @@ eos_soundspeed_mhd(Real& c, Real as, Real ca, Real bd) {
   c = std::sqrt(c);
 }
 
-AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-void
-qflux(Real* qflx, Real* flx, Real* q_zone) {
-
-  // Calculate the C to P Jacobian applied to the fluxes
-
-  // this is step 10 in the paper, just after Eq. 48
-
-  // this implements dW/dU . qflux, where dW/dU is the Jacobian of
-  // the primitive quantities (W) with respect to conserved quantities (U)
-
-  qflx[QRHO] = flx[URHO];
-  qflx[QU] = (flx[UMX] - flx[URHO] * q_zone[QU]) / q_zone[QRHO];
-  qflx[QV] = (flx[UMY] - flx[URHO] * q_zone[QV]) / q_zone[QRHO];
-  qflx[QW] = (flx[UMZ] - flx[URHO] * q_zone[QW]) / q_zone[QRHO];
-
-  for (int n = 0; n < NumSpec; n++) {
-    qflx[QFS+n] = (flx[UFS+n] - flx[URHO] * q_zone[QFS+n]) / q_zone[QRHO];
-  }
-
-  eos_t eos_state;
-  eos_state.rho = q_zone[QRHO];
-  eos_state.p = q_zone[QPRES];
-  eos_state.T = 100.0_rt; // dummy initial guess
-  for (int n = 0; n < NumSpec; n++) {
-    eos_state.xn[n] = q_zone[QFS+n];
-  }
-#if NAUX_NET > 0
-  for (int n = 0; n < NumAux; n++) {
-    eos_state.aux[n] = q_zone[QFX+n];
-  }
-#endif
-
-  eos(eos_input_rp, eos_state);
-
-  Real dedrho = eos_state.dedr - eos_state.dedT * eos_state.dpdr * 1.0_rt / eos_state.dpdT;
-  Real dedp = eos_state.dedT * 1.0_rt / eos_state.dpdT;
-
-  qflx[QPRES] = (-q_zone[QMAGX] * flx[UMAGX] - q_zone[QMAGY] * flx[UMAGY] -
-                  q_zone[QMAGZ] * flx[UMAGZ] + flx[UEDEN] -
-                 flx[UMX] * q_zone[QU] - flx[UMY] * q_zone[QV] - flx[UMZ] * q_zone[QW] +
-                 flx[URHO] * (0.5_rt * (q_zone[QU] * q_zone[QU] +
-                                        q_zone[QV] * q_zone[QV] +
-                                        q_zone[QW] * q_zone[QW]) -
-                              eos_state.e - q_zone[QRHO] * dedrho)) /
-    (dedp * q_zone[QRHO]);
-
-  qflx[QMAGX] = flx[UMAGX];
-  qflx[QMAGY] = flx[UMAGY];
-  qflx[QMAGZ] = flx[UMAGZ];
-
-  qflx[QTEMP] = 0.0_rt;
-}
-
 
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
 void

Source/mhd/mhd_util.cpp

@@ -124,30 +124,35 @@ Castro::prim_half(const Box& bx,
   [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k)
   {
 
-    Real divF[NUM_STATE+3];
-    Real divF_q[NQ];
-    Real q_zone[NQ];
+    Array1D<Real, 0, NQ-1> q_zone;
+    Array1D<Real, 0, NUM_STATE+2> U_zone;
+
+    // we come in with q, so we need to compute this to U
+    for (int n = 0; n < NQ; n++) {
+      q_zone(n) = q_arr(i,j,k,n);
+    }
+
+    Real gam1;
+    PToC(q_zone, U_zone, gam1);
+
+    // now do the conservative update of U
 
     for (int n = 0; n < NUM_STATE+3; n++) {
-      divF[n] = (flxx(i+1,j,k,n) - flxx(i,j,k,n)) / dx[0];
+      U_zone(n) -= 0.5 * dt * (flxx(i+1,j,k,n) - flxx(i,j,k,n)) / dx[0];
 #if AMREX_SPACEDIM >= 2
-      divF[n] += (flxy(i,j+1,k,n) - flxy(i,j,k,n)) / dx[1];
+      U_zone(n) -= 0.5 * dt * (flxy(i,j+1,k,n) - flxy(i,j,k,n)) / dx[1];
 #endif
 #if AMREX_SPACEDIM == 3
-      divF[n] += (flxz(i,j,k+1,n) - flxz(i,j,k,n)) / dx[2];
+      U_zone(n) -= 0.5 * dt * (flxz(i,j,k+1,n) - flxz(i,j,k,n)) / dx[2];
 #endif
     }
 
-    // that is a flux of conserved variables -- transform it to primitive
-    for (int n = 0; n < NQ; n++) {
-      q_zone[n] = q_arr(i,j,k,n);
-    }
-
-    qflux(divF_q, divF, q_zone);
+    // now convert this to q
+    ConsToPrim(q_zone.arr, U_zone.arr);
 
     // Right below eq. 48
     for (int n = 0; n < NQ; n++) {
-      q2D(i,j,k,n) = q_arr(i,j,k,n) - 0.5_rt * dt * divF_q[n];
+      q2D(i,j,k,n) = q_zone(n);
     }
   });
 }