Revert "implement full well-balancing in PPM (#2945)"

This reverts commit db92cba.

Docs/source/hse.rst

@@ -74,12 +74,6 @@ then done on the parabola, again working with :math:`\tilde{p}`.
 Finally, the parabola values are updated to include the hydrostatic
 pressure.
 
-.. index:: castro.ppm_well_balanced
-
-We can do better with PPM, and only use the perturbational pressure,
-$\tilde{p}$, in the characteristic tracing and then add back the
-hydrostatic pressure to the interface afterwards.  This is done via
-``castro.ppm_well_balanced=1``.
 
 Fully fourth-order method
 -------------------------
@@ -171,3 +165,5 @@ test the different HSE approaches.  This sets up a 1-d X-ray burst
 atmosphere (based on the ``flame_wave`` setup).  Richardson
 extrapolation can be used to measure the convergence rate (or just
 look at how the peak velocity changes).
+
+

Exec/gravity_tests/hse_convergence/README.md

@@ -31,13 +31,10 @@ To run this problem, use one of the convergence scripts:
 
   * convergence_sdc.sh :
 
-    this uses the `TRUE_SDC` integration, with the following variations:
-    1. SDC-2 + PLM and reflecting BCs
-    2. SDC-2 + PPM and reflecting BCs
-    3. SDC-2 + PLM with HSE BCs
-    4. SDC-2 + PPM with HSE BCs
-    5. SDC-4 + reflect
+    this uses the TRUE_SDC integration, first with SDC-2 + PLM and
+    reflecting BCs, the SDC-2 + PPM and reflecting BCs, then the same
+    but HSE BCs, and finally SDC-4 + reflect
 
     These tests show that the PLM + reflect (which uses the
-    well-balanced `use_pslope`) and the SDC-4 + reflect give the lowest
+    well-balanced use_pslope) and the SDC-4 + reflect give the lowest
     errors and expected (or better) convergence.

Exec/gravity_tests/hse_convergence/convergence_ppm.sh

@@ -102,30 +102,3 @@ ${EXEC} inputs.ppm.512 ${RUNPARAMS} >& 512.out
 pfile=`ls -t | grep -i hse_512_plt | head -1`
 fextrema.gnu.ex -v magvel ${pfile} | grep -i magvel >> ${ofile}
 
-
-## ppm + reflect + well-balanced
-
-ofile=ppm-reflect-wellbalanced.converge.out
-
-RUNPARAMS="
-castro.lo_bc=3
-castro.hi_bc=3
-castro.ppm_well_balanced=1
-"""
-
-${EXEC} inputs.ppm.64 ${RUNPARAMS} >& 64.out
-pfile=`ls -t | grep -i hse_64_plt | head -1`
-fextrema.gnu.ex -v magvel ${pfile} | grep -i magvel > ${ofile}
-
-${EXEC} inputs.ppm.128 ${RUNPARAMS} >& 128.out
-pfile=`ls -t | grep -i hse_128_plt | head -1`
-fextrema.gnu.ex -v magvel ${pfile} | grep -i magvel >> ${ofile}
-
-${EXEC} inputs.ppm.256 ${RUNPARAMS} >& 256.out
-pfile=`ls -t | grep -i hse_256_plt | head -1`
-fextrema.gnu.ex -v magvel ${pfile} | grep -i magvel >> ${ofile}
-
-${EXEC} inputs.ppm.512 ${RUNPARAMS} >& 512.out
-pfile=`ls -t | grep -i hse_512_plt | head -1`
-fextrema.gnu.ex -v magvel ${pfile} | grep -i magvel >> ${ofile}
-

Exec/gravity_tests/hse_convergence/convergence_sdc.sh

@@ -67,15 +67,14 @@ pfile=`ls -t | grep -i hse_512_plt | head -1`
 fextrema.gnu.ex -v magvel ${pfile} | grep -i magvel >> ${ofile}
 
 
-## sdc-2 + HSE
+## sdc-2 + ppm
 
-ofile=sdc2.converge.out
+ofile=sdc2-ppm.converge.out
 
 RUNPARAMS="
 castro.time_integration_method=2
 castro.sdc_order=2
-castro.ppm_type=0
-castro.use_pslope=1
+castro.ppm_type=1
 castro.limit_fourth_order=1
 castro.use_reconstructed_gamma1=1
 """
@@ -97,14 +96,15 @@ pfile=`ls -t | grep -i hse_512_plt | head -1`
 fextrema.gnu.ex -v magvel ${pfile} | grep -i magvel >> ${ofile}
 
 
-## sdc-2 + ppm
+## sdc-2 + HSE
 
-ofile=sdc2-ppm.converge.out
+ofile=sdc2.converge.out
 
 RUNPARAMS="
 castro.time_integration_method=2
 castro.sdc_order=2
