Rckirby/explicit

demos/mixed_wave/README

@@ -1 +1,4 @@
-This directory contains two demos for the mixed form of the wave equation using the TimeStepper and AdaptiveTimeStepper classes.  The point is that if one uses a symplectic RK method, then energy will be conserved up to roundoff error and solver tolerances for linear problems.
+This directory contains two demos for the mixed form of the wave equation using the TimeStepper class.
+We have one demo doing a fully implicit Gauss-Legendre method that exactly conserves energy but requires a stage-coupled linear solve at each time step.
+We compare this to the explicit time stepping widget with pseudo-energy-conserving methods, which require many more stages, but these can be computed independently by only inverting a mass matrix.
+

demos/mixed_wave/demo_RTwave_PEP.py.rst

@@ -0,0 +1,109 @@
+Solving the Mixed Wave Equation: Energy conservation
+====================================================
+
+Let :math:`\Omega` be the unit square with boundary :math:`\Gamma`.  We write
+the wave equation as a first-order system of PDE:
+
+.. math::
+
+   u_t + \nabla p & = 0
+
+   p_t + \nabla \cdot u & = 0
+
+together with homogeneous Dirichlet boundary conditions
+
+.. math::
+
+   p = 0 \quad \textrm{on}\ \Gamma
+
+In this form, at each time, :math:`u` is a vector-valued function in the Sobolev space :math:`H(\mathrm{div})` and `p` is a scalar-valued function.  If we select appropriate test functions :math:`v` and :math:`w`, then we can arrive at the weak form
+
+.. math::
+
+   (u_t, v) - (p, \nabla \cdot v) & = 0
+
+   (p_t, w) + (\nabla \cdot u, w) & = 0
+
+Note that in mixed formulations, the Dirichlet boundary condition is weakly
+enforced via integration by parts rather than strongly in the definition of
+the approximating space.
+
+In this example, we will use the next-to-lowest order Raviart-Thomas elements
+for the velocity variable :math:`u` and discontinuous piecewise linear
+polynomials for the scalar variable :math:`p`.
+
+Here is some typical Firedrake boilerplate and the construction of a simple
+mesh and the approximating spaces::
+
+  from firedrake import *
+  from irksome import Dt, MeshConstant, TimeStepper, PEPRK
+
+  N = 10
+
+  msh = UnitSquareMesh(N, N)
+  V = FunctionSpace(msh, "RT", 2)
+  W = FunctionSpace(msh, "DG", 1)
+  Z = V*W
+
+Now we can build the initial condition, which has zero velocity and a sinusoidal displacement::
+
+  x, y = SpatialCoordinate(msh)
+  up0 = project(as_vector([0, 0, sin(pi*x)*sin(pi*y)]), Z)
+  u0, p0 = split(up0)
+
+
+We build the variational form in UFL::
+
+  v, w = TestFunctions(Z)
+  F = inner(Dt(u0), v)*dx + inner(div(u0), w) * dx + inner(Dt(p0), w)*dx - inner(p0, div(v)) * dx
+
+Energy conservation is an important principle of the wave equation, and we can
+test how well the spatial discretization conserves energy by creating a
+UFL expression and evaluating it at each time step::
+
+  E = 0.5 * (inner(u0, u0)*dx + inner(p0, p0)*dx)
+
+The time and time step variables::
+
+  MC = MeshConstant(msh)
+  t = MC.Constant(0.0)
+  dt = MC.Constant(0.2/N)
+
+The PEP RK methods of de Leon, Ketcheson, and Ranoch offer explicit RK methods
+that preserve the energy up to a given order in step size.  They have more
+stages than classical explicit methods but have much better energy conservation.::
+
+  butcher_tableau = PEPRK(4, 2, 5)
+
+We'll use a simple iterative method to invert the mass matrix for each stage::
+
+  params = {"snes_type": "ksponly",
+            "ksp_type": "cg",
+            "pc_type": "icc"}
+
+  stepper = TimeStepper(F, butcher_tableau, t, dt, up0,
+                        stage_type="explicit",
+                        solver_parameters=params)
+
+
+And, as with the heat equation, our time-stepping logic is quite simple.  At each time step, we print out the energy in the system::
+
+  print("Time    Energy")
+  print("==============")
+
+  while (float(t) < 1.0):
+      if float(t) + float(dt) > 1.0:
+          dt.assign(1.0 - float(t))
+
+      stepper.advance()
+
+      t.assign(float(t) + float(dt))
+      print("{0:1.1e} {1:5e}".format(float(t), assemble(E)))
+
+If all is well with the world, the energy will be nearly identical (up
+to roundoff error) at each time step because the PEP methods conserve
+energy to quite high order.  The reader can compare this to the mixed
+wave demo using Gauss-Legendre methods (which exactly conserve energy up
+to roundoff and solver tolerances.
+
+

docs/source/index.rst

@@ -71,6 +71,13 @@ We now have support for DIRKs:
 
    demos/demo_heat_dirk.py
 
+and for explicit schemes:
+
+.. toctree::
+   :maxdepth: 1
+
+   demos/demo_RTwave_PEP.py
+
 Or check out an IMEX-type method for the monodomain equations:
 
 .. toctree::

irksome/__init__.py

@@ -1,22 +1,23 @@
-from .ButcherTableaux import Alexander  # noqa: F401
+from .ButcherTableaux import Alexander      # noqa: F401
 from .ButcherTableaux import BackwardEuler  # noqa: F401
 from .ButcherTableaux import GaussLegendre  # noqa: F401
-from .ButcherTableaux import LobattoIIIA  # noqa: F401
-from .ButcherTableaux import LobattoIIIC  # noqa: F401
+from .ButcherTableaux import LobattoIIIA    # noqa: F401
+from .ButcherTableaux import LobattoIIIC    # noqa: F401
 from .ButcherTableaux import PareschiRusso  # noqa: F401
-from .ButcherTableaux import QinZhang  # noqa: F401
-from .ButcherTableaux import RadauIIA  # noqa: F401
-from .ButcherTableaux import WSODIRK432  # noqa: F401
-from .ButcherTableaux import WSODIRK433  # noqa: F401
-from .ButcherTableaux import WSODIRK643  # noqa: F401
-from .ButcherTableaux import WSODIRK744  # noqa: F401
-from .ButcherTableaux import WSODIRK1254  # noqa: F401
-from .ButcherTableaux import WSODIRK1255  # noqa: F401
-from .deriv import Dt  # noqa: F401
-from .dirk_stepper import DIRKTimeStepper  # noqa: F401
-from .getForm import getForm  # noqa: F401
-from .imex import RadauIIAIMEXMethod  # noqa: F401
-from .pc import RanaBase, RanaDU, RanaLD  # noqa: F401
-from .stage import StageValueTimeStepper  # noqa: F401
-from .stepper import TimeStepper  # noqa: F401
-from .tools import MeshConstant  # noqa: F401
+from .ButcherTableaux import QinZhang       # noqa: F401
+from .ButcherTableaux import RadauIIA       # noqa: F401
+from .ButcherTableaux import WSODIRK432     # noqa: F401
+from .ButcherTableaux import WSODIRK433     # noqa: F401
+from .ButcherTableaux import WSODIRK643     # noqa: F401
+from .ButcherTableaux import WSODIRK744     # noqa: F401
+from .ButcherTableaux import WSODIRK1254    # noqa: F401
+from .ButcherTableaux import WSODIRK1255    # noqa: F401
+from .pep_explicit_rk import PEPRK          # noqa: F401
+from .deriv import Dt                       # noqa: F401
+from .dirk_stepper import DIRKTimeStepper   # noqa: F401
+from .getForm import getForm                # noqa: F401
+from .imex import RadauIIAIMEXMethod        # noqa: F401
+from .pc import RanaBase, RanaDU, RanaLD    # noqa: F401
+from .stage import StageValueTimeStepper    # noqa: F401
+from .stepper import TimeStepper            # noqa: F401
+from .tools import MeshConstant             # noqa: F401

irksome/explicit_stepper.py

@@ -0,0 +1,21 @@
+from .dirk_stepper import DIRKTimeStepper
+
+
+# We can reuse the DIRK stepper to do one-stage at a time, but since we're
+# just solving a mass matrix at each time step we can optimize to
+# never rebuild the jacobian or preconditioner.
+class ExplicitTimeStepper(DIRKTimeStepper):
+    def __init__(self, F, butcher_tableau, t, dt, u0, bcs=None,
+                 solver_parameters=None,
+                 appctx=None):
+        assert butcher_tableau.is_explicit
+        # we just have one mass matrix we're reusing for each time step and
+        # each stage, so we can nudge this along
+        solver_parameters["snes_lag_jacobian_persists"] = "true"
+        solver_parameters["snes_lag_jacobian"] = -2
+        solver_parameters["snes_lag_preconditioner_persists"] = "true"
+        solver_parameters["snes_lag_preconditioner"] = -2
+        super(ExplicitTimeStepper, self).__init__(
+            F, butcher_tableau, t, dt, u0, bcs=bcs,
+            solver_parameters=solver_parameters, appctx=appctx,
+            nullspace=None)

irksome/pep_explicit_rk.py

@@ -0,0 +1,70 @@
+from .ButcherTableaux import ButcherTableau
+import numpy as np
+
+
+pep425A = np.array(
+    [[0, 0, 0, 0],
+     [1/10, 0, 0, 0],
+     [-35816/35721, 56795/35721, 0, 0],
+     [11994761/5328000, -11002961/4420800, 215846127/181744000, 0]])
+pep425b = np.array([-17/222, 6250/15657, 5250987/10382126, 4000/23307])
+pep425c = np.array([0, 1/10, 19/20, 37/63])
+
+pep526A = np.array(
+    [[0, 0, 0, 0, 0],
+     [0.193445628056365, 0, 0, 0, 0],
+     [-0.090431947690469, 0.646659568003039, 0, 0, 0],
+     [-0.059239621354435, 0.598571867726670, -0.010476084304794, 0, 0],
+     [0.173154586278662, 0.043637751980064, 0.949323298732961, -0.262838451019868, 0]])
+pep526b = np.array([0.054828314201395, 0.310080077556546, 0.531276882919990, -0.135494569336049, 0.239309294658118])
+pep526c = [0, 0.193445628056365, 0.55622762031257, 0.528856162067441, 0.9032771859718189]
+
+pep636A = np.array(
+    [[0, 0, 0, 0, 0, 0],
+     [0.12316523079127038, 0, 0, 0, 0, 0],
+     [-0.53348119048187126, 1.1200645707708279, 0, 0, 0, 0],
+     [0.35987162974687092, -0.17675778446586507, 0.7331973326225617, 0, 0, 0],
+     [0.015700424346522388, 0.02862938097533644, -0.014047147149911631, -0.015653338246176568, 0, 0],
+     [-1.9608805853984794, -0.82154709029385564, -0.0033631561953843502, 0.046367461001250457, 2.782035718578454, 0]])
+pep636b = np.array([0.78642719559722885, 0.69510370728230297, 0.42190724518033551, 0.21262030193155254, -0.70167978222250704, -0.41437866776891263])
+pep636c = np.array([0, 0.12316523079127044, 0.58658338028895673, 0.91631117790356775, 0.014629319925770667, 0.042612347691984923])
+
+pep746A = np.array(
+    [[0, 0, 0, 0, 0, 0, 0],
+     [-0.10731260966924323, 0, 0, 0, 0, 0, 0],
+     [0.14772934954602848, -0.12537555684690285, 0, 0, 0, 0, 0],
+     [0.7016079790308741, -0.75094597518803941, 0.76631666070124027, 0, 0, 0, 0],
+     [-0.8967481787471202, -0.43795858531068965, 1.7727346351832869, 0.1706052810617312, 0, 0, 0],
+     [1.6243872270239892, -0.69700589895015241, -0.3861309831750398, -0.032848941899304235, 0.30227620385295728, 0, 0],
+     [-0.32463926305048885, -0.3480143346241919, 1.3500419757109139, 0.039096802121597336, -0.17851883247877129, 0.010142489530892661, 0]])
+pep746b = np.array([-0.69203318482299292, 0.0074442860308153933, 0.93216717844052677, -1.159431111205361, 0.27787978605406632, 0.93890392164164138, 0.69506912386130404])
+pep746c = np.array([0, -0.10731260966924323, 0.022353792699125609, 0.71697866454407488, 0.60863315218720804, 0.81067760685245005, 0.54810883720995185])
+
+pep756A = np.array(
+    [[0, 0, 0, 0, 0, 0, 0],
+     [0.34288981581855521, 0, 0, 0, 0, 0, 0],
+     [0.16800230418143236, 0.1262987524809161, 0, 0, 0, 0, 0],
+     [0.4326925567104672, -0.24221982610439177, 0.15241708521248304, 0, 0, 0, 0],
+     [0.019843989305203335, 0.20330206481276515, -0.3494376489494413, 0.09780248603799992, 0, 0, 0],
+     [3.5441758455721732, 9.884560134482289, -3.7993663287883006, -6.07804112569088, -2.820029405964353, 0, 0],
+     [-16.625817935606782, -49.999620978741511, 22.3661445506308, 30.50526767511958, 13.408435545803448, 1.3455911427944685, 0]])
+pep756b = np.array([0.15881394125505754, 3.390357323579911e-13, 0.4109696726168125, -1.6409254928717294e-13, -0.056173857997504642, 0.40542999348169673, 0.08096025064376304])
+pep756c = np.array([0, 0.34288981581855521, 0.2943010566234846, 0.34288981581855849, -0.028489108793472939, 0.73129911961092908, 1.0000000000000007])
+
+pepdict = {
+    (4, 2, 5): (pep425A, pep425b, pep425c),
+    (5, 2, 6): (pep526A, pep526b, pep526c),
+    (6, 3, 6): (pep636A, pep636b, pep636c),
+    (7, 4, 6): (pep746A, pep746b, pep746c),
+    (7, 5, 6): (pep756A, pep756b, pep756c)
+}
+
+
+class PEPRK(ButcherTableau):
+    def __init__(self, ns, order, peporder):
+        try:
+            A, b, c = pepdict[ns, order, peporder]
+        except KeyError:
+            raise NotImplementedError("No PEP method for that combination of stages, order and pseudo-energy preserving order")
+        self.peporder = peporder
+        super(PEPRK, self).__init__(A, b, None, c, order)

irksome/stepper.py

@@ -3,6 +3,7 @@
 from firedrake import NonlinearVariationalSolver as NLVS
 from firedrake.dmhooks import pop_parent, push_parent
 from .dirk_stepper import DIRKTimeStepper
+from .explicit_stepper import ExplicitTimeStepper
 from .getForm import AI, getForm
 from .stage import StageValueTimeStepper
 from .imex import RadauIIAIMEXMethod
@@ -55,6 +56,7 @@ def TimeStepper(F, butcher_tableau, t, dt, u0, **kwargs):
         "value": ["stage_type", "bcs", "nullspace", "solver_parameters",
                   "update_solver_parameters", "appctx", "splitting"],
         "dirk": ["stage_type", "bcs", "nullspace", "solver_parameters", "appctx"],
+        "explicit": ["stage_type", "bcs", "solver_parameters", "appctx"],
         "imex": ["Fexp", "stage_type", "bcs", "nullspace",
                  "it_solver_parameters", "prop_solver_parameters",
                  "splitting", "appctx",
@@ -97,6 +99,13 @@ def TimeStepper(F, butcher_tableau, t, dt, u0, **kwargs):
         return DIRKTimeStepper(
             F, butcher_tableau, t, dt, u0, bcs,
             solver_parameters, appctx, nullspace)
+    elif stage_type == "explicit":
+        bcs = kwargs.get("bcs")
+        appctx = kwargs.get("appctx")
+        solver_parameters = kwargs.get("solver_parameters")
+        return ExplicitTimeStepper(
+            F, butcher_tableau, t, dt, u0, bcs,
+            solver_parameters, appctx)
     elif stage_type == "imex":
         Fexp = kwargs.get("Fexp")
         assert Fexp is not None

tests/test_split.py

@@ -157,8 +157,8 @@ def Ffull(z, test):
         "snes_type": "ksponly",
         "ksp_type": "preonly",
         "pc_type": "lu",
-        "pc_factor_mat_solver_type": "pastix",
-        "pc_factor_shift_type": "inblocks"}
+        "pc_factor_mat_solver_type": "mumps"
+    }
 
     F_full = Ffull(z_imp, test_z)
     F_imp = Fimp(z_split, test_z)

