add some protection against negatives in 4th order

pyro/compressible_fv4/fluxes.py

@@ -57,6 +57,21 @@ def fluxes(myd, rp, ivars):
     U_cc[:, :, ivars.iymom] = myd.to_centers("y-momentum")
     U_cc[:, :, ivars.iener] = myd.to_centers("energy")
 
+    # the mask will be 1 in any zone where the density or energy
+    # is unphysical do to the conversion from averages to centers
+
+    rhoe = U_cc[..., ivars.iener] - 0.5 * (U_cc[..., ivars.ixmom]**2 +
+                                           U_cc[..., ivars.iymom]**2) / U_cc[..., ivars.idens]
+
+    mask = myg.scratch_array(dtype=np.uint8)
+    mask[:, :] = np.where(np.logical_or(U_cc[:, :, ivars.idens] < 0, rhoe < 0),
+                          1, 0)
+
+    U_cc[..., ivars.idens] = np.where(mask == 1, U_avg[..., ivars.idens], U_cc[..., ivars.idens])
+    U_cc[..., ivars.ixmom] = np.where(mask == 1, U_avg[..., ivars.ixmom], U_cc[..., ivars.ixmom])
+    U_cc[..., ivars.iymom] = np.where(mask == 1, U_avg[..., ivars.iymom], U_cc[..., ivars.iymom])
+    U_cc[..., ivars.iener] = np.where(mask == 1, U_avg[..., ivars.iener], U_cc[..., ivars.iener])
+
     # compute the primitive variables of both the cell-center and averages
     q_bar = comp.cons_to_prim(U_avg, gamma, ivars, myd.grid)
     q_cc = comp.cons_to_prim(U_cc, gamma, ivars, myd.grid)
@@ -66,6 +81,15 @@ def fluxes(myd, rp, ivars):
     for n in range(ivars.nq):
         q_avg.v(n=n, buf=3)[:, :] = q_cc.v(n=n, buf=3) + myg.dx**2/24.0 * q_bar.lap(n=n, buf=3)
 
+    # enforce positivity
+    q_avg.v(n=ivars.irho, buf=3)[:, :] = np.where(q_avg.v(n=ivars.irho, buf=3) > 0,
+                                                  q_avg.v(n=ivars.irho, buf=3),
+                                                  q_cc.v(n=ivars.irho, buf=3))
+
+    q_avg.v(n=ivars.ip, buf=3)[:, :] = np.where(q_avg.v(n=ivars.ip, buf=3) > 0,
+                                                q_avg.v(n=ivars.ip, buf=3),
+                                                q_cc.v(n=ivars.ip, buf=3))
+
     # flattening -- there is a single flattening coefficient (xi) for all directions
     use_flattening = rp.get_param("compressible.use_flattening")
 

pyro/mesh/fv.py

@@ -3,6 +3,8 @@
 
 """
 
+import numpy as np
+
 from pyro.mesh.patch import CellCenterData2d
 
 
@@ -13,14 +15,18 @@ class FV2d(CellCenterData2d):
 
     """
 
-    def to_centers(self, name):
+    def to_centers(self, name, is_positive=False):
         """ convert variable name from an average to cell-centers """
 
         a = self.get_var(name)
         c = self.grid.scratch_array()
         ng = self.grid.ng
         c[:, :] = a[:, :]
         c.v(buf=ng-1)[:, :] = a.v(buf=ng-1) - self.grid.dx**2*a.lap(buf=ng-1)/24.0
+
+        if is_positive:
+            c.v(buf=ng-1)[:, :] = np.where(c.v(buf=ng-1) >= 0.0, c.v(buf=ng-1), a.v(buf=ng-1))
+
         return c
 
     def from_centers(self, name):

pyro/mesh/patch.py

@@ -146,15 +146,15 @@ def __init__(self, nx, ny, *, ng=1,
         self.xr2d = ArrayIndexer(d=xr2d, grid=self)
         self.yr2d = ArrayIndexer(d=yr2d, grid=self)
 
-    def scratch_array(self, *, nvar=1):
+    def scratch_array(self, *, nvar=1, dtype=np.float64):
         """
         return a standard numpy array dimensioned to have the size
         and number of ghostcells as the parent grid
         """
         if nvar == 1:
-            _tmp = np.zeros((self.qx, self.qy), dtype=np.float64)
+            _tmp = np.zeros((self.qx, self.qy), dtype=dtype)
         else:
-            _tmp = np.zeros((self.qx, self.qy, nvar), dtype=np.float64)
+            _tmp = np.zeros((self.qx, self.qy, nvar), dtype=dtype)
         return ArrayIndexer(d=_tmp, grid=self)
 
     def coarse_like(self, N):