HIGH-ORDER CONSERVATIVE AND ACCURATELY DISSIPATIVE NUMERICAL INTEGRATORS VIA AUXILIARY VARIABLES

Boris Andrews, Patrick Farrell

16.JUL.2024 (arXiv)

[…] we propose an approach for the construction of timestepping schemes that preserve dissipation laws and conserve multiple general invariants, via finite elements in time and the systematic introduction of auxiliary variables. […] We [devise] novel arbitrary-order schemes that conserve to machine precision all known invariants of Hamiltonian ODEs […] and arbitrary-order schemes for the compressible Navier–Stokes equations that conserve mass, momentum, and energy, and provably possess non-decreasing entropy.

FULL ABSTRACT

Numerical methods for the simulation of transient systems with structure-preserving properties are known to exhibit greater accuracy and physical reliability, in particular over long durations. These schemes are often built on powerful geometric ideas for broad classes of problems, such as Hamiltonian or reversible systems. However, there remain difficulties in devising timestepping schemes that conserve non-quadratic invariants or dissipation laws. In this work, we propose an approach for the construction of timestepping schemes that preserve dissipation laws and conserve multiple general invariants, via finite elements in time and the systematic introduction of auxiliary variables. The approach generalises several existing ideas in the literature, including Gauss methods, the framework of Cohen & Hairer, and the energy- and helicity-conserving scheme of Rebholz. We demonstrate the ideas by devising novel arbitrary-order schemes that conserve to machine precision all known invariants of Hamiltonian ODEs, including the Kepler and Kovalevskaya problems, and arbitrary-order schemes for the compressible Navier-Stokes equations that conserve mass, momentum, and energy, and provably possess non-decreasing entropy.

While the results of this work are more general, I would like to provide some exposition for it through the lens of:

Hamiltonian systems
symplectic integrators
the Benjamin–Bona–Mahony (BBM) equation

I find these results be informative, curious, motivating, and (not least) cool!

Symplectic integrators are frequently lauded for their “energy-conserving properties”. Their status as the gold standard for simulating Hamiltonian systems is often put down to this.
Yet, this belief is not entirely accurate.

Symplecticity enhances the collective behaviour of a group of simulations, which is beneficial! However, it does not guarantee energy conservation.
In fact, as noted by Ge and Marsden (1988):

Symplectic integrators cannot* conserve energy.
(*in general)

ge_marsden_quote

This limitation is evident in the Benjamin–Bona–Mahony (BBM) equation, a model for phenomena including long water waves.
Solutions to the BBM equation conserve energy, \(\int(\frac{1}{2}u^2 + \frac{1}{6}u^3)\), contributing to their stability and persistence over time.

waves_off_timor_sea

Simulating the BBM equations using the 2-stage Gauss method, a symplectic integrator, we observe a gradual decline in the simulated energy.
This decline manifests as artificial, unphysical oscillations in the solution.

FULL INTEGRATOR SPECIFICATIONS

The full BBM equation is \[ \dot{u} - \dot{u}_{xx} = - u_x - uu_x. \] Define the domain \(\Omega \coloneqq (-50, 50)\). Up to projection, take initial conditions to be a soliton of speed \(\frac{1 + \sqrt{5}}{2}\), \[ u(0) = \frac{3\sqrt{5} - 3}{2}{\rm sech}\!\left(\!\left(\sqrt{5} - 1\right)x\right)^2, \] and assume periodic boundary conditions. Define \(U\) to be the space of Hermite finite elements of uniform width 2, and take a uniform timestep \(1\).
Take the following semi-discrete variational formulation: find \(u \in C^1(\mathbb{R}_+; U)\) such that for all \(v \in U\), \[ \int_\Omega(\dot{u}v + \dot{u}_xv_x) = \int_\Omega\!\left(u + \frac{1}{2}u^2\right)v_x, \] at all times \(t \in \mathbb{R}_+\). Discretise this in time using the 2-stage Gauss method.
The above video is at time \(20000\) with a moving camera of speed \(\frac{1 + \sqrt{5}}{2}\).

In our preprint, Patrick Farrell and I propose a framework to modify numerical time discretisations to preserve conservation laws exactly. We achieve this through:

