Merge branch 'main' into one-sided-fourth-order

docs/source/index.rst

@@ -72,6 +72,7 @@ new ideas.
 
    compressible-rt-compare.ipynb
    adding_a_problem_jupyter.ipynb
+   rt_average.ipynb
    advection-error.ipynb
    compressible-convergence.ipynb
 

docs/source/output.rst

@@ -42,20 +42,41 @@ Reading and plotting manually
 pyro output data can be read using the :func:`util.io_pyro.read <pyro.util.io_pyro.read>` method. The following
 sequence (done in a python session) reads in stored data (from the
 compressible Sedov problem) and plots data falling on a line in the x
-direction through the y-center of the domain (note: this will include
-the ghost cells).
+direction through the y-center of the domain.  The return value of
+``read`` is a ``Simulation`` object.
+
 
 .. code-block:: python
 
    import matplotlib.pyplot as plt
    import pyro.util.io_pyro as io
-   sim = io.read("sedov_unsplit_0000.h5")
+
+   sim = io.read("sedov_unsplit_0290.h5")
    dens = sim.cc_data.get_var("density")
-   plt.plot(dens.g.x, dens[:,dens.g.ny//2])
-   plt.show()
+
+   fig, ax = plt.subplots()
+   ax.plot(dens.g.x, dens[:,dens.g.qy//2])
+   ax.grid()
 
 .. image:: manual_plot.png
    :align: center
 
-Note: this includes the ghost cells, by default, seen as the small
-regions of zeros on the left and right.
+.. note::
+
+   This includes the ghost cells, by default, seen as the small
+   regions of zeros on the left and right.  The total number of cells,
+   including ghost cells in the y-direction is ``qy``, which is why
+   we use that in our slice.
+
+If we wanted to exclude the ghost cells, then we could use the ``.v()`` method
+on the density array to exclude the ghost cells, and then manually index ``g.x``
+to just include the valid part of the domain:
+
+.. code:: python
+
+   ax.plot(dens.g.x[g.ilo:g.ihi+1], dens.v()[:, dens.g.ny//2])
+
+.. note::
+
+   In this case, we are using ``ny`` since that is the width of the domain
+   excluding ghost cells.

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/problems/convection.py

@@ -42,6 +42,9 @@ def init_data(my_data, rp):
     ymom[:, :] = 0.0
     dens[:, :] = dens_cutoff
 
+    # create a seeded random number generator
+    rng = np.random.default_rng(12345)
+
     # set the density to be stratified in the y-direction
     myg = my_data.grid
 
@@ -75,7 +78,7 @@ def init_data(my_data, rp):
     ener[:, :] = p[:, :]/(gamma - 1.0)
 
     # pairs of random numbers between [-1, 1]
-    vel_pert = 2.0 * np.random.random_sample((myg.qx, myg.qy, 2)) - 1
+    vel_pert = 2.0 * rng.random(size=(myg.qx, myg.qy, 2)) - 1
 
     cs = np.sqrt(gamma * p / dens)
 

pyro/compressible/problems/inputs.rt_multimode

@@ -0,0 +1,33 @@
+# simple inputs files for the four-corner problem.
+
+[driver]
+max_steps = 10000
+tmax = 3.0
+
+
+[io]
+basename = rt_
+n_out = 100
+
+
+[mesh]
+nx = 64
+ny = 192
+xmax = 1.0
+ymax = 3.0
+
+xlboundary = periodic
+xrboundary = periodic
+
+ylboundary = hse
+yrboundary = hse
+
+
+[rt_multimode]
+amp = 0.25
+nmodes = 12
+
+[compressible]
+grav = -1.0
+
+limiter = 2

pyro/compressible/problems/rt_multimode.py

@@ -0,0 +1,85 @@
+"""A multi-mode Rayleigh-Taylor instability."""
+
+import numpy as np
+
+from pyro.util import msg
+
+DEFAULT_INPUTS = "inputs.rt_multimode"
+
+PROBLEM_PARAMS = {"rt_multimode.dens1": 1.0,
+                  "rt_multimode.dens2": 2.0,
+                  "rt_multimode.amp": 1.0,
+                  "rt_multimode.sigma": 0.1,
+                  "rt_multimode.nmodes": 10,
+                  "rt_multimode.p0": 10.0}
+
+
+def init_data(my_data, rp):
+    """ initialize the rt problem """
+
+    # see the random number generator
+    rng = np.random.default_rng(12345)
+
+    if rp.get_param("driver.verbose"):
+        msg.bold("initializing the rt problem...")
+
+    # get the density, momenta, and energy as separate variables
+    dens = my_data.get_var("density")
+    xmom = my_data.get_var("x-momentum")
+    ymom = my_data.get_var("y-momentum")
+    ener = my_data.get_var("energy")
+
+    gamma = rp.get_param("eos.gamma")
+
+    grav = rp.get_param("compressible.grav")
+
+    dens1 = rp.get_param("rt_multimode.dens1")
+    dens2 = rp.get_param("rt_multimode.dens2")
+    p0 = rp.get_param("rt_multimode.p0")
+    amp = rp.get_param("rt_multimode.amp")
+    sigma = rp.get_param("rt_multimode.sigma")
+    nmodes = rp.get_param("rt_multimode.nmodes")
+
+    # initialize the components, remember, that ener here is
+    # rho*eint + 0.5*rho*v**2, where eint is the specific
+    # internal energy (erg/g)
+    xmom[:, :] = 0.0
+    ymom[:, :] = 0.0
+    dens[:, :] = 0.0
+
+    # set the density to be stratified in the y-direction
+    myg = my_data.grid
+
+    ycenter = 0.5*(myg.ymin + myg.ymax)
+
+    p = myg.scratch_array()
+
+    j = myg.jlo
+    while j <= myg.jhi:
+        if myg.y[j] < ycenter:
+            dens[:, j] = dens1
+            p[:, j] = p0 + dens1*grav*myg.y[j]
+
+        else:
+            dens[:, j] = dens2
+            p[:, j] = p0 + dens1*grav*ycenter + dens2*grav*(myg.y[j] - ycenter)
+
+        j += 1
+
+    # add multiple modes to the vertical velocity
+    L = myg.xmax - myg.xmin
+    for k in range(1, nmodes+1):
+        phase = rng.random() * 2 * np.pi
+        mode_amp = amp * rng.random()
+        ymom[:, :] += (mode_amp * np.cos(2.0 * np.pi * k*myg.x2d / L + phase) *
+                       np.exp(-(myg.y2d - ycenter)**2 / sigma**2))
+    ymom /= nmodes
+    ymom *= dens
+
+    # set the energy (P = cs2*dens)
+    ener[:, :] = p[:, :]/(gamma - 1.0) + \
+        0.5*(xmom[:, :]**2 + ymom[:, :]**2)/dens[:, :]
+
+
+def finalize():
+    """ print out any information to the user at the end of the run """

pyro/compressible/simulation.py

@@ -65,6 +65,9 @@ def cons_to_prim(U, gamma, ivars, myg):
                   out=np.zeros_like(U[:, :, ivars.iener]),
                   where=(U[:, :, ivars.idens] != 0.0))
 
+    assert e.v().min() > 0.0
+    assert q.v(n=ivars.irho).min() > 0.0
+
     q[:, :, ivars.ip] = eos.pres(gamma, q[:, :, ivars.irho], e)
 
     if ivars.naux > 0:
@@ -156,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
@@ -392,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