[wip] Add small e protection

pyro/compressible/simulation.py

@@ -68,6 +68,13 @@ def cons_to_prim(U, gamma, ivars, myg):
     e_min = e.v().min()
     rho_min = q.v(n=ivars.irho).min()
 
+    if e_min < 0:
+        eidx = np.asarray(e < 0).nonzero()
+        i_idx = eidx[0]
+        j_jdx = eidx[1]
+        for i, j in zip(i_idx, j_jdx):
+            print(f" e < 0: {i}, {j}, {e[i, j]}")
+
     assert e_min > 0.0 and rho_min > 0.0, f"invalid state, min(rho) = {rho_min}, min(e) = {e_min}"
 
     q[:, :, ivars.ip] = eos.pres(gamma, q[:, :, ivars.irho], e)
@@ -452,8 +459,19 @@ def evolve(self):
     def clean_state(self, U):
         """enforce minimum density and eint on the conserved state U"""
 
-        U.v(n=self.ivars.idens)[:, :] = np.maximum(U.v(n=self.ivars.idens),
-                                                   self.rp.get_param("compressible.small_dens"))
+        U[..., self.ivars.idens] = np.maximum(U[..., self.ivars.idens],
+                                              self.rp.get_param("compressible.small_dens"))
+
+        if self.small_eint > 0:
+
+            ekin = 0.5 * (U[..., self.ivars.ixmom]**2 +
+                          U[..., self.ivars.iymom]**2) / U[..., self.ivars.idens]
+
+            rhoe = U[..., self.ivars.iener] - ekin
+
+            U[..., self.ivars.iener] = np.where(rhoe < self.small_eint,
+                                                U[..., self.ivars.idens] * self.small_eint + ekin,
+                                                U[..., self.ivars.iener])
 
     def dovis(self):
         """

pyro/compressible_fv4/_defaults

@@ -24,6 +24,7 @@ grav = 0.0                ; gravitational acceleration (in y-direction)
 riemann = CGF
 
 small_dens = -1.e200      ; minimum allowed density
+small_eint_factor = -1.e200  ; multiplicative factor on the initial minimum eint to use for limiting e
 
 [sponge]
 do_sponge = 0             ; do we include a sponge source term

pyro/compressible_fv4/simulation.py

@@ -1,6 +1,6 @@
 import pyro.compressible_fv4.fluxes as flx
 from pyro import compressible_rk
-from pyro.compressible import get_external_sources, get_sponge_factor
+from pyro.compressible import cons_to_prim, get_external_sources, get_sponge_factor
 from pyro.mesh import fv
 
 
@@ -73,3 +73,12 @@ def preevolve(self):
         # we just initialized cell-centers, but we need to store averages
         for var in self.cc_data.names:
             self.cc_data.from_centers(var)
+
+        efac = self.rp.get_param("compressible.small_eint_factor")
+
+        print(type(self.cc_data))
+
+        q = cons_to_prim(self.cc_data.data, self.rp.get_param("eos.gamma"), self.ivars, self.cc_data.grid)
+
+        self.small_eint = efac * q.v(n=self.ivars.ip).min()
+        print("small_eint = ", self.small_eint)

pyro/compressible_sdc/_defaults

@@ -21,6 +21,8 @@ riemann = CGF
 
 small_dens = -1.e200      ; minimum allowed density
 
+small_eint_factor = -1.e200  ; multiplicative factor on the initial minimum eint to use for limiting e
+
 [sponge]
 do_sponge = 0             ; do we include a sponge source term