diff --git a/Source/hydro/Castro_ctu_hydro.cpp b/Source/hydro/Castro_ctu_hydro.cpp index 9d37a0a7bf..5365ac854e 100644 --- a/Source/hydro/Castro_ctu_hydro.cpp +++ b/Source/hydro/Castro_ctu_hydro.cpp @@ -569,9 +569,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 @@ -613,9 +613,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 @@ -697,9 +697,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))); @@ -725,9 +725,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] @@ -773,9 +773,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))); @@ -804,9 +804,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] @@ -852,9 +852,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))); @@ -883,9 +883,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 @@ -946,9 +946,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 SIMPLIFIED_SDC #ifdef REACTIONS @@ -1025,9 +1025,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 SIMPLIFIED_SDC #ifdef REACTIONS @@ -1106,9 +1106,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 SIMPLIFIED_SDC #ifdef REACTIONS diff --git a/Source/hydro/trans.cpp b/Source/hydro/trans.cpp index 724b8e2a66..98f600d88c 100644 --- a/Source/hydro/trans.cpp +++ b/Source/hydro/trans.cpp @@ -154,135 +154,73 @@ 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] = 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] = 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); - -#ifdef RADIATION - Real lambda[NGROUPS]; - Real ergp[NGROUPS]; - Real ergm[NGROUPS]; - - 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); - } -#endif + // TODO: how does the hybrid stuff fit in here? - // 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 + // now do the transverse flux update #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; -#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); + const Real volinv = 1.0_rt / vol(il,jl,kl); #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); -#endif + for (int n = 0; n < NUM_STATE; n++) { -#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]; - } + if (n == UTEMP) { + continue; + } - if (radiation::fspace_advection_type == 1 && radiation::comoving) { - for (int g = 0; g < NGROUPS; ++g) { - Real eddf = Edd_factor(lambda[g]); - Real f1 = 0.5_rt * (1.0_rt - eddf); - der[g] = cdtdx * uav * f1 * (ergp[g] - ergm[g]); +#ifdef SHOCK_VAR + if (n == USHK) { + continue; } - } - else if (radiation::fspace_advection_type == 2) { +#endif + #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]); - Real f1 = 0.5_rt * (1.0_rt - eddf); - der[g] = -hdt * f1 * 0.5_rt * (ergp[g] + ergm[g]) * divu; - } + 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 - for (int g = 0; g < NGROUPS; g++) { - Real eddf = Edd_factor(lambda[g]); - 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]; - } - } + U_int[n] += - cdtdx * (flux_t(ir,jr,kr,n) - flux_t(il,jl,kl,n)); #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 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 AMREX_SPACEDIM == 2 + Real du = area_t(ir,jr,kr) * ugp - area_t(il,jl,kl) * ugm; +#else + 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 // @@ -291,144 +229,86 @@ 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 + // 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; - for (int ipassive = 0; ipassive < npassive; ++ipassive) { - int nqp = qpassmap(ipassive); - qo_arr(i,j,k,nqp) = q_arr(i,j,k,nqp); - } -#ifdef RADIATION - for (int g = 0; g < NGROUPS; ++g) { - ernewn[g] = ern[g]; - } -#endif + if (reset_density == 1 && U_int[URHO] < 0.0_rt) { 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; + + 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; + } else { + qo_arr(i,j,k,QREINT) = U_int[UEINT]; } - // If (rho e) is negative by this point, - // set it back to the original interface state, - // which turns off the transverse correction. + qo_arr(i,j,k,QREINT) = amrex::max(qo_arr(i,j,k,QREINT), small_dens * small_ener); - if (qo_arr(i,j,k,QREINT) <= 0.0_rt) { - qo_arr(i,j,k,QREINT) = q_arr(i,j,k,QREINT); + // 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; } - // Pretend QREINT has been fixed and transverse_use_eos != 1. - // If we are wrong, we will fix it later. + eos_rep_t eos_state; - // 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)); + 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) * rhoinv; + 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 - qo_arr(i,j,k,QPRES) = amrex::max(pnewn, small_p); - } - 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); -#ifdef RADIATION - for (int g = 0; g < NGROUPS; ++g) { - qo_arr(i,j,k,QRAD + g) = ernewn[g]; - } + qo_arr(i,j,k,QTEMP) = eos_state.T; + qo_arr(i,j,k,QPRES) = amrex::max(eos_state.p, small_pres); - 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]; - } + } else { - 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 + for (int n = 0; n < NQ; n++) { + qo_arr(i,j,k,n) = q_arr(i,j,k,n); + } + } }); } @@ -597,25 +477,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 - for (int ipassive = 0; ipassive < npassive; ++ipassive) { + 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] = 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] = q_arr(i,j,k,QREINT); + + 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 @@ -625,231 +528,86 @@ 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 (radiation::fspace_advection_type == 1 && radiation::comoving) { - for (int g = 0; g < NGROUPS; ++g) { - Real eddf = Edd_factor(lambda[g]); - 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 (radiation::fspace_advection_type == 2) { - for (int g = 0; g < NGROUPS; ++g) { - Real eddf = Edd_factor(lambda[g]); - 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 + if (! reset_state) { - // 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; - for (int ipassive = 0; ipassive < npassive; ++ipassive) { - int nqp = qpassmap(ipassive); - qo_arr(i,j,k,nqp) = q_arr(i,j,k,nqp); - } -#ifdef RADIATION - for (int g = 0; g < NGROUPS; ++g) { - ernewn[g] = ern[g]; + // 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; + } else { + qo_arr(i,j,k,QREINT) = U_int[UEINT]; } -#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 = 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) * rhoinv; + 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) = amrex::max(eos_state.p, small_pres); -#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 });