Merge branch 'master' into smaclachlan/add_dirk_imex

demos/bbm/demo_bbm_galerkin.py.rst

@@ -0,0 +1,112 @@
+Solving the Benjamin-Bona-Mahony equation with Galerkin-in-Time
+===============================================================
+
+This demo solves the Benjamin-Bona-Mahony equation:
+
+.. math::
+
+   u_t + u_x + u u_x - u_{txx} = 0
+
+typically posed on :math:`\mathbb{R}` or a bounded interval with periodic
+boundaries.
+
+It is interesting because it is nonlinear (:math:`u u_x`) and also a Sobolev-type equation, with spatial derivatives acting on a time derivative.  We can obtain a weak form in the standard way:
+
+.. math::
+
+   (u_t, v) + (u_x, v) + (u u_x, v) + (u_{tx}, v_x) = 0
+
+BBM is known to have a Hamiltonian structure, and there are three canonical polynomial invariants:
+
+.. math::
+
+   I_1 & = \int u \, dx
+
+   I_2 & = \int u^2 + (u_x)^2 \, dx
+
+   I_3 & = \int (u_x)^2 + \tfrac{1}{3} u^3 \, dx
+
+We are mainly interested in accuracy and in conserving these quantities reasonably well.
+
+
+Firedrake imports::
+
+  from firedrake import *
+  from irksome import Dt, GalerkinTimeStepper, MeshConstant
+
+This function seems to be left out of UFL, but solitary wave solutions for BBM need it::
+
+  def sech(x):
+      return 2 / (exp(x) + exp(-x))
+
+Set up problem parameters, etc::
+
+  N = 1000
+  L = 100
+  h = L / N
+  msh = PeriodicIntervalMesh(N, L)
+
+  MC = MeshConstant(msh)
+  t = MC.Constant(0)
+  dt = MC.Constant(10*h)
+
+  x, = SpatialCoordinate(msh)
+
+Here is the true solution, which is right-moving solitary wave (but not a soliton)::
+
+  c = Constant(0.5)
+
+  center = 30.0
+  delta = -c * center
+
+  uexact = 3 * c**2 / (1-c**2) * sech(0.5 * (c * x - c * t / (1 - c ** 2) + delta))**2
+
+Create the approximating space and project true solution::
+
+  V = FunctionSpace(msh, "CG", 1)
+  u = project(uexact, V)
+  v = TestFunction(V)
+
+  F = (inner(Dt(u), v) * dx
+       + inner(u.dx(0), v) * dx
+       + inner(u * u.dx(0), v) * dx
+       + inner((Dt(u)).dx(0), v.dx(0)) * dx)
+
+For a 1d problem, we don't worry much about efficient solvers.::
+
+  params = {"mat_type": "aij",
+            "ksp_type": "preonly",
+            "pc_type": "lu"}
+
+The Galerkin-in-Time approach should have symplectic properties::
+
+  stepper = GalerkinTimeStepper(F, 2, t, dt, u,
+                        solver_parameters=params)
+
+UFL for the mathematical invariants and containers to track them over time::
+
+  I1 = u * dx
+  I2 = (u**2 + (u.dx(0))**2) * dx
+  I3 = ((u.dx(0))**2 - u**3 / 3) * dx
+
+  I1s = []
+  I2s = []
+  I3s = []
+
+
+Time-stepping loop, keeping track of :math:`I` values::
+
+  tfinal = 18.0
+  while (float(t) < tfinal):
+      if float(t) + float(dt) > tfinal:
+          dt.assign(tfinal - float(t))
+      stepper.advance()
+
+      I1s.append(assemble(I1))
+      I2s.append(assemble(I2))
+      I3s.append(assemble(I3))
+
+      print('%.15f %.15f %.15f %.15f' % (float(t), I1s[-1], I2s[-1], I3s[-1]))
+      t.assign(float(t) + float(dt))
+
+  print(errornorm(uexact, u) / norm(uexact))

