From c2ecf50590221b08a2e1e15567b3781eb9892b5a Mon Sep 17 00:00:00 2001
From: Michael Zingale <michael.zingale@stonybrook.edu>
Date: Sat, 25 Jan 2025 12:01:49 -0500
Subject: [PATCH] implement a clean_state routine for compressible (#329)

this will enforce small_dens
---
 pyro/compressible/_defaults         | 2 ++
 pyro/compressible/simulation.py     | 8 ++++++++
 pyro/compressible_fv4/_defaults     | 1 +
 pyro/compressible_fv4/simulation.py | 2 ++
 pyro/compressible_rk/_defaults      | 1 +
 pyro/compressible_rk/simulation.py  | 2 ++
 pyro/compressible_sdc/_defaults     | 1 +
 7 files changed, 17 insertions(+)

diff --git a/pyro/compressible/_defaults b/pyro/compressible/_defaults
index 850d3d729..37c2c9b0a 100644
--- a/pyro/compressible/_defaults
+++ b/pyro/compressible/_defaults
@@ -21,6 +21,8 @@ grav = 0.0                ; gravitational acceleration (in y-direction)
 
 riemann = HLLC            ; HLLC or CGF
 
+small_dens = -1.e200      ; minimum allowed density
+small_eint = -1.e200      ; minimum allowed specific internal energy
 
 [sponge]
 do_sponge = 0             ; do we include a sponge source term
diff --git a/pyro/compressible/simulation.py b/pyro/compressible/simulation.py
index 645d8277a..5f03f68e4 100644
--- a/pyro/compressible/simulation.py
+++ b/pyro/compressible/simulation.py
@@ -293,6 +293,8 @@ def evolve(self):
         timestep dt.
         """
 
+        self.clean_state(self.cc_data.data)
+
         tm_evolve = self.tc.timer("evolve")
         tm_evolve.begin()
 
@@ -437,6 +439,12 @@ def evolve(self):
 
         tm_evolve.end()
 
+    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"))
+
     def dovis(self):
         """
         Do runtime visualization.
diff --git a/pyro/compressible_fv4/_defaults b/pyro/compressible_fv4/_defaults
index 18d1c3707..ae1c2979c 100644
--- a/pyro/compressible_fv4/_defaults
+++ b/pyro/compressible_fv4/_defaults
@@ -23,6 +23,7 @@ grav = 0.0                ; gravitational acceleration (in y-direction)
 
 riemann = CGF
 
+small_dens = -1.e200      ; minimum allowed density
 
 [sponge]
 do_sponge = 0             ; do we include a sponge source term
diff --git a/pyro/compressible_fv4/simulation.py b/pyro/compressible_fv4/simulation.py
index 11c9a4aa5..2fff3541d 100644
--- a/pyro/compressible_fv4/simulation.py
+++ b/pyro/compressible_fv4/simulation.py
@@ -19,6 +19,8 @@ def substep(self, myd):
         compute the advective source term for the given state
         """
 
+        self.clean_state(myd.data)
+
         myg = myd.grid
 
         # compute the source terms -- we need to do this first
diff --git a/pyro/compressible_rk/_defaults b/pyro/compressible_rk/_defaults
index 8f1aa0467..ea0aa9814 100644
--- a/pyro/compressible_rk/_defaults
+++ b/pyro/compressible_rk/_defaults
@@ -25,6 +25,7 @@ riemann = HLLC            ; HLLC or CGF
 
 well_balanced = 0         ; use a well-balanced scheme to keep the model in hydrostatic equilibrium
 
+small_dens = -1.e200      ; minimum allowed density
 
 [sponge]
 do_sponge = 0             ; do we include a sponge source term
diff --git a/pyro/compressible_rk/simulation.py b/pyro/compressible_rk/simulation.py
index 7378fd478..2a28b931e 100644
--- a/pyro/compressible_rk/simulation.py
+++ b/pyro/compressible_rk/simulation.py
@@ -15,6 +15,8 @@ def substep(self, myd):
         conservative state defined as part of myd
         """
 
+        self.clean_state(myd.data)
+
         myg = myd.grid
 
         # source terms -- note: this dt is the entire dt, not the
diff --git a/pyro/compressible_sdc/_defaults b/pyro/compressible_sdc/_defaults
index 00740ac78..11568dd35 100644
--- a/pyro/compressible_sdc/_defaults
+++ b/pyro/compressible_sdc/_defaults
@@ -19,6 +19,7 @@ grav = 0.0                ; gravitational acceleration (in y-direction)
 
 riemann = CGF
 
+small_dens = -1.e200      ; minimum allowed density
 
 [sponge]
 do_sponge = 0             ; do we include a sponge source term