High-Order Locally A-Stable Implicit Schemes for Linear ODEs

Abstract

Accurate simulations of wave propagation in complex media like Earth subsurface can be performed with a reasonable computational burden by using hybrid meshes stuffing fine and coarse cells. Locally implicit time discretizations are then of great interest. They indeed allow using unconditionally stable Schemes in the regions of computational domain covered by small cells. The receivable values of the time step are then increased which reduces the computational costs while limiting the dispersion effects. In this work we construct a method that combines optimized explicit Schemes and implicit Schemes to form locally implicit schemes for linear ODEs, including in particular semi-discretized wave problems that are considered herein for numerical experiments. Both the explicit and implicit schemes used are one-step methods constructed using their stability function. The stability function of the explicit schemes are computed by maximizing the time step that can be chosen. The implicit schemes used are unconditionally stable and do not necessary require the same number of stages as the explicit schemes. The performance assessment we provide shows a very good level of accuracy for locally implicit schemes. It also shows that a locally implicit scheme is a good compromise between purely explicit and purely implicit schemes in terms of computational time and memory usage.

Remarks on the Implementation of the HDG Formulation of the Acoustic Wave Equation.

For the first order acoustic wave Eq. (24), the local variational formulation set on an element K in the HDG formulation is given as

$$\begin{aligned} \left\{ \begin{aligned}&\frac{d}{dt} \int _K \rho \, u \, \varphi _i \, dx + \int _K v \cdot \nabla \varphi _i \, dx - \int _{\partial K} v \cdot n \, \varphi _i \, dx + \int _{\partial K} \tau (u-\lambda ) \, \varphi _i \, dx = 0 \\&\frac{d}{dt} \int _K \mu ^{-1} v \cdot \psi _i \, dx - \int _K \nabla u \cdot \psi _i + \int _{\partial K} (u-\lambda ) \psi _i \cdot n \, dx = 0 \\&\int _{\partial K} (-v \cdot n + \tau (u- \lambda ) ) \, q \, dx = 0 \end{aligned} \right. \end{aligned}$$

(25)

where u, v are volume unknowns (discontinuous) and $\lambda$ a surface unknown. $\varphi _i, \psi _i, q$ are the test functions associated with u, v and $\lambda$. The penalization parameter $\tau$ is equal to $\sqrt{\rho \, \mu }$. The last equation is modified when the boundary of K meets the boundary of the computational domain. For instance, is a first-order absorbing boundary condition is set

$$\begin{aligned} v \cdot n + \sqrt{\rho \, \mu } \, u = f_A . \end{aligned}$$

the last equation of (25) becomes:

$$\begin{aligned} \int _{\varGamma _N} (-v \cdot n + \tau (u - 2 \lambda )) \, q \, dx = -\int _{\varGamma _N} f_A \, q \, dx. \end{aligned}$$

The HDG formulation provides locally for each element K the semi-discrete system

$$\begin{aligned} \left\{ \begin{aligned}&M_u \dfrac{dU}{dt} + K_u V + C_u \varLambda = 0 \\&M_v \dfrac{dV}{dt} + K_v U + C_v \varLambda = 0 \\&C_\lambda \varLambda + C_u^T U + C_v^T V = 0 \\ \end{aligned} \right. \end{aligned}$$

(26)

$C_\lambda$ is multiplied by two for an edge with an absorbing boundary condition. When the ODE (20) (which is set in the fine region) is solved with an implicit scheme, a linear system of the form

$$\begin{aligned} \beta Y + A P Y = F \end{aligned}$$

has to be solved for close dofs (with $Y=(U, V)$). For elements located in the fine region, we have

$$\begin{aligned} \left\{ \begin{aligned}&\beta M_u U + K_u V + C_u \varLambda = F_u \\&\beta M_v V + K_v U + C_v \varLambda = F_v \\&C_\lambda \varLambda + C_u^T U + C_v^T V = 0 \\ \end{aligned} \right. \end{aligned}$$

(27)

Unknowns U and V are eliminated element-wise to obtain a local system in $\varLambda$

$$\begin{aligned} \left[ C_\lambda - \left( C_u^T C_v^T \right) \left( \begin{array}{cc} \beta M_u &{} K_u \\ K_v &{} \beta M_v \end{array} \right) ^{-1} \left( \begin{array}{c} C_u \\ C_v \end{array} \right) \right] \varLambda = - \left( C_u^T C_v^T \right) \left( \begin{array}{cc} \beta M_u &{} K_u \\ K_v &{} \beta M_v \end{array} \right) ^{-1} \left( \begin{array}{c} F_u \\ F_v \end{array} \right) . \end{aligned}$$

This equation has to be assembled with all other elements to obtain the final system solved by $\varLambda$. For adjacent elements (located on the coarse region), we have:

$$\begin{aligned} \left\{ \begin{aligned}&\beta M_u U + C_u \varLambda = F_u \\&\beta M_v V + C_v \varLambda = F_v \\&C_\lambda \varLambda = 0 \\ \end{aligned} \right. \end{aligned}$$

(28)

where only unknowns on edges of the fine region are concerned for $\varLambda$. The equation to be assembled with other elements is given as

$$\begin{aligned} C_\lambda \varLambda = 0. \end{aligned}$$

As a result, when the linear system is assembled for $\varLambda$, only unknowns located on edges (faces in 3-D) of the fine region are involved. It is actually equivalent to add the contribution $C_\lambda \varLambda = 0$ or impose a fictitious homogeneous absorbing boundary condition on edges located at the interface between the fine and coarse region. Once $\varLambda$ has been computed on the edges of the fine region ($\varLambda$ is null for other edges), the unknown U and V are reconstructed element-wise, e.g. in adjacent elements:

$$\begin{aligned} U= & {} \dfrac{1}{\beta } M_u^{-1} \left( F_u - C_u \varLambda \right) \\ V= & {} \dfrac{1}{\beta } M_v^{-1} \left( F_v - C_v \varLambda \right) \end{aligned}$$

For quadrilateral/hexahedral elements, the mass matrices $M_u$ and $M_v$ are diagonal, and the elimination/reconstruction of U and V can be lead efficiently as detailed in [21].