add some protection against negatives in 4th order

pyro/compressible_fv4/fluxes.py

@@ -52,10 +52,10 @@ def fluxes(myd, rp, ivars):
     # convert U from cell-averages to cell-centers
     U_cc = np.zeros_like(U_avg)
 
-    U_cc[:, :, ivars.idens] = myd.to_centers("density")
+    U_cc[:, :, ivars.idens] = myd.to_centers("density", is_positive=True)
     U_cc[:, :, ivars.ixmom] = myd.to_centers("x-momentum")
     U_cc[:, :, ivars.iymom] = myd.to_centers("y-momentum")
-    U_cc[:, :, ivars.iener] = myd.to_centers("energy")
+    U_cc[:, :, ivars.iener] = myd.to_centers("energy", is_positive=True)
 
     # compute the primitive variables of both the cell-center and averages
     q_bar = comp.cons_to_prim(U_avg, gamma, ivars, myd.grid)
@@ -66,6 +66,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/pyro_sim.py

@@ -30,6 +30,8 @@
                  "lm_atm",
                  "swe"]
 
+import warnings
+warnings.filterwarnings("error")
 
 class Pyro:
     """