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 27edf68e4..850d3d729 100644 --- a/pyro/compressible/_defaults +++ b/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 diff --git a/pyro/compressible/simulation.py b/pyro/compressible/simulation.py index 993f10b8f..dd9e65e71 100644 --- a/pyro/compressible/simulation.py +++ b/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) 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