Solving the Mixed Wave Equation: Energy conservation

Let \(\Omega\) be the unit square with boundary \(\Gamma\). We write the wave equation as a first-order system of PDE:

\[ \begin{align}\begin{aligned}u_t + \nabla p & = 0\\p_t + \nabla \cdot u & = 0\end{aligned}\end{align} \]

together with homogeneous Dirichlet boundary conditions

\[p = 0 \quad \textrm{on}\ \Gamma\]

In this form, at each time, \(u\) is a vector-valued function in the Soboleve space \(H(\mathrm{div})\) and p is a scalar-valued function. If we select appropriate test functions \(v\) and \(w\), then we can arrive at the weak form

\[ \begin{align}\begin{aligned}(u_t, v) - (p, \nabla \cdot v) & = 0\\(p_t, w) + (\nabla \cdot u, w) & = 0\end{aligned}\end{align} \]

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 \(u\) and discontinuous piecewise linear polynomials for the scalar variable \(p\).

Here is some typical Firedrake boilerplate and the construction of a simple mesh and the approximating spaces:

from firedrake import *
from irksome import GaussLegendre, Dt, MeshConstant, TimeStepper

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(1.0/N)

The two-stage Gauss-Legendre method is, like all instances of that family, A-stable and symplectic. This gives us a fourth order method in time, although our spatial accuracy is of lower order. Feel free to experiment!:

butcher_tableau = GaussLegendre(2)

Like the heat equation demo, we are just using a direct method to solve the system at each time step:

params = {"mat_type": "aij",
          "snes_type": "ksponly",
          "ksp_type": "preonly",
          "pc_type": "lu"}

stepper = TimeStepper(F, butcher_tableau, t, dt, up0,
                      solver_parameters=params)

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 GL methods are symplectic and applied to a linear Hamiltonian system. As an exercise, the reader should edit this code to use other RK methods. In particular, LobattoIIIC and BackwardEuler or other Radau methods as well as the symplectic DIRK QinZhang.