Geophysical simulations with dynamic adaptivity

My work as an ASP postdoc is to participate to the development of a new generation of models for the simulation of oceanic and atmospheric flows. The goal is to develop a set of numerical techniques and demonstrate their potential for solving multiscale geophysical flows in an efficient manner.

I currently focus on the Shallow Water Equations that are discretized using the Discontinuous Galerkin method. This method takes advantage of the potential of unstructured meshes to increase the resolution exactly when and where it is needed, giving rise to multi-scale/physics numerical simulations. With this approach, a wider range of space- and time-scales of motion can be represented without having recourse to the brute-force approach consisting in enhancing the resolution homogeneously.

The other key aspect of the model is dynamic adaptation: the mesh and the order of interpolation are adapted dynamically during runtime to minimize the discretization error and capture the key aspects of the flows. The mesh is refined using the Adaptive Mesh Refinement (AMR) technique: it is modified by splitting elements recursively into finer sub-meshes, until the desired level of refinement is reached. A local conservative projection of the data is performed when the order of interpolation changes or when the mesh is adapted to transfer the solution form the old representation to the new adapted representation.

In order to take advantage of parallel computing, the mesh is partitioned into sub-domains attributed to each processors. The high level of locality of the method (few communications between elements lying on different processors) makes it very efficient on large parallel computers. Load balancing is performed to equilibrate the work of each processor that can be modified by adaptation steps.

Application to a global tsunami

The model has been applied to compute the propagation of the 2010 Chilean tsunami through the global ocean. The bottom topography of the sea is very steep, and the water depth can vary from 9000 m to 0 m in a single element. This application is then a good testcase to check the stability of the model when applied to tough realistic configurations.

As long as the tsunami propagates through the Pacific Ocean, the order of interpolation and the mesh are adapted to track the wave. The computational power is used effectively by concentrating the load where it is needed (front of the wave). The resolution remains however high in the shallow areas along the earthquake initial uplift, where the propagation of the wave is slow due to the very low depth.

Propagation of the wave. UP: Free-surface elevation [m]. DOWN: State of the mesh with order of interpolation.

The elevation of the free-surface has been compared with the DART data from the NOAA Center for Tsunami Research at eight different stations. It is seen that the model estimates accurately the time at which the tsunami reaches the different stations. The amplitude of the waves is also well predicted, except for the stations 5 and 8. Those stations being located in very shallow areas with an irregular an bathymetry and several small islands, the resolution of the model is probably not sufficient to reproduce accurately the height of the waves.

Free-surface elevation: model data (blue) and DART data (red). The different plot boxes are aligned with the timeline to indicate the time at which the tsunami reaches the different stations (t=0 at the moment of the initial earthquake).