Finite elements in time
The systematic introduction of auxiliary variables

TECHNICAL NOTE ON SYMMETRY

N.B. Unlike other techniques for the construction of energy-conserving schemes (e.g. projection methods) this approach preserves the symmetry of the initial timestepping scheme, another property that is crucial for realistic simulations.

Applying our framework to Hamiltonian systems (incl. the BBM equation) we derive a general Hamiltonian integrator with exact energy conservation. Simulating the BBM equations using the modified 2-stage Gauss method, we observe exact energy conservation in the simulation. This avoids the artificial oscillations and provides far more qualitatively accurate results.
(Note: The video should appear stationary. It’s loading fine; this is the correct solution behaviour!)

FULL INTEGRATOR SPECIFICATIONS

Over a timestep \(T_n = [t_n, t_{n+1}]\), define the space-time finite element space \[ X_n \coloneqq \{u \in P_2(T_n; U) : u(t_n) \text{ satisfies known initial data}\}, \] where \(P_s(T_n; U)\) is the space of degree-\(s\) polynomials from \(T_n\) to \(U\). Note that \(\dot{X}_n = P_1(T_n; U)\).
Take the following fully discrete variational formulation over \(T_n\): find \((u, \tilde{w}) \in X_n \times \dot{X}_n\) such that for all \((v, \tilde{v}) \in \dot{X}_n \times \dot{X}_n\), \[ \int_\Omega(\dot{u}v + \dot{u}_xv_x) = \int_\Omega(\tilde{w}v_x + \tilde{w}_xv_{xx}), \\ \int_\Omega(\tilde{w}\tilde{v} + \tilde{w}_x\tilde{v}_x) = \int_\Omega\left(u + \frac{1}{2}u^2\right)\tilde{v}_x. \] The introduced auxiliary variable \(\tilde{w}\) as introduced by our framework approximates \((1 - \partial_x^2)^{-1}\!\left[u + \frac{1}{2}u^2\right]\!\). We can see this scheme conserves \(H\) over \(T_n\) by taking \((v, \tilde{v}) = (\tilde{w}, \dot{u})\).
Again, the above video is at time \(20000\) with a moving camera of speed \(\frac{1 + \sqrt{5}}{2}\).

Crucially however, our framework extends beyond Hamiltonian systems, and beyond conservation laws. For instance, we use it to develop numerical schemes for the compressible Navier–Stokes equations that:

Conserve mass, momentum, and energy
Increase total entropy

FULL INTEGRATOR SPECIFICATIONS

Honestly, this is way too complicated for a little box on a website, sorry! Check out Section 5 of the preprint for the full details.

Further applications of the framework can be found in our preprint, and we are actively pursuing many more at the moment!
The framework is general and powerful. If you are investigating any type of transient system, we hope our work can provide a simple approach for generating more physically realistic simulations.

We would both gladly discuss it further!

Co-authors

Patrick Farrell

Affiliation | University of Oxford
Joint publications | High-order conservative and accurately dissipative numerical integrators via auxiliary variables, High-order structure-preserving discretisation for ideal MHD arising in the Parker problem (Upcoming), An augmented Lagrangian preconditioner for natural convection at high Reynolds number (Upcoming), Conservative integrators exhibit greater stability than symplectic integrators on the Toda lattice (Upcoming)

Talks

2025

Invited talk, Brown Unversity
Numerical Mathematics & Scientific Computing seminar, Rice University

2024

⬆️ UPCOMING | PAST ⬇️

External ("tiny desk") seminar, Rice University
Computing Division technical meeting, UKAEA
Firedrake User Meeting 2024, University of Oxford
PDEsoft, University of Cambridge
Finite Element Fair, University College London (UCL)
Exploiting Algebraic and Geometric Structure in Time-integration Methods workshop, University of Pisa
UKAEA PhD student engagement day, UKAEA
Junior Applied Mathematics Seminar, University of Warwick

2023

ICIAM 2023, Waseda University
Numerical analysis group internal seminar, University of Oxford
Junior Applied Mathematics Seminar, University of Oxford
Met Office presentation, University of Oxford