ABSTRACT. With the growing complexity of modern computing hardware, the composition of appropriate domain-specific abstractions is becoming increasingly important to achieving the vaunted `three Ps' of performance, portability and productivity. Here, we present preliminary work on pyop3, a new domain-specific language embedded in Python for automating mesh-based simulations. The primary contribution of pyop3 is the introduction of a hierarchical data layout language that associates topological information with the locations of the degrees-of-freedom. This additional information raises the level of abstraction and allows pyop3 to perform advanced reasoning about mesh data layouts. In particular, we hope to provide Firedrake with the tools it needs to handle, among other things, p-adaptivity and unstructured hexahedral meshes.

ABSTRACT. Memory managed languages such as Python and Julia are being used increasingly often on high performance computers (HPC) to drive extremely large, pre-exascale simulations. Advantages of working in a managed language include the flexibility to change discretisations and solvers for better hardware utilisation, as well as being more productive for scientists and engineers. However, one drawback is that the memory management of such a language can create havoc when attempting to clean up distributed objects.

When running Firedrake in parallel, it is possible for Python's garbage collector to cause a deadlock when attempting to clean up distributed PETSc data structures that require collective destruction. Turning off Python's garbage collection is a poor workaround at best and a catastrophic memory leak at worst.

We outline an algorithm for the safe parallel destruction of distributed objects that can be used in any managed language and assess its impact on performance using different HPC facilities using a reference implementation in PETSc/petsc4py. Furthermore, we discuss some of the new very large scale simulations written using Firedrake that these improvements allow us to run.

ABSTRACT. NETGEN is an automatic mesh generator, led by Joachim Schoberl. In this talk, we introduce the recently developed NETGEN-Firedrake integration. This integration allows the user to define the geometry of their problem in the same Python environment used to implement the desired FEM solver with Firedrake; the CAD and UFL can live in the same place. Some of the features to be presented include describing 2D and 3D geometries using NETGEN's Constructive Solid Geometry features; converting NETGEN meshes to Firedrake meshes; and adaptive mesh refinement that conforms to the CAD geometry. We also describe ongoing work on parallelising the integration between Firedrake and NETGEN.

ABSTRACT. One important aspect of scientific and engineering simulations is scalable input/output of solutions. In this work we enhanced HDF5 input/output capabilities of Firedrake and PETSc to allow for saving/loading finite element solutions in association with the mesh of the computational domain efficiently in parallel. Particularly, our new implementation allows for using different number of MPI processes for saving and for loading. Transparent interfaces for mesh extrusion and timestepping problems are also provided.

ABSTRACT. We present numerical simulations for an idealized coupled ocean-atmosphere climate model. Our climate model belongs to the class of intermediate coupled models which are much simpler than the coupled general circulation models of the ocean-atmosphere system but still allow to study the fundamental aspects of ocean-atmosphere interactions. Our model couples an atmosphere system, described by the compressible two-dimensional (2D) Navier-Stokes equations and an advection-diffusion equation for temperature, to an ocean system, given by 2D incompressible Navier-Stokes equations and an advection-diffusion equation for temperature. Finite element method is used to discretize the system of PDEs representing the climate model on a 2D periodic domain and the discrete model is solved using Firedrake, which is an efficient automated finite element method library. To ensure the accuracy of simulation results of the coupled model, we have carried out detailed numerical investigation of its atmosphere and ocean components separately and tested our codes against some benchmark problems available in the literature. Our final goal is to incorporate stochasticity into the coupled ocean-atmosphere model following the Hasselmann's paradigm and use the model to study key features of climate phenomena like El-Nino Southern Oscillation (ENSO) which occur in the tropical eastern Pacific Ocean as a result of strong ocean-atmosphere interactions. This will be the subject of our future work.

ABSTRACT. Accurate modelling of basal melting beneath ice shelves is key to reducing the uncertainty in forecasts of ice-shelf stability and, thus, the Antarctic contribution to sea level rise. However, the lack of flexibility inherent to traditional ocean models can pose problems. Obtaining accurate melt estimates requires capturing the turbulent exchange of momentum, heat and salt at the ice-ocean interface, which may be modulated by the competing effects of stratification and basal slope. There are still significant uncertainties surrounding the trade-off between the simplicity of the melt parameterisation and the processes that need to be resolved by the numerical ocean model near the boundary. Real ice-shelf cavity geometries are complicated. Bathymetric valleys are common and provide pathways for warm circumpolar deep water. The ice base is marked by channels, crevasses and terraces. These features will affect the boundary flow, with an added complication that melting plays a role in their formation. It is very difficult to model such flow regimes using a traditional ocean model not only because of the resolution constraints imposed by inflexible grids, but also due to the inbuilt assumptions of large aspect ratio processes and domains that may be violated when flow occurs past these features. To model these effects requires the use of an ocean model that contains a flexible, unstructured mesh, is applicable at a range of length scales and, crucially, is still valid when the vertical-to-horizontal grid aspect ratio approaches order one. We are developing such a model for simulating flow under ice shelf cavities using the Firedrake finite element framework. We present our 3d simulations of ISOMIP+ experiments alongside simulations using the MITgcm ocean model, a commonly used z-layer (constant vertical resolution) model and investigate the trade-offs between spurious numerical mixing and melt rate convergence. We also show a new MMS test that includes melt boundary conditions. Finally, we will highlight some investigations into the dynamics of ocean flow within basal crevasses beneath ice shelves.

ABSTRACT. We have employed Firedrake to simulate extreme water waves, defined as waves with heights of more than twice the significant wave height. While these waves tend to be unusually large, unpredictable and rare, we have focused on deterministic modelling of such waves [1]. One theory about extreme waves, using analytical solutions of the Kadomtsev-Petviashvilli (KP) equation, reveals how three-line solitons can interact to produce waves with up to nine times the amplitude of the far-field solitons. We have seeded variationally consistent finite-element discretisations of Benney-Luke equations (BLE) and potential-flow equations (PFE) with KP-solutions (well prior to reaching the maximum) as initial conditions, in order to assess this maximum amplification in more realistic models. Another theory of extreme waves uses the (nonlinear) coalescence of dispersive waves. Numerical challenges we have hitherto overcome, to be presented, with Firedrake have been the following [1,2]: p-refinement studies to deal with the narrow wave peak, numerical conservation in space and time of the wave amplitude, validation against wave-tank experiments, dealing with non-autonomous aspects within a variational framework, and solving the hyperbolic-elliptic PFE by using CGFEM in the horizontal and higher-order p-refinement over one element in the vertical. Numerical challenges to be discussed and overcome are the following: solver optimization such that PFE-solvers can possibly start to match the many (commercial) PFE solvers, (r-)mesh-refinement in an x-periodic domain, validation against wave-tank experiments driven by wave flaps (leading to strong non-autonomous aspects), and using “fderivative” to automatically generate efficient versions of the algebraically complicated weak formulations from (time-discrete) variational principles of the wave dynamics, including elimination of the need for p-refinement of the PFE over one element in the vertical.

[1] Cho et al. 2022: Numerical experiments on extreme waves through oblique–soliton interactions. https://link.springer.com/article/10.1007/s42286-022-00059-3

[2] Gidel. Lu et al. 2022: Variational and numerical modelling strategies for cost-effective simulations of driven free-surface waves. https://eartharxiv.org/repository/view/3411/

ABSTRACT. This talk will show some recent work I've been doing on mesh generation. Many algorithms in computational geometry are built as a sequence of local transformations to the topology of the mesh. The simplest of these is flipping the edge between two triangles for computing Delaunay triangulations. The transformations can become much more involved, especially in 3D. Implementing them is famously annoying and error-prone using either of the common data structures used to represent mesh topologies in applications. Here I'll show how representing the mesh topology using the boundary operators on a chain complex makes it easy to implement one of the most challenging transformations: splitting a union of polytopes (in arbitrary dimension!) on a vertex.