Search

In this paper we present a new high order semi-implicit DG scheme on two-dimensional staggered triangular meshes applied to different nonlinear systems of hyperbolic conservation laws such as advection-diffusion models, incompressible Navier-Stokes equations and natural convection problems. While the temperature and pressure field are defined on a triangular main grid, the velocity field is defined on a quadrilateral edge-based staggered mesh. A semi-implicit time discretization is proposed, which separates slow and fast time scales by treating them explicitly and implicitly, respectively. The nonlinear convection terms are evolved explicitly using a semi-Lagrangian approach, whereas we consider an implicit discretization for the diffusion terms and the pressure contribution. High-order of accuracy in time is achieved using a new flexible and general framework of IMplicit-EXplicit (IMEX) Runge-Kutta schemes specifically designed to operate with semi-Lagrangian methods. To improve the efficiency in the computation of the DG divergence operator and the mass matrix, we propose to approximate the numerical solution with a less regular polynomial space on the edge-based mesh, which is defined on two sub-triangles that split the staggered quadrilateral elements. Due to the implicit treatment of the fast scale terms, the resulting numerical scheme is unconditionally stable for the considered governing equations. Contrarily to a genuinely space-time discontinuous-Galerkin scheme, the IMEX discretization permits to preserve the symmetry and the positive semi-definiteness of the arising linear system for the pressure that can be solved at the aid of an efficient matrix-free implementation of the conjugate gradient method. We present several convergence results, including nonlinear transport and density currents, up to third order of accuracy in both space and time.

The study of moving particles (e.g. molecules, virus, vesicles, organelles, or whole cells) is crucial to decipher a plethora of cellular mechanisms within physiological and pathological conditions. Powerful live-imaging approaches enable life scientists to capture particle movements at different scale from cells to single molecules, that are collected in a series of frames. However, although these events can be captured, an accurate quantitative analysis of live-imaging experiments still remains a challenge. Two main approaches are currently used to study particle kinematics: kymographs, which are graphical representation of spatial motion over time, and single particle tracking (SPT) followed by linear linking. Both kymograph and SPT apply a space-time approximation in quantifying particle kinematics, considering the velocity constant either over several frames or between consecutive frames, respectively. Thus, both approaches intrinsically limit the analysis of complex motions with rapid changes in velocity. Therefore, we design, implement and validate a novel reconstruction algorithm aiming at supporting tracking particle trafficking analysis with mathematical foundations. Our method is based on polynomial reconstruction of 4D (3D+time) particle trajectories, enabling to assess particle instantaneous velocity and acceleration, at any time, over the entire trajectory. Here, the new algorithm is compared to state-of-the-art SPT followed by linear linking, demonstrating an increased accuracy in quantifying particle kinematics. Our approach is directly derived from the governing equations of motion, thus it arises from physical principles and, as such, it is a versatile and reliable numerical method for accurate particle kinematics analysis which can be applied to any live-imaging experiment where the space-time coordinates can be retrieved.

We propose a novel numerical method for the solution of the shallow water equations in different regimes of the Froude number making use of general polygonal meshes. The fluxes of the governing equations are split such that advection and acoustic-gravity sub-systems are derived, hence separating slow and fast phenomena. This splitting allows the nonlinear convective fluxes to be discretized explicitly in time, while retaining an implicit time marching for the acoustic-gravity terms. Consequently, the novel schemes are particularly well suited in the low Froude limit of the model, since no numerical viscosity is added in the implicit solver. Besides, stability follows from a milder CFL condition which is based only on the advection speed and not on the celerity. High order time accuracy is achieved using the family of semi-implicit IMEX Runge-Kutta schemes, while high order in space is granted relying on two discretizations: (i) a cell-centered finite volume (FV) scheme for the nonlinear convective contribution on the polygonal cells; (ii) a staggered discontinuous Galerkin (DG) scheme for the solution of the linear system associated to the implicit discretization of the pressure sub-system. Therefore, three different meshes are used, namely a polygonal Voronoi mesh, a triangular subgrid and a staggered quadrilateral subgrid. The novel schemes are proved to be Asymptotic Preserving (AP), hence a consistent discretization of the limit model is retrieved for vanishing Froude numbers, which is the given by the so-called "lake at rest" equations. Furthermore, the novel methods are well-balanced by construction, and this property is also demonstrated. Accuracy and robustness are then validated against a set of benchmark test cases with Froude numbers ranging in the interval $\Fr \approx [10^{-6};5]$, hence showing that multiple time scales can be handled by the novel methods.

