ABSTRACT. The classical no-slip (or homogeneous Dirichlet) boundary condition for fluid flow is not appropriate in many situations, as fluids will often slip at solid walls. In this work we propose a theoretical framework that is able to capture a wide variety of slip models, including linear (Navier) slip, non-linear (and possibly non-smooth) slip described by monotone graphs, non-monotone slip, and dynamic (time-dependent) slip. A finite element scheme is proposed, in which the non-penetrability condition at the wall is enforced with a Nitsche formulation. Numerical experiments implemented in firedrake will also be presented.

ABSTRACT. Diffusion-driven multicomponent flows can be modelled well using the Onsager-Stefan-Maxwell (OSM) equations. We start by introducing a new weak formulation and discretisation for the OSM equations, and discuss how this can be preconditioned using geometric multigrid with Vanka patch smoothers. Time permitting, I will also discuss how this can be extended to the Stokes-Onsager-Stefan-Maxwell (SOSM) equations in which the bulk velocity is modelled using Stokes flow.

ABSTRACT. I have got totally fed up with splitting methods for atmosphere and ocean models and have been exploring the use of ASM smoothers/preconditioners for the whole monolithic system, inspired by the efforts of the Irksome team. The idea is to use a single unsplit implicit timestepping method for the full coupled system of velocity, density, temperature etc. The timestepping method becomes simpler at the expense of having to do something more sophisticated in the iterative solver.

Working entirely experimentally and unaided by theory, with ASM smoothers we have good multigrid behaviour for mixed Hdiv/L2 elliptic systems that arise in the linearisation of models around a state of rest. This is lost if we try to use them for the full Jacobian linearised about the current state arising in e.g. Crank-Nicholson time discretisation of a nonlinear geophysical fluid models. However, we don't care, because in those models we always want to keep the Courant number fixed, so dt decreases proportionally to dx. When we scale like that, we return to dx independent convergence rates.

In my talk I will briefly explain what ASM smoothers are, for the uninitiated. Then I will go through various examples I have been playing with recently: the rotating shallow water equations on the sphere, and the vertical slice compressible Euler equations (in work recently published in Cotter and Shipton (2023)). Then, it turns out that it is possible to apply the same approach to models with moisture parameterisations, where the monolithic system is enlarged with variables for water vapour, cloud water, rain, and other species of H20, so I will explain that as well. Finally if there is time I will discuss my ideas to use variational inequality solvers for moisture parameterisations, so that the constraints of maximum humidity (which depends on temperature and pressure etc) can be seamlessly incorporated into a monolithic implicit timestepping method instead of fixing up values at the end of each step as is currently done in weather and climate models.

ABSTRACT. Quasivariational inequalities (QVIs) are ubiquitous but, in particular, arise in PDE-constrained optimization in cases where the constraint set depends on the solution itself. In obstacle-type QVIs, this manifests as an obstacle that bends according to the state of the system. QVIs are notoriously hard to analyze, especially in the infinite-dimensional setting, and developing fast solvers posed in infinite dimensions has proven particularly challenging. As such most solvers in the literature rely on fixed point algorithms which can be slow to converge. In this talk, we introduce the first semismooth Newton method, posed in a Banach space setting, for such problems. We will see that the solver enjoys favourable properties such as local superlinear convergence and mesh independence.

ABSTRACT. The Grad–Shafranov equation (GSE) for axisymmetric magnetohydrodynamic (MHD) equilibrium is a nonlinear, scalar PDE which in principle can have zero, one or more non-trivial solutions. This is the equation used to find plasma equilibria in tokamak fusion experiments. The conditions for the existence of multiple solutions have been little explored in the literature so far. We develop a simple analytic model to calculate multiple solutions in the large aspect ratio limit. We then compare to numerical results using Firedrake to solve the GSE and DefCon, which implements the deflated continuation method, to find multiple solutions of the equation. The analytic model is surprisingly accurate in calculating multiple solutions of the GSE for given boundary conditions and the two methods agree well in limiting cases. We then examine the effect of plasma shaping and aspect ratio on the multiple solutions numerically and show that shaping generally does not alter the number of solutions. We discuss future work in the area of plasma equilibrium modelling.

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.

