From 4a074a698b51e531e2b9c21365e337297c347b26 Mon Sep 17 00:00:00 2001
From: Michael Zingale <michael.zingale@stonybrook.edu>
Date: Thu, 9 Jan 2025 15:49:30 -0500
Subject: [PATCH] add a sponge to the compressible solvers for the CTU solver,
 the sponge update is implicit in time for the RK/SDC solvers, it is an
 explicit source this damps the velocity field over a timescale
 "sponge_timescale" The sponge is essentially the same as in Castro

---
 pyro/compressible/__init__.py       | 2 +-
 pyro/compressible/_defaults         | 1 +
 pyro/compressible_fv4/_defaults     | 8 ++++++++
 pyro/compressible_fv4/simulation.py | 9 ++++++++-
 pyro/compressible_rk/_defaults      | 6 ++++++
 pyro/compressible_rk/simulation.py  | 7 +++++++
 pyro/compressible_sdc/_defaults     | 8 ++++++++
 7 files changed, 39 insertions(+), 2 deletions(-)

diff --git a/pyro/compressible/__init__.py b/pyro/compressible/__init__.py
index d78b69b9b..9aacc0ca0 100644
--- a/pyro/compressible/__init__.py
+++ b/pyro/compressible/__init__.py
@@ -7,4 +7,4 @@
 __all__ = ["simulation"]
 
 from .simulation import (Simulation, Variables, cons_to_prim,
-                         get_external_sources, prim_to_cons)
+                         get_external_sources, get_sponge_factor, prim_to_cons)
diff --git a/pyro/compressible/_defaults b/pyro/compressible/_defaults
index 148426b8b..850d3d729 100644
--- a/pyro/compressible/_defaults
+++ b/pyro/compressible/_defaults
@@ -21,6 +21,7 @@ grav = 0.0                ; gravitational acceleration (in y-direction)
 
 riemann = HLLC            ; HLLC or CGF
 
+
 [sponge]
 do_sponge = 0             ; do we include a sponge source term
 
diff --git a/pyro/compressible_fv4/_defaults b/pyro/compressible_fv4/_defaults
index c9f936e0d..18d1c3707 100644
--- a/pyro/compressible_fv4/_defaults
+++ b/pyro/compressible_fv4/_defaults
@@ -24,4 +24,12 @@ grav = 0.0                ; gravitational acceleration (in y-direction)
 riemann = CGF
 
 
+[sponge]
+do_sponge = 0             ; do we include a sponge source term
+
+sponge_rho_begin = 1.e-2         ; density below which to begin the sponge
+sponge_rho_full = 1.e-3          ; density below which the sponge is fully enabled
+sponge_timescale = 1.e-2         ; the timescale over which the sponge should act
+
+
 
diff --git a/pyro/compressible_fv4/simulation.py b/pyro/compressible_fv4/simulation.py
index df662ae4a..0097832b6 100644
--- a/pyro/compressible_fv4/simulation.py
+++ b/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
+from pyro.compressible import get_external_sources, get_sponge_factor
 from pyro.mesh import fv, integration
 
 
@@ -48,6 +48,13 @@ def substep(self, myd):
                (flux_x.v(n=n) - flux_x.ip(1, n=n))/myg.dx + \
                (flux_y.v(n=n) - flux_y.jp(1, n=n))/myg.dy + S.v(n=n)
 
+        # finally, add the sponge source, if desired
+        if self.rp.get_param("sponge.do_sponge"):
+            kappa_f = get_sponge_factor(myd.data, self.ivars, self.rp, myg)
+
+            k.v(n=self.ivars.ixmom)[:, :] -= kappa_f.v() * myd.data.v(n=self.ivars.ixmom)
+            k.v(n=self.ivars.iymom)[:, :] -= kappa_f.v() * myd.data.v(n=self.ivars.iymom)
+
         return k
 
     def preevolve(self):
diff --git a/pyro/compressible_rk/_defaults b/pyro/compressible_rk/_defaults
index 293a4ec6b..8f1aa0467 100644
--- a/pyro/compressible_rk/_defaults
+++ b/pyro/compressible_rk/_defaults
@@ -26,3 +26,9 @@ riemann = HLLC            ; HLLC or CGF
 well_balanced = 0         ; use a well-balanced scheme to keep the model in hydrostatic equilibrium
 
 
+[sponge]
+do_sponge = 0             ; do we include a sponge source term
+
+sponge_rho_begin = 1.e-2         ; density below which to begin the sponge
+sponge_rho_full = 1.e-3          ; density below which the sponge is fully enabled
+sponge_timescale = 1.e-2         ; the timescale over which the sponge should act
diff --git a/pyro/compressible_rk/simulation.py b/pyro/compressible_rk/simulation.py
index f0306403a..7378fd478 100644
--- a/pyro/compressible_rk/simulation.py
+++ b/pyro/compressible_rk/simulation.py
@@ -33,6 +33,13 @@ def substep(self, myd):
                (flux_x.v(n=n) - flux_x.ip(1, n=n))/myg.dx + \
                (flux_y.v(n=n) - flux_y.jp(1, n=n))/myg.dy + S.v(n=n)
 
+        # finally, add the sponge source, if desired
+        if self.rp.get_param("sponge.do_sponge"):
+            kappa_f = compressible.get_sponge_factor(myd.data, self.ivars, self.rp, myg)
+
+            k.v(n=self.ivars.ixmom)[:, :] -= kappa_f.v() * myd.data.v(n=self.ivars.ixmom)
+            k.v(n=self.ivars.iymom)[:, :] -= kappa_f.v() * myd.data.v(n=self.ivars.iymom)
+
         return k
 
     def method_compute_timestep(self):
diff --git a/pyro/compressible_sdc/_defaults b/pyro/compressible_sdc/_defaults
index 438576e9d..1fb4e9e13 100644
--- a/pyro/compressible_sdc/_defaults
+++ b/pyro/compressible_sdc/_defaults
@@ -22,3 +22,11 @@ temporal_method = RK4     ; integration method (see mesh/integration.py)
 grav = 0.0                ; gravitational acceleration (in y-direction)
 
 riemann = CGF
+
+
+[sponge]
+do_sponge = 0             ; do we include a sponge source term
+
+sponge_rho_begin = 1.e-2         ; density below which to begin the sponge
+sponge_rho_full = 1.e-3          ; density below which the sponge is fully enabled
+sponge_timescale = 1.e-2         ; the timescale over which the sponge should act