-castro.ppm_type=1
+castro.ppm_type=0
+castro.use_pslope=1
 castro.limit_fourth_order=1
 castro.use_reconstructed_gamma1=1
 """
@@ -126,6 +126,7 @@ pfile=`ls -t | grep -i hse_512_plt | head -1`
 fextrema.gnu.ex -v magvel ${pfile} | grep -i magvel >> ${ofile}
 
 
+
 ## sdc-4 + reflect
 
 ofile=sdc4-reflect.converge.out

Exec/science/flame_wave/ci-benchmarks/job_info_params.txt

@@ -64,7 +64,6 @@
     castro.dual_energy_eta1 = 1
     castro.dual_energy_eta2 = 0.0001
     castro.use_pslope = 0
-    castro.ppm_well_balanced = 0
     castro.pslope_cutoff_density = -1e+20
     castro.limit_fluxes_on_small_dens = 0
     castro.speed_limit = 0

Source/driver/_cpp_parameters

@@ -161,10 +161,6 @@ dual_energy_eta2             Real          1.0e-4
 # does well with HSE
 use_pslope                   bool           0
 
-# for PPM, do we only use the perturbational pressure in the characteristic
-# tracing?  This is more indepth than the simple `use_pslope` approach.
-ppm_well_balanced            bool           0
-
 # if we are using pslope, below what density to we turn off the well-balanced
 # reconstruction?
 pslope_cutoff_density        Real          -1.e20

Source/hydro/ppm.H

@@ -1,18 +1,16 @@
 
 #include <cmath>
-
-#include <AMReX_REAL.H>
 #include <reconstruction.H>
 
 #ifndef PPM_H
 #define PPM_H
 
-using namespace amrex::literals;
+using namespace amrex;
 using namespace reconstruction;
 
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
 int
-check_trace_source(amrex::Array4<amrex::Real const> const& srcQ, const int idir,
+check_trace_source(Array4<Real const> const& srcQ, const int idir,
                    const int i, const int j, const int k, const int ncomp) {
 
     int do_trace = 0;
@@ -55,29 +53,29 @@ check_trace_source(amrex::Array4<amrex::Real const> const& srcQ, const int idir,
 ///
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
 void
-ppm_reconstruct(const amrex::Real* s,
-                const amrex::Real flatn, amrex::Real& sm, amrex::Real& sp) {
+ppm_reconstruct(const Real* s,
+                const Real flatn, Real& sm, Real& sp) {
 
 
     if (ppm_do_limiting) {
 
         // Compute van Leer slopes
 
-        amrex::Real dsl = 2.0_rt * (s[im1] - s[im2]);
-        amrex::Real dsr = 2.0_rt * (s[i0] - s[im1]);
+        Real dsl = 2.0_rt * (s[im1] - s[im2]);
+        Real dsr = 2.0_rt * (s[i0] - s[im1]);
 
-        amrex::Real dsvl_l = 0.0_rt;
+        Real dsvl_l = 0.0_rt;
         if (dsl*dsr > 0.0_rt) {
-          amrex::Real dsc = 0.5_rt * (s[i0] - s[im2]);
+          Real dsc = 0.5_rt * (s[i0] - s[im2]);
           dsvl_l = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc), std::abs(dsl), std::abs(dsr));
         }
 
         dsl = 2.0_rt * (s[i0] - s[im1]);
         dsr = 2.0_rt * (s[ip1] - s[i0]);
 
-        amrex::Real dsvl_r = 0.0_rt;
+        Real dsvl_r = 0.0_rt;
         if (dsl*dsr > 0.0_rt) {
-          amrex::Real dsc = 0.5_rt * (s[ip1] - s[im1]);
+          Real dsc = 0.5_rt * (s[ip1] - s[im1]);
           dsvl_r = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc), std::abs(dsl),std::abs(dsr));
         }
 
@@ -96,17 +94,17 @@ ppm_reconstruct(const amrex::Real* s,
 
         dsvl_l = 0.0_rt;
         if (dsl*dsr > 0.0_rt) {
-            amrex::Real dsc = 0.5_rt * (s[ip1] - s[im1]);
-            dsvl_l = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc), std::abs(dsl), std::abs(dsr));
+          Real dsc = 0.5_rt * (s[ip1] - s[im1]);
+          dsvl_l = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc), std::abs(dsl), std::abs(dsr));
         }
 
         dsl = 2.0_rt * (s[ip1] - s[i0]);
         dsr = 2.0_rt * (s[ip2] - s[ip1]);
 
         dsvl_r = 0.0_rt;
         if (dsl*dsr > 0.0_rt) {
-            amrex::Real dsc = 0.5_rt * (s[ip2] - s[i0]);
-            dsvl_r = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc), std::abs(dsl), std::abs(dsr));
+          Real dsc = 0.5_rt * (s[ip2] - s[i0]);
+          dsvl_r = std::copysign(1.0_rt, dsc) * amrex::min(std::abs(dsc), std::abs(dsl), std::abs(dsr));
         }
 
         // Interpolate s to edges
@@ -155,7 +153,7 @@ ppm_reconstruct(const amrex::Real* s,
 /// the gravitational force by only limiting on the wave-generating pressure.
 /// This uses the standard PPM limiters described in Colella & Woodward (1984)
 ///
-/// @param rho         Real[nslp] giving the density in zones i-2, i-1, i, i+1, i+2
+/// @param rho          Real[nslp] giving the density in zones i-2, i-1, i, i+1, i+2
 /// @param p           Real[nslp] giving the pressure in zones i-2, i-1, i, i+1, i+2
 /// @param src         Real[nslp] the source in the velocity equation (e.g. g) in zones
 ///                    i-2, i-1, i, i+2, i+2
@@ -166,25 +164,23 @@ ppm_reconstruct(const amrex::Real* s,
 /// @param sm     The value of the parabola on the left edge of the zone
 /// @param sp     The value of the parabola on the right edge of the zone
 ///
-/// @return a bool indicating whether HSE correction was done
-///
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-bool
-ppm_reconstruct_pslope(const amrex::Real* rho, const amrex::Real* p, const amrex::Real* src, const amrex::Real flatn,
-                       const amrex::Real dx,
-                       amrex::Real& sm, amrex::Real& sp) {
+void
+ppm_reconstruct_pslope(const Real* rho, const Real* p, const Real* src, const Real flatn,
+                       const Real dx,
+                       Real& sm, Real& sp) {
 
-    amrex::Real tp[nslp];
+    Real tp[nslp];
 
     // compute the hydrostatic pressure in each zone center starting with i0
 
-    amrex::Real p0_hse = p[i0];
+    Real p0_hse = p[i0];
 
-    amrex::Real pp1_hse = p0_hse + 0.25_rt*dx * (rho[i0] + rho[ip1]) * (src[i0] + src[ip1]);
-    amrex::Real pp2_hse = pp1_hse + 0.25_rt*dx * (rho[ip1] + rho[ip2]) * (src[ip1] + src[ip2]);
+    Real pp1_hse = p0_hse + 0.25_rt*dx * (rho[i0] + rho[ip1]) * (src[i0] + src[ip1]);
+    Real pp2_hse = pp1_hse + 0.25_rt*dx * (rho[ip1] + rho[ip2]) * (src[ip1] + src[ip2]);
 
-    amrex::Real pm1_hse = p0_hse - 0.25_rt*dx * (rho[i0] + rho[im1]) * (src[i0] + src[im1]);
-    amrex::Real pm2_hse = pm1_hse - 0.25_rt*dx * (rho[im1] + rho[im2]) * (src[im1] + src[im2]);
+    Real pm1_hse = p0_hse - 0.25_rt*dx * (rho[i0] + rho[im1]) * (src[i0] + src[im1]);
+    Real pm2_hse = pm1_hse - 0.25_rt*dx * (rho[im1] + rho[im2]) * (src[im1] + src[im2]);
 
     // this only makes sense if the pressures are positive
     bool ptest = p0_hse < 0.0 ||
@@ -194,9 +190,7 @@ ppm_reconstruct_pslope(const amrex::Real* rho, const amrex::Real* p, const amrex
                  pm2_hse < 0.0;
 
 
-    bool do_hse = !ptest && rho[i0] >= pslope_cutoff_density;
-
-    if (do_hse) {
+    if (!ptest && rho[i0] >= pslope_cutoff_density) {
 
         // subtract off the hydrostatic pressure
 
@@ -224,16 +218,12 @@ ppm_reconstruct_pslope(const amrex::Real* rho, const amrex::Real* p, const amrex
 
     // now correct sm and sp to be back to the full pressure by
     // adding in the hydrostatic pressure at the interface
-    // if we are doing the full well-balanced method, then we
-    // don't do this until after the characteristic tracing
 
-    if (do_hse && !castro::ppm_well_balanced) {
+    if (!ptest && rho[i0] >= pslope_cutoff_density) {
         sp += p[i0] + 0.5 * dx * rho[i0] * src[i0];
         sm += p[i0] - 0.5 * dx * rho[i0] * src[i0];
     }
 
-    return do_hse;
-
 }
 
 
@@ -254,14 +244,13 @@ ppm_reconstruct_pslope(const amrex::Real* rho, const amrex::Real* p, const amrex
 ///
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
 void
-ppm_int_profile(const amrex::Real sm, const amrex::Real sp, const amrex::Real sc,
-                const amrex::Real u, const amrex::Real c, const amrex::Real dtdx,
-                amrex::Real* Ip, amrex::Real* Im) {
+ppm_int_profile(const Real sm, const Real sp, const Real sc,
+                const Real u, const Real c, const Real dtdx,
+                Real* Ip, Real* Im) {
 
   // Integrate the parabolic profile to the edge of the cell.
 
   // compute x-component of Ip and Im
-
   Real s6 = 6.0_rt * sc - 3.0_rt * (sm + sp);
 
   // Ip/m is the integral under the parabola for the extent