tests/test_wave_energy.py

@@ -1,17 +1,16 @@
-import pytest
 import numpy as np
-from firedrake import (inner, dx, UnitIntervalMesh, FunctionSpace,
-                       assemble, TestFunctions, SpatialCoordinate,
-                       project, as_vector, sin, pi, split)
-
-from irksome import Dt, MeshConstant, TimeStepper, GaussLegendre
+import pytest
+from firedrake import (FunctionSpace, SpatialCoordinate, TestFunctions,
+                       UnitIntervalMesh, as_vector, assemble, dx, inner, pi,
+                       project, sin, split)
+from irksome import PEPRK, Dt, GaussLegendre, MeshConstant, TimeStepper
 from irksome.tools import AI, IA
 
 # test the energy conservation of the 1d wave equation in mixed form
 # various time steppers.
 
 
-def wave(n, deg, butcher_tableau, splitting=AI):
+def wave(n, deg, alpha, butcher_tableau, **kwargs):
     N = 2**n
     msh = UnitIntervalMesh(N)
 
@@ -28,7 +27,7 @@ def wave(n, deg, butcher_tableau, splitting=AI):
 
     MC = MeshConstant(msh)
     t = MC.Constant(0.0)
-    dt = MC.Constant(2.0 / N)
+    dt = MC.Constant(alpha / N)
 
     up = project(as_vector([0, sin(pi*x)]), Z)
     u, p = split(up)
