diff --git a/Source/CAMR.cpp b/Source/CAMR.cpp index 84da8c5..f359f98 100644 --- a/Source/CAMR.cpp +++ b/Source/CAMR.cpp @@ -183,6 +183,14 @@ CAMR::read_params() } } + // Check on PLM slope order -- the slope routines only distinguish + // between 1 (piecewise constant), 2 (MC-limited centered slope) and + // 4 (fourth-order blend), so any other value would silently be + // treated as if it were 2. + if (plm_iorder != 1 && plm_iorder != 2 && plm_iorder != 4) { + amrex::Error("CAMR::plm_iorder must be 1, 2 or 4"); + } + #ifdef AMREX_USE_EB // We only support PLM (not PPM) with using EB if (do_mol == 0 && ppm_type != 0) { diff --git a/Source/Hydro/Godunov/Godunov_utils_2D.H b/Source/Hydro/Godunov/Godunov_utils_2D.H index 923d786..3c82204 100644 --- a/Source/Hydro/Godunov/Godunov_utils_2D.H +++ b/Source/Hydro/Godunov/Godunov_utils_2D.H @@ -90,6 +90,16 @@ hydro_transd ( dAu = ugp - ugm; } + // If the transverse update would make the density negative then we reset the + // state below. We must test for that here as well so that the passive + // variables are reset along with the rest of the state -- otherwise they + // would be divided by the negative pre-reset density. These predicates must + // stay identical to the ones used in the QP / QM blocks below. + const bool reset_state_r = (transverse_reset_density == 1 && + qnormp(iv , QRHO) - flxrho < amrex::Real(0.0)); + const bool reset_state_l = (transverse_reset_density == 1 && + qnormm(ivpn, QRHO) - flxrho < amrex::Real(0.0)); + // **************************************************************************** // Update passive variables // **************************************************************************** @@ -102,19 +112,27 @@ hydro_transd ( amrex::Real rr, rrnew, compo; if (lo_face_not_covered) { - rr = qnormp(iv, QRHO); - rrnew = rr - flxrho; - compo = rr * qnormp(iv, nqp) - compn; - qp(iv, nqp) = compo / rrnew + hdt * srcpass; + if (reset_state_r) { + qp(iv, nqp) = qnormp(iv, nqp) + hdt * srcpass; + } else { + rr = qnormp(iv, QRHO); + rrnew = rr - flxrho; + compo = rr * qnormp(iv, nqp) - compn; + qp(iv, nqp) = compo / rrnew + hdt * srcpass; + } } else { qp(iv, nqp) = qnormp(iv, nqp); } if (hi_face_not_covered) { - rr = qnormm(ivpn, QRHO); - rrnew = rr - flxrho; - compo = rr * qnormm(ivpn, nqp) - compn; - qm(ivpn, nqp) = compo / rrnew + hdt * srcpass; + if (reset_state_l) { + qm(ivpn, nqp) = qnormm(ivpn, nqp) + hdt * srcpass; + } else { + rr = qnormm(ivpn, QRHO); + rrnew = rr - flxrho; + compo = rr * qnormm(ivpn, nqp) - compn; + qm(ivpn, nqp) = compo / rrnew + hdt * srcpass; + } } else { qm(ivpn, nqp) = qnormm(ivpn, nqp); } @@ -143,13 +161,12 @@ hydro_transd ( amrex::Real rvnewr = rvr - flxv; amrex::Real renewr = rer - flxe; - reset_state = false; - if (transverse_reset_density == 1 && rrnewr < 0.) { + reset_state = reset_state_r; + if (reset_state) { rrnewr = rrr; runewr = rur; rvnewr = rvr; renewr = rer; - reset_state = true; } // Convert back to primitive @@ -200,13 +217,12 @@ hydro_transd ( amrex::Real rvnewl = rvl - flxv; amrex::Real renewl = rel - flxe; - reset_state=false; - if (transverse_reset_density == 1 && rrnewl < 0.) { + reset_state = reset_state_l; + if (reset_state) { rrnewl = rrl; runewl = rul; rvnewl = rvl; renewl = rel; - reset_state = true; } qm(ivpn, QRHO) = rrnewl + hdt * srcr; diff --git a/Source/Hydro/Godunov/Godunov_utils_3D.H b/Source/Hydro/Godunov/Godunov_utils_3D.H index f0e0fd5..ddb54e6 100644 --- a/Source/Hydro/Godunov/Godunov_utils_3D.H +++ b/Source/Hydro/Godunov/Godunov_utils_3D.H @@ -78,6 +78,16 @@ hydro_transdo ( const amrex::Real flxe = no_cov_face ? cdtdx * (flxx(ivpt, UEDEN) - flxx(iv, UEDEN)) : amrex::Real(0.0); const amrex::Real c = qa(iv, QGAMC); + // If the transverse update would make the density negative then we reset the + // state below. We must test for that here as well so that the passive + // variables are reset along with the rest of the state -- otherwise they + // would be divided by the negative pre-reset density. These predicates must + // stay identical to the ones used in the QP / QM blocks below. + const bool reset_state_r = (transverse_reset_density == 1 && + qnormp(iv , QRHO) - flxrho < amrex::Real(0.0)); + const bool reset_state_l = (transverse_reset_density == 1 && + qnormm(ivpn, QRHO) - flxrho < amrex::Real(0.0)); + // Update passive variables for (int ipass = 0; ipass < NPASSIVE; ++ipass) { const int n = pmap.upassMap[ipass]; @@ -86,7 +96,7 @@ hydro_transdo ( amrex::Real rr, rrnew, compo; - if (lo_face_not_covered) + if (lo_face_not_covered && !reset_state_r) { rr = qnormp(iv, QRHO); rrnew = rr - flxrho; @@ -96,7 +106,7 @@ hydro_transdo ( qp(iv, nqp) = qnormp(iv, nqp); } - if (hi_face_not_covered) + if (hi_face_not_covered && !reset_state_l) { rr = qnormm(ivpn, QRHO); rrnew = rr - flxrho; @@ -143,15 +153,14 @@ hydro_transdo ( amrex::Real rwnewr = rwr - flxw; amrex::Real renewr = rer - flxe; - bool reset_state = false; + const bool reset_state = reset_state_r; - if (transverse_reset_density == 1 && rrnewr < 0.) { + if (reset_state) { rrnewr = rrr; runewr = rur; rvnewr = rvr; rwnewr = rwr; renewr = rer; - reset_state = true; } // Convert back to primitive @@ -201,14 +210,13 @@ hydro_transdo ( amrex::Real rwnewl = rwl - flxw; amrex::Real renewl = rel - flxe; - bool reset_state = false; - if (transverse_reset_density == 1 && rrnewl < 0.) { + const bool reset_state = reset_state_l; + if (reset_state) { rrnewl = rrl; runewl = rul; rvnewl = rvl; rwnewl = rwl; renewl = rel; - reset_state = true; } qm(ivpn, QRHO) = rrnewl; @@ -328,6 +336,14 @@ hydro_transdd ( amrex::Real rrnewl = rrl - flxrho; amrex::Real rrnewr = rrr - flxrho; + // If the transverse update would make the density negative then we reset the + // state below. We must test for that here as well so that the passive + // variables are reset along with the rest of the state -- otherwise they + // would be divided by the negative pre-reset density. These predicates must + // stay identical to the ones used in the qp / qm blocks below. + const bool reset_state_r = (transverse_reset_density == 1 && rrnewr < 0.); + const bool reset_state_l = (transverse_reset_density == 1 && rrnewl < 0.); + for (int ipass = 0; ipass < NPASSIVE; ++ipass) { const int n = pmap.upassMap[ipass]; const int nqp = pmap.qpassMap[ipass]; @@ -341,8 +357,12 @@ hydro_transdd ( // qp if (lo_face_not_covered) { - compo = rrr * qnormp(iv, nqp) - compn; - qp(iv, nqp) = compo / rrnewr + hdt * srcpass; + if (reset_state_r) { + qp(iv, nqp) = qnormp(iv, nqp) + hdt * srcpass; + } else { + compo = rrr * qnormp(iv, nqp) - compn; + qp(iv, nqp) = compo / rrnewr + hdt * srcpass; + } } else { qp(iv, nqp) = qnormp(iv, nqp); } @@ -351,8 +371,12 @@ hydro_transdd ( if (hi_face_not_covered) { - compo = rrl * qnormm(ivpn, nqp) - compn; - qm(ivpn, nqp) = compo / rrnewl + hdt * srcpass; + if (reset_state_l) { + qm(ivpn, nqp) = qnormm(ivpn, nqp) + hdt * srcpass; + } else { + compo = rrl * qnormm(ivpn, nqp) - compn; + qm(ivpn, nqp) = compo / rrnewl + hdt * srcpass; + } } else { qm(ivpn, nqp) = qnormm(ivpn, nqp); } @@ -397,15 +421,14 @@ hydro_transdd ( amrex::Real rwnewr = rwr - flxw; amrex::Real renewr = rer - flxe; - bool reset_state = false; + const bool reset_state = reset_state_r; - if (transverse_reset_density == 1 && rrnewr < 0.) { + if (reset_state) { rrnewr = rrr; runewr = rur; rvnewr = rvr; rwnewr = rwr; renewr = rer; - reset_state = true; } qp(iv, QRHO) = rrnewr + hdt * srcrho; @@ -445,14 +468,13 @@ hydro_transdd ( amrex::Real rwnewl = rwl - flxw; amrex::Real renewl = rel - flxe; - bool reset_state = false; - if (transverse_reset_density == 1 && rrnewl < 0.) { + const bool reset_state = reset_state_l; + if (reset_state) { rrnewl = rrl; runewl = rul; rvnewl = rvl; rwnewl = rwl; renewl = rel; - reset_state = true; } qm(ivpn, QRHO) = rrnewl + hdt * srcrho; diff --git a/Source/Hydro/Godunov/PLM.H b/Source/Hydro/Godunov/PLM.H index 7c8fccf..66d3b11 100644 --- a/Source/Hydro/Godunov/PLM.H +++ b/Source/Hydro/Godunov/PLM.H @@ -56,7 +56,7 @@ plm_slope ( dfm = dsgn * std::min(dlim, std::abs(dcen)); dlft = qp - qc; - drgt = (order == 4) ? qp2 - qp : 0.0; + drgt = qp2 - qp; dcen = 0.5 * (dlft + drgt); dsgn = amrex::Math::copysign(1.0, dcen); dlim = (dlft * drgt >= 0.0) @@ -72,7 +72,9 @@ plm_slope ( dlim = (dlft * drgt >= 0.0) ? 2.0 * std::min(std::abs(dlft), std::abs(drgt)) : 0.0; - dtemp = 4.0 / 3.0 * dcen - 1.0 / 6.0 * (dfp + dfm); + // The fourth-order blend only applies if we computed dfm/dfp above; + // for order == 2 this is just the MC-limited centered slope + dtemp = (order == 4) ? 4.0 / 3.0 * dcen - 1.0 / 6.0 * (dfp + dfm) : dcen; return flat * dsgn * std::min(dlim, std::abs(dtemp)); } @@ -155,7 +157,9 @@ plm_pslope ( dlim = (dlft * drgt >= 0.0) ? 2.0 * std::min(std::abs(dlft), std::abs(drgt)) : 0.0; - dtemp = 4.0 / 3.0 * dcen - 1.0 / 6.0 * (dfp + dfm); + // The fourth-order blend only applies if we computed dfm/dfp above; + // for order == 2 this is just the MC-limited centered slope + dtemp = (order == 4) ? 4.0 / 3.0 * dcen - 1.0 / 6.0 * (dfp + dfm) : dcen; amrex::Real dp = flat * dsgn * std::min(dlim, std::abs(dtemp)); diff --git a/Source/Hydro/Godunov/PPM.cpp b/Source/Hydro/Godunov/PPM.cpp index 2490623..52cccb4 100644 --- a/Source/Hydro/Godunov/PPM.cpp +++ b/Source/Hydro/Godunov/PPM.cpp @@ -395,7 +395,7 @@ trace_ppm( qm(ivp1, QPRES) = p_ref + (alphap + alpham) * csq_ref; qm(ivp1, QRHO) = std::max( qm(ivp1, QRHO), small_dens); - qm(ivp1, QPRES) = std::max( qm(ivp1, QPRES), small_dens); + qm(ivp1, QPRES) = std::max( qm(ivp1, QPRES), small_pres); // transverse velocities qm(ivp1, QUT) = Ip[QUT][1] /*+ hdt * Ip_src[QUT][1]*/;