This article is devoted to the construction of a new class of semi-Lagrangian (SL) schemes with implicit-explicit (IMEX) Runge-Kutta (RK) time stepping for PDEs involving multiple space-time scales. The semi-Lagrangian (SL) approach fully couples the space and time discretization, thus making the use of RK strategies particularly difficult to be combined with. First, a simple scalar advection-diffusion equation is considered as a prototype PDE for the development of a high order formulation of the semi-Lagrangian IMEX algorithms. The advection part of the PDE is discretized explicitly at the aid of a SL technique, while an implicit discretization is employed for the diffusion terms. Second, the SL-IMEX approach is extended to deal with hyperbolic systems with multiple scales, including balance laws, that involve shock waves and other discontinuities. A novel SL technique is proposed, which is based on the integration of the governing equations over the space-time control volume which arises from the motion of each grid point. High order of accuracy is ensured by the usage of IMEX RK schemes combined with a Cauchy-Kowalevskaya procedure that provides a predictor solution within each space-time element. The one-dimensional shallow water equations (SWE) are chosen to validate the new conservative SL-IMEX schemes, where convection and pressure fluxes are treated explicitly and implicitly, respectively. The asymptotic-preserving (AP) property of the novel schemes is also studied considering a relaxation PDE system for the SWE. A large suite of convergence studies for both the non-conservative and the conservative version of the novel class of methods demonstrates that the formal order of accuracy is achieved and numerical evidences about the conservation property are shown. The AP property for the corresponding relaxation system is also investigated.

Earthquake fault zones are more complex, both geometrically and rheologically, than an idealised infinitely thin plane embedded in linear elastic material. To incorporate nonlinear material behaviour, natural complexities, and multi-physics coupling within and outside of fault zones, here we present a first-order hyperbolic and thermodynamically compatible mathematical model for a continuum in a gravitational field which provides a unified description of nonlinear elasto-plasticity, material damage and of viscous Newtonian flows with phase transition between solid and liquid phases. The fault geometry and secondary cracks are described via a scalar function $\xi \in [0,1]$ that indicates the local level of material damage. The model also permits the representation of arbitrarily complex geometries via a diffuse interface approach based on the solid volume fraction function $\alpha \in [0,1]$. Neither of the two scalar fields $\xi$ and $\alpha$ needs to be mesh-aligned, allowing thus faults and cracks with complex topology and the use of adaptive Cartesian meshes (AMR). The model shares common features with phase-field approaches but substantially extends them. We show a wide range of numerical applications that are relevant for dynamic earthquake rupture in fault zones, including the co-seismic generation of secondary off-fault shear cracks, tensile rock fracture in the Brazilian disc test, as well as a natural convection problem in molten rock-like material., Comment: Supplementary Materials are at the end of the main manuscript

A new high order accurate staggered semi-implicit space-time discontinuous Galerkin (DG) method is presented for the simulation of viscous incompressible flows on unstructured triangular grids in two space dimensions. The staggered DG scheme defines the discrete pressure on the primal triangular mesh, while the discrete velocity is defined on a staggered edge-based dual quadrilateral mesh. In this paper, a new pair of equal-order-interpolation velocity-pressure finite elements is proposed. On the primary triangular mesh (the pressure elements) the basis functions are piecewise polynomials of degree $N$ and are allowed to jump on the boundaries of each triangle. On the dual mesh instead (the velocity elements), the basis functions consist in the union of piecewise polynomials of degree $N$ on the two subtriangles that compose each quadrilateral and are allowed to jump only on the dual element boundaries, while they are continuous inside. In other words, the basis functions on the dual mesh are built by continuous finite elements on the subtriangles. This choice allows the construction of an efficient, quadrature-free and memory saving algorithm. In our coupled space-time pressure correction formulation for the incompressible Navier-Stokes equations, arbitrary high order of accuracy in time is achieved through the use of time-dependent test and basis functions, in combination with simple and efficient Picard iterations. Several numerical tests on classical benchmarks confirm that the proposed method outperforms existing staggered semi-implicit space-time DG schemes, not only from a computer memory point of view, but also concerning the computational time., Comment: 42 pages, 22 figures. Submitted to "Communications on Applied Mathematics and Computation"