@@ -42,7 +41,7 @@ def wave(n, deg, butcher_tableau, splitting=AI):
 
     stepper = TimeStepper(F, butcher_tableau, t, dt, up,
                           solver_parameters=params,
-                          splitting=splitting)
+                          **kwargs)
 
     energies = []
 
@@ -62,6 +61,19 @@ def wave(n, deg, butcher_tableau, splitting=AI):
                          [(1, i) for i in (1, 2)]
                          + [(2, i) for i in (2, 3)])
 def test_wave_eq(deg, N, time_stages, splitting):
-    energy = wave(N, deg, GaussLegendre(time_stages), splitting)
-    print(energy)
+    kwargs = {"splitting": splitting}
+    energy = wave(N, deg, 2.0, GaussLegendre(time_stages), **kwargs)
+    assert np.allclose(energy[1:], energy[:-1])
+
+
+@pytest.mark.parametrize('N', [2**j for j in range(2, 3)])
+@pytest.mark.parametrize('deg', (1, 2))
+@pytest.mark.parametrize('pep', ((4, 2, 5),
+                                 (5, 2, 6),
+                                 (6, 3, 6),
+                                 (7, 4, 6),
+                                 (7, 5, 6)))
+def test_wave_eq_peprk(deg, N, pep):
+    kwargs = {"stage_type": "explicit"}
+    energy = wave(N, deg, 0.3, PEPRK(*pep), **kwargs)
     assert np.allclose(energy[1:], energy[:-1])