ABSTRACT. Double-precision arithmetic has become standard in numerical codes. However, in real applications, uncertainties and model error are many orders of magnitude larger than the machine epsilon, and most of the precision is unnecessary. Notably, several operational weather prediction models now successfully use single-precision in some or all of their code. In an academic context, people have investigated using arbitrary levels of precision, not just single- and double-, even though language and hardware support is currently lacking. There are therefore a number of software emulation tools which emulate the effect of reduced-precision arithmetic. These typically work by declaring new types, or by extending the compiler. Here, I explore using one of these tools with the assembly kernels generated by Firedrake's form compilation toolchain, the effects on a low-resolution shallow-water simulation, and some simple ideas for increasing the resilience of finite element codes to reductions in numerical precision.

ABSTRACT. We are exploring ways of overcoming the timestep restrictions due to oscillatory stiffness in geophysical fluid models by introducing parallelism into the timestepping algorithm, following the ideas of Wingate and Haut (2014). One part of doing this is by replacing the full model with a phase-averaged version which filters oscillations, with the averaging being done by a discrete average over a set of phase values which can be calculated in parallel. The averaging error can then be corrected if necessary with Parareal or other time-parallel algorithms. The phase averaging requires application of the exponential of the wave operator (which has purely imaginary eigenvalues). One way of doing this is by using rational approximations, which result in solving many parallel systems of the form (a + ib)u + Lu = u_0, with large |b|. I will discuss the challenges of solving these systems and some preconditioner approaches for them, and will discuss progress towards a parallel phase-averaged model of rotating shallow water on the sphere.

ABSTRACT. Domain decomposition has been an effective method for parallelising many computational problems. However its scalability is limited, as the communication overhead increases the more the domain is sub-divided. Making effective use of future computing architectures will require methods for exploiting very high degrees of parallelism which will require looking beyond domain decomposition.

Fortunately it is possible to exploit parallelism in the time dimension, as well as in spatial dimensions. One approach to time parallelism is the use of rational approximations to exponential integrators (REXI). REXI enables a calculation to take a large time step by using sum of rational terms which can be computed in parallel. This talk will present initial results from applying the REXI method to the shallow water equations in Firedrake looking at accuracy and performance considerations.

ABSTRACT. We consider metric-based mesh adaptation methods for steady-state PDEs, solved using Firedrake. In this work, a number of adaptation methods are implemented in Firedrake, each seeking meshes upon which a scalar quantity of interest (QoI) may be accurately approximated.

Through the QoI we define adjoint equations, which convey information relating its sensitivities to aspects of the PDE solution. A goal-oriented mesh adaptation approach is considered, based on dual weighted residual type error estimation. strategy. Metrics giving rise to isotropic and anisotropic meshes are considered, each of which are able to achieve the same relative error in approximating the QoI as with uniform refinement, but using fewer elements.

For validation purposes, we compare QoI values resulting from these approaches against analytical values which may be extracted for a particular advection-diffusion based test case. Potential applications in desalination plant outfall modelling are discussed.

ABSTRACT. In this session we introduce a recently added feature of Firedrake to apply boundary conditions in equation form using `EquationBC` class. Designed to be a generalisation of DirichletBC class for applying Dirichlet boundary conditions, EquationBC class allows users to force degrees of freedom on boundary to satisfy some partial differential equations. This generalisation becomes useful when, e.g., solving problems that have distinct physics in the domain and on the boundary. We will illustrate the use of EquationBC using simple examples.

ABSTRACT. In many practical cases, PDEs are not enough to accurately describe the physical problem of interest. For example, it may be necessary to include features not represented in the differential equations, or to include closures for unresolved scales. While the PDEs themselves are often based on fundamental physical laws, these additional terms are often more heuristic in nature. They may be represented by a neural network and learned entirely from observed data, or they might be based on looser physical arguments about the system.

Here we present extensions to UFL and Firedrake which enable the inclusion of arbitrary external operators. These operators are not limited to vector calculus expressions written in UFL, but can instead be anything that can be computed from the state of the system. We also present the differentiation and assembly mechanisms that enable variational problems involving these new operators to be solved. We illustrate the new capabilities with two new classes of operator:

- Pointwise solve operator : This operator enables the inclusion of pointwise nonlinear relationships. This facilitates the modelling of implicit constitutive relations, i.e. nonlinear solids or fluids where the stress is computed by solving a pointwise nonlinear equation every time the operator is evaluated.

- Neural Network : The pointwise neural network operator provides values by evaluating a neural network whose inputs are other state variables. This presents the opportunity to create new models of phenomena by learning from experimental or computed data. This may be applicable in a wide range of areas where the theoretical basis for models is less strong. Potential examples include turbulence closures and regularisation terms in inverse problems.