We are concerned with the numerical solution of a unified first order hyperbolic formulation of continuum mechanics that originates from the work of Godunov, Peshkov and Romenski (GPR model) and which is an extension of nonlinear hyperelasticity that is able to describe simultaneously nonlinear elasto-plastic solids at large strain, as well as viscous and ideal fluids. The proposed governing PDE system also contains the effect of heat conduction and can be shown to be symmetric and thermodynamically compatible. In this paper we extend the GPR model to the simulation of nonlinear dynamic rupture processes and material fatigue effects, by adding a new scalar variable to the governing PDE system. This extra parameter describes the material damage and is governed by an advection-reaction equation, where the stiff and highly nonlinear reaction mechanisms depend on the ratio of the local von Mises stress to the yield stress of the material. The stiff reaction mechanisms are integrated in time via an efficient exponential time integrator. Due to the multiple space-time scales, the model is solved on space-time adaptive Cartesian meshes using high order discontinuous Galerkin finite element schemes with a posteriori subcell finite volume limiting. A key feature of our new model is the use of a twofold diffuse interface approach that allows the cracks to form anywhere and at any time, independently of the chosen computational grid, without requiring that the geometry of the rupture fault be known a priori. We furthermore make use of a scalar volume fraction function that indicates whether a given point is inside the solid or outside, allowing the description of solids of arbitrarily complex shape.

We propose a new family of high order staggered semi-implicit discontinuous Galerkin (DG) methods for the simulation of natural convection problems. Assuming small temperature fluctuations, the Boussinesq approximation is valid and the flow can simply be modeled by the incompressible Navier-Stokes equations coupled with a transport equation for the temperature and a buoyancy source term in the momentum equation. Our numerical scheme is developed starting from the work presented in [TD14], in which the spatial domain is discretized using a face-based staggered unstructured mesh. For the computation of the advection and diffusion terms, two different algorithms are presented: i) a purely Eulerian explicit upwind-type scheme and ii) a semi-Lagrangian approach. The first methodology leads to a conservative scheme whose major drawback is the time step restriction imposed by the CFL stability condition. On the contrary, computational efficiency can be notably improved relying on a semi-Lagrangian approach. This method leads to an unconditionally stable scheme if the diffusive terms are discretized implicitly. Once the advection and diffusion contributions have been computed, the pressure Poisson equation is solved and the velocity is updated. As a second model for the computation of buoyancy-driven flows, we also consider the full compressible Navier-Stokes equations. The staggered semi-implicit DG method first proposed in [TD17] for all Mach number flows is properly extended to account for the gravity source terms arising in the momentum and energy conservation laws. The validity and the robustness of our novel class of staggered semi-implicit DG schemes is assessed at the aid of several classical benchmark problems, showing in all cases a good agreement with available numerical reference data. Finally, a detailed comparison between the incompressible and the compressible solver is presented.

ExaHyPE ("An Exascale Hyperbolic PDE Engine") is a software engine for solving systems of first-order hyperbolic partial differential equations (PDEs). Hyperbolic PDEs are typically derived from the conservation laws of physics and are useful in a wide range of application areas. Applications powered by ExaHyPE can be run on a student's laptop, but are also able to exploit thousands of processor cores on state-of-the-art supercomputers. The engine is able to dynamically increase the accuracy of the simulation using adaptive mesh refinement where required. Due to the robustness and shock capturing abilities of ExaHyPE's numerical methods, users of the engine can simulate linear and non-linear hyperbolic PDEs with very high accuracy. Users can tailor the engine to their particular PDE by specifying evolved quantities, fluxes, and source terms. A complete simulation code for a new hyperbolic PDE can often be realised within a few hours - a task that, traditionally, can take weeks, months, often years for researchers starting from scratch. In this paper, we showcase ExaHyPE's workflow and capabilities through real-world scenarios from our two main application areas: seismology and astrophysics.

In this paper we discuss a new and very efficient implementation of high order accurate ADER discontinuous Galerkin (ADER-DG) finite element schemes on modern massively parallel supercomputers. The numerical methods apply to a very broad class of nonlinear systems of hyperbolic partial differential equations. ADER-DG schemes are by construction communication avoiding and cache blocking and are furthermore very well-suited for vectorization, so that they appear to be a good candidate for the future generation of exascale supercomputers. We introduce the numerical algorithm and show some applications to a set of hyperbolic equations with increasing level of complexity, ranging from the compressible Euler equations over the equations of linear elasticity and the unified Godunov-Peshkov-Romenski (GPR) model of continuum mechanics to general relativistic magnetohydrodynamics (GRMHD) and the Einstein field equations of general relativity. We present strong scaling results of the new ADER-DG schemes up to 180,000 CPU cores. To our knowledge, these are the largest runs ever carried out with high order ADER-DG schemes for nonlinear hyperbolic PDE systems. We also provide a detailed performance comparison with traditional Runge-Kutta DG schemes., Comment: 26 pages, 10 figures, 2 tables

