[WIP] Setting the PPM advection problem

docs/source/advection_basics.rst

@@ -19,7 +19,7 @@ with different spatial and temporal integration schemes.
 :py:mod:`pyro.advection` implements the directionally unsplit corner
 transport upwind algorithm :cite:`colella:1990` with piecewise linear reconstruction.
 This is an overall second-order accurate method, with timesteps restricted
-by 
+by
 
 .. math::
 
@@ -29,6 +29,23 @@ The parameters for this solver are:
 
 .. include:: advection_defaults.inc
 
+``advection_ppm`` solver
+------------------------
+
+:py:mod:`pyro.advection_ppm` applies a piecewise parabolic reconstruction (PPM) in the directionally
+unsplit corner transport upwind algorithm (CTU) :cite:`colella:1984`. This reconstruction is an extension
+of a higher-order Gudunov's method, designed to capture steeper representation of discontinuities, with
+emphasis on contact discontinuities.
+
+This is an overall second-order accurate method, with timesteps restricted by:
+
+.. math::
+
+  \Delta t < \min \left \{ \frac{\Delta x}{|u|}, \frac{\Delta y}{|v|} \right \}
+
+The parameters for this solver are:
+
+.. include:: advection_ppm_defaults.inc
 
 ``advection_fv4`` solver
 ------------------------
@@ -39,8 +56,8 @@ method with RK4 time integration, following the ideas in
 method-of-lines integration, and as such has a slightly more restrictive
 timestep:
 
-.. math:: 
- 
+.. math::
+
   \Delta t \lesssim \left [ \frac{|u|}{\Delta x} + \frac{|v|}{\Delta y} \right ]^{-1}
 
 The main complexity comes from needing to average the flux over the

docs/source/advection_ppm_defaults.inc

@@ -0,0 +1,27 @@
+* section: [advection]
+
+  +----------------------------------+----------------+----------------------------------------------------+
+  | option                           | value          | description                                        |
+  +==================================+================+====================================================+
+  | u                                | ``1.0``        | advective velocity in x-direction                  |
+  +----------------------------------+----------------+----------------------------------------------------+
+  | v                                | ``1.0``        | advective velocity in y-direction                  |
+  +----------------------------------+----------------+----------------------------------------------------+
+
+* section: [driver]
+
+  +----------------------------------+----------------+----------------------------------------------------+
+  | option                           | value          | description                                        |
+  +==================================+================+====================================================+
+  | cfl                              | ``0.8``        | advective CFL number                               |
+  +----------------------------------+----------------+----------------------------------------------------+
+
+* section: [particles]
+
+  +----------------------------------+----------------+----------------------------------------------------+
+  | option                           | value          | description                                        |
+  +==================================+================+====================================================+
+  | do_particles                     | ``0``          |                                                    |
+  +----------------------------------+----------------+----------------------------------------------------+
+  | particle_generator               | ``grid``       |                                                    |
+  +----------------------------------+----------------+----------------------------------------------------+

docs/source/refs.bib

@@ -12,6 +12,19 @@ @ARTICLE{colella:1990
   adsnote = {Provided by the SAO/NASA Astrophysics Data System}
 }
 
