use one-sided PPM stencils at reflecting BCs

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

@@ -49,6 +49,7 @@
     castro.hybrid_hydro = 0
     castro.ppm_type = 1
     castro.ppm_do_limiting = 1
+    castro.ppm_one_sided_stencils = 0
     castro.mhd_limit_characteristic = 1
     castro.ppm_temp_fix = 0
     castro.plm_iorder = 2

Source/driver/_cpp_parameters

@@ -89,6 +89,9 @@ ppm_type                     int           1
 # do we limit the ppm parabola?
 ppm_do_limiting              bool           1
 
+# do we use one-sided stencils at reflecting boundaries?
+ppm_one_sided_stencils       bool           0
+
 # For MHD + PLM, do we limit on characteristic or primitive variables
 mhd_limit_characteristic     bool           1
 

Source/hydro/Castro_mol.cpp

@@ -55,7 +55,7 @@ Castro::mol_plm_reconstruct(const Box& bx,
     Real s[nslp];
     Real flat = flatn_arr(i,j,k);
 
-    load_stencil(q_arr, idir, i, j, k, n, s);
+    load_stencil(q_arr, reconstruction::Centering::zone_centered, idir, i, j, k, n, s);
 
     // normal velocity?
     bool vtest = n == QU+idir;
@@ -83,9 +83,9 @@ Castro::mol_plm_reconstruct(const Box& bx,
                                     (idir == 1 && j == domhi[1]) ||
                                     (idir == 2 && k == domhi[2]));
 
-      load_stencil(q_arr, idir, i, j, k, QPRES, s);
-      load_stencil(q_arr, idir, i, j, k, QRHO, trho);
-      load_stencil(src_q_arr, idir, i, j, k, QU+idir, src);
+      load_stencil(q_arr, reconstruction::Centering::zone_centered, idir, i, j, k, QPRES, s);
+      load_stencil(q_arr, reconstruction::Centering::zone_centered, idir, i, j, k, QRHO, trho);
+      load_stencil(src_q_arr, reconstruction::Centering::zone_centered, idir, i, j, k, QU+idir, src);
 
       Real dp = dq(i,j,k,QPRES);
       pslope(trho, s, src, flat, lo_bc_test, hi_bc_test, dx[idir], dp);
@@ -164,8 +164,9 @@ Castro::mol_ppm_reconstruct(const Box& bx,
     Real sm;
     Real sp;
 
-    load_stencil(q_arr, idir, i, j, k, n, s);
-    ppm_reconstruct(s, flat, sm, sp);
+    load_stencil(q_arr, reconstruction::Centering::zone_centered, idir, i, j, k, n, s);
+    ppm_reconstruct(s, reconstruction::Centering::zone_centered,
+                    flat, sm, sp);
 
     if (idir == 0) {
       // right state at i-1/2
@@ -351,7 +352,7 @@ Castro::compute_flux_from_q(const Box& bx,
   //
 
   int iu, iv1, iv2;
-  int im1, im2, im3;
+  int imx1, imx2, imx3;
 
   auto coord = geom.Coord();
   auto mom_check = mom_flux_has_p(idir, idir, coord);
@@ -360,25 +361,25 @@ Castro::compute_flux_from_q(const Box& bx,
     iu = QU;
     iv1 = QV;
     iv2 = QW;
-    im1 = UMX;
-    im2 = UMY;
-    im3 = UMZ;
+    imx1 = UMX;
+    imx2 = UMY;
+    imx3 = UMZ;
 
   } else if (idir == 1) {
     iu = QV;
     iv1 = QU;
     iv2 = QW;
-    im1 = UMY;
-    im2 = UMX;
-    im3 = UMZ;
+    imx1 = UMY;
+    imx2 = UMX;
+    imx3 = UMZ;
 
   } else {
     iu = QW;
     iv1 = QU;
     iv2 = QV;
-    im1 = UMZ;
-    im2 = UMX;
-    im3 = UMY;
+    imx1 = UMZ;
+    imx2 = UMX;
+    imx3 = UMY;
   }
 
 #ifdef HYBRID_MOMENTUM
@@ -395,12 +396,12 @@ Castro::compute_flux_from_q(const Box& bx,
     // Compute fluxes, order as conserved state (not q)
     F(i,j,k,URHO) = qint(i,j,k,QRHO)*u_adv;
 
-    F(i,j,k,im1) = F(i,j,k,URHO)*qint(i,j,k,iu);
+    F(i,j,k,imx1) = F(i,j,k,URHO)*qint(i,j,k,iu);
     if (mom_check) {
-      F(i,j,k,im1) += qint(i,j,k,QPRES);
+      F(i,j,k,imx1) += qint(i,j,k,QPRES);
     }
-    F(i,j,k,im2) = F(i,j,k,URHO)*qint(i,j,k,iv1);
-    F(i,j,k,im3) = F(i,j,k,URHO)*qint(i,j,k,iv2);
+    F(i,j,k,imx2) = F(i,j,k,URHO)*qint(i,j,k,iv1);
+    F(i,j,k,imx3) = F(i,j,k,URHO)*qint(i,j,k,iv2);
 
     Real rhoetot = rhoeint + 0.5_rt * qint(i,j,k,QRHO)*
       (qint(i,j,k,iu)*qint(i,j,k,iu) +