demos/mixed_wave/demo_RTwave_PEP.py.rst

@@ -1,5 +1,5 @@
-Solving the Mixed Wave Equation: Energy conservation
-====================================================
+Solving the Mixed Wave Equation: Energy conservation and explicit RK
+====================================================================
 
 Let :math:`\Omega` be the unit square with boundary :math:`\Gamma`.  We write
 the wave equation as a first-order system of PDE:

demos/mixed_wave/demo_RTwave_galerkin.py.rst

@@ -0,0 +1,130 @@
+Solving the Mixed Wave Equation: Energy conservation, Multigrid, Galerkin-in-Time
+=================================================================================
+
+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 Soboleve 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.  We are going to use a multigrid preconditioner for each timestep, so we create a MeshHierarchy as well::
+
+  from firedrake import *
+  from irksome import Dt, MeshConstant, GalerkinTimeStepper
+
+  N = 10
+
+  base = UnitSquareMesh(N, N)
+  mh = MeshHierarchy(base, 2)
+  msh = mh[-1]
+  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(1.0/N)
+
+Here, we experiment with a multigrid preconditioner for the CG(2)-in-time discretization::
+
+  mgparams = {"mat_type": "aij",
+              "snes_type": "ksponly",
+              "ksp_type": "fgmres",
+	      "ksp_rtol": 1e-8,
+              "pc_type": "mg",
+              "mg_levels": {
+                  "ksp_type": "chebyshev",
+                  "ksp_max_it": 2,
+                  "ksp_convergence_test": "skip",
+                  "pc_type": "python",
+                  "pc_python_type": "firedrake.PatchPC",
+                  "patch": {
+                      "pc_patch": {
+                          "save_operators": True,
+                          "partition_of_unity": False,
+                          "construct_type": "star",
+                          "construct_dim": 0,
+                          "sub_mat_type": "seqdense",
+                          "dense_inverse": True,
+                          "precompute_element_tensors": None},
+                      "sub": {
+                          "ksp_type": "preonly",
+                          "pc_type": "lu"}}},
+              "mg_coarse": {
+                  "pc_type": "lu",
+                  "pc_factor_mat_solver_type": "mumps"}
+              }
+  
+
+  stepper = GalerkinTimeStepper(F, 2, t, dt, up0,
+                                solver_parameters=mgparams)
+
+
+And, as with the heat equation, our time-stepping logic is quite simple.  At easch 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 Galerkin-in-time methods
+are symplectic and applied to a linear Hamiltonian system.
+
+We can also confirm that the multigrid preconditioner is effective, by computing the average number of linear iterations per time-step::
+
+  (steps, nl_its, linear_its) = stepper.solver_stats()
+  print(f"The average number of multigrid iterations per time-step is {linear_its/steps}.")

demos/preconditioning/demo_heat_mg_dg.py.rst