+@ARTICLE{colella:1984,
+  title = {The Piecewise Parabolic Method (PPM) for gas-dynamical simulations},
+  journal = {Journal of Computational Physics},
+  volume = {54},
+  number = {1},
+  pages = {174-201},
+  year = {1984},
+  issn = {0021-9991},
+  doi = {https://doi.org/10.1016/0021-9991(84)90143-8},
+  url = {https://www.sciencedirect.com/science/article/pii/0021999184901438},
+  author = {Phillip Colella and Paul R Woodward},
+}
+
 @ARTICLE{mccorquodalecolella,
    author = {{McCorquodale}, P. and {Colella}, P.},
   title = "{A high-order finite-volume method for conservation laws on

pyro/advection_ppm/__init__.py

@@ -0,0 +1,10 @@
+"""
+The Piecewise Parabolic reconstruction Method (PPM) + CTU unsplit dimensional implementation. This problem is a second-order
+advection solver, based on the Colella & Woodward (1984), Colella (1990), Colella & Miller (2002) and Colella & Sekora (2002)
+papers descriptions. The advantage of choosing the PPM parabolic reconstruction instead of the Piecewise Linear Reconstruction
+(PLR), is the increased resolution of interfaces at contact discontinuities.
+
+"""
+
+__all__ = ['simulation']
+from .simulation import Simulation

pyro/advection_ppm/_defaults

@@ -0,0 +1,14 @@
+[driver]
+cfl       = 0.8           ; advective CFL number
+
+
+[advection]
+u = 1.0      ; advective velocity in x-direction
+v = 1.0      ; advective velocity in y-direction
+
+limiter = 1  ; limiter (0 = none, 1 = 2nd order, 2 = 4th order)
+
+
+[particles]
+do_particles = 0
+particle_generator = grid

pyro/advection_ppm/interface.py

@@ -0,0 +1,59 @@
+from pyro.mesh.reconstruction import ppm_reconstruction
+
+
+def ppm_interface(a, myg, rp, dt):
+
+    # get the advection velocities
+    u = rp.get_param("advection.u")
+    v = rp.get_param("advection.v")
+
+    limiter = rp.get_param("advection.limiter")
+
+    cx = u*dt/myg.dx
+    cy = v*dt/myg.dy
+
+    # --------------------------------------------------------------------------
+    # monotonized central differences
+    # --------------------------------------------------------------------------
+
+    delta_ax = myg.scratch_array()
+    a6x = myg.scratch_array()
+    delta_ay = myg.scratch_array()
+    a6y = myg.scratch_array()
+
+    # upwind
+    a_x = myg.scratch_array()
+
+    ar, al = ppm_reconstruction(a, myg, idir=1, limiter=limiter)
+    delta_ax.v(buf=1)[:, :] = al.v(buf=1) - ar.v(buf=1)
+    a6x.v(buf=1)[:, :] = 6.0 * (a.v(buf=1) - 0.5*(ar.v(buf=1) + al.v(buf=1)))
+
+    if u < 0:
+        cx = -cx
+        # a_x[i,j] = a[i,j] - 0.5*(1.0 + cx)*ldelta_a[i,j]
+        a_x.v(buf=1)[:, :] = ar.v(buf=1) + 0.5 * cx * (delta_ax.v(buf=1) +
+                                (1 - 2.0 * cx / 3.0) * a6x.v(buf=1))
+    else:
+        # a_x[i,j] = a[i-1,j] + 0.5*(1.0 - cx)*ldelta_a[i-1,j]
+        a_x.v(buf=1)[:, :] = al.ip(-1, buf=1) - 0.5 * cx * (delta_ax.ip(-1, buf=1) -
+                                (1 - 2.0 * cx / 3.0) * a6x.ip(-1, buf=1))
+
+    # y-direction
+    a_y = myg.scratch_array()
+
+    ar, al = ppm_reconstruction(a, myg, idir=2, limiter=limiter)
+    delta_ay.v(buf=1)[:, :] = al.v(buf=1) - ar.v(buf=1)
+    a6y.v(buf=1)[:, :] = 6.0 * (a.v(buf=1) - 0.5*(ar.v(buf=1) + al.v(buf=1)))
+
+    # upwind
+    if v < 0:
+        cy = -cy
+        # a_y[i,j] = a[i,j] - 0.5*(1.0 + cy)*ldelta_a[i,j]
+        a_y.v(buf=1)[:, :] = ar.v(buf=1) + 0.5 * cy * (delta_ay.v(buf=1) +
+                                (1 - 2.0 * cy / 3.0) * a6y.v(buf=1))
+    else:
+        # a_y[i,j] = a[i,j-1] + 0.5*(1.0 - cy)*ldelta_a[i,j-1]
+        a_y.v(buf=1)[:, :] = al.jp(-1, buf=1) - 0.5 * cy * (delta_ay.jp(-1, buf=1) -
+                                (1 - 2.0 * cy / 3.0) * a6y.jp(-1, buf=1))
+
+    return u, v, a_x, a_y

pyro/advection_ppm/problems/__init__.py

@@ -0,0 +1 @@
+__all__ = ['smooth', 'tophat']

pyro/advection_ppm/problems/inputs.smooth

@@ -0,0 +1 @@
+../../advection/problems/inputs.smooth

pyro/advection_ppm/problems/inputs.tophat

@@ -0,0 +1 @@
+../../advection/problems/inputs.tophat

pyro/advection_ppm/problems/smooth.py

@@ -0,0 +1 @@
+../../advection/problems/smooth.py

pyro/advection_ppm/problems/tophat.py

@@ -0,0 +1 @@
+../../advection/problems/tophat.py

pyro/advection_ppm/simulation.py

@@ -0,0 +1,42 @@
+import pyro.advection.advective_fluxes as flx
+from pyro import advection
+from pyro.advection_ppm.interface import ppm_interface
+
+
+class Simulation(advection.Simulation):
+
+    def evolve(self):
+        """
+        Evolve the linear advection equation through one timestep.  We only
+        consider the "density" variable in the CellCenterData2d object that
+        is part of the Simulation.
+        """
+
+        tm_evolve = self.tc.timer("evolve")
+        tm_evolve.begin()
+
+        dtdx = self.dt/self.cc_data.grid.dx
+        dtdy = self.dt/self.cc_data.grid.dy
+
+        Fx, Fy = flx.unsplit_fluxes(self.cc_data, self.rp, self.dt, "density", ppm_interface)
+
+        dens = self.cc_data.get_var("density")
+
+        dens.v()[:, :] = dens.v() + dtdx*(Fx.v() - Fx.ip(1)) + \
+                                    dtdy*(Fy.v() - Fy.jp(1))
+
+        if self.particles is not None:
+            myg = self.cc_data.grid
+            u = self.rp.get_param("advection.u")
+            v = self.rp.get_param("advection.v")
+
+            u2d = myg.scratch_array() + u
+            v2d = myg.scratch_array() + v
+
+            self.particles.update_particles(self.dt, u2d, v2d)
+
+        # increment the time
+        self.cc_data.t += self.dt
+        self.n += 1
+
+        tm_evolve.end()

pyro/advection_ppm/tests/test_advection_ppm.py

@@ -0,0 +1,37 @@
+import pyro.advection_nonuniform.simulation as sn
+from pyro.util import runparams
+
+
+class TestSimulation:
+    @classmethod
+    def setup_class(cls):
+        """ this is run once for each class before any tests """
+
+    @classmethod
+    def teardown_class(cls):
+        """ this is run once for each class after all tests """
+
+    def setup_method(self):
+        """ this is run before each test """
+        self.rp = runparams.RuntimeParameters()
+
+        self.rp.params["mesh.nx"] = 8
+        self.rp.params["mesh.ny"] = 8
+        self.rp.params["particles.do_particles"] = 0
+
+        self.sim = sn.Simulation("advection_ppm", "test", self.rp)
+        self.sim.initialize()
+
+    def teardown_method(self):
+        """ this is run after each test """
+        self.rp = None
+        self.sim = None
+
+    def test_initializationst(self):
+        dens = self.sim.cc_data.get_var("density")
+        u = self.sim.cc_data.get_var("x-velocity")
+        v = self.sim.cc_data.get_var("y-velocity")
+
+        assert dens.min() == 1.0 and dens.max() == 1.0
+        assert u.min() == 1.0 and u.max() == 1.0
+        assert v.min() == 1.0 and v.max() == 1.0