implement a sponge term in compressible (#313)

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

@@ -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)

pyro/compressible/_defaults

@@ -21,6 +21,14 @@ grav = 0.0                ; gravitational acceleration (in y-direction)
 
 riemann = HLLC            ; HLLC or 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
+
 [particles]
 do_particles = 0
 particle_generator = grid

pyro/compressible/simulation.py

@@ -159,6 +159,29 @@ def get_external_sources(t, dt, U, ivars, rp, myg, *, U_old=None, problem_source
     return S
 
 
+def get_sponge_factor(U, ivars, rp, myg):
+    """compute the sponge factor, f / tau, that goes into a
+    sponge damping term of the form S = - (f / tau) (rho U)"""
+
+    rho = U[:, :, ivars.idens]
+    rho_begin = rp.get_param("sponge.sponge_rho_begin")
+    rho_full = rp.get_param("sponge.sponge_rho_full")
+
+    assert rho_begin > rho_full
+
+    f = myg.scratch_array()
+
+    f[:, :] = np.where(rho > rho_begin,
+                       0.0,
+                       np.where(rho < rho_full,
+                                1.0,
+                                0.5 * (1.0 - np.cos(np.pi * (rho - rho_begin) /
+                                                    (rho_full - rho_begin)))))
+
+    tau = rp.get_param("sponge.sponge_timescale")
+    return f / tau
+
+
 class Simulation(NullSimulation):
     """The main simulation class for the corner transport upwind
     compressible hydrodynamics solver
@@ -395,6 +418,14 @@ def evolve(self):
             var = self.cc_data.get_var_by_index(n)
             var.v()[:, :] += 0.5 * self.dt * (S_new.v(n=n) - S_old.v(n=n))
 
+        # finally, do the sponge, if desired -- this is formulated as an
+        # implicit update to the velocity
+        if self.rp.get_param("sponge.do_sponge"):
+            kappa_f = get_sponge_factor(self.cc_data.data, self.ivars, self.rp, myg)
+
+            self.cc_data.data[:, :, self.ivars.ixmom] /= (1.0 + self.dt * kappa_f)
+            self.cc_data.data[:, :, self.ivars.iymom] /= (1.0 + self.dt * kappa_f)
+
         if self.particles is not None:
             self.particles.update_particles(self.dt)
 

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
+
+
 

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):

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

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):

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