Solving the Benjamin-Bona-Mahony equation

This demo solves the Benjamin-Bona-Mahony equation:

\[u_t + u_x + u u_x - u_{txx} = 0\]

typically posed on \(\mathbb{R}\) or a bounded interval with periodic boundaries.

It is interesting because it is nonlinear (\(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:

\[(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:

\[ \begin{align}\begin{aligned}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\end{aligned}\end{align} \]

We are mainly interested in accuracy and in conserving these quantities reasonably well.

Firedrake imports:

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

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)

The symplectic Gauss-Legendre methods are of interest here:

butcher_tableau = GaussLegendre(2)

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"}

stepper = TimeStepper(F, butcher_tableau, 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 \(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))