Hamiltonian PDE demos

demos/bbm_auxiliary/demo_bbm_aux.py.rst

@@ -1,5 +1,5 @@
-Hamiltonian-structure-preserving implementation of the Benjamin-Bona-Mahoney equation
-=====================================================================================
+Mixed Hamiltonian-preserving formulation of the Benjamin-Bona-Mahoney equation
+==============================================================================
 
 This demo solves the Benjamin-Bona-Mahony equation:
 
@@ -15,45 +15,47 @@ BBM is known to have a Hamiltonian structure, and there are several canonical po
 
    I_1 & = \int u \, \mathrm{d}x
 
-   I_2 & = \int u^2 + (u_x)^2 \, \mathrm{d}x
+   I_2 & = \int \frac{u^2 + (u_x)^2}{2} \, \mathrm{d}x
 
    I_3 & = \int \frac{u^2}{2} + \frac{u^3}{6} \, \mathrm{d}x
 
 The BBM invariants are the total momentum :math:`I_1`, the :math:`H^1`-energy
-norm :math:`I_2`, and the Hamiltonian :math:`I_3`.  
-The Hamiltonian variational formulation reads
+:math:`I_2`, and the Hamiltonian :math:`I_3`.  
+
+Standard Gauss-Legendre and continuous Petrov-Galerkin (cPG) methods conserve
+the first two invariants exactly (up to roundoff and solver tolerances).  They
+do quite well, but are inexact for the cubic one.  In this demo, we consider a
+mixed Hamiltonian formulation that preserves the third invariant at the expense
+of the second. The mixed formulation solves
 
 .. math::
 
-   (\partial_t u + \partial_x \tilde{w}_H, v)_{H^1} & = 0
+   (\partial_t u, v)_{H^1} & = (\tilde{w}_H, \partial_x v)_{L^2}
 
-   (\tilde{w}_H, v_H)_{H^1} & = \langle \frac{\delta I_3}{\delta u}, v_H \rangle 
+   (\tilde{w}_H, v_H)_{L^2} & = \langle \frac{\delta I_3}{\delta u}, v_H \rangle 
 
-For all test functions :math:`v, v_H` in a suitable function space.
-The numerical scheme in this demo introduces
-the :math:`H^1`-Riesz representative of the Fréchet derivative of the
-Hamiltonian :math:`\frac{\delta I_3}{\delta u}` 
-as the auxiliary variable :math:`\tilde{w}_H`.
+for all test functions :math:`v, v_H` in a suitable function space :math:`V \times V`.
+In this demo we choose to discretize :math:`V` with :math:`C^0`
+Lagrange elements by introducing the auxiliary variable :math:`\tilde{w}_H \in V`
+that holds the :math:`L^2`-Riesz representative of the Fréchet derivative of the
+Hamiltonian :math:`\frac{\delta I_3}{\delta u}`.
 
-Standard Gauss-Legendre and continuous Petrov-Galerkin (cPG) methods conserve
-the first two invariants exactly (up to roundoff and solver tolerances).  They
-do quite well, but are inexact for the cubic one.  Here, we consider the
-reformulation in Andrews and Farrell, "Enforcing conservation laws and dissipation
-inequalities numerically via auxiliary variables" (`arXiv:2407.11904 <https://arxiv.org/abs/2407.11904>`_, to appear
-in SIAM J. Scientific Computing) that preserves the third invariant at
-the expense of the second.  This method has an auxiliary variable in the system
-and requires a continuously differentiable spatial discretization (1D Hermite
-elements in this case).  The time discretization puts the main unknown in a
-continuous space and the auxiliary variable in a discontinuous one.  See
-equation (7.17) of Boris Andrews' thesis for the particular formulation.
+.. note::
+
+   Here, we consider the framework in Andrews and Farrell, "Enforcing conservation laws and dissipation
+   inequalities numerically via auxiliary variables" (`arXiv:2407.11904 <https://arxiv.org/abs/2407.11904>`_, to appear
+   in SIAM J. Scientific Computing).  Their method has an auxiliary variable in the system
+   and requires a continuously differentiable spatial discretization (1D Hermite
+   elements in their case).  The time discretization puts the main unknown in a
+   continuous space and the auxiliary variable in a discontinuous one. See
+   equation (7.17) of Boris Andrews' thesis for their particular formulation.
 
 
 Firedrake, Irksome, and other imports::
 
   from firedrake import (Constant, Function, FunctionSpace,
       PeriodicIntervalMesh, SpatialCoordinate, TestFunction, TrialFunction,
-      assemble, derivative, dx, errornorm, exp, grad, inner,
-      interpolate, norm, plot, project, replace, solve, split
+      assemble, derivative, dx, exp, grad, inner, plot, replace, solve, split
   )
 
   from irksome import Dt, GalerkinTimeStepper, TimeQuadratureLabel
@@ -85,13 +87,13 @@ Next, we define the domain and the exact solution ::
   uexact = 3 * c**2 / (1-c**2) \
       * sech(0.5 * (c * x - c * t / (1 - c ** 2) + delta))**2
 
-This sets up the function space for the unknown :math:`u` and
+This sets up the mixed function space for the unknown :math:`u` and
 auxiliary variable :math:`\tilde{w}_H`::
 
