[WIP] add clipping to the species interface states for true SDC

Source/hydro/Castro_hydro.H

@@ -795,6 +795,14 @@
                                amrex::Array4<amrex::Real> const& qm,
                                amrex::Array4<amrex::Real> const& qp);
 
+///
+/// Clip the mass fractions to lie between [0, 1] and renormalize
+/// so they sum to 1.  This operates on interfaces.
+///
+   void enforce_species_sum(const amrex::Box& bx,
+                            amrex::Array4<amrex::Real> const& qm,
+                            amrex::Array4<amrex::Real> const& qp);
+
 
     void scale_flux(const amrex::Box& bx,
 #if AMREX_SPACEDIM == 1

Source/hydro/Castro_mol.cpp

@@ -95,44 +95,47 @@ Castro::mol_plm_reconstruct(const Box& bx,
   {
 
 
-   // this is a loop over zones.  For each slope in the zone, fill the
-   // two adjacent edge states (e.g., the right state at i-1/2 and the
-   // left state at i+1/2
+      // this is a loop over zones.  For each slope in the zone, fill the
+      // two adjacent edge states (e.g., the right state at i-1/2 and the
+      // left state at i+1/2
 
-   if (idir == 0) {
+      if (idir == 0) {
 
-     // left state at i+1/2 interface
-     qm(i+1,j,k,n) = q_arr(i,j,k,n) + 0.5_rt*dq(i,j,k,n);
+          // left state at i+1/2 interface
+          qm(i+1,j,k,n) = q_arr(i,j,k,n) + 0.5_rt*dq(i,j,k,n);
 
-     // right state at i-1/2 interface
-     qp(i,j,k,n) = q_arr(i,j,k,n) - 0.5_rt*dq(i,j,k,n);
+          // right state at i-1/2 interface
+          qp(i,j,k,n) = q_arr(i,j,k,n) - 0.5_rt*dq(i,j,k,n);
 
 
 #if AMREX_SPACEDIM >= 2
-   } else if (idir == 1) {
+      } else if (idir == 1) {
 
-     // left state at j+1/2 interface
-     qm(i,j+1,k,n) = q_arr(i,j,k,n) + 0.5_rt*dq(i,j,k,n);
+          // left state at j+1/2 interface
+          qm(i,j+1,k,n) = q_arr(i,j,k,n) + 0.5_rt*dq(i,j,k,n);
 
-     // right state at j-1/2 interface
-     qp(i,j,k,n) = q_arr(i,j,k,n) - 0.5_rt*dq(i,j,k,n);
+          // right state at j-1/2 interface
+          qp(i,j,k,n) = q_arr(i,j,k,n) - 0.5_rt*dq(i,j,k,n);
 #endif
 
 #if AMREX_SPACEDIM == 3
-   } else {
+      } else {
 
-     // left state at k+1/2 interface
-     qm(i,j,k+1,n) = q_arr(i,j,k,n) + 0.5_rt*dq(i,j,k,n);
+          // left state at k+1/2 interface
+          qm(i,j,k+1,n) = q_arr(i,j,k,n) + 0.5_rt*dq(i,j,k,n);
 
-     // right state at k-1/2 interface
-     qp(i,j,k,n) = q_arr(i,j,k,n) - 0.5_rt*dq(i,j,k,n);
+          // right state at k-1/2 interface
+          qp(i,j,k,n) = q_arr(i,j,k,n) - 0.5_rt*dq(i,j,k,n);
 
 #endif
-   }
+      }
 
   });
 
 
+  // normalize species on interfaces
+  enforce_species_sum(bx, qm, qp);
+
   // special care for reflecting BCs
 
   // we have to do this after the loops above, since here we will
@@ -188,6 +191,10 @@ Castro::mol_ppm_reconstruct(const Box& bx,
 
   });
 
+
+  // normalize species on interfaces
+  enforce_species_sum(bx, qm, qp);
+
   // special care for reflecting BCs
 
   // we have to do this after the loops above, since here we will

Source/hydro/Castro_mol_hydro.cpp

@@ -257,6 +257,10 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update)
                      qm_arr, qp_arr);
             }
 
+            // enforce that the species sum to 1 on interfaces
+
+            enforce_species_sum(ibx[idir], qm_arr, qp_arr);
+
             // get the face-averaged state and flux, <q> and F(<q>),
             // in the idir direction by solving the Riemann problem
             // operate on ibx[idir]

Source/hydro/advection_util.cpp

@@ -745,3 +745,34 @@ Castro::enforce_reflect_states(const Box& bx, const int idir,
 
     }
 }
+
+
+void
+Castro::enforce_species_sum(const Box& bx,
+                            Array4<Real> const& qm,
+                            Array4<Real> const& qp) {
+
+    // this is a loop over interfaces
+
+    amrex::ParallelFor(bx,
+    [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+    {
+        Real sum_xm{0.0};
+        Real sum_xp{0.0};
+
+        for (int n = 0; n < NumSpec; ++n) {
+            qm(i,j,k,QFS+n) = std::max(0.0, std::min(1.0, qm(i,j,k,QFS+n)));
+            sum_xm += qm(i,j,k,QFS+n);
+
+            qp(i,j,k,QFS+n) = std::max(0.0, std::min(1.0, qp(i,j,k,QFS+n)));
+            sum_xp += qp(i,j,k,QFS+n);
+        }
+
+        for (int n = 0; n < NumSpec; ++n) {
+            qm(i,j,k,QFS+n) /= sum_xm;
+            qp(i,j,k,QFS+n) /= sum_xp;
+        }
+
+    });
+
+}