ABSTRACT. We provide an iterative solution approach for the indefinite Helmholtz equation discretised using finite elements, based upon a Hermitian Skew-Hermitian Splitting (HSS) iteration applied to the shifted operator, and prove that the iteration is k- and mesh-robust when O(k) HSS iterations are performed. The HSS iterations involve solving a shifted operator that is suitable for approximation by multigrid using standard smoothers and transfer operators, and hence we can use O(N) parallel processors in a high performance computing implementation, where N is the total number of degrees of freedom. We argue that the algorithm converges in O(k) wallclock time when within the range of scalability of the multigrid. We provide numerical results verifying our proofs and demonstrating this claim, establishing a method that can make use of large scale high performance computing systems.

ABSTRACT. Although many problems in science and engineering are modelled by well-established PDEs, they often involve unknown or incomplete relationships, such as material constitutive laws or thermal response, that limit accuracy and generality. Existing surrogate-modelling approaches directly approximate PDE solutions but remain tied to a specific geometry, boundary conditions, and set of physical constraints. To address these limitations, we introduce a fully differentiable finite element-based machine learning (FEBML) framework that embeds trainable operators for unknown physics within a state-of-the-art, general FEM solver, enabling true end-to-end differentiation. At its core, FEBML represents each unknown operator as an encode-process-decode pipeline over finite-element degrees of freedom: field values are projected to nodal coefficients, transformed by a neural network, and then lifted back to a continuous FE function, ensuring the learned physics respects the variational structure. We demonstrate its versatility by recovering nonlinear stress-strain laws from laboratory tests, applying the learned model to a new mechanical scenario without retraining, and identifying temperature-dependent conductivity in transient heat flow.

ABSTRACT. 4D variational data assimilation (4DVar) can significantly improve the accuracy of numerical weather prediction by calculating the atmospheric state which best matches real world observations. However, the Gauss-Newton optimisation method used is computationally expensive, requiring the action of the adjoint, tangent linear, and Hessian models. This also makes 4DVar time-consuming to implement, so 4DVar research is often restricted to simple models.

We present new functionality being developed in Firedrake for automatically constructing 4DVar systems. Both strong constraint and weak constraint 4DVar are implemented, but we focus on weak constraint methods which allow for parallel-in-time evaluation of the forward and adjoint models. We show various preconditioners implemented using PETSc/TAO, including a recent saddle point reformulation which is amenable to Schur complement preconditioning, and demonstrate the flexibility for configuring the solvers.

This new functionality will enable exploration of 4DVar methods on a wide range of PDE and preconditioner combinations.

ABSTRACT. Fully implicit timestepping methods have several potential advantages for atmospheric and oceanic simulations. First, being unconditionally stable, they degrade more gracefully as the Courant number increases, typically requiring more solver iterations rather than suddenly blowing up. Second, particular choices of implicit timestepping methods can extend energy conservation properties of spatial discretisations to the fully discrete method. Third, these methods avoid issues related to splitting errors that can occur in some situations, and avoid the complexities of splitting methods.

Fully implicit timestepping methods have had limited application in geophysical fluid dynamics due to challenges of finding suitable iterative solvers, since the coupled treatment of advection prevents the standard elimination techniques. However, overlapping Additive Schwarz methods, as introduced for geophysical fluid dynamics by Cotter and Shipton (2023), provide a robust, scalable iterative approach for solving the monolithic coupled system for all fields and Runge-Kutta stages.

In this study, we investigate this approach applied to the rotating shallow water equations, facilitated by the Irksome package (Farrell et al, 2021), which provides automated code generation for implicit Runge-Kutta methods. We compare various schemes in terms of accuracy and efficiency using an implicit/explicit splitting method, namely the ARK2 scheme of Giraldo et al (2013), as a benchmark. This provides an initial look at whether implicit Runge-Kutta methods can be viable for atmosphere and ocean simulation.

