From 3a5d173f7da5b6e5ab80bf2216a7ddd309626af7 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Thu, 24 Jun 2021 21:28:16 -0400 Subject: [PATCH 1/5] work in progress --- Source/hydro/trans.cpp | 296 +++++++++++------------------------------ 1 file changed, 77 insertions(+), 219 deletions(-) diff --git a/Source/hydro/trans.cpp b/Source/hydro/trans.cpp index 26b587a195..c3b02571bc 100644 --- a/Source/hydro/trans.cpp +++ b/Source/hydro/trans.cpp @@ -103,13 +103,6 @@ Castro::actual_trans_single(const Box& bx, bool reset_rhoe = transverse_reset_rhoe; Real small_p = small_pres; -#ifdef RADIATION - int fspace_t = Radiation::fspace_advection_type; - int comov = Radiation::comoving; - int limiter = Radiation::limiter; - int closure = Radiation::closure; -#endif - amrex::ParallelFor(bx, [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) { @@ -161,135 +154,63 @@ Castro::actual_trans_single(const Box& bx, kr += d; } - // Update all of the passively-advected quantities with the - // transverse term and convert back to the primitive quantity. + // construct the conserved variables + + Real U_int[NUM_STATE]; + + U_int[URHO] = q_arr(i,j,k,QRHO); + U_int[UMX] = U_int[URHO] * q_arr(i,j,k,QU); + U_int[UMY] = U_int[URHO] * q_arr(i,j,k,QV); + U_int[UMZ] = U_int[URHO] * q_arr(i,j,k,QW); + U_int[UEDEN] = U_int[URHO] * q_arr(i,j,k,QREINT) + + 0.5_rt * (U_int[UMX] * U_int[UMX] + + U_int[UMY] * U_int[UMY] + + U_int[UMZ] * U_int[UMZ]) / U_int[URHO]; + U_int[UEINT] = U_int[URHO] * q_arr(i,j,k,QREINT); -#if AMREX_SPACEDIM == 2 - const Real volinv = 1.0_rt / vol(il,jl,kl); -#endif for (int ipassive = 0; ipassive < npassive; ipassive++) { int n = upassmap(ipassive); int nqp = qpassmap(ipassive); -#if AMREX_SPACEDIM == 2 - Real rrnew = q_arr(i,j,k,QRHO) - hdt * (area_t(ir,jr,kr) * flux_t(ir,jr,kr,URHO) - - area_t(il,jl,kl) * flux_t(il,jl,kl,URHO)) * volinv; - Real compu = q_arr(i,j,k,QRHO) * q_arr(i,j,k,nqp) - hdt * (area_t(ir,jr,kr) * flux_t(ir,jr,kr,n) - - area_t(il,jl,kl) * flux_t(il,jl,kl,n)) * volinv; - qo_arr(i,j,k,nqp) = compu / rrnew; -#else - Real rrnew = q_arr(i,j,k,QRHO) - cdtdx * (flux_t(ir,jr,kr,URHO) - flux_t(il,jl,kl,URHO)); - Real compu = q_arr(i,j,k,QRHO) * q_arr(i,j,k,nqp) - cdtdx * (flux_t(ir,jr,kr,n) - flux_t(il,jl,kl,n)); - qo_arr(i,j,k,nqp) = compu / rrnew; -#endif + U_int[n] = U_int[URHO] * q_arr(i,j,k,nqp); } - Real pgp = q_t(ir,jr,kr,GDPRES); - Real pgm = q_t(il,jl,kl,GDPRES); - Real ugp = q_t(ir,jr,kr,GDU+idir_t); - Real ugm = q_t(il,jl,kl,GDU+idir_t); + // TODO: how does the hybrid stuff fit in here? -#ifdef RADIATION - Real lambda[NGROUPS]; - Real ergp[NGROUPS]; - Real ergm[NGROUPS]; + // now do the transverse flux update - for (int g = 0; g < NGROUPS; ++g) { - lambda[g] = qaux_arr(il,jl,kl,QLAMS+g); - ergp[g] = q_t(ir,jr,kr,GDERADS+g); - ergm[g] = q_t(il,jl,kl,GDERADS+g); - } +#if AMREX_SPACEDIM == 2 + const Real volinv = 1.0_rt / vol(il,jl,kl); #endif - // we need to augment our conserved system with either a p - // equation or gammae (if we have ppm_predict_gammae = 1) to - // be able to deal with the general EOS + for (int n = 0; n < NUM_STATE; n++) { #if AMREX_SPACEDIM == 2 - Real dup = area_t(ir,jr,kr) * pgp * ugp - area_t(il,jl,kl) * pgm * ugm; - Real du = area_t(ir,jr,kr) * ugp - area_t(il,jl,kl) * ugm; + U_int[n] += - hdt * (area_t(ir,jr,kr) * flux_t(ir,jr,kr,URHO) - + area_t(il,jl,kl) * flux_t(il,jl,kl,URHO)) * volinv; #else - Real dup = pgp * ugp - pgm * ugm; - Real du = ugp - ugm; -#endif - Real pav = 0.5_rt * (pgp + pgm); -#ifdef RADIATION - Real uav = 0.5_rt * (ugp + ugm); -#endif - - // this is the gas gamma_1 -#ifdef RADIATION - Real gamc = qaux_arr(il,jl,kl,QGAMCG); -#else - Real gamc = qaux_arr(il,jl,kl,QGAMC); + U_int[n] += - cdtdx * (flux_t(ir,jr,kr,n) - flux_t(il,jl,kl,n)); #endif -#ifdef RADIATION - Real lamge[NGROUPS]; - Real luge[NGROUPS]; - Real der[NGROUPS]; + } - Real dmom = 0.0_rt; - Real dre = 0.0_rt; - for (int g = 0; g < NGROUPS; ++g) { - lamge[g] = lambda[g] * (ergp[g] - ergm[g]); - dmom += -cdtdx * lamge[g]; - luge[g] = uav * lamge[g]; - dre += -cdtdx * luge[g]; - } + Real pgp = q_t(ir,jr,kr,GDPRES); + Real pgm = q_t(il,jl,kl,GDPRES); + Real ugp = q_t(ir,jr,kr,GDU+idir_t); + Real ugm = q_t(il,jl,kl,GDU+idir_t); - if (fspace_t == 1 && comov) { - for (int g = 0; g < NGROUPS; ++g) { - Real eddf = Edd_factor(lambda[g], limiter, closure); - Real f1 = 0.5_rt * (1.0_rt - eddf); - der[g] = cdtdx * uav * f1 * (ergp[g] - ergm[g]); - } - } - else if (fspace_t == 2) { #if AMREX_SPACEDIM == 2 - Real divu = (area_t(ir,jr,kr) * ugp - area_t(il,jl,kl) * ugm) * volinv; - for (int g = 0; g < NGROUPS; g++) { - Real eddf = Edd_factor(lambda[g], limiter, closure); - Real f1 = 0.5_rt * (1.0_rt - eddf); - der[g] = -hdt * f1 * 0.5_rt * (ergp[g] + ergm[g]) * divu; - } + Real du = area_t(ir,jr,kr) * ugp - area_t(il,jl,kl) * ugm; #else - for (int g = 0; g < NGROUPS; g++) { - Real eddf = Edd_factor(lambda[g], limiter, closure); - Real f1 = 0.5_rt * (1.0_rt - eddf); - der[g] = cdtdx * f1 * 0.5_rt * (ergp[g] + ergm[g]) * (ugm - ugp); - } -#endif - } - else { // mixed frame - for (int g = 0; g < NGROUPS; g++) { - der[g] = cdtdx * luge[g]; - } - } -#endif - - // Convert to conservation form - Real rrn = q_arr(i,j,k,QRHO); - Real run = rrn * q_arr(i,j,k,QU); - Real rvn = rrn * q_arr(i,j,k,QV); - Real rwn = rrn * q_arr(i,j,k,QW); - Real ekenn = 0.5_rt * rrn * (q_arr(i,j,k,QU) * q_arr(i,j,k,QU) + q_arr(i,j,k,QV) * q_arr(i,j,k,QV) + q_arr(i,j,k,QW) * q_arr(i,j,k,QW)); - Real ren = q_arr(i,j,k,QREINT) + ekenn; -#ifdef RADIATION - Real ern[NGROUPS]; - for (int g = 0; g < NGROUPS; ++g) { - ern[g] = q_arr(i,j,k,QRAD+g); - } + Real du = ugp - ugm; #endif + Real pav = 0.5_rt * (pgp + pgm); #if AMREX_SPACEDIM == 2 - // Add transverse predictor - Real rrnewn = rrn - hdt * (area_t(ir,jr,kr) * flux_t(ir,jr,kr,URHO) - - area_t(il,jl,kl) * flux_t(il,jl,kl,URHO)) * volinv; - // Note that pressure may be treated specially here, depending on - // the geometry. Our y-interface equation for (rho u) is: + // For the x-momentum, we need to consider the grad p term. + // Our y-interface equation for (rho u) is: // // d(rho u)/dt + d(rho u v)/dy = - 1/r d(r rho u u)/dr - dp/dr // @@ -298,140 +219,77 @@ Castro::actual_trans_single(const Box& bx, // are no area factors. For this geometry, we do not // include p in our definition of the flux in the // x-direction, for we need to fix this now. - Real runewn = run - hdt * (area_t(ir,jr,kr) * flux_t(ir,jr,kr,UMX) - - area_t(il,jl,kl) * flux_t(il,jl,kl,UMX)) * volinv; - if (idir_t == 0 && !mom_flux_has_p(0, idir_t, coord)) { - runewn = runewn - cdtdx * (pgp - pgm); - } - Real rvnewn = rvn - hdt * (area_t(ir,jr,kr) * flux_t(ir,jr,kr,UMY) - - area_t(il,jl,kl) * flux_t(il,jl,kl,UMY)) * volinv; - Real rwnewn = rwn - hdt * (area_t(ir,jr,kr) * flux_t(ir,jr,kr,UMZ) - - area_t(il,jl,kl) * flux_t(il,jl,kl,UMZ)) * volinv; - Real renewn = ren - hdt * (area_t(ir,jr,kr) * flux_t(ir,jr,kr,UEDEN) - - area_t(il,jl,kl) * flux_t(il,jl,kl,UEDEN)) * volinv; - -#ifdef RADIATION - if (idir_t == 0) { - Real lamge_sum = 0.0_rt; - for (int g = 0; g < NGROUPS; ++g) - lamge_sum = lamge_sum + lamge[g]; - - runewn = runewn - 0.5_rt * hdt * (area_t(ir,jr,kr) + area_t(il,jl,kl)) * lamge_sum * volinv; - } - else { - rvnewn = rvnewn + dmom; - } - renewn = renewn + dre; - - Real ernewn[NGROUPS]; - for (int g = 0; g < NGROUPS; ++g) { - ernewn[g] = ern[g] - hdt * (area_t(ir,jr,kr) * rflux_t(ir,jr,kr,g) - - area_t(il,jl,kl) * rflux_t(il,jl,kl,g)) * volinv + der[g]; + if (idir_t == 0 && !mom_flux_has_p(0, idir_t, coord)) { + U_int[UMX] += - cdtdx * (pgp - pgm); } #endif + // For the dual energy approach, we need to correct the + // internal energy equation with the p dV term + +#if AMREX_SPACEDIM == 2 + U_int[UEINT] += - hdt * pav * du * volinv; #else - // Add transverse predictor - Real rrnewn = rrn - cdtdx * (flux_t(ir,jr,kr,URHO) - flux_t(il,jl,kl,URHO)); - Real runewn = run - cdtdx * (flux_t(ir,jr,kr,UMX) - flux_t(il,jl,kl,UMX)); - Real rvnewn = rvn - cdtdx * (flux_t(ir,jr,kr,UMY) - flux_t(il,jl,kl,UMY)); - Real rwnewn = rwn - cdtdx * (flux_t(ir,jr,kr,UMZ) - flux_t(il,jl,kl,UMZ)); - Real renewn = ren - cdtdx * (flux_t(ir,jr,kr,UEDEN) - flux_t(il,jl,kl,UEDEN)); -#ifdef RADIATION - runewn = runewn + dmom; - renewn = renewn + dre; - Real ernewn[NGROUPS]; - for (int g = 0; g < NGROUPS; ++g) { - ernewn[g] = ern[g] - cdtdx * (rflux_t(ir,jr,kr,g) - rflux_t(il,jl,kl,g)) + der[g]; - } -#endif + U_int[UEINT] += - cdtdx * pav * du; #endif // Reset to original value if adding transverse terms made density negative + bool reset_state = false; if (reset_density == 1 && rrnewn < 0.0_rt) { - rrnewn = rrn; - runewn = run; - rvnewn = rvn; - rwnewn = rwn; - renewn = ren; -#ifdef RADIATION - for (int g = 0; g < NGROUPS; ++g) { - ernewn[g] = ern[g]; - } -#endif reset_state = true; } - // Convert back to primitive form - qo_arr(i,j,k,QRHO) = rrnewn; - Real rhoinv = 1.0_rt / rrnewn; - qo_arr(i,j,k,QU) = runewn * rhoinv; - qo_arr(i,j,k,QV) = rvnewn * rhoinv; - qo_arr(i,j,k,QW) = rwnewn * rhoinv; + if (! reset_state) { - // note: we run the risk of (rho e) being negative here - Real rhoekenn = 0.5_rt * (runewn * runewn + rvnewn * rvnewn + rwnewn * rwnewn) * rhoinv; - qo_arr(i,j,k,QREINT) = renewn - rhoekenn; + // Convert back to primitive form - if (!reset_state) { - // do the transverse terms for p, gamma, and rhoe, as necessary + qo_arr(i,j,k,QRHO) = U_int[URHO]; + Real rhoinv = 1.0_rt / U_int[URHO]; - if (reset_rhoe == 1 && qo_arr(i,j,k,QREINT) <= 0.0_rt) { - // If it is negative, reset the internal energy by - // using the discretized expression for updating (rho e). -#if AMREX_SPACEDIM == 2 - qo_arr(i,j,k,QREINT) = q_arr(i,j,k,QREINT) - hdt * (area_t(ir,jr,kr) * flux_t(ir,jr,kr,UEINT) - - area_t(il,jl,kl) * flux_t(il,jl,kl,UEINT) + - pav * du) * volinv; -#else - qo_arr(i,j,k,QREINT) = q_arr(i,j,k,QREINT) - cdtdx * (flux_t(ir,jr,kr,UEINT) - flux_t(il,jl,kl,UEINT) + pav * du); -#endif - } + qo_arr(i,j,k,QU) = U_int[UMX] * rhoinv; + qo_arr(i,j,k,QV) = U_int[UMY] * rhoinv; + qo_arr(i,j,k,QW) = U_int[UMZ] * rhoinv; - // If (rho e) is negative by this point, - // set it back to the original interface state, - // which turns off the transverse correction. + Real kineng = 0.5_rt * qo_arr(i,j,k,QRHO) * + (qo_arr(i,j,k,QU) * qo_arr(i,j,k,QU) + + qo_arr(i,j,k,QV) * qo_arr(i,j,k,QV) + + qo_arr(i,j,k,QW) * qo_arr(i,j,k,QW)); - if (qo_arr(i,j,k,QREINT) <= 0.0_rt) { - qo_arr(i,j,k,QREINT) = q_arr(i,j,k,QREINT); + if ((U_int[UEDEN] - kineng) > castro::dual_energy_eta1 * U_int[UEDEN]) { + qo_arr(i,j,k,QREINT) = (U_int[UEDEN] - kineng) * rhoinv; + } else { + qo_arr(i,j,k,QREINT) = U_int[UEINT] * rhoinv; } - // Pretend QREINT has been fixed and transverse_use_eos != 1. - // If we are wrong, we will fix it later. + // Load passively advected quatities into q + for (int ipassive = 0; ipassive < npassive; ipassive++) { + int n = upassmap(ipassive); + int iq = qpassmap(ipassive); + qo_arr(i,j,k,iq) = U_int[n] * rhoinv; + } - // Add the transverse term to the p evolution eq here. -#if AMREX_SPACEDIM == 2 - // the divergences here, dup and du, already have area factors - Real pnewn = q_arr(i,j,k,QPRES) - hdt * (dup + pav * du * (gamc - 1.0_rt)) * volinv; -#else - Real pnewn = q_arr(i,j,k,QPRES) - cdtdx * (dup + pav * du * (gamc - 1.0_rt)); -#endif - qo_arr(i,j,k,QPRES) = amrex::max(pnewn, small_p); + eos_rep_t eos_state; - } - else { - qo_arr(i,j,k,QPRES) = q_arr(i,j,k,QPRES); - qo_arr(i,j,k,QREINT) = q_arr(i,j,k,QREINT); - } + eos_state.T = qo_arr(i,j,k,QTEMP); + eos_state.rho = qo_arr(i,j,k,QRHO); + eos_state.e = qo_arr(i,j,k,QREINT); + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = qo_arr(i,j,k,QFS+n); + } +#if NAUX_NET > 0 + for (int n = 0; n < NumAux; n++) { + eos_state.aux[n] = qo_arr(i,j,k,QFX+n); + } +#endif -#ifdef RADIATION - for (int g = 0; g < NGROUPS; ++g) { - qo_arr(i,j,k,QRAD + g) = ernewn[g]; - } + eos(eos_input_re, eos_state); - qo_arr(i,j,k,QPTOT) = qo_arr(i,j,k,QPRES); - for (int g = 0; g < NGROUPS; ++g) { - qo_arr(i,j,k,QPTOT) += lambda[g] * ernewn[g]; - } + qo_arr(i,j,k,QTEMP) = eos_state.T; + qo_arr(i,j,k,QPRES) = eos_state.p; - qo_arr(i,j,k,QREITOT) = qo_arr(i,j,k,QREINT); - for (int g = 0; g < NGROUPS; ++g) { - qo_arr(i,j,k,QREITOT) += qo_arr(i,j,k,QRAD + g); } -#endif - }); } From c3ee56cbaa4eb1e98a13880e4fb5b237a52488cb Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 25 Jun 2021 10:30:07 -0400 Subject: [PATCH 2/5] first cut --- Source/hydro/Castro_ctu_hydro.cpp | 44 ++--- Source/hydro/trans.cpp | 316 ++++++++++-------------------- 2 files changed, 130 insertions(+), 230 deletions(-) diff --git a/Source/hydro/Castro_ctu_hydro.cpp b/Source/hydro/Castro_ctu_hydro.cpp index 157f9187e7..9aa816598b 100644 --- a/Source/hydro/Castro_ctu_hydro.cpp +++ b/Source/hydro/Castro_ctu_hydro.cpp @@ -629,9 +629,9 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) vol_arr, hdt, hdtdy); - reset_edge_state_thermo(xbx, ql.array()); + //reset_edge_state_thermo(xbx, ql.array()); - reset_edge_state_thermo(xbx, qr.array()); + //reset_edge_state_thermo(xbx, qr.array()); // solve the final Riemann problem axross the x-interfaces @@ -671,9 +671,9 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) vol_arr, hdt, hdtdx); - reset_edge_state_thermo(ybx, ql.array()); + //reset_edge_state_thermo(ybx, ql.array()); - reset_edge_state_thermo(ybx, qr.array()); + //reset_edge_state_thermo(ybx, qr.array()); // solve the final Riemann problem axross the y-interfaces @@ -754,9 +754,9 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) qgdnvtmp1_arr, hdt, cdtdx); - reset_edge_state_thermo(tyxbx, qmyx.array()); + //reset_edge_state_thermo(tyxbx, qmyx.array()); - reset_edge_state_thermo(tyxbx, qpyx.array()); + //reset_edge_state_thermo(tyxbx, qpyx.array()); // [lo(1), lo(2)-1, lo(3)], [hi(1), hi(2)+1, hi(3)+1] const Box& tzxbx = amrex::grow(zbx, IntVect(AMREX_D_DECL(0,1,0))); @@ -782,9 +782,9 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) qgdnvtmp1_arr, hdt, cdtdx); - reset_edge_state_thermo(tzxbx, qmzx.array()); + //reset_edge_state_thermo(tzxbx, qmzx.array()); - reset_edge_state_thermo(tzxbx, qpzx.array()); + //reset_edge_state_thermo(tzxbx, qpzx.array()); // compute F^y // [lo(1)-1, lo(2), lo(3)-1], [hi(1)+1, hi(2)+1, hi(3)+1] @@ -830,9 +830,9 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) qgdnvtmp1_arr, hdt, cdtdy); - reset_edge_state_thermo(txybx, qmxy.array()); + //reset_edge_state_thermo(txybx, qmxy.array()); - reset_edge_state_thermo(txybx, qpxy.array()); + //reset_edge_state_thermo(txybx, qpxy.array()); // [lo(1)-1, lo(2), lo(3)], [hi(1)+1, hi(2), lo(3)+1] const Box& tzybx = amrex::grow(zbx, IntVect(AMREX_D_DECL(1,0,0))); @@ -861,9 +861,9 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) qgdnvtmp1_arr, hdt, cdtdy); - reset_edge_state_thermo(tzybx, qmzy.array()); + //reset_edge_state_thermo(tzybx, qmzy.array()); - reset_edge_state_thermo(tzybx, qpzy.array()); + //reset_edge_state_thermo(tzybx, qpzy.array()); // compute F^z // [lo(1)-1, lo(2)-1, lo(3)], [hi(1)+1, hi(2)+1, hi(3)+1] @@ -909,9 +909,9 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) qgdnvtmp1_arr, hdt, cdtdz); - reset_edge_state_thermo(txzbx, qmxz.array()); + //reset_edge_state_thermo(txzbx, qmxz.array()); - reset_edge_state_thermo(txzbx, qpxz.array()); + //reset_edge_state_thermo(txzbx, qpxz.array()); // [lo(1)-1, lo(2), lo(3)], [hi(1)+1, hi(2)+1, lo(3)] const Box& tyzbx = amrex::grow(ybx, IntVect(AMREX_D_DECL(1,0,0))); @@ -940,9 +940,9 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) qgdnvtmp1_arr, hdt, cdtdz); - reset_edge_state_thermo(tyzbx, qmyz.array()); + //reset_edge_state_thermo(tyzbx, qmyz.array()); - reset_edge_state_thermo(tyzbx, qpyz.array()); + //reset_edge_state_thermo(tyzbx, qpyz.array()); // we now have q?zx, q?yx, q?zy, q?xy, q?yz, q?xz @@ -1003,9 +1003,9 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) qgdnvtmp2_arr, hdtdy, hdtdz); - reset_edge_state_thermo(xbx, ql.array()); + //reset_edge_state_thermo(xbx, ql.array()); - reset_edge_state_thermo(xbx, qr.array()); + //reset_edge_state_thermo(xbx, qr.array()); #ifdef PRIM_SPECIES_HAVE_SOURCES add_species_source_to_states(xbx, 0, dt, @@ -1081,9 +1081,9 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) qgdnvtmp1_arr, hdtdx, hdtdz); - reset_edge_state_thermo(ybx, ql.array()); + //reset_edge_state_thermo(ybx, ql.array()); - reset_edge_state_thermo(ybx, qr.array()); + //reset_edge_state_thermo(ybx, qr.array()); #ifdef PRIM_SPECIES_HAVE_SOURCES add_species_source_to_states(ybx, 1, dt, @@ -1161,9 +1161,9 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) qgdnvtmp2_arr, hdtdx, hdtdy); - reset_edge_state_thermo(zbx, ql.array()); + //reset_edge_state_thermo(zbx, ql.array()); - reset_edge_state_thermo(zbx, qr.array()); + //reset_edge_state_thermo(zbx, qr.array()); #ifdef PRIM_SPECIES_HAVE_SOURCES add_species_source_to_states(zbx, 2, dt, diff --git a/Source/hydro/trans.cpp b/Source/hydro/trans.cpp index c3b02571bc..b9a72ed466 100644 --- a/Source/hydro/trans.cpp +++ b/Source/hydro/trans.cpp @@ -185,6 +185,16 @@ Castro::actual_trans_single(const Box& bx, for (int n = 0; n < NUM_STATE; n++) { + if (n == UTEMP) { + continue; + } + +#ifdef SHOCK_VAR + if (n == USHK) { + continue; + } +#endif + #if AMREX_SPACEDIM == 2 U_int[n] += - hdt * (area_t(ir,jr,kr) * flux_t(ir,jr,kr,URHO) - area_t(il,jl,kl) * flux_t(il,jl,kl,URHO)) * volinv; @@ -234,10 +244,11 @@ Castro::actual_trans_single(const Box& bx, U_int[UEINT] += - cdtdx * pav * du; #endif - // Reset to original value if adding transverse terms made density negative + // Reset to original value if adding transverse terms made + // density negative bool reset_state = false; - if (reset_density == 1 && rrnewn < 0.0_rt) { + if (reset_density == 1 && U_int[URHO] < 0.0_rt) { reset_state = true; } @@ -289,6 +300,12 @@ Castro::actual_trans_single(const Box& bx, qo_arr(i,j,k,QTEMP) = eos_state.T; qo_arr(i,j,k,QPRES) = eos_state.p; + } else { + + for (int n = 0; n < NQ; n++) { + qo_arr(i,j,k,n) = q_arr(i,j,k,n); + } + } }); @@ -465,25 +482,48 @@ Castro::actual_trans_final(const Box& bx, } - // Update all of the passively-advected quantities with the - // transverse terms and convert back to the primitive quantity. + // construct the conserved variables + + Real U_int[NUM_STATE]; + + U_int[URHO] = q_arr(i,j,k,QRHO); + U_int[UMX] = U_int[URHO] * q_arr(i,j,k,QU); + U_int[UMY] = U_int[URHO] * q_arr(i,j,k,QV); + U_int[UMZ] = U_int[URHO] * q_arr(i,j,k,QW); + U_int[UEDEN] = U_int[URHO] * q_arr(i,j,k,QREINT) + + 0.5_rt * (U_int[UMX] * U_int[UMX] + + U_int[UMY] * U_int[UMY] + + U_int[UMZ] * U_int[UMZ]) / U_int[URHO]; + U_int[UEINT] = U_int[URHO] * q_arr(i,j,k,QREINT); - for (int ipassive = 0; ipassive < npassive; ++ipassive) { + for (int ipassive = 0; ipassive < npassive; ipassive++) { int n = upassmap(ipassive); int nqp = qpassmap(ipassive); - Real rrn = q_arr(i,j,k,QRHO); - Real compn = rrn * q_arr(i,j,k,nqp); - Real rrnewn = rrn - cdtdx_t1 * (flux_t1(ir_t1,jr_t1,kr_t1,URHO) - - flux_t1(il_t1,jl_t1,kl_t1,URHO)) - - cdtdx_t2 * (flux_t2(ir_t2,jr_t2,kr_t2,URHO) - - flux_t2(il_t2,jl_t2,kl_t2,URHO)); - Real compnn = compn - cdtdx_t1 * (flux_t1(ir_t1,jr_t1,kr_t1,n) - - flux_t1(il_t1,jl_t1,kl_t1,n)) - - cdtdx_t2 * (flux_t2(ir_t2,jr_t2,kr_t2,n) - - flux_t2(il_t2,jl_t2,kl_t2,n)); - - qo_arr(i,j,k,nqp) = compnn / rrnewn; + U_int[n] = U_int[URHO] * q_arr(i,j,k,nqp); + } + + // TODO: how does the hybrid stuff fit in here? + + // now do the transverse flux update + + for (int n = 0; n < NUM_STATE; n++) { + + if (n == UTEMP) { + continue; + } + +#ifdef SHOCK_VAR + if (n == USHK) { + continue; + } +#endif + + U_int[n] += - cdtdx_t1 * (flux_t1(ir_t1,jr_t1,kr_t1,n) - + flux_t1(il_t1,jl_t1,kl_t1,n)) + - cdtdx_t2 * (flux_t2(ir_t2,jr_t2,kr_t2,n) - + flux_t2(il_t2,jl_t2,kl_t2,n)); + } // Add the transverse differences to the normal states for the @@ -493,227 +533,87 @@ Castro::actual_trans_final(const Box& bx, Real pgt1m = q_t1(il_t1,jl_t1,kl_t1,GDPRES); Real ugt1p = q_t1(ir_t1,jr_t1,kr_t1,GDU+idir_t1); Real ugt1m = q_t1(il_t1,jl_t1,kl_t1,GDU+idir_t1); -#ifdef RADIATION - Real ergt1p[NGROUPS]; - Real ergt1m[NGROUPS]; - for (int g = 0; g < NGROUPS; ++g) { - ergt1p[g] = q_t1(ir_t1,jr_t1,kr_t1,GDERADS+g); - ergt1m[g] = q_t1(il_t1,jl_t1,kl_t1,GDERADS+g); - } -#endif Real pgt2p = q_t2(ir_t2,jr_t2,kr_t2,GDPRES); Real pgt2m = q_t2(il_t2,jl_t2,kl_t2,GDPRES); Real ugt2p = q_t2(ir_t2,jr_t2,kr_t2,GDU+idir_t2); Real ugt2m = q_t2(il_t2,jl_t2,kl_t2,GDU+idir_t2); -#ifdef RADIATION - Real ergt2p[NGROUPS]; - Real ergt2m[NGROUPS]; - for (int g = 0; g < NGROUPS; ++g) { - ergt2p[g] = q_t2(ir_t2,jr_t2,kr_t2,GDERADS+g); - ergt2m[g] = q_t2(il_t2,jl_t2,kl_t2,GDERADS+g); - } -#endif - Real dupt1 = pgt1p * ugt1p - pgt1m * ugt1m; Real pt1av = 0.5_rt * (pgt1p + pgt1m); Real dut1 = ugt1p - ugt1m; -#ifdef RADIATION - Real pt1new = cdtdx_t1 * (dupt1 + pt1av * dut1 * (qaux_arr(iln,jln,kln,QGAMCG) - 1.0_rt)); -#else - Real pt1new = cdtdx_t1 * (dupt1 + pt1av * dut1 * (qaux_arr(iln,jln,kln,QGAMC) - 1.0_rt)); -#endif - Real dupt2 = pgt2p * ugt2p - pgt2m * ugt2m; Real pt2av = 0.5_rt * (pgt2p + pgt2m); Real dut2 = ugt2p - ugt2m; -#ifdef RADIATION - Real pt2new = cdtdx_t2 * (dupt2 + pt2av * dut2 * (qaux_arr(iln,jln,kln,QGAMCG) - 1.0_rt)); -#else - Real pt2new = cdtdx_t2 * (dupt2 + pt2av * dut2 * (qaux_arr(iln,jln,kln,QGAMC) - 1.0_rt)); -#endif -#ifdef RADIATION - Real lambda[NGROUPS]; - Real lget1[NGROUPS]; - Real lget2[NGROUPS]; - Real dmt1 = 0.0_rt; - Real dmt2 = 0.0_rt; - Real luget1[NGROUPS]; - Real luget2[NGROUPS]; - Real dre = 0.0_rt; - - for (int g = 0; g < NGROUPS; ++g) { - lambda[g] = qaux_arr(iln,jln,kln,QLAMS+g); - lget1[g] = lambda[g] * (ergt1p[g] - ergt1m[g]); - lget2[g] = lambda[g] * (ergt2p[g] - ergt2m[g]); - dmt1 -= cdtdx_t1 * lget1[g]; - dmt2 -= cdtdx_t2 * lget2[g]; - luget1[g] = 0.5_rt * (ugt1p + ugt1m) * lget1[g]; - luget2[g] = 0.5_rt * (ugt2p + ugt2m) * lget2[g]; - dre -= cdtdx_t1 * luget1[g] + cdtdx_t2 * luget2[g]; - } + // For the dual energy approach, we need to correct the + // internal energy equation with the p dV term - Real der[NGROUPS]; + U_int[UEINT] += - cdtdx_t1 * pt1av * dut1 - cdtdx_t2 * pt2av * dut2; - if (fspace_t == 1 && comov) { - for (int g = 0; g < NGROUPS; ++g) { - Real eddf = Edd_factor(lambda[g], limiter, closure); - Real f1 = 0.5_rt * (1.0_rt - eddf); - der[g] = f1 * (cdtdx_t1 * 0.5_rt * (ugt1p + ugt1m) * (ergt1p[g] - ergt1m[g]) + - cdtdx_t2 * 0.5_rt * (ugt2p + ugt2m) * (ergt2p[g] - ergt2m[g])); - } - } - else if (fspace_t == 2) { - for (int g = 0; g < NGROUPS; ++g) { - Real eddf = Edd_factor(lambda[g], limiter, closure); - Real f1 = 0.5_rt * (1.0_rt - eddf); - der[g] = f1 * (cdtdx_t1 * 0.5_rt * (ergt1p[g] + ergt1m[g]) * (ugt1m - ugt1p) + - cdtdx_t2 * 0.5_rt * (ergt2p[g] + ergt2m[g]) * (ugt2m - ugt2p)); - } - } - else { // mixed frame - for (int g = 0; g < NGROUPS; ++g) { - der[g] = cdtdx_t1 * luget1[g] + cdtdx_t2 * luget2[g]; - } - } -#endif + // Reset to original value if adding transverse terms made + // density negative - // Convert to conservation form - Real rrn = q_arr(i,j,k,QRHO); - Real run = rrn * q_arr(i,j,k,QU); - Real rvn = rrn * q_arr(i,j,k,QV); - Real rwn = rrn * q_arr(i,j,k,QW); - Real ekenn = 0.5_rt * rrn * (q_arr(i,j,k,QU) * q_arr(i,j,k,QU) + q_arr(i,j,k,QV) * q_arr(i,j,k,QV) + q_arr(i,j,k,QW) * q_arr(i,j,k,QW)); - Real ren = q_arr(i,j,k,QREINT) + ekenn; -#ifdef RADIATION - Real ern[NGROUPS]; - for (int g = 0; g < NGROUPS; ++g) { - ern[g] = q_arr(i,j,k,QRAD+g); + bool reset_state = false; + if (reset_density == 1 && U_int[URHO] < 0.0_rt) { + reset_state = true; } -#endif - // Add transverse predictor - Real rrnewn = rrn - cdtdx_t1 * (flux_t1(ir_t1,jr_t1,kr_t1,URHO) - - flux_t1(il_t1,jl_t1,kl_t1,URHO)) - - cdtdx_t2 * (flux_t2(ir_t2,jr_t2,kr_t2,URHO) - - flux_t2(il_t2,jl_t2,kl_t2,URHO)); - Real runewn = run - cdtdx_t1 * (flux_t1(ir_t1,jr_t1,kr_t1,UMX) - - flux_t1(il_t1,jl_t1,kl_t1,UMX)) - - cdtdx_t2 * (flux_t2(ir_t2,jr_t2,kr_t2,UMX) - - flux_t2(il_t2,jl_t2,kl_t2,UMX)); - Real rvnewn = rvn - cdtdx_t1 * (flux_t1(ir_t1,jr_t1,kr_t1,UMY) - - flux_t1(il_t1,jl_t1,kl_t1,UMY)) - - cdtdx_t2 * (flux_t2(ir_t2,jr_t2,kr_t2,UMY) - - flux_t2(il_t2,jl_t2,kl_t2,UMY)); - Real rwnewn = rwn - cdtdx_t1 * (flux_t1(ir_t1,jr_t1,kr_t1,UMZ) - - flux_t1(il_t1,jl_t1,kl_t1,UMZ)) - - cdtdx_t2 * (flux_t2(ir_t2,jr_t2,kr_t2,UMZ) - - flux_t2(il_t2,jl_t2,kl_t2,UMZ)); - Real renewn = ren - cdtdx_t1 * (flux_t1(ir_t1,jr_t1,kr_t1,UEDEN) - - flux_t1(il_t1,jl_t1,kl_t1,UEDEN)) - - cdtdx_t2 * (flux_t2(ir_t2,jr_t2,kr_t2,UEDEN) - - flux_t2(il_t2,jl_t2,kl_t2,UEDEN)); -#ifdef RADIATION - if (idir_n == 0) { - rvnewn = rvnewn + dmt1; - rwnewn = rwnewn + dmt2; - } - else if (idir_n == 1) { - runewn = runewn + dmt1; - rwnewn = rwnewn + dmt2; - } - else { - runewn = runewn + dmt1; - rvnewn = rvnewn + dmt2; - } - renewn = renewn + dre; - - Real ernewn[NGROUPS]; - for (int g = 0; g < NGROUPS; ++g) { - ernewn[g] = ern[g] - cdtdx_t1 * (rflux_t1(ir_t1,jr_t1,kr_t1,g) - - rflux_t1(il_t1,jl_t1,kl_t1,g)) - - cdtdx_t2 * (rflux_t2(ir_t2,jr_t2,kr_t2,g) - - rflux_t2(il_t2,jl_t2,kl_t2,g)) - + der[g]; - } -#endif - // Reset to original value if adding transverse terms - // made density negative - bool reset_state = false; - if (reset_density == 1 && rrnewn < 0.0_rt) { - rrnewn = rrn; - runewn = run; - rvnewn = rvn; - rwnewn = rwn; - renewn = ren; -#ifdef RADIATION - for (int g = 0; g < NGROUPS; ++g) { - ernewn[g] = ern[g]; + if (! reset_state) { + + // Convert back to primitive form + + qo_arr(i,j,k,QRHO) = U_int[URHO]; + Real rhoinv = 1.0_rt / U_int[URHO]; + + qo_arr(i,j,k,QU) = U_int[UMX] * rhoinv; + qo_arr(i,j,k,QV) = U_int[UMY] * rhoinv; + qo_arr(i,j,k,QW) = U_int[UMZ] * rhoinv; + + Real kineng = 0.5_rt * qo_arr(i,j,k,QRHO) * + (qo_arr(i,j,k,QU) * qo_arr(i,j,k,QU) + + qo_arr(i,j,k,QV) * qo_arr(i,j,k,QV) + + qo_arr(i,j,k,QW) * qo_arr(i,j,k,QW)); + + if ((U_int[UEDEN] - kineng) > castro::dual_energy_eta1 * U_int[UEDEN]) { + qo_arr(i,j,k,QREINT) = (U_int[UEDEN] - kineng) * rhoinv; + } else { + qo_arr(i,j,k,QREINT) = U_int[UEINT] * rhoinv; } -#endif - reset_state = true; - } - qo_arr(i,j,k,QRHO) = rrnewn; - qo_arr(i,j,k,QU) = runewn / rrnewn; - qo_arr(i,j,k,QV) = rvnewn / rrnewn; - qo_arr(i,j,k,QW) = rwnewn / rrnewn; - - // note: we run the risk of (rho e) being negative here - Real rhoekenn = 0.5_rt * (runewn * runewn + rvnewn * rvnewn + rwnewn * rwnewn) / rrnewn; - qo_arr(i,j,k,QREINT) = renewn - rhoekenn; - - if (!reset_state) { - if (reset_rhoe == 1 && qo_arr(i,j,k,QREINT) <= 0.0_rt) { - // If it is negative, reset the internal energy by - // using the discretized expression for updating - // (rho e). - qo_arr(i,j,k,QREINT) = q_arr(i,j,k,QREINT) - - cdtdx_t1 * (flux_t1(ir_t1,jr_t1,kr_t1,UEINT) - - flux_t1(il_t1,jl_t1,kl_t1,UEINT) + pt1av * dut1) - - cdtdx_t2 * (flux_t2(ir_t2,jr_t2,kr_t2,UEINT) - - flux_t2(il_t2,jl_t2,kl_t2,UEINT) + pt2av * dut2); + // Load passively advected quatities into q + for (int ipassive = 0; ipassive < npassive; ipassive++) { + int n = upassmap(ipassive); + int iq = qpassmap(ipassive); + qo_arr(i,j,k,iq) = U_int[n] * rhoinv; } - // If (rho e) is negative by this point, - // set it back to the original interface state, - // which turns off the transverse correction. + eos_rep_t eos_state; - if (qo_arr(i,j,k,QREINT) <= 0.0_rt) { - qo_arr(i,j,k,QREINT) = q_arr(i,j,k,QREINT); + eos_state.T = qo_arr(i,j,k,QTEMP); + eos_state.rho = qo_arr(i,j,k,QRHO); + eos_state.e = qo_arr(i,j,k,QREINT); + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = qo_arr(i,j,k,QFS+n); } +#if NAUX_NET > 0 + for (int n = 0; n < NumAux; n++) { + eos_state.aux[n] = qo_arr(i,j,k,QFX+n); + } +#endif - // Pretend QREINT has been fixed and transverse_use_eos != 1. - // If we are wrong, we will fix it later. - - // add the transverse term to the p evolution eq here - Real pnewn = q_arr(i,j,k,QPRES) - pt1new - pt2new; - qo_arr(i,j,k,QPRES) = pnewn; - } - else { - qo_arr(i,j,k,QPRES) = q_arr(i,j,k,QPRES); - qo_arr(i,j,k,QREINT) = q_arr(i,j,k,QREINT); - } + eos(eos_input_re, eos_state); - qo_arr(i,j,k,QPRES) = amrex::max(qo_arr(i,j,k,QPRES), small_p); + qo_arr(i,j,k,QTEMP) = eos_state.T; + qo_arr(i,j,k,QPRES) = eos_state.p; -#ifdef RADIATION - for (int g = 0; g < NGROUPS; ++g) { - qo_arr(i,j,k,QRAD+g) = ernewn[g]; - } + } else { - qo_arr(i,j,k,QPTOT) = qo_arr(i,j,k,QPRES); - for (int g = 0; g < NGROUPS; ++g) { - qo_arr(i,j,k,QPTOT) += lambda[g] * ernewn[g]; - } + for (int n = 0; n < NQ; n++) { + qo_arr(i,j,k,n) = q_arr(i,j,k,n); + } - qo_arr(i,j,k,QREITOT) = qo_arr(i,j,k,QREINT); - for (int g = 0; g < NGROUPS; ++g) { - qo_arr(i,j,k,QREITOT) += qo_arr(i,j,k,QRAD+g); } -#endif }); From b1de06d42c54e778c2970fed4281ae1281165e74 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 25 Jun 2021 11:30:42 -0400 Subject: [PATCH 3/5] fix undefined T --- Source/hydro/trans.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/hydro/trans.cpp b/Source/hydro/trans.cpp index b9a72ed466..e9883909af 100644 --- a/Source/hydro/trans.cpp +++ b/Source/hydro/trans.cpp @@ -283,7 +283,7 @@ Castro::actual_trans_single(const Box& bx, eos_rep_t eos_state; - eos_state.T = qo_arr(i,j,k,QTEMP); + eos_state.T = T_guess; // the input T may not be valid eos_state.rho = qo_arr(i,j,k,QRHO); eos_state.e = qo_arr(i,j,k,QREINT); for (int n = 0; n < NumSpec; n++) { @@ -590,7 +590,7 @@ Castro::actual_trans_final(const Box& bx, eos_rep_t eos_state; - eos_state.T = qo_arr(i,j,k,QTEMP); + eos_state.T = T_guess; // the input T may not be valid eos_state.rho = qo_arr(i,j,k,QRHO); eos_state.e = qo_arr(i,j,k,QREINT); for (int n = 0; n < NumSpec; n++) { From 1a32a8e61ab783c0dd0b4bcc3e1861610e00f5b9 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 25 Jun 2021 13:33:28 -0400 Subject: [PATCH 4/5] fix internal energyt --- Source/hydro/trans.cpp | 26 ++++++++++++++------------ 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/Source/hydro/trans.cpp b/Source/hydro/trans.cpp index e9883909af..9b876c2986 100644 --- a/Source/hydro/trans.cpp +++ b/Source/hydro/trans.cpp @@ -162,11 +162,11 @@ Castro::actual_trans_single(const Box& bx, U_int[UMX] = U_int[URHO] * q_arr(i,j,k,QU); U_int[UMY] = U_int[URHO] * q_arr(i,j,k,QV); U_int[UMZ] = U_int[URHO] * q_arr(i,j,k,QW); - U_int[UEDEN] = U_int[URHO] * q_arr(i,j,k,QREINT) + + U_int[UEDEN] = q_arr(i,j,k,QREINT) + 0.5_rt * (U_int[UMX] * U_int[UMX] + U_int[UMY] * U_int[UMY] + U_int[UMZ] * U_int[UMZ]) / U_int[URHO]; - U_int[UEINT] = U_int[URHO] * q_arr(i,j,k,QREINT); + U_int[UEINT] = q_arr(i,j,k,QREINT); for (int ipassive = 0; ipassive < npassive; ipassive++) { int n = upassmap(ipassive); @@ -269,11 +269,13 @@ Castro::actual_trans_single(const Box& bx, qo_arr(i,j,k,QW) * qo_arr(i,j,k,QW)); if ((U_int[UEDEN] - kineng) > castro::dual_energy_eta1 * U_int[UEDEN]) { - qo_arr(i,j,k,QREINT) = (U_int[UEDEN] - kineng) * rhoinv; + qo_arr(i,j,k,QREINT) = U_int[UEDEN] - kineng; } else { - qo_arr(i,j,k,QREINT) = U_int[UEINT] * rhoinv; + qo_arr(i,j,k,QREINT) = U_int[UEINT]; } + qo_arr(i,j,k,QREINT) = amrex::max(qo_arr(i,j,k,QREINT), small_dens * small_ener); + // Load passively advected quatities into q for (int ipassive = 0; ipassive < npassive; ipassive++) { int n = upassmap(ipassive); @@ -285,7 +287,7 @@ Castro::actual_trans_single(const Box& bx, eos_state.T = T_guess; // the input T may not be valid eos_state.rho = qo_arr(i,j,k,QRHO); - eos_state.e = qo_arr(i,j,k,QREINT); + eos_state.e = qo_arr(i,j,k,QREINT) * rhoinv; for (int n = 0; n < NumSpec; n++) { eos_state.xn[n] = qo_arr(i,j,k,QFS+n); } @@ -298,7 +300,7 @@ Castro::actual_trans_single(const Box& bx, eos(eos_input_re, eos_state); qo_arr(i,j,k,QTEMP) = eos_state.T; - qo_arr(i,j,k,QPRES) = eos_state.p; + qo_arr(i,j,k,QPRES) = amrex::max(eos_state.p, small_pres); } else { @@ -490,11 +492,11 @@ Castro::actual_trans_final(const Box& bx, U_int[UMX] = U_int[URHO] * q_arr(i,j,k,QU); U_int[UMY] = U_int[URHO] * q_arr(i,j,k,QV); U_int[UMZ] = U_int[URHO] * q_arr(i,j,k,QW); - U_int[UEDEN] = U_int[URHO] * q_arr(i,j,k,QREINT) + + U_int[UEDEN] = q_arr(i,j,k,QREINT) + 0.5_rt * (U_int[UMX] * U_int[UMX] + U_int[UMY] * U_int[UMY] + U_int[UMZ] * U_int[UMZ]) / U_int[URHO]; - U_int[UEINT] = U_int[URHO] * q_arr(i,j,k,QREINT); + U_int[UEINT] = q_arr(i,j,k,QREINT); for (int ipassive = 0; ipassive < npassive; ipassive++) { int n = upassmap(ipassive); @@ -576,9 +578,9 @@ Castro::actual_trans_final(const Box& bx, qo_arr(i,j,k,QW) * qo_arr(i,j,k,QW)); if ((U_int[UEDEN] - kineng) > castro::dual_energy_eta1 * U_int[UEDEN]) { - qo_arr(i,j,k,QREINT) = (U_int[UEDEN] - kineng) * rhoinv; + qo_arr(i,j,k,QREINT) = U_int[UEDEN] - kineng; } else { - qo_arr(i,j,k,QREINT) = U_int[UEINT] * rhoinv; + qo_arr(i,j,k,QREINT) = U_int[UEINT]; } // Load passively advected quatities into q @@ -592,7 +594,7 @@ Castro::actual_trans_final(const Box& bx, eos_state.T = T_guess; // the input T may not be valid eos_state.rho = qo_arr(i,j,k,QRHO); - eos_state.e = qo_arr(i,j,k,QREINT); + eos_state.e = qo_arr(i,j,k,QREINT) * rhoinv; for (int n = 0; n < NumSpec; n++) { eos_state.xn[n] = qo_arr(i,j,k,QFS+n); } @@ -605,7 +607,7 @@ Castro::actual_trans_final(const Box& bx, eos(eos_input_re, eos_state); qo_arr(i,j,k,QTEMP) = eos_state.T; - qo_arr(i,j,k,QPRES) = eos_state.p; + qo_arr(i,j,k,QPRES) = amrex::max(eos_state.p, small_pres); } else { From 91667fffebfe44b2742b3708a5a57a046732802a Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 25 Jun 2021 13:53:52 -0400 Subject: [PATCH 5/5] fix 2-d update --- Source/hydro/trans.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/hydro/trans.cpp b/Source/hydro/trans.cpp index 9b876c2986..1ccb5a6205 100644 --- a/Source/hydro/trans.cpp +++ b/Source/hydro/trans.cpp @@ -196,8 +196,8 @@ Castro::actual_trans_single(const Box& bx, #endif #if AMREX_SPACEDIM == 2 - U_int[n] += - hdt * (area_t(ir,jr,kr) * flux_t(ir,jr,kr,URHO) - - area_t(il,jl,kl) * flux_t(il,jl,kl,URHO)) * volinv; + U_int[n] += - hdt * (area_t(ir,jr,kr) * flux_t(ir,jr,kr,n) - + area_t(il,jl,kl) * flux_t(il,jl,kl,n)) * volinv; #else U_int[n] += - cdtdx * (flux_t(ir,jr,kr,n) - flux_t(il,jl,kl,n)); #endif