@@ -0,0 +1,123 @@
+Solving the heat equation with monolithic multigrid and Discontinuous Galerkin-in-Time
+======================================================================================
+
+This reprise of the heat equation demo uses a monolithic multigrid
+algorithm suggested by Patrick Farrell to perform time advancement,
+but now applies it to the Discontinuous Galerkin-in-Time discretization.
+
+We consider the heat equation on :math:`\Omega = [0,10]
+\times [0,10]`, with boundary :math:`\Gamma`: giving rise to the weak form
+
+.. math::
+
+   (u_t, v) + (\nabla u, \nabla v) = (f, v)
+
+This demo implements an example used by Solin with a particular choice
+of :math:`f` given below
+
+We perform similar imports and setup as before::
+
+  from firedrake import *
+  from irksome import DiscontinuousGalerkinTimeStepper, Dt, MeshConstant
+  from ufl.algorithms.ad import expand_derivatives
+
+
+However, we need to set up a mesh hierarchy to enable geometric multigrid
+within Firedrake::
+
+  N = 128
+  x0 = 0.0
+  x1 = 10.0
+  y0 = 0.0
+  y1 = 10.0
+
+  from math import log
+  coarseN = 8  # size of coarse grid
+  nrefs = log(N/coarseN, 2)
+  assert nrefs == int(nrefs)
+  nrefs = int(nrefs)
+  base = RectangleMesh(coarseN, coarseN, x1, y1)
+  mh = MeshHierarchy(base, nrefs)
+  msh = mh[-1]
+
+From here, setting up the function space, manufactured solution, etc,
+are just as for the regular heat equation demo::
+
+  V = FunctionSpace(msh, "CG", 1)
+
+  MC = MeshConstant(msh)
+  dt = MC.Constant(10.0 / N)
+  t = MC.Constant(0.0)
+
+  x, y = SpatialCoordinate(msh)
+  S = Constant(2.0)
+  C = Constant(1000.0)
+  B = (x-Constant(x0))*(x-Constant(x1))*(y-Constant(y0))*(y-Constant(y1))/C
+  R = (x * x + y * y) ** 0.5
+  uexact = B * atan(t)*(pi / 2.0 - atan(S * (R - t)))
+  rhs = expand_derivatives(diff(uexact, t)) - div(grad(uexact))
+
+  u = Function(V)
+  u.interpolate(uexact)
+  v = TestFunction(V)
+  F = inner(Dt(u), v)*dx + inner(grad(u), grad(v))*dx - inner(rhs, v)*dx
+  bc = DirichletBC(V, 0, "on_boundary")
+
+And now for the solver parameters.  Note that we are solving a
+block-wise system with all stages coupled together.  This performs a
+monolithic multigrid with pointwise block Jacobi preconditioning::
+
+  mgparams = {"mat_type": "aij",
+              "snes_type": "ksponly",
+              "ksp_type": "gmres",
+              "ksp_monitor_true_residual": None,
+              "pc_type": "mg",
+              "mg_levels": {
+                  "ksp_type": "chebyshev",
+                  "ksp_max_it": 1,
+                  "ksp_convergence_test": "skip",
+                  "pc_type": "python",
+                  "pc_python_type": "firedrake.PatchPC",
+                  "patch": {
+                      "pc_patch": {
+                          "save_operators": True,
+                          "partition_of_unity": False,
+                          "construct_type": "star",
+                          "construct_dim": 0,
+                          "sub_mat_type": "seqdense",
+                          "dense_inverse": True,
+                          "precompute_element_tensors": None},
+                      "sub": {
+                          "ksp_type": "preonly",
+                          "pc_type": "lu"}}},
+              "mg_coarse": {
+                  "pc_type": "lu",
+                  "pc_factor_mat_solver_type": "mumps"}
+              }
+  
+These solver parameters work just fine in the :class:`.DiscontinuousGalerkinTimeStepper`::
+
+  stepper = DiscontinuousGalerkinTimeStepper(F, 2, t, dt, u, bcs=bc,
+                                    solver_parameters=mgparams)
+
+And we can advance the solution in time in typical fashion::
+
+  while (float(t) < 1.0):
+      if (float(t) + float(dt) > 1.0):
+          dt.assign(1.0 - float(t))
+      stepper.advance()
+      print(float(t), flush=True)
+      t.assign(float(t) + float(dt))
+
+After the solve, we can retrieve some statistics about the solver::
+
+  steps, nonlinear_its, linear_its = stepper.solver_stats()
+
+  print("Total number of timesteps was %d" % (steps))
+  print("Average number of nonlinear iterations per timestep was %.2f" % (nonlinear_its/steps))
+  print("Average number of linear iterations per timestep was %.2f" % (linear_its/steps))
+
+Finally, we print out the relative :math:`L^2` error::
+
+  print()
+  print(norm(u-uexact)/norm(uexact))