Dr Duru's research interest encompasses all components of the development of numerical algorithms and simulation tools for partial differential equations (PDEs):

• Convergent  volume discretisation in complex geometries using high order methods (finite difference methods, discontinuous Galerkin methods, finite element and spectral methods).

• Mathematical and numerical analysis of (hyperbolic) PDEs.

• Development and analysis of perfectly matched layers (PML)  – to absorb outgoing waves at domain boundaries.

• Efficient implementation of scalable simulation codes on modern, heterogeneous many-core HPC platforms.

• Numerical simulations of real-world wave propagation problems.

Some ongoing and future projects are listed below for prospective research/PhD students and postdocs.

High fidelity and extreme-scale methods for multi-scale wave propagation problems

Multi-physics/multi-scale problems, simultaneously coupling different physical phenomena at different scales in space and time, occur in many applications. One example is simulations of coupled elastic deformation and magma flows during volcanic eruptions, and the associated waves, as well as earthquake source process, coupling frictional failure on the fault to wave propagation off the fault, in the elastic medium. To effectively simulate such problems, advanced methods and hybrid solvers coupling different numerical methods with dynamic adaptivity are essential, as well as the availability of sufficient computational resources.

The main aims of this topic include:

• Develop provably stable hybrid solvers that couple DG methods to SBP FD methods for efficient numerical simulation of multi-scale and/or multi-physics wave propagation problems.

• Provide a rigorous analysis of a stable and high order accurate coupling with nontrivial boundary conditions for a class of DG, finite volume and FD approximations with efficient explicit or semi-explicit time discretisation.

• Implement the numerical methods in a scalable framework such as ExaHyPE or WaveQLab.

Perfectly matched layer (PML), absorbing layers and absorbing boundary conditions

The main features of propagating waves is that they can propagate long distances relative to their characteristic length-scale, the wavelength. For numerical simulations, it is precisely this essential feature of waves, the radiation of waves to far field, that leads to the greatest difficulties. The perfectly matched layer (PML) is a method which, when stable, can provide a domain truncation which is convergent with increasing layer width/damping. The difficulties in using the PML are primarily associated with stability, which can be present at the continuous level or be triggered by numerical approximations.  The mathematical and numerical analysis of the PML for hyperbolic wave propagation problems has been area of active research. It is now possible to construct stable and high order accurate numerical wave solvers using  summation-by-parts finite difference methods, continuous or discontinuous Galerkin finite element methods.

The main aim of this topic is to answer some of the following important questions:

• How can PML be adapted to the presence of surface and interface waves?

• Can we derive a stable PML for general heterogeneous isotropic or anisotropic elastic media?

• Can we derive a stable and effective absorbing layer for nonlinear problems, with examples from nonlinear aero-acoustics and general relativity?

Dynamic earthquake rupture simulations 

Earthquakes are catastrophic natural events that can take a heavy toll on our societies, through life and property losses. By improving our understanding on how earthquakes nucleate and propagate, we may make progress towards mitigating their negative impact. Insights from modelling and simulations may help in the development of better seismic hazard models as well as in the construction of more efficient earthquake early warning systems. However, earthquakes present us with some of the most difficult mathematical, numerical and computational challenges. One main objective under this topic is to help seismologists correctly model earthquake source processes to better quantify and constrain seismic hazard at frequencies relevant to engineering applications.

This topic will address some of these important mathematical and numerical challenges:

• Prove well/ill-posedness of the IBVP for a class of modern friction laws for crustal rocks.

• Derive stable numerical methods to efficiently solve quasi-static elasticity coupled to nonlinear friction laws and enable efficient simulation of earthquake sequences.

• Derive stable numerical methods to efficiently couple quasi-static elasticity to dynamic elasticity subject to nonlinear friction laws.

Space-time adaptive DG for Einstein’s equations in second order form.

Second order hyperbolic systems often describe problems where wave propagation is dominant. Typical examples are the acoustic wave equation, elastic wave equation, and the Einstein’s equation of general relativity. However, many solvers for wave equations and the Einstein’s equations are designed for first order systems. In particular multi-domain spectral methods which are increasingly becoming attractive because they are optimal, in terms of efficiency and accuracy, are commonly implemented as first-order systems. That is, the system of second order hyperbolic equations

are first reduced to a system of first order hyperbolic PDEs before numerical approximations are introduced. The main reason being that the theory and numerical methods to solve hyperbolic PDEs are well developed for first order hyperbolic systems, and less developed for second order systems.

There are many advantages to solving the equations in second order form. The reduction of the equation to first-order hyperbolic form, has the disadvantage of introducing auxiliary variables with their constraints and boundary conditions. For example, in the harmonic description of general relativity, the Einstein’s equations are a system of 10 curved space second order wave equations, while the corresponding reduction to first order systems will involve about 60 equations. The reduction to first-order form is also less attractive from a computational point of view considering the efficiency and accuracy of numerical approximation.

It has been a long held ambition of the computational relativity community to develop the theory and and provably stable numerical techniques for second order hyperbolic systems such that Einstein’s equations can be solved efficiently. This has proven to be an incredibly difficult task. In particular, it appears to be more difficult to guarantee stability for naturally second order systems than for first order reductions of them.

The overall goal in this project is to design a new class of spectrally accurate and robust numerical algorithms to solve Einstein’s equations in second order form, and simulate relativistic astrophysical wave phenomena at an unprecedented precision never achieved before.