diff --git a/demos/mixed_wave/README b/demos/mixed_wave/README index 0318ea6e..e3ad6ff4 100644 --- a/demos/mixed_wave/README +++ b/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. \ No newline at end of file +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. + diff --git a/demos/mixed_wave/demo_RTwave_PEP.py.rst b/demos/mixed_wave/demo_RTwave_PEP.py.rst new file mode 100644 index 00000000..60f9258c --- /dev/null +++ b/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. + + diff --git a/docs/source/index.rst b/docs/source/index.rst index 541f0176..2e6fe8d6 100644 --- a/docs/source/index.rst +++ b/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:: diff --git a/irksome/__init__.py b/irksome/__init__.py index 5fe129ac..46433bdd 100644 --- a/irksome/__init__.py +++ b/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 diff --git a/irksome/explicit_stepper.py b/irksome/explicit_stepper.py new file mode 100644 index 00000000..065c2a94 --- /dev/null +++ b/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) diff --git a/irksome/pep_explicit_rk.py b/irksome/pep_explicit_rk.py new file mode 100644 index 00000000..e3e396d0 --- /dev/null +++ b/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) diff --git a/irksome/stepper.py b/irksome/stepper.py index 799128fc..8152cb9e 100644 --- a/irksome/stepper.py +++ b/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 diff --git a/tests/test_split.py b/tests/test_split.py index 1c49c427..83e07141 100644 --- a/tests/test_split.py +++ b/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) diff --git a/tests/test_wave_energy.py b/tests/test_wave_energy.py index 13132aae..1735420d 100644 --- a/tests/test_wave_energy.py +++ b/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])