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

Boris Andrews | Patrick Farrell

16.JUL.2024 (arXiv)

CHECK OUT ON 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. However, there remain difficulties in devising geometric numerical integrators that preserve dissipation laws and conserve non-quadratic invariants.

In this work, we propose a framework for the construction of timestepping schemes that preserve dissipation laws and conserve multiple general invariants. The framework employs finite elements in time and systematically introduces auxiliary variables; it extends to arbitrary order in time. We demonstrate the ideas by devising novel integrators that conserve (to machine precision) all known invariants of general conservative ODEs, energy-conserving finite-element discretisations of general Hamiltonian PDEs, and finite-element schemes for the compressible Navier-Stokes equations that conserve mass, momentum, and energy, and provably possess non-decreasing entropy. The approach generalises and unifies several existing ideas in the literature, including Gauss methods, the framework of Cohen & Hairer, and the energy- and helicity-conserving scheme of Rebholz.

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)

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.

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.

CHECK OUT ON ARXIV!

We would both gladly discuss it further!

• boris.andrews@maths.ox.ac.uk
• patrick.farrell@maths.ox.ac.uk

RELATED WORKS

For a neat application of these ideas to a problem in magnetic relaxation that really highlights their importance, check out my subsequent work with Mingdong He, Patrick Farrell & Kaibo Hu, on structure-preserving integrators for the magneto-frictional equations.

RELATED OPEN PROBLEMS

Problem	Reward
Conservative integrators for PDEs with arbitrarily many invariants	★★★★★
Dissipative integrators for ODEs with arbitrarily many dissipated quantities	★★★★★
Discrete maximum principles via the auxiliary variable framework	★★★★★
Connections between the latent variable proximal point algorithm and my auxiliary variable framework	★★★★☆
Extension of the auxiliary variable framework to SDEs	★★★★☆
Extension of the auxiliary variable framework to general adiabatic invariants	★★★☆☆
Superconvergence of the auxiliary variable framework	★★★☆☆
Extension of the auxiliary variable framework to time discretisations beyond collocation RK methods	★★★☆☆
Application of the auxiliary variable framework to a viscoelastic fluid system	★★☆☆☆
Application of the auxiliary variable framework to delay DEs	★★☆☆☆
Conservative integrators for Hamiltonian systems in Lie groups	★★☆☆☆
Stable reduced order models from the auxiliary variable framework	★☆☆☆☆
Structure-preserving integrators for compressible MHD	★☆☆☆☆

CO-AUTHORS

Patrick Farrell

• Affiliation | University of Oxford
• Joint publications | High-order conservative and accurately dissipative numerical integrators via auxiliary variables / Topology-preserving discretization for the magneto-frictional equations arising in the Parker conjecture / An augmented Lagrangian preconditioner for natural convection at high Reynolds number (In preparation)

TALKS

2025

• Biennial Numerical Analysis Conference, University of Strathclyde
• Numerical Analysis Group internal seminar, University of Oxford
⬆️ UPCOMING | PAST ⬇️
• Numerical Mathematics & Scientific Computing seminar, Rice University
• SIAM CSE, Fort Worth, Texas
• Scientific Computing seminar, Brown Unversity

2024

• External ("tiny desk") seminar, Rice University
• Computing Division technical meeting, UKAEA
• Firedrake User Meeting, 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, Waseda University
• Numerical Analysis Group internal seminar, University of Oxford
• Junior Applied Mathematics Seminar, University of Oxford
• Met Office presentation, University of Oxford