ABSTRACT. We introduce a family of finite element methods, known as recovered finite element methods (R-FEM), developed for the time-harmonic Maxwell equations. These methods utilize recovery operators that map discontinuous Lagrange spaces to H(curl)-conforming spaces. In this talk, we motivate and investigate the stability of R-FEM solutions at sufficiently large wavenumbers and explore the potential applications of R-FEM in various contexts.

ABSTRACT. Discretisation methods that are of high order in space and time are highly desirable when solving the equations of fluid dynamics since they show superlinear convergence and make efficient use of modern computer hardware.

We develop an efficient timestepping method for the incompressible Euler equations based on combining implicit-explicit (IMEX) time-integrators, high-order Hybridised Discontinuous Galerkin (HDG) discretisations and splitting methods. The two computational bottlenecks are the calculation of a tentative velocity and the solution of a velocity-pressure system to enforce incompressibility. The second solve is preconditioned with a non-nested hybridised multigrid preconditioner as proposed by Cockburn, Dubois, Gopalakrishnan and Tan for the Schur-complement system that arises from static condensation to eliminate the pressure/velocity unknowns in favour of variables on the mesh-facets. The tentative velocity solve is preconditioned with ILU0. With this choice of preconditioners, the number of solver iterations depends only weakly on the polynomial degree and the grid-spacing; in other words, empirically the method is approximately h- and p-robust. As a consequence, the cost of a single timestep is roughly proportional to the total number of unknowns. At high order this allows the computation of highly accurate solutions at a low computational cost.

The code has been implemented in Firedrake, using static condensation technology and PETSc to construct a fairly complex solver/preconditioner. Numerical results are presented for the simulation of a two-dimensional Taylor-Green vertex.

ABSTRACT. The viscous flow simulations utilizing Stokes equations is used extensively in the study of the deformation of the earth crust and mantle and flow of the glaciers. Using proper boundary conditions in these simulations is rarely discussed but is a source of difficulty, sometimes leading to unrealistic models in order to avoid boundary effects such as concentration of the stresses for poorly chosen boundaries due to excessive constraining of the velocity field. The existence of pressure nullspace in the Stokes equations is well known and well studied. The pressure nullspace arises when only Dirichlet boundary conditions exist. Pressure nullspace can be tackled with the elimination of the constant pressure mode from the Krylov solution space. When models have free top surfaces and open boundaries at the bottom of the domain in order to allow for inflow or outflow of the material, part of the boundary becomes an interface. This eliminates the pressure nullspace. However a new nullspace in the velocity field exists for this configuration when the velocity field perpendicular to the direction of compression or extension is not restrained at the boundary. A common practice is to remove the nullspace using a Lagrangian penalty method in conjunction with local mass matrix inversion which shifts the eigenvalues of the Picard matrix away from zero. We show that the weak form of the Stokes flow for the problems involving free surfaces should take into account the variation in the domain of the flow field. With the application of Reynolds transport theorem when taking directional Gateaux derivative, integral terms appear on the boundary facets that have the effect of making the problem well posed. Therefore Schur complement can be applied successfully without resorting to regularization such as penalty methods. We utilize Firedrake to easily and effectively implement these boundary integrals. As a testament to Firedrake’s versatility, the modifications required are terse and readily implemented. A benchmark problem involving Rayleigh-Taylor instability is tested against the analytic linearly perturbed solution. Also a slightly modified nonlinear benchmark problem from the geodynamics community is solved with the implications on the choice of boundary conditions in order to create realistic models.

