From 5a01e72a157c6870afab6a146b9feaa1b7c62c05 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sat, 4 Jul 2020 21:01:39 -0400 Subject: [PATCH] change how we do MHD prim_half to use our existing ConsToPrim routine --- Source/mhd/mhd_util.H | 52 ----------------------------------------- Source/mhd/mhd_util.cpp | 31 +++++++++++++----------- 2 files changed, 18 insertions(+), 65 deletions(-) diff --git a/Source/mhd/mhd_util.H b/Source/mhd/mhd_util.H index 65f4764cc8..742bee526a 100644 --- a/Source/mhd/mhd_util.H +++ b/Source/mhd/mhd_util.H @@ -21,58 +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]; - } - for (int n = 0; n < NumAux; n++) { - eos_state.aux[n] = q_zone[QFX+n]; - } - - 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 diff --git a/Source/mhd/mhd_util.cpp b/Source/mhd/mhd_util.cpp index ad9b787f4b..1ded06ff02 100644 --- a/Source/mhd/mhd_util.cpp +++ b/Source/mhd/mhd_util.cpp @@ -236,30 +236,35 @@ Castro::prim_half(const Box& bx, [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) noexcept { - Real divF[NUM_STATE+3]; - Real divF_q[NQ]; - Real q_zone[NQ]; + Array1D q_zone; + Array1D 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); } }); }