ABSTRACT. The Geodynamic ADjoint Optimisation PlaTform (G-ADOPT) aims to develop and support a transformational new computational modelling framework for inverse geodynamics. It builds on several recent breakthroughs: (i) an ongoing surge in accessible observational datasets that constrain the structure, dynamics, and evolution of Earth’s mantle; (ii) advances in inversion methods, using sophisticated adjoint techniques, that provide a mechanism for fusing these observations with mantle models; and (iii) two novel software libraries, Firedrake and dolfin-adjoint, which when combined, provide a state-of-the-art finite element platform, with a fully-automated adjoint system, that offers a radical new approach for rigorously integrating geoscientific data (and uncertainties) with multi-resolution, time-dependent geodynamical models, through high-performance computing.

In the initial stage of this project we have developed a suite of models of increasing complexity, based on well established benchmarks, that tests and validates the forward modelling capabilities for this project. The results of this have just appeared in the Firedrake special issue of GMD: https://doi.org/10.5194/gmd-15-5127-2022. At present, we are extending the series of geodynamical examples presented in this paper, to facilitate multi-material simulations, which are used extensively by the geodynamical modelling community. Concurrently, we are exploring the use of new iterative solution strategies for the Stokes system, to improve computational efficiency: a major benefit of Firedrake is access to the wide variety of solution algorithms and preconditioning strategies provided by the PETSc library in combination with the ideas of, and infrastructure for composable block preconditioning techniques, scalable augmented Lagrangian and multigrid approaches developed in the Firedrake community.

The focus of current development, however, is the validation, verification and optimisation of adjoint models, under different physical approximations. As part of this effort, checkpointing functionality is being developed, both to address the memory footprint of adjoint calculations through disk writing strategies in dolfin-adjoint, as well as for restart functionality in the outer, optimisation loop through improvements in the PyRoL2 python interface to the Rapid Optimization Library.

Finally, the design of our platform has been guided by the desire to be able to rapidly transfer the developed technologies to other related fields in geodynamics. As such we will be extending the framework to enable the modelling of Glacial Isostatic Adjustment (GIA), the response of Earth’s surface and global sea level to melting polar ice sheets, in a new project associated with the Australian Centre of Excellence for Antarctic Science (ACEAS).

ABSTRACT. We develop a structure preserving mesh-adaptive finite elements method to simulate Maxwell’s eigenvalue problem, leveraging a subdivision-based construction of k-form basis functions. Subdivision is an algorithm that iteratively refines the mesh and obtains refined functions on the mesh by taking weighted linear combinations. To achieve efficient computations, we use the subdivision library OpenSubdiv to carry out the refinement and interpolate the resulting k-form into the lowest order k-form spaces over the refined mesh. We implemented an interface that makes these interpolants available in FEniCS and intend to do so for Firedrake. The resulting finite element spaces of subdivision k-forms preserve a finite-dimensional de Rham complex under subdivision; their elements are refinable, i.e., every coarse basis function exactly decomposes into fine basis functions using the weights given by the subdivision scheme. For this reason, basis functions constructed by subdivision open up an avenue towards adaptive, de Rham complex preserving simulations through local refinement. As a test case, we study the Maxwell eigenvalue problem and compare the convergence rates of the eigenvalues with standard H(curl) finite elements.

ABSTRACT. The Onsager framework for linear irreversible thermodynamics provides a thermodynamically consistent model of mass transport in a phase consisting of multiple species, via the Stefan–Maxwell equations, but a complete description of the overall transport problem necessitates also solving the momentum equations for the flow velocity of the medium. We derive a novel nonlinear variational formulation of this coupling, called the (Navier–)Stokes–Onsager–Stefan–Maxwell system, which governs molecular diffusion and convection within a non-ideal, single-phase fluid composed of multiple species, in the regime of low Reynolds number in the steady state. We propose an appropriate Picard linearisation posed in a novel Sobolev space relating to the diffusional driving forces, and prove convergence of a structure-preserving finite element discretisation. This represents some of the first rigorous numerics for the coupling of multicomponent molecular diffusion with compressible convective flow. The broad applicability of our theory is illustrated with simulations of the centrifugal separation of noble gases and the microfluidic mixing of hydrocarbons.