ABSTRACT. We have devised a novel wave-to-wire wave-energy model. It consists of water-wave motion in a wave-enhancing contraction, constrained (vertical) buoy motion therein and an electro-magnetic generator. The water waves are modelled as potential flow. Buoy and generator are modelled by 4 to 8 ordinary differential equations. Both equality and inequality constraints have been used to couple wave and buoy motion. Submodel by submodel, I will show the ongoing build-up to a complete (nonlinear) finite element wave-to-wire model within Firedrake by using (time-discrete) variational principles augmented with the nontrivial damping of the electric circuits and energy-harvesting load.

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. In this talk we consider the Stokes–Onsager–Stefan–Maxwell (SOSM) equations, which model the flow of concentrated mixtures of distinct chemical species in a common thermodynamic phase. The equations account for both the diffusive interactions between chemical species and the bulk convection. Our aim is to develop computationally efficient high-order finite element schemes that discretize these nonlinear equations in two and three spatial dimensions. Because the SOSM equations relate many unknown variables (e.g. the bulk and species velocities, pressure, concentrations, chemical potentials, etc.), this is a difficult task. In particular, there are many choices of which variables should be explicitly solved for in the formulation, and it is not clear which discrete finite element function spaces should be employed. To tackle this challenge, we derive a novel weak formulation of the SOSM problem in which the species mass fluxes are treated as unknowns. We show that this new formulation naturally leads to a large class of high-order finite element discretizations that are straightforward to implement and have desirable linear-algebraic properties. From a theoretical standpoint, we are able to prove that when applied to a linearized version of the SOSM problem, the proposed finite element schemes are convergent. We present preliminary simulations on the microfluidic mixing of hydrocarbons.

ABSTRACT. We discuss integration of the submesh feature in Firedrake that enables one to extract parts of an existing mesh and make them work together to solve finite element problems. We focus on cell submeshes, in which topological dimensions of the submeshes are the same as that of the parent mesh. We briefly discuss the implementation and the user interface, and then show interesting applications including fluid-structure interaction benchmark problem [Turek, S., Hron, J. (2006). Proposal for Numerical Benchmarking of Fluid-Structure Interaction between an Elastic Object and Laminar Incompressible Flow. In: Bungartz, HJ., Schäfer, M. (eds) Fluid-Structure Interaction. Lecture Notes in Computational Science and Engineering, vol 53. Springer, Berlin, Heidelberg].

ABSTRACT. We present ngsPETSc, an interface between the NETGEN mesher, the finite element library

NGSolve and the Portable, Extensible Toolkit for Scientific computation (PETSc). ngsPETSc

interface is written in Python and takes advantage of the Python bindings of Netgen/NGSolve and petsc4py. The key components of ngsPETSc are the interface between NETGEN and PETSc DMPlex and the interface between NGSolve and PETSc PC, KSP and SNES.

In particular, thanks to the first component it is possible to export mesh generated from geometries described via constructive solid geometry (CSG) using OpenCASCADE as PETSc DMPlex objects. Based on this component new features have been added to Firedrake. In particular, now Firedrake supports linear and higher-order meshes that conform to the constructive solid geometry. Furthermore, ngsPETSc allows the construction of mesh hierarchies for multigrid solvers in two dimensions. Lastly, different mesh splits such as Alfeld and Powell-Sabin splits are supported in Firedrake.

The PETSc PC and KSP interface allows the use of PETSc solvers in NGSolve. In particular,

the PETSc PC interface also allows to use of the PETSc preconditioners directly in the linear algebra solvers implemented in NGSolve. Lastly, the PETSc SNES interface gives NGSolve users access to the full range of non-linear solvers present in PETSc. In particular, PETSc SNES can be used to solve complex non-linear problems and to solve the optimisation problems related to the energy of a particular system.

ABSTRACT. The accuracy of a numerical solution is strongly dependent on the underlying computational mesh. The goal of mesh adaptation is to improve the mesh-based discretization of the solution space to achieve better accuracy and stability of the numerical solution at a comparably diminished computational cost. $r$-adaptivity, or mesh movement, seeks to optimise the mesh through the movement of the mesh nodes only, without altering the underlying topology. Alternatively, $h$-adaptivity optimizes the mesh through re-gridding operations, which in two dimensions include node insertion, edge splitting, node deletion, edge collapse, and edge swap. Many unsteady partial differential equations (PDEs) at different points in the evolution of the numerical solution field can benefit from applying either $r$-adaptivity or $h$-adaptivity methods towards increasing robustness, flexibility, and computational cost reductions. We explore the strengths and limitations of alternating these methods utilizing a suite of well-known 2D benchmark PDE test cases, such as flow past a cylinder and a lid-driven cavity, in the context of metric-based, anisotropic mesh adaptation and focusing on issues of mesh quality, error measures, and computational efficiency.