In most classical approaches of computational geophysics for seismic wave propagation problems, complex surface topography is either accounted for by boundary-fitted unstructured meshes, or, where possible, by mapping the complex computational domain from physical space to a topologically simple domain in a reference coordinate system. In this paper we propose a completely different strategy. We address the problem of geometrically complex free surface boundary conditions with a novel diffuse interface method on adaptive Cartesian meshes that consists in the introduction of a characteristic function $ 0\leq\alpha\leq 1$ which identifies the location of the solid medium and the surrounding air and thus implicitly defines the location of the free surface boundary. Our new approach completely avoids the problem of mesh generation, since all that is needed for the definition of the complex surface topography is to set a scalar color function to unity inside the regions covered by the solid and to zero outside. An analysis of the eigenvalues of the PDE system shows that the complexity of the geometry has no influence on the admissible time step size due to the CFL condition. The model reduces to the classical linear elasticity equations inside the solid medium where the gradients of $\alpha$ are zero, while in the diffuse interface zone at the free surface boundary the governing PDE system becomes nonlinear. We can prove that the solution of the Riemann problem with arbitrary data and a jump in $\alpha$ from unity to zero yields a Godunov-state at the interface that satisfies the free-surface boundary condition exactly. In order to reduce numerical dissipation, we use high order DG finite element schemes on adaptive AMR grids together with a high resolution shock capturing subcell finite volume (FV) limiter in the diffuse interface region.

In this paper we propose a new high order accurate space-time DG finite element scheme for the solution of the linear elastic wave equations in first order velocity-stress formulation in two and three-space dimensions on staggered unstructured triangular and tetrahedral meshes. The method reaches arbitrary high order of accuracy in both space and time via the use of space-time basis and test functions. Within the staggered mesh formulation, we define the discrete velocity field in the control volumes of a primary mesh, while the discrete stress tensor is defined on a face-based staggered dual mesh. The space-time DG formulation leads to an implicit scheme that requires the solution of a linear system for the unknown degrees of freedom at the new time level. The number of unknowns is reduced at the aid of the Schur complement, so that in the end only a linear system for the degrees of freedom of the velocity field needs to be solved. Thanks to the use of a spatially staggered mesh, the stencil of the final velocity system involves only the element and its direct neighbors. The resulting linear system can be efficiently solved via matrix-free iterative methods. The chosen discretization and the linear nature of the governing PDE system lead to an unconditionally stable scheme, which allows large time steps even for low quality meshes that contain sliver elements. The fully discrete staggered space-time DG method is proven to be energy stable for any order of accuracy, for any mesh and for any time step size. For the particular case of Crank-Nicolson time discretization and homogeneous material, the final velocity system can be proven to be symmetric and positive definite and in this case the scheme is also exactly energy preserving. The new scheme is applied to several test problems in two and three space dimensions, providing also a comparison with high order explicit ADER-DG schemes

In this paper we present a novel pressure-based semi-implicit finite volume solver for the equations of compressible ideal, viscous and resistive magnetohydrodynamics (MHD). The new method is conservative for mass, momentum and total energy and in multiple space dimensions it is constructed in such a way as to respect the divergence-free condition of the magnetic field exactly, also in the presence of resistive effects. This is possible via the use of multi-dimensional Riemann solvers on an appropriately staggered grid for the time evolution of the magnetic field and a double curl formulation of the resistive terms. The new semi-implicit method for the MHD equations proposed here discretizes all terms related to the pressure in the momentum equation and the total energy equation implicitly, making again use of a properly staggered grid for pressure and velocity. The time step of the scheme is restricted by a CFL condition based only on the fluid velocity and the Alfv\'en wave speed and is not based on the speed of the magnetosonic waves. Our new method is particularly well-suited for low Mach number flows and for the incompressible limit of the MHD equations, for which it is well-known that explicit density-based Godunov-type finite volume solvers become increasingly inefficient and inaccurate due to the increasingly stringent CFL condition and the wrong scaling of the numerical viscosity in the incompressible limit. We show a relevant MHD test problem in the low Mach number regime where the new semi-implicit algorithm is a factor of 50 faster than a traditional explicit finite volume method, which is a very significant gain in terms of computational efficiency. However, our numerical results confirm that our new method performs well also for classical MHD test cases with strong shocks. In this sense our new scheme is a true all Mach number flow solver., Comment: 26 pages, 12 figures,1 table