ABSTRACT. The development of morphodynamic models to simulate sediment transport accurately is a challenging process that is becoming crucial with sea level rise and the potential increase in strength and frequency of storms due to a changing climate. These are complex models given the non-linear nature of the sediment transport problem. We implement a new depth-averaged coupled hydrodynamic and sediment transport model within the coastal ocean model Thetis, that is built using the code generating framework Firedrake which facilitates code flexibility and optimisation benefits. To the best of our knowledge, this represents the first full morphodynamic model including both bedload and suspended sediment transport which uses a discontinuous Galerkin based finite element discretisation.

We implement new functionalities within Thetis extending its existing capacity to model scalar transport to a capacity to model suspended sediment transport and incorporating within Thetis options to model bedload transport and bedlevel changes. We apply our model to problems with non-cohesive sediment and account for effects of gravity and helical flow by adding slope gradient terms and parametrising secondary currents. For validation purposes and to demonstrate the model capability, we present results from test cases of a migrating trench and a meandering channel comparing against experimental data and the widely-used model Telemac-Mascaret.

ABSTRACT. In this talk, we present a composable preconditioner interface for multigrid methods using a hierarchy of non-nested function spaces. In particular, we present an N+2 multigrid method with applications in hybridizable mixed and hybridizable discontinuous Galerkin (HDG) discretizations. This multigrid procedure starts with a two-level approach, with the finest problem defined on the space of discontinuous traces. The coarse grid solver is another N-level multigrid method on a different function space for a field which is spectrally equivalent to the trace variables. Our implementation uses the Python preconditioner interfaces provided by the Firedrake finite element library.

ABSTRACT. Authors: Jack Betteridge, Colin Cotter, Thomas Gibson, Ivan Graham, Lawrence Mitchell, Eike Müller

For problems in Numerical Weather Prediction (NWP), time to solution is a critical factor. Semi-implicit time-stepping methods can speed up geophysical fluid dynamics simulations by taking larger time-steps than explicit methods. This is possible because they treat the fast (but physically less important) waves implicitly, and the time-step size is not restricted by the CFL condition for these waves. One disadvantage of this method is that an expensive linear solve must be performed at every time step, however, using an effective preconditioner for an iterative method significantly reduces the computational cost of this solve, making a semi-implicit scheme faster overall.

Higher-order Discontinuous Galerkin (DG) methods are known for having high arithmetic intensity making them well suited for modern HPC hardware, but are difficult to precondition due to the large number of coupled degrees of freedom. This coupling arises since the numerical flux introduces off diagonal artificial diffusion terms. By using a hybridised DG method we can eliminate the original coupling and instead couple the equations to a smaller global system on the trace space, which is easier to precondition. This is achieved by considering the numerical flux variables which only lie on the facets of the mesh. Recent work by Kang, Giraldo and Bui-Thanh[1] solves the resultant system directly. However, this becomes impractical for high resolution simulations. We build on this work by solving the resultant system using a non-nested geometric multigrid technique instead[2].

We discretise and solve the non-linear shallow water equations, an important model system in geophysical fluid dynamics, and demonstrate the effectiveness of the multigrid preconditioner for a semi-implicit IMEX time-stepper. The method is implemented in the SLATE language, which is part of the Firedrake project.

References:

[1] Kang, Shinhoo and Giraldo, Francis X and Bui-Thanh, Tan. IMEX HDG-DG: a coupled implicit hybridized discontinuous Galerkin (HDG) and explicit discontinuous Galerkin (DG) approach for shallow water systems arXiv preprint arXiv:1711.02751, 2017

[2] Cockburn, Bernardo and Dubois, Olivier and Gopalakrishnan, Jay and Tan, Shuguang. Multigrid for an HDG method IMA Journal of Numerical Analysis 34(4):1386–1425, 2014

ABSTRACT. The Kolmogorov Superposition Theorem shows that you can exactly represent any continuous function as the composition of one-dimensional functions. The basis for the inner functions is universal in that it need not depend on the function $f$ being represented. The outer basis does depend on $f$, but appears to be purely interpolatory, and thus is a perfect fit for finite elements. We present preliminary approximation results using this framework with adaptive linear elements, and present a plan for incorporating approximations of this type into Firedrake.

ABSTRACT. Effective relaxation methods are necessary for good multigrid convergence. For many equations, standard Jacobi and Gau\ss--Seidel are inadequate, and more sophisticated space decompositions are required. Examples include problems with semidefinite terms or saddle point structure.

In this talk, I will present PCPATCH, which offers a uniform and flexible interface for the topological specification of space decompositions for multigrid relaxation methods.

Implemented in PETSc, and accessed through Firedrake's Python preconditioners, PCPATCH facilitates access to a wide range of relaxation schemes by simply changing solver options. More complex patches can be built by writing a short piece of Python code. I will illustrate the utility of this approach with some examples.