ABSTRACT. We present a structure-preserving space-time finite element discretisation of the rotating shallow water equations via Poisson bracket. We develop a formulation which conserves energy for high-order time integrators, works as a standard space–time scheme. By preserving the underlying Hamiltonian structure, this geometric integrator ensures long-term stability for complex geophysical flows.

ABSTRACT. The typical energy estimate for the Navier--Stokes equations provides a bound for the gradient of the velocity; energy-stable numerical methods that preserve this estimate preserve this bound. However, the bound scales with the Reynolds number (Re) causing solutions to be numerically unstable (i.e. exhibit spurious oscillations) on under-resolved meshes. The dissipation of enstrophy on the other hand provides, in the transient 2D case, a bound for the gradient that is independent of Re.

We propose a finite-element integrator for the Navier-Stokes equations that preserves the evolution of both the energy and enstrophy, implying gradient bounds that are, in the 2D case, independent of Re. Our scheme is a mixed velocity-vorticity discretisation, making use of a discrete Stokes complex. While we introduce an auxiliary vorticity in the discretisation, the energy- and enstrophy-stability results both hold on the primal variable, the velocity; our scheme thus exhibits greater numerical stability at large Re than traditional methods.

We conclude with a demonstration of numerical results, and a discussion of the existence and uniqueness of solutions.

ABSTRACT. The ability to create submeshes out of existing meshes in Firedrake is a recent development <a href="https://pdesoft.org/talks/seamless-integration-of-the-submesh-feature-in-firedrake/?utm_source=chatgpt.com">[1]</a>.

The feature is particularly relevant for ongoing work around the modelling of a novel wave-energy device, in which a subsection of the hydrodnamic domain is in contact with a buoy, over which an equation is required to be solved.

The Submesh functionality is first investigated using simple 1D and 2D toy problems, before being implemented as part of a wider solution methodology for the device. It is hoped that utilising Submesh will enable basis functions with degree greater than 1 to be used - something that the current methodology does not handle properly.

ABSTRACT. Tightly coupled multi-domain problems are crucial for modelling complex physical systems. These problems can be formulated as mixed systems, and admit a natural block matrix structure. The off-diagonal blocks will contain the interactions between the domains, and these are calculated by cross-mesh interpolation. A barrier for implementing programming abstractions for these problems in Firedrake has been the difficulty in assembling matrices representing cross-mesh interpolation.

This talk presents how cross-mesh interpolation is implemented in Firedrake, how the interpolation matrix is constructed, the goals we're aiming for in two-way coupling, and possible application areas.

ABSTRACT. We investigate the coupling of nonlinear water-wave motion to buoy dynamics in the presence of an inequality constraint. Building on augmented Lagrangian variational principles (VPs) developed by Burman and others, we impose constraints of the form G(q)≥0, where q are system variables, through a Lagrange multiplier L. The strict Kuhn–Karush–Tucker (KKT) conditions {L·G=0,G(q)≥0,L≤0} are replaced by those smooth approximations of the involved function F(cG(q)−L)=max(cG(q)−L,0), with constant c>0, which allow explicit computation of the multiplier L as (part of) a force. Our approach combines: (a) an Average Vector Field (AVF) time-stepping method, extended to water-wave systems with an auxiliary field, enforcing energy conservation in the discrete system; (b) a (novel?) smooth relation L(G) that regularises the KKT conditions by approximating the solution G=0 with L≤0 and G>0 with L=0 in the (L,G)-plane, but leading to an implicit definition of F(cG(q)−L). This framework has been implemented in Firedrake, leading to improved benchmark problems: (i) a point particle under gravity bouncing off a rigid table, (ii) a particle moving in a rectangular (“billiard”) domain, and (iii) forced (Boussinesq-type) nonlinear waves in a horizontal channel heaving a buoy in a wave-enhancing contraction. The latter, finite-element, model supports design of a prototype wave-energy device for enhanced energy capture. More generally, this work aims to develop analytical and computational tools for finite-element coupling of nonlinear wave dynamics in fluid-structure interactions, here exemplified by the vertical (heave) motion of the buoy.