-  space_deg = 3
+  space_deg = 1
   time_deg = 1
 
-  V = FunctionSpace(msh, "Hermite", space_deg)
+  V = FunctionSpace(msh, "CG", space_deg)
   Z = V * V
 
 We next define the BBM invariants. Again, the discrete formulation preserves 
@@ -105,41 +107,41 @@ but :math:`I_2` is preserved up to a bounded oscillation ::
       return u * dx
 
   def I2(u):
-      return h1inner(u, u) * dx
+      return 0.5 * h1inner(u, u) * dx
 
   def I3(u):
       return (u**2 / 2 + u**3 / 6) * dx
 
-We project the initial condition on :math:`u`, but we also need a consistent
+We project the initial condition on :math:`u` in the :math:`H^1` norm, but we also need a consistent
 initial condition for the auxiliary variable.  We need to find :math:`\tilde{w}_H \in V` such that
 
 .. math::
 
-   (\tilde{w}_H, v)_{H^1} = \langle \frac{\delta I_3}{\delta u}, v \rangle \text{ for all } v \in V
+   (\tilde{w}_H, v)_{L^2} = \langle \frac{\delta I_3}{\delta u}, v \rangle \text{ for all } v \in V
 
 ::
 
-  uwH = Function(Z)
-  u0, wH0 = uwH.subfunctions
+  uw = Function(Z)
+  u0, w0 = uw.subfunctions
 
   v = TestFunction(V)
   w = TrialFunction(V)
-  a = h1inner(w, v) * dx
-  dHdu = derivative(I3(u0), u0, v)
 
-  solve(a == h1inner(uexact, v)*dx, u0)
-  solve(a == dHdu, wH0)
+  solve(h1inner(w, v)*dx == h1inner(uexact, v)*dx, u0)
+
+  dHdu = derivative(I3(u0), u0, v)
+  solve(inner(w, v)*dx == dHdu, w0)
 
 Visualize the initial condition::
 
   fig, axes = plt.subplots(1)
-  plot(Function(FunctionSpace(msh, "CG", 1)).interpolate(u0), axes=axes)
+  plot(u0, axes=axes)
   axes.set_title("Initial condition")
   axes.set_xlabel("x")
   axes.set_ylabel("u")
-  plt.savefig("bbm_init.png")
+  plt.savefig("bbm_aux_init.png")
 
-.. figure:: bbm_init.png
+.. figure:: bbm_aux_init.png
    :align: center  
 
 Create time quadrature labels::
@@ -153,10 +155,9 @@ Create time quadrature labels::
 This tags several of the terms with a low-order time integration scheme,
 but forces a higher-order method on the nonlinear term::
 
-  u, wH = split(uwH)
+  u, w = split(uw)
   v, vH = split(TestFunction(Z))
-
-  Flow = h1inner(Dt(u) + wH.dx(0), v) * dx + h1inner(wH, vH) * dx
+  Flow = h1inner(Dt(u), v) * dx - inner(w, v.dx(0))*dx + inner(w, vH)*dx
   Fhigh = replace(dHdu, {u0: u})
 
   F = Llow(Flow) - Lhigh(Fhigh(vH))
@@ -165,8 +166,10 @@ This sets up the cPG time stepper.  There are two fields in the unknown, we
 indicate the second one is an auxiliary and hence to be discretized in the DG
 test space instead by passing the `aux_indices` keyword::
 
-  stepper = GalerkinTimeStepper(
-      F, time_deg, t, dt, uwH, aux_indices=[1])
+  sparams = {"snes_atol": 0, "snes_rtol": 1E-14}
+  stepper = GalerkinTimeStepper(F, time_deg, t, dt, uw,
+                                aux_indices=[1],
+                                solver_parameters=sparams)
 
 UFL expressions for the invariants, which we are going to track as we go
 through time steps::
@@ -201,7 +204,7 @@ Visualize invariant preservation::
   axes.set_xlabel("Time")
   axes.set_ylabel("I(t)")
   axes.legend()
-  plt.savefig("invariants.png")
+  plt.savefig("bbm_aux_invariants.png")
   axes.clear()
 
   for i in (0, 1, 2):
@@ -210,22 +213,22 @@ Visualize invariant preservation::
   axes.set_xlabel("Time")
   axes.set_ylabel("|1-I/I(0)|")  
   axes.legend()  
-  plt.savefig("invariant_errors.png")
+  plt.savefig("bbm_aux_errors.png")
 
-.. figure:: invariants.png
+.. figure:: bbm_aux_invariants.png
    :align: center
 
-.. figure:: invariant_errors.png
+.. figure:: bbm_aux_errors.png
    :align: center
 
 Visualize the solution at final time step::
 
   axes.clear()
-  plot(Function(FunctionSpace(msh, "CG", 1)).interpolate(u0), axes=axes)
+  plot(u0, axes=axes)
   axes.set_title(f"Solution at time {tfinal}")
   axes.set_xlabel("x")
   axes.set_ylabel("u")  
-  plt.savefig("bbm_final.png") 
+  plt.savefig("bbm_aux_final.png") 
 
-.. figure:: bbm_final.png
+.. figure:: bbm_aux_final.png
    :align: center