ABSTRACT. In optimization problems with control variables living in infinite dimensional spaces, the choice of inner product on the control space affects the definition of the gradient, which is the direction of steepest descent with respect to the norm induced by the inner product.

Shape optimization is not an exception, with the control variable living in an infinite dimensional space of vector fields. In this talk, we will study the impact of endowing this space of vector fields with an H^2 inner product on performing shape optimization.

ABSTRACT. UFL provides an excellent description of variational problems, capturing a wide range of finite element discretizations. However, because it lacks abstractions for time dependence, users have typically hand-written their own low-order time-stepping loops. We claim this is a major mismatch, and an opportunity for innovation.

Irksome provides a way around this -- it allows description of time-dependent problems and generates UFL for problems to be solved for each time step of a Runge-Kutta method. Working in terms of general Butcher tableaux, even fully implicit RK methods are supported. While these are traditionally eschewed in practice, we will advocate for their adoption in practice since i) they can be automatically generated ii) they typically have the "best" theoretical properties for a given application and iii) recent advances in preconditioners, plus Firedrake's solver infrastructure, means the solutions can be obtained very efficiently.

We will describe many of options now available in Irksome and show some numerical examples

ABSTRACT. I will start by describing some of the ongoing work within the Dedalus project, an effort that is comparable to Firedrake in many ways. In particular, I will discuss some of the design choices and application motivations that initiated the project over ten years ago. I will also mention some lessons we’ve learned along the way.

The Dedalus project maintains a similar philosophy to Firedrake in we seek to reduce barriers to flexible, accurate and efficient scientific computing. However, Dedalus achieves this aim partially through one major constraint: we focus on geometric domains with significant symmetry and regularity. While ostensibly limiting, we have often found many surprisingly suitable and pragmatic workarounds. However, the upside is that highly symmetric geometries allow powerful global spectral methods with calculus founded on group representation theory. To this end, I will briefly discuss a small fraction of the wonderful mathematics at the core of the Dedalus framework.

ABSTRACT. In this talk we are presenting a new approach for compatible finite element discretisations for atmospheric flows on a terrain following mesh. In classical compatible finite element discretisations, the H(div)-velocity space involves the application of Piola transforms when mapping from a reference element to the physical element in order to guarantee normal continuity. In the case of a terrain following mesh, this causes an undesired coupling of the horizontal and vertical velocity components. We are proposing a new finite element space, that drops the Piola transform. For solving the equations we introduce a hybridisable formulation with trace variables supported on horizontal cell faces in order to enforce the normal continuity of the velocity in the solution. Alongside the discrete formulation for various fluid equations we discuss solver approaches that are compatible with them and present our latest numerical results.

ABSTRACT. The Boussinesq approximation describes a fluid that is driven by buoyancy forces arising from differences in temperature. For the velocity field often times no-slip or free-slip boundary conditions are used. Here Navier-slip boundary conditions, which interpolate between these two cases, are used. In this talk we first mention the mathematical properties of the solutions, including the heat transport scaling in terms of the Rayleigh number, describing the strength of the buoyancy force. We then discuss how to simulate the problem with Firedrake, which of the previously derived properties are visible in the simulations and how to solve the obstacles arising from the boundary conditions.

ABSTRACT. We introduce a stochastic representation of the rotating shallow water equations and a corresponding structure preserving finite element discretization in Firedrake. The stochastic flow model follows from using a stochastic transport principle and a decomposition of the fluid flow into a large-scale component and a noise term that models the unresolved flow components. Similarly to the deterministic case, this stochastic model (denoted as modeling under location uncertainty (LU)) conserves the global energy of any realization. Consequently, it permits us to generate an ensemble of physically relevant random simulations with a good trade-off between the representation of the model error and the ensemble's spread. Applying a compatible finite element discretisation of the deterministic part of the equations combined with a standard weak finite element discretization of the stochastic terms, the resulting stochastic scheme preserves (spatially) the total energy. To address the enstrophy accumulation at the grid scale, we applied an anticipated potential vorticity method (AVPM) to stabilize the stochastic scheme. Using this setup, we compare different realizations of noise parametrizations in the context of geophysical flow phenomena and study potential pathways to fully energy preserving stochastic discretizations.

ABSTRACT. Standard finite element discretisations of the vector transport equation with lowest-order HDiv elements have a low order of accuracy. This talk will present two methods that do have higher-order accuracy when the equation is solved on curved 2D manifolds.

The first method uses a recovered finite element method, in which the transported field is reconstructed in a higher-order function space to be transported. The second scheme applies a mixed finite element formulation to the vector transport equation, simultaneously solving for the transported field and its vorticity. An approach to stabilising this mixed vector-vorticity formulation is presented that uses a Streamline Upwind Petrov-Galerkin (SUPG) method.

The accuracy of the methods will be illustrated through a new convergence test for vector-valued fields on the surface of a sphere.

ABSTRACT. The challenge with solving the indefinite Helmholtz equation is that standard approaches like GMRES with multigrid preconditioning do not work well. One way of achieving robust k-independent convergence is by using the shifted Helmholtz equation as preconditioner. In this talk, I will present shifted Hermitian/skew-Hermitian Splitting (HSS) preconditioners for the primal and mixed formulations of this problem and show some numerical examples. The result is a preconditioning strategy which is robust in mesh resolution and suitable in parallel.

ABSTRACT. Spectral coarse spaces are instrumental for the efficiency of domain decomposition methods making them one of the most efficient preconditioners for solving linear systems. PCHPDDM, a preconditioning package accessible through PETSc, provides a variety of spectral coarse spaces suitable for a wide range of matrices, especially those arising from the discretisation of PDEs. Most of these coarse spaces can be constructed in a completely algebraic way — no information is required on the problem behind them. In this talk, I will give an overview of PCHPDDM and its available functionalities.

ABSTRACT. In many applications of interest - such as magnetic confinement fusion - strongly anisotropic diffusion operators have to be solved with a high degree of accuracy in order to avoid excessive pollution across the direction of anisotropy. One way of doing so is by formulating the operator as a mixed system, with an auxilary variable denoting the diffused field's derivative in the direction of anisotropy.

Unfortunately, the highly anisotropic nature of the problem makes it difficult to be solved with standard multigrid methods. In this talk, we present a novel solver technique for resolving this issue for problems where all field lines defining the direction of anisotropy are open. By exploiting the fact that the auxiliary variable setup splits the diffusion formulation into two transport operators, we apply AIR - an algebraic multigrid solver well adapted for transport problems - to the mixed formulation. Here we highlight the issues and the formulation of the anisotropic diffusion operator, identify the novel solver technique and demonstrate the solver's effectiveness on some test problems implemented in Firedrake.

ABSTRACT. The Riesz maps of the $L^2$ de Rham complex frequently arise as subproblems in the construction of fast preconditioners for more complicated problems. In this work we present multigrid solvers for high-order finite element discretizations of these Riesz maps with optimal complexity in polynomial degree, i.e. with the same time and space complexity as sum-factorized operator application.

The key idea of our approach is to build new finite elements for each space in the de Rham complex with orthogonality properties in both the $L^2$ and $H(\mathrm{d})$ inner products ($\mathrm{d} \in \{\mathrm{grad}, \mathrm{curl}, \mathrm{div}\})$ on the reference hexahedron. The resulting sparsity enables the fast solution of the patch problems arising in the Arnold--Falk--Winther and Hiptmair space decompositions, in the separable case. In the non-separable case, the method can be applied to a spectrally-equivalent auxiliary operator.

With exact Cholesky factorizations of the sparse patch problems, the application complexity is optimal but the setup costs and storage are not. We overcome this with the use of incomplete Cholesky factorizations with carefully specified sparsity patterns arising from static condensation. This yields multigrid relaxations with computational complexity and storage that are both optimal in the polynomial degree.