Add SSPRK43 implemented by THG (#863)

Co-authored-by: Mike Campbell <mtcampbe@illinois.edu>

mirgecom/integrators/__init__.py

@@ -26,6 +26,7 @@
 
 from .explicit_rk import rk4_step                          # noqa: F401
 from .lsrk import euler_step, lsrk54_step, lsrk144_step    # noqa: F401
+from .ssprk import ssprk43_step                            # noqa: F401
 
 __doc__ = """
 .. automodule:: mirgecom.integrators.explicit_rk

mirgecom/integrators/ssprk.py

@@ -0,0 +1,40 @@
+"""Timestepping routines for strong-stability preserving Runge-Kutta methods.
+
+.. autofunction:: ssprk43_step
+"""
+
+__copyright__ = """
+Copyright (C) 2022 University of Illinois Board of Trustees
+"""
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+
+def ssprk43_step(state, t, dt, rhs):
+    """Take one step using an explicit 4-stage, 3rd-order, SSPRK method."""
+
+    def rhs_update(t, y):
+        return y + dt*rhs(t, y)
+
+    y1 = 1/2*state + 1/2*rhs_update(t, state)
+    y2 = 1/2*y1 + 1/2*rhs_update(t + dt/2, y1)
+    y3 = 2/3*state + 1/6*y2 + 1/6*rhs_update(t + dt, y2)
+    y4 = 1/2*y3 + 1/2*rhs_update(t + dt/2, y3)
+
+    return y4

test/test_time_integrators.py

@@ -32,10 +32,10 @@
     pytest_generate_tests_for_pyopencl_array_context
     as pytest_generate_tests)
 
-from mirgecom.integrators import (euler_step,
-                                  lsrk54_step,
-                                  lsrk144_step,
-                                  rk4_step)
+from mirgecom.integrators import (
+    euler_step, lsrk54_step, lsrk144_step,
+    rk4_step, ssprk43_step
+)
 
 logger = logging.getLogger(__name__)
 
@@ -44,7 +44,8 @@
                          [(euler_step, 1),
                           (lsrk54_step, 4),
                           (lsrk144_step, 4),
-                          (rk4_step, 4)])
+                          (rk4_step, 4),
+                          (ssprk43_step, 3)])
 @pytest.mark.parametrize("local_dt", [True, False])
 def test_integrator_order(integrator, method_order, local_dt):
     """Test that time integrators have correct order."""
@@ -87,7 +88,8 @@ def rhs(t, state):
                          [(euler_step, 1),
                           (lsrk54_step, 4),
                           (lsrk144_step, 4),
-                          (rk4_step, 4)])
+                          (rk4_step, 4),
+                          (ssprk43_step, 3)])
 @pytest.mark.parametrize("local_dt", [True, False])
 def test_state_advancer(integrator, method_order, local_dt):
     """Test that time integrators have correct order."""