Convergence analysis of energy conserving explicit local time-stepping methods for the wave equation

TL;DR

This paper provides a rigorous convergence analysis of an explicit local time-stepping method based on leap-frog for the wave equation, enabling efficient simulations on locally refined meshes with complex geometries.

Contribution

It offers the first rigorous proof of convergence for the fully-discrete LTS-LF method combined with FEM, demonstrating its effectiveness in handling corner singularities.

Findings

Proves convergence of the LTS-LF Galerkin FEM method.

Shows effectiveness in complex geometries with corner singularities.

Enables efficient wave simulations on locally refined meshes.

Abstract

Local adaptivity and mesh refinement are key to the efficient simulation of wave phenomena in heterogeneous media or complex geometry. Locally refined meshes, however, dictate a small time-step everywhere with a crippling effect on any explicit time-marching method. In [18] a leap-frog (LF) based explicit local time-stepping (LTS) method was proposed, which overcomes the severe bottleneck due to a few small elements by taking small time-steps in the locally refined region and larger steps elsewhere. Here a rigorous convergence proof is presented for the fully-discrete LTS-LF method when combined with a standard conforming finite element method (FEM) in space. Numerical results further illustrate the usefulness of the LTS-LF Galerkin FEM in the presence of corner singularities.

Equations

$$a\left(u,v\right)=a\left(v,u\right)\forall u,v\in V$$

$$\left|a\left(u,v\right)\right|\le C_{\mathrm{cont}}{\Vert u\Vert }_{H^1\left(\Omega \right)}{\Vert v\Vert }_{H^1\left(\Omega \right)}\forall u,v\in V$$

$$a\left(u,u\right)\ge c_{\mathrm{coer}}{\Vert u\Vert }_{H^1\left(\Omega \right)}^2\forall u\in V.$$

$$\left(\ddot{u},w\right)+a\left(u,w\right)=F\left(w\right)\forall w\in V,t>0$$

$$u\left(0\right)=u_0\text{and}\dot{u}\left(0\right)=v_0.$$

$$\begin{array}{rl}{u}_{tt}-\nabla \cdot (c^2\nabla u)& =f\text{mbox}in\Omega \times (0,T),\\ u& =0\text{mbox}on{\Gamma }_D\times (0,T),\\ \frac{\partial u}{\partial \nu }& =0\text{mbox}on{\Gamma }_N\times (0,T),\\ u{|}_{t=0}& =u_0\text{mbox}in\Omega ,\\ u_t{|}_{t=0}& =u_0\text{mbox}in\Omega .\\ & \end{array}$$

$$h_{\tau }:=\mathrm{diam}\tau \text{and}h:=\underset{\tau \in \mathcal{T}}{\max }{h}_{\tau }\text{and}{h}_{\min }:=\underset{\tau \in \mathcal{T}}{\min }{h}_{\tau }$$

$$\gamma\left(\mathcal{T}\right):=\max_{\tau}\left\{\begin{array}[c]{ll}\max\left\{\frac{h_{\tau}}{h_{t}}:t\in\mathcal{T}:t\cap\tau\neq\emptyset\right\}&d=1,\\ \dfrac{h_{\tau}}{\rho_{\tau}}&d=2,3,\end{array}\right.$$

$$C_{\mathrm{qu}}:=\frac{h}{h_{\min }}.$$

$$S_{\mathcal{T}}^m:=\left\{u\in C^0\left(\Omega \right)\mid \forall \tau \in \mathcal{T}:{u|}_{\tau }\in {\mathbb{P}}_m\right\},$$

$$\hat{\tau }:=\left\{\mathbf{x}={\left(x_i\right)}_{i=1}^d\in {\mathbb{R}}_{\ge 0}^d:\sum _{i=1}^d{x}_i\le 1\right\}$$

$${\Sigma }_{\tau }^m:=\left\{{\varphi }_{\tau }\left(z\right):z\in {\hat{\Sigma }}^m\right\}.$$

$${\Sigma }_{\mathcal{T}}^m:={\bigcup }_{\tau \in \mathcal{T}}{\Sigma }_{\tau }^m.$$

$$b_{z,m}\in S_{\mathcal{T}}^{m}\quad\text{and\quad}\forall z^{\prime}\in\Sigma_{\mathcal{T}}^{m}\text{ it holds }b_{z,m}\left(z^{\prime}\right)=\left\{\begin{array}[c]{ll}1&z=z^{\prime},\\ 0&\text{otherwise.}\end{array}\right.$$

$$P_{\Sigma }\mathbf{u}=\sum _{z\in \Sigma }{u}_z{b}_{z,m}\text{and}\left({\mathbf{R}}_{\Sigma }v\right)={\left({\int }_{\Omega }v{b}_{z,m}\right)}_{z\in \Sigma }.$$

$${\mathbf{M}}_{\Sigma }:={\left({\int }_{\Omega }{b}_{z,m}{b}_{z^{\prime },m}\right)}_{z,z^{\prime }\in \Sigma }\text{.}$$

$$\left(u,v\right)_{\mathcal{T}}:=\sum_{\tau\in\mathcal{T}}\left|\tau\right|\sum_{\mathbf{z}\in\Sigma_{\tau}^{m}}u_{z}v_{z}=\left\langle\mathbf{D}_{\Sigma_{\mathcal{T}}^{m}}\mathbf{u},\mathbf{v}\right\rangle\text{\quad with\quad}\left\{\begin{array}[c]{l}\mathbf{D}_{\Sigma_{\mathcal{T}}^{m}}=\operatorname*{diag}\left[d_{z}:z\in\Sigma_{\mathcal{T}}^{m}\right],\\ d_{z}:=\left|\operatorname*{supp}b_{z,m}\right|,\end{array}\right.$$

$${\Vert u\Vert }_{\mathcal{T}}:={\left(u,u\right)}_{\mathcal{T}}^{1/2}.$$

$$d_{\mathcal{N},z}:=\left\{\begin{array}[c]{ll}d_{z}&z\in\mathcal{N},\\ 0&z\in\Sigma_{\mathcal{T}}^{m}\backslash\mathcal{N}.\end{array}\right.$$

$$R_{\mathcal{N}}={\mathbf{R}}^{-1}{\mathbf{D}}_{\mathcal{N}}{P}^{-1}.$$

$$\left(R_{\mathcal{N}}P\mathbf{u},P\mathbf{v}\right)=⟨{\mathbf{D}}_{\mathcal{N}}\mathbf{u},\mathbf{v}⟩=\sum _{z\in \mathcal{N}}{d}_z{u}_z{v}_z.$$

$$\mathrm{supp}\left(R_{\mathcal{N}}u\right)\subset {\Omega }_{\mathcal{N}}:={\bigcup }_{\begin{array}{c}\tau \in \mathcal{T}\\ \tau \cap \mathcal{N}\ne \varnothing \end{array}}\tau .$$

$$V_{\mathcal{T}}^m:=S_{\mathcal{T}}^m\cap V\text{.}$$

$${⟨Au,v⟩}_{V^{\prime }\times V}$$

$$\left(A_{S}u,v\right)$$

$$\left({\ddot{u}}_S,v\right)+a\left(u_S,v\right)=F\left(v\right)\forall v\in S,t>0$$

$$\left.\begin{array}[c]{c}\left(u_{S}\left(0\right),w\right)=\left(u_{0},w\right)\\ \\ \left(\dot{u}_{S}\left(0\right),w\right)=\left(v_{0},w\right)\end{array}\right\}\quad\forall w\in S.$$

$$\left(f_S,w\right)=F\left(w\right)\forall w\in S\text{and}{f}_S^{\left(n\right)}:=f_S\left(t_n\right)\text{,}$$

$$f_{S,m}^{\left(n\right)}:=f_S\left(t_n+\frac{m}{p}\Delta t\right),\text{mbox}with{f}_{S,m}^{\left(n\right)}=P{\mathbf{f}}_{s,m}^{\left(n\right)},0\le m\le p.$$

$$\mathbf{w}={\mathbf{M}}^{-1}\left(\left(\mathbf{M}-{\mathbf{D}}_{\mathcal{N}}\right){\mathbf{f}}_S^{\left(n\right)}-\mathbf{A}\left(\mathbf{I}-{\mathbf{M}}^{-1}{\mathbf{D}}_{\mathcal{N}}\right){\mathbf{u}}_S^{\left(n\right)}\right).$$

$${\tilde{\mathbf{u}}}_{S,1}^{\left(n\right)}={\tilde{\mathbf{u}}}_{S,0}^{\left(n\right)}+\frac{1}{2}{\left(\frac{\Delta t}{p}\right)}^2\left({\mathbf{w}+\mathbf{M}}^{-1}\left({\mathbf{D}}_{\mathcal{N}}{\mathbf{f}}_S^{\left(n\right)}-{\mathbf{AM}}^{-1}{\mathbf{D}}_{\mathcal{N}}{\tilde{\mathbf{u}}}_{S,0}^{\left(n\right)}\right)\right).$$

Full text

Keywords: wave propagation, finite element methods, explicit time integration, leap-frog method, error analysis, convergence theory

AMS-Classification: 65M12, 65M20, 65M60, 65L06, 65L20

1 Introduction

Efficient numerical methods are crucial for the simulation of time-dependent acoustic, electromagnetic or elastic wave phenomena. Finite element methods (FEM), in particular, easily accommodate varying mesh sizes or polynomial degrees. Hence, they are remarkably effective and widely used for the spatial discretization in heterogeneous media or complex geometry. However, as spatial discretizations become increasingly accurate and flexible, the need for more sophisticated time-integration methods for the resulting systems of ordinary differential equations (ODE) becomes all the more apparent.

Today’s standard use of local adaptivity and mesh refinement causes a severe bottleneck for any standard explicit time integration. Even if the refined region consists of only a few small elements, those smallest elements will impose a tiny time-step everywhere for stability reasons. To overcome that geometry induced stiffness, various local time integration strategies were devised in recent years. Typically the mesh is partitioned into a “coarse” part, where most of the elements are located, and a “fine” part, which contains the remaining few smallest elements. Inside the “coarse” part, standard explicit methods are used for time integration. Inside the “fine” part, local time-stepping (LTS) methods either use implicit or explicit time integration.

Locally implicit methods are based on implicit-explicit (IMEX) approaches commonly used in CFD for operator splitting [2, 31]. They require the solution of a linear system inside the refined region at every time-step, which becomes increasingly expensive (and ill-conditioned) as the mesh size decreases [33]. Alternatively, exponential Adams methods [29] apply the matrix exponential locally in the fine part while reducing to the underlying Adams-Bashforth scheme elsewhere.

Locally implicit or exponential time integrators typically use the same time-step everywhere but apply different methods in the ”fine” and the ”coarse” part. In contrast, explicit LTS methods typically use the same method everywhere but take smaller time-steps inside the “fine” region [24]; hence, they remain fully explicit. Since the finite-difference based adaptive mesh refinement (AMR) method by Berger and Oliger [5], various explicit LTS were proposed in the context of discontinuous Galerkin (DG) FEM, which permit a different time-step inside each individual element [23, 35, 21, 46, 14, 15]. In [16] multiple time-stepping algorithms were presented which allow any choice of explicit Adams type or predictor-corrector scheme for the integration of the coarse region and any choice of ODE solver for the integration of the fine part. High-order explicit LTS methods for wave propagation were derived in [26, 27, 25] starting either from Leap-Frog, Adams-Bashforth or Runge-Kutta methods.

In [11, 4, 13], Collino et al. proposed a first energy conserving LTS method for the wave equation which was analyzed in [12, 32]. This second-order method conserves a discrete energy and thereby guarantees stability, but it requires at every time-step the solution of a linear system at the interface between the fine and the coarser elements; hence, it is not fully explicit. A fully explicit second-order LTS method was proposed for Maxwell’s equations by Piperno  [41] and further developed in [20, 37]. In [36, 42], the high-order energy conserving explicit LTS method proposed in [18] was successfully applied to 3D seismic wave propagation on a large-scale parallel computer architecture.

Despite the many different explicit LTS methods that were proposed and successfully used for wave propagation in recent years, a rigorous fully discrete space-time convergence theory is still lacking. In fact, convergence has been proved only for the method of Collino et al. [12, 11, 32] and very recently for the locally implicit method for Maxwell’s equations by Verwer [47, 17, 30], neither fully explicit. Indeed, the difficulty in proving convergence of fully explicit LTS methods is twofold. On the one hand, classical proofs of convergence [22, 3] always assume standard time discretizations, while proofs for multirate schemes (in the ODE literature) are always restricted to the finite-dimensional case. Hence, standard convergence analysis cannot be easily extended to LTS methods for partial differential equations. On the other hand, when explicit LTS schemes are reformulated as perturbed one-step schemes, they involve products of differential and restriction operators, which do not commute and seem to inevitably lead to a loss of regularity.

Our paper is structured as follows. In Section 2, we consider a general second-order wave equation and introduce (the notation for) conforming finite element spaces on simplicial meshes with local polynomial order $m$. Next, we define finite-dimensional restriction operators to the ”fine” grid and formulate the leap-frog (LF) based LTS method from [18] in a Galerkin conforming finite element setting. In Section 3, we prove continuity and coercivity estimates for the LTS operator that are robust with respect to the number of local time-steps $p$, provided a genuine CFL condition is satisfied. Here, new estimates on the coefficients that appear when rewriting the LTS-LF scheme in ”leap-frog manner” play a key-role – see Appendix. Those estimates pave the way for the stability estimate of the time iteration operator, for which we then prove a stability bound independently of $p$. In doing so, the truncation errors are estimated through standard Taylor arguments for the leap-frog method. Due to the local restriction, however, a judicious splitting of the iteration operator and its inverse is required to avoid negative powers of $h$ via inverse inequalities. By combining our analysis of the semi-discrete formulation, which takes into account the effect of local time-stepping, with classical error estimates [3], we eventually obtain optimal convergence rates explicit with respect to the time step $\Delta t$, the mesh size $h$, the right-hand side, the initial data and the final time $T$, which hold uniformly with respect to the number of local time-steps $p$. Finally, in Section 4, we report on some numerical experiments inside an L-shaped domain. By applying the LTS method in the locally refined region near the re-entrant corner, we obtain a significant speedup over a standard leap-frog method with a small time-step everywhere.

2 Galerkin Discretization with Leap-Frog Based Local Time-Stepping

2.1 The Wave Equation

Let $\Omega\subset\mathbb{R}^{d}$ be a Lipschitz domain and $L^{2}\left(\Omega\right)$ denote the space of square integrable, real-valued functions with scalar product denoted by $\left(\cdot,\cdot\right)$ and corresponding norm by $\left\|\cdot\right\|=\left(\cdot,\cdot\right)^{1/2}$. Next, let $H^{1}\left(\Omega\right)$ denote the standard Sobolev space of all square integrable, real-valued functions whose first (weak) derivatives are also square integrable; as usual, $H^{1}\left(\Omega\right)$ is equipped with the norm $\left\|u\right\|_{H^{1}\left(\Omega\right)}=(\left\|u\right\|^{2}+\left\|\nabla u\right\|^{2})^{1/2}$.

We now let $V\subset H^{1}\left(\Omega\right)$ denote a closed subspace of $H^{1}\left(\Omega\right)$, such as $V=H^{1}\left(\Omega\right)$ or $V=H_{0}^{1}\left(\Omega\right)$, and consider a bilinear form $a:V\times V\rightarrow\mathbb{R}$ which is symmetric, continuous, and coercive:

[TABLE]

[TABLE]

and

[TABLE]

and

[TABLE]

For given $u_{0}\in V,v_{0}\in L^{2}\left(\Omega\right)$ and $F:\left[0,T\right]\rightarrow V^{\prime}$, we consider the wave equation: Find $u:\left[0,T\right]\rightarrow V$ such that

[TABLE]

with initial conditions

[TABLE]

It is well known that (2)–(3) is well-posed for sufficiently regular $u_{0}$, $v_{0}$ and $F$ [34]. In fact, the weak solution $u$ can be shown to be continuous in time, that is, $u\in C^{0}(0,T;V),\dot{u}\in C^{0}(0,T;L^{2}\left(\Omega\right))$ – see [[34], Chapter III, Theorems 8.1 and 8.2] for details – which implies that the initial conditions (3) are well defined.

Example 1

The classical second-order wave equation in strong form is given by

[TABLE]

In this case, we have $V:=H_{D}^{1}\left(\Omega\right):=\left\{w\in H^{1}\left(\Omega\right):\left.w\right|_{\Gamma_{D}}=0\right\}$; the bilinear form is given by $a\left(u,v\right):=\left(c^{2}\nabla u,\nabla u\right)$ and the right-hand side by $F\left(w\right)=\left(f,w\right)$ for all $w\in V$.

2.2 Galerkin Finite Element Discretization

For the semi-discretization in space, we employ the Galerkin finite element method and we first have to introduce some notation. We assume for the spatial dimension $d\in\left\{1,2,3\right\}$ and that the bounded Lipschitz domain $\Omega\subset\mathbb{R}^{d}$ is an interval for $d=1$, a polygonal domain for $d=2$, and a polyhedral domain for $d=3$. Let $\mathcal{T}:=\left\{\tau_{i}:1\leq i\leq N_{\mathcal{T}}\right\}$ denote a conforming (i.e.: no hanging nodes), simplicial finite element mesh for $\Omega$. Let

[TABLE]

and denote by $\rho_{\tau}$ the diameter of the largest inscribed ball in $\tau$. As a convention, the simplices $\tau\in\mathcal{T}$ are closed sets. The shape regularity constant $\gamma$ of the mesh $\mathcal{T}$ is defined by

[TABLE]

and the quasi-uniformity constant by

[TABLE]

For $m\in\mathbb{N}$, we define the continuous, piecewise polynomial finite element space by

[TABLE]

where $\mathbb{P}_{m}$ is the space to $d$-variate polynomials of maximal total degree $m$. The definition of a Lagrangian nodal basis is standard and employs the concept of a reference element. Let

[TABLE]

denote the reference element. For $\tau\in\mathcal{T}$, let $\phi_{\tau}:\widehat{\tau}\rightarrow\tau$ denote an affine pullback. For $m\geq 1$, we denote by $\hat{\Sigma}^{m}$ a set of nodal points in $\hat{\tau}$ unisolvent on $\mathbb{P}_{m}$, which allow to impose continuity across simplex faces. The nodal points on a simplex $\tau\in\mathcal{T}$ are then given by lifting those of the reference element:

[TABLE]

The set of global nodal points is given by

[TABLE]

A Lagrange basis for $S_{\mathcal{T}}^{m}$ is given by $\left(b_{z,m}\right)_{z\in\Sigma_{\mathcal{T}}^{m}}$ via the conditions

[TABLE]

For a subset $\Sigma\subset\Sigma_{\mathcal{T}}^{m}$, we define a prolongation map $P_{\Sigma}:\mathbb{R}^{\Sigma}\rightarrow S_{\mathcal{T}}^{m}$ and a restriction map $\mathbf{R}_{\Sigma}:S_{\mathcal{T}}^{m}\rightarrow\mathbb{R}^{\Sigma}$ by

[TABLE]

The mass matrix, $\mathbf{M}_{\Sigma}$, is given by

[TABLE]

If $\Sigma=\Sigma_{\mathcal{T}}^{m}$ holds, we write $P,\mathbf{R}$, $\mathbf{M}$ short for $P_{\Sigma},\mathbf{R}_{\Sigma}$, $\mathbf{M}_{\Sigma}$.

Remark 2

Since $\mathbf{M}_{\Sigma}=\mathbf{R}_{\Sigma}P_{\Sigma}$, we also have $P_{\Sigma}^{-1}=\mathbf{M}_{\Sigma}^{-1}\mathbf{R}_{\Sigma}$.

The matrix $\mathbf{M}_{\Sigma}$ is the matrix representation of the $L^{2}$-scalar product with respect to the basis $\left(b_{z,m}\right)_{z\in\Sigma}$. We introduce a diagonally weighted, mesh dependent Euclidean scalar product which is equivalent to the bilinear form $\left\langle\mathbf{u},\mathbf{M}_{\Sigma}\mathbf{v}\right\rangle$ (cf. Lemma 8), where $\left\langle\cdot,\cdot\right\rangle$ denotes the Euclidean scalar product on $\mathbb{R}^{\Sigma}$.

For $u=P\mathbf{u}$ and $v=P\mathbf{v}$ with $\mathbf{u}=\left(u_{z}\right)_{z\in\Sigma_{\mathcal{T}}^{m}}$ and $\mathbf{v}=\left(v_{z}\right)_{z\in\Sigma_{\mathcal{T}}^{m}}$ we set

[TABLE]

where, for a measurable set $\omega\subset\mathbb{R}^{d}$, we denote by $\left|\omega\right|$ its $d$-dimensional volume. The norm is given by

[TABLE]

For later use, we define a localized version of $\mathbf{D}_{\Sigma_{\mathcal{T}}^{m}}$. Let $\mathcal{N}\subset\Sigma_{\mathcal{T}}^{m}$ and define the diagonal matrix $\mathbf{D}_{\mathcal{N}}=\operatorname{diag}\left[d_{\mathcal{N},z}:z\in\Sigma_{\mathcal{T}}^{m}\right]$ by

[TABLE]

We define the fine grid restriction operator $R_{\mathcal{N}}:S_{\mathcal{T}}^{m}\rightarrow S_{\mathcal{T}}^{m}$ by

[TABLE]

Remark 3

Note that the diagonal matrix $\mathbf{D}_{\mathcal{N}}$ corresponds to the matrix representation of $R_{\mathcal{N}}$:

[TABLE]

For the support of $R_{\mathcal{N}}u$ it holds

[TABLE]

The operator $R_{\mathcal{N}}$ is symmetric positive semi-definite, which follows from $d_{z}\geq 0$ and the symmetry of the right-hand side in (6).

We define conforming subspaces of $V$ by

[TABLE]

Notation 4

We write $S$ short for $V_{\mathcal{T}}^{m}$ if no confusion is possible. Since $S=S_{\mathcal{T}}^{m}\cap V,$ we may assume that there is a subset $\Sigma_{S}\subset\Sigma_{\mathcal{T}}^{m}$ such that $S=\operatorname{span}\left\{b_{z,m}:z\in\Sigma_{S}\right\}$.

The operators associated to the continuous and discrete bilinear form are the linear mappings $A:V\rightarrow V^{\prime}$ and $A_{S}:S\rightarrow S$ defined by

[TABLE]

Here $\left\langle\cdot,\cdot\right\rangle_{V^{\prime}\times V}$ is the continuous extension of the $L^{2}\left(\Omega\right)$ scalar product to the dual pairing $\left\langle\cdot,\cdot\right\rangle_{V^{\prime}\times V}$.

Example 5

If homogeneous Dirichlet boundary conditions are imposed for the wave equation we have $V:=H_{0}^{1}\left(\Omega\right):=\left\{u\in H^{1}\left(\Omega\right)\mid\left.u\right|_{\partial\Omega}=0\right\}$. The nodal points $\Sigma_{\mathcal{T}}^{1}$ for the $\mathbb{P}_{1}$ finite element space are the inner triangle vertices and $b_{z,1}$ is the usual continuous, piecewise affine basis function for the nodal point $z$.

The semi-discrete wave equation then is given by: find $u_{S}:\left[0,T\right]\rightarrow S$ such that

[TABLE]

[TABLE]

with initial conditions

[TABLE]

2.3 Discrete LTS-Galerkin FE Formulation

Starting from the leap-frog based local time-stepping LTS-LF scheme from [18], we now present the fully discrete space-time Galerkin FE formulation. First we let the (global) time-step $\Delta t=T/N$ and denote by $u_{S}^{\left(n\right)}=P\mathbf{u}_{S}^{\left(n\right)}$ the FE approximation at time $t_{n}=n\Delta t$ for the corresponding coefficient vector (nodal values) $\mathbf{u}_{S}^{\left(n\right)}\in\mathbb{R}^{\Sigma}$ . Similarly we define the right-hand sides $f_{S}:\left[0,T\right]\rightarrow S$ and $f_{S}^{\left(n\right)}\in S$ by

[TABLE]

where again $f_{S}^{\left(n\right)}=P\,\mathbf{f}_{S}^{\left(n\right)}$ with corresponding coefficients $\mathbf{f}_{S}^{\left(n\right)}\in\mathbb{R}^{\Sigma}$.

Given the numerical solution at times $t_{n-1}$ and $t_{n}$, the LTS-LF method then computes the numerical solution of (7) at $t_{n+1}$ by using a smaller time-step $\Delta\tau=\Delta t/p$ inside the regions of local refinement; here, $p\geq 2$ denotes the ”coarse” to ”fine” mesh size ratio. Clearly, if the maximal velocity in the coarse and the fine regions differ significantly, the choice of $p$ should also reflect that variation and instead denote the local CFL number ratio. In the ”fine” region, the right-hand side is also evaluated at the intermediate times $t_{n+\frac{m}{p}}=t_{n}+m\Delta\tau$ and we let

[TABLE]

In Algorithm 6, we list the full second-order LTS-LF Algorithm ([18], [26, Alg. 1]) for the sake of completeness. All computations in Steps 2 and 3 that involve the right-hand side $\mathbf{f}_{S,m}^{\left(n\right)}$ or the stiffness matrix $\mathbf{A}$ only affect those degrees of freedom inside the region of local refinement or directly adjacent to it. The successive updates of the coarse unknowns involving $\mathbf{w}$ during sub-steps reduce to a single standard LF step of size $\Delta t$ and, in fact, can be replaced by it. In that sense, Algorithm 6 yields a local time-stepping method. We remark that higher order LTS-LF methods of arbitrarily high (even) accuracy were derived and implemented in [18].

Algorithm 6

LTS-LF Galerkin FE Algorithm

Set $\mathbf{\tilde{u}}_{S,0}^{\left(n\right)}:=\mathbf{u}_{S}^{\left(n\right)}$ and compute $\mathbf{w}$ as

[TABLE] 2. 2.

Compute

[TABLE] 3. 3.

For $m=1,\ldots,p-1$, compute

[TABLE] 4. 4.

Compute

[TABLE]

Like the standard leap-frog method (without local time-stepping), the LTS-LF Algorithm requires in principle the solution of a linear system involving ${\mathbf{M}}$ at every time-step. Although the mass matrix is sparse, positive definite, and well-conditioned so that solving linear systems with this matrix is relatively cheap, this computational effort is commonly avoided by using either mass-lumping techniques [10, 38], spectral elements [7, 9] or discontinuous Galerkin finite elements [1, 28]. The resulting LTS-LF scheme is then fully explicit.

In [18], the above LTS-LF Algorithm was rewritten in “leap-frog manner” by introducing the perturbed bilinear form $a_{p}:S\times S\rightarrow\mathbb{R}$:

[TABLE]

with associated operator

[TABLE]

Here the constants $\alpha_{j}^{m}$, $j=1,\dots,m-1$ are recursively defined for $m\geq 2$ by

[TABLE]

Then the LTS-LF scheme (Algorithm 6) is equivalent to

[TABLE]

Neither the equivalent formulation (12) nor the constants $\alpha_{j}^{m}$ are ever used in practice but only for the purpose of analysis; in fact, the constants $\alpha_{j}^{m}$ do not appear in Algorithm 1.

Remark 7

*In (12) the term $a\left(u_{0},w\right)$ in the third equation could be replaced by $a_{p}\left(u_{0},w\right)$ which allows for local time-stepping already during the very first time-step. In that case, the analysis below also applies but requires a minor change, namely, replacing $A_{S}$ by $A_{S,p}$ in (74) and (75). This modification neither affects the stability nor the convergence rate of the overall LTS-LF scheme. *

3 Stability and Convergence Analysis

3.1 Estimates of the Bilinearform

The following equivalence of the continuous $L^{2}\left(\Omega\right)$- and mesh-dependent norm is well known.

Lemma 8

$\left\|\cdot\right\|_{\mathcal{T}}$* and $\left\|\cdot\right\|$ are equivalent norms on $S_{\mathcal{T}}^{m}$. The constants $c_{\operatorname*{eq}}$, $C_{\operatorname*{eq}}$ in the equivalence estimates*

[TABLE]

only depend on the polynomial degree $m$ and the shape regularity constant $\gamma\left(\mathcal{T}\right)$.

It is also well known that the functions in $S_{\mathcal{T}}^{m}$ satisfy an inverse inequality (for a proof we refer, e.g., [8, (3.2.33) with $m=1$, $q=r=2$, $l=0$, $n=d$.]111There is a misprint in this reference: $m-1$ should be replaced by $m-\ell$, see also [6, (4.5.3) Lemma].).

Lemma 9

There exists a constant $C_{\operatorname*{inv}}>0$, which only depends on $\gamma\left(\mathcal{T}\right)$ and $m$, such that for all $\tau\in\mathcal{T}$

[TABLE]

The global versions of the inverse inequality involves also the quasi-uniformity constant

[TABLE]

for all $u\in S_{\mathcal{T}}^{m}$.

In the next step, we will estimate $\left\|A_{S}u\right\|$ in terms of $\left\|u\right\|_{H^{1}\left(\Omega\right)}$.

Lemma 10

It holds

[TABLE]

Proof. Since $A_{S}$ is a self-adjoint, positive operator there exists an orthonormal system $\left(\eta_{\nu}\right)_{\nu=1}^{M}$ such that

[TABLE]

and

[TABLE]

where $M:=\dim S$. Hence, every function $v\in S$ has a representation

[TABLE]

For $s\in\mathbb{R}$ we define the norm on $S$

[TABLE]

It is obvious that for all $v\in S$, it holds

[TABLE]

Note that

[TABLE]

We assume that the eigenvalues $\lambda_{\nu}$ are ordered increasingly. From Lemma 9 we conclude that

[TABLE]

holds. Hence,

[TABLE]

Next, we will estimate the bilinear form $a_{p}\left(\cdot,\cdot\right)$.

Lemma 11

The operator $R_{\mathcal{N}}$ as in (5) has bounded $L^{2}\left(\Omega\right)$ norm:

[TABLE]

For $u\in\mathcal{S}_{\mathcal{T}}^{m}$ it holds

[TABLE]

Proof. Let $u=P\mathbf{u}$ and $v=P\mathbf{v}$ with $\mathbf{u}=\left(u_{z}\right)_{z\in\Sigma_{\mathcal{T}}^{m}}$, $\mathbf{v}=\left(v_{z}\right)_{z\in\Sigma_{\mathcal{T}}^{m}}$. We employ

[TABLE]

Hence

[TABLE]

For the second estimate we employ (15) and (14) to obtain

[TABLE]

for all $u\in\mathcal{S}_{\mathcal{T}}^{m}$.

Lemma 12

Let the bilinear form $a\left(\cdot,\cdot\right)$ satisfy (1) and let the CFL condition

[TABLE]

*hold.

Then, the bilinear form $a_{p}\left(\cdot,\cdot\right)$ is continuous,*

[TABLE]

*with *

[TABLE]

and symmetric, $a_{p}\left(u,v\right)=a_{p}\left(v,u\right)$ for all $u,v\in S.$ Moreover, for any $f\in L^{2}\left(\Omega\right)$, the problem: Find $u\in S$ such that

[TABLE]

has a unique solution, which satisfies

[TABLE]

Remark 13

In (19) the condition on the time-step $\Delta t$ implies that $\Delta t$ is essentially proportional to $h$ and inversely proportional to $\sqrt{C_{\operatorname*{cont}}}$, as $c_{\operatorname*{coer}}\leq C_{\operatorname*{cont}}$. Hence (19) corresponds to a genuine CFL condition since $\sqrt{C_{\operatorname*{cont}}}$ usually corresponds to the maximal (physical) wave speed.

Proof of Lemma 12. If $p=1$, the two bilinear forms $a_{p}$ and $a$ coincide and the result trivially follows. Thus, we now assume that $p\geq 2$.

**a) Continuity. **Let $u,v\in S$ and

[TABLE]

Then, by definition of $a_{p}$ and continuity of $a$, we have

[TABLE]

By applying the triangle inequality to (21) we obtain

[TABLE]

From (1), it follows that

[TABLE]

Hence,

[TABLE]

with

[TABLE]

The operator $A_{S}^{1/2}R_{\mathcal{N}}A_{S}^{1/2}$ is self-adjoint with respect to the $L^{2}\left(\Omega\right)$ scalar product and positive semi-definite. It is well-known that under these conditions we have

[TABLE]

From (17) we conclude that the spectrum $\sigma\left(A_{S}^{1/2}R_{\mathcal{N}}A_{S}^{1/2}\right)$ is contained in the interval $\left[0,\frac{C_{\operatorname*{cont}}}{c_{\operatorname*{eq}}^{2}}\left(1+\frac{C_{\operatorname*{inv}}^{2}C_{\operatorname*{qu}}^{2}}{h^{2}}\right)\right]$ so that

[TABLE]

with $\kappa$ as in (20). The CFL condition (19), together with the continuity and the coercivity of $a$ and $p\geq 2$, implies $\kappa\in\left[0,4p^{2}\right]$. Thus, Lemma 19 (Appendix) implies

[TABLE]

which we insert in (22) to obtain

[TABLE]

**b) Symmetry. **This follows since $A_{S}$, $R_{\mathcal{N}}$ are self-adjoint with respect to the $L^{2}\left(\Omega\right)$ scalar product.

**c) Coercivity. **Note that the problem: Find $u\in S$ such that

[TABLE]

can be solved in two steps: Find $w\in S$ such that

[TABLE]

Then $u$ is the solution of

[TABLE]

By the similar arguments as in the first part of this proof, one concludes that the CFL-condition (19) implies

[TABLE]

so that

[TABLE]

The well-posedness of problem (24) follows from the Lax-Milgram lemma as well as the estimate

[TABLE]

Corollary 14

The bilinear form $a_{p}\left(u,v\right)$ is symmetric, continuous and coercive. Hence, there exists an $L^{2}\left(\Omega\right)$-orthonormal eigensystem $\left(\lambda_{S,p,k},\eta_{S,p,k}\right)_{k=1}^{M}$ for $a_{p}\left(\cdot,\cdot\right)$, i.e.,

[TABLE]

with real and positive eigenvalues $\lambda_{S,p,k}>0$. Let the CFL condition (19) be satisfied. Then, the smallest and largest eigenvalue satisfy

[TABLE]

Proof. We start with the smallest eigenvalue. It holds

[TABLE]

with $\kappa$ as in (20). Hence,

[TABLE]

The CFL condition (19) implies

[TABLE]

[TABLE]

which yields the lower bound on the smallest eigenvalue $\lambda_{p}^{\min}$.

For the largest eigenvalue $\lambda_{p}^{\max}$, we get by using the CFL condition and (14) that

[TABLE]

from which the upper bound on $\lambda_{p}^{\max}$ follows.

Corollary 15

Let the assumptions of Lemma 12 be satisfied. Then

[TABLE]

uniformly in $p$.

Proof. We write

[TABLE]

Note that for all $w\in S$ it holds

[TABLE]

Since $R_{\mathcal{N}}$ is symmetric, positive semi-definite (see Remark 3), we infer from (16) that $\left\|R_{\mathcal{N}}^{1/2}v\right\|\leq c_{\operatorname*{eq}}^{-1}\left\|v\right\|$ holds for all $v\in S$. From Lemmas 9 and 10 we obtain for all $v\in S$

[TABLE]

Thus, we argue as for (22) and get

[TABLE]

with

[TABLE]

From Lemma 19 we conclude that $C_{p}^{\prime}\leq(p^{2}-1)/12\leq p^{2}/12$ so that (19) implies

[TABLE]

Thus, we have proved

[TABLE]

From (1c) we conclude that

[TABLE]

which together with (27) leads to the assertion.

3.2 Error equation and estimates

To derive a priori error estimates for the LTS/FE-Galerkin solution of (12), we first introduce the new function

[TABLE]

and rewrite (12) as a one-step method

[TABLE]

The elimination of $v_{S}^{\left(n+1/2\right)}$ from the second equation by using the first one leads to the operator equation

[TABLE]

[TABLE]

with $A_{S,p}$ as in (10), $f_{S}^{\left(n\right)}$ as in (8), and

[TABLE]

Next, we will derive a recursion for the error

[TABLE]

where $u$ is the solution of (2)-(3) and $v$ the solution of the corresponding first-order formulation: Find $u,v:[0,T]\rightarrow V$ such that

[TABLE]

and initial conditions $u(0)=u_{0}$ and $v(0)=v_{0}.$

To split the error we introduce the first-order formulation of the semi-discrete problem (7). Find $u_{S},v_{S}:\left[0,T\right]\rightarrow S$ such that

[TABLE]

Hence, we may write $\mathbf{e}^{\left(n+1\right)}:=\left(e_{v}^{\left(n+\frac{1}{2}\right)},e_{u}^{\left(n+1\right)}\right)^{\intercal}=\mathbf{e}_{S}^{\left(n+1\right)}+\mathbf{e}_{S,\Delta t}^{\left(n+1\right)}$ with

[TABLE]

We first investigate the error $\mathbf{e}_{S,\Delta t}^{\left(n+1\right)}$ and introduce

[TABLE]

[TABLE]

These equations can be written in the form

[TABLE]

By subtracting the first equation in (29) from (43) and the second equation in (29) from (44) we obtain

[TABLE]

Eliminating the term $e_{v,S,\Delta t}^{\left(n+1/2\right)}$ in the second equation by using the first one yields

[TABLE]

We rewrite it in operator form by using the operator $\mathfrak{S}$ as in (30)

[TABLE]

with

[TABLE]

This recursion can be resolved

[TABLE]

Let I_{S}^{2\times 2}:=\left[\begin{array}[c]{cc}I_{S}&0\\ 0&I_{S}\end{array}\right] and observe that

[TABLE]

and

[TABLE]

We introduce

[TABLE]

and the differences

[TABLE]

and use (3.2) to rewrite the error representation (3.2) as

[TABLE]

3.2.1 Stability

As usual, the convergence analysis can be split into an estimate for the stability of the iteration operator $\mathfrak{S}$ (corresponding to a homogeneous right-hand side) and a consistency estimate. We begin with the analysis of the stability.

Theorem 16** (Stability)**

Let the CFL condition (19) be satisfied. Then the leap-frog scheme (12) is stable

[TABLE]

where $C_{0}$ is independent of $n$, $\Delta t$, $h$, and $T$.

Proof. We choose the eigensystem as introduced in Corollary 14 and expand

[TABLE]

Inserting this into the recursion \left(\begin{array}[c]{c}v_{S}^{\left(n+1/2\right)}\\ u_{S}^{\left(n+1\right)}\end{array}\right)=\mathfrak{S}\left(\begin{array}[c]{c}v_{S}^{\left(n-1/2\right)}\\ u_{S}^{\left(n\right)}\end{array}\right) leads to a recursion for the coefficients $\beta_{S,p,k}^{\left(n+1/2\right)}$, $\chi_{S,p,k}^{\left(n+1\right)}$:

[TABLE]

with

[TABLE]

The eigenvalues of $\mathbf{S}_{p}$ are given by

[TABLE]

The CFL condition (19) implies $\left(\Delta t\right)^{2}\lambda_{p}^{\max}<4$ so that the eigenvalues are different and $\mathbf{S}_{p}$ is diagonalizable. From [45, Satz (6.9.2)(2)] we conclude that there is a norm $\left|\kern-1.00006pt\left|\kern-1.00006pt\left|\cdot\right|\kern-1.00006pt\right|\kern-1.00006pt\right|$ in $\mathbb{R}^{2}$ such that the associated matrix norm $\left|\kern-1.00006pt\left|\kern-1.00006pt\left|\mathbf{S}_{p}\right|\kern-1.00006pt\right|\kern-1.00006pt\right|$ is bounded from above by the spectral radius:

[TABLE]

Hence

[TABLE]

Since all norms in $\mathbb{R}^{2}$ are equivalent there exists a constant $C$ such that

[TABLE]

The eigenfunctions $\eta_{S,p,k}$ are chosen to be an orthonormal system in $L^{2}\left(\Omega\right)$ so that

[TABLE]

which shows the $L^{2}\left(\Omega\right)$-stability of the method.

3.2.2 Error Estimates

In this section we first estimate the discrete error $e_{u,S,\Delta t}^{\left(n+1\right)}$. Standard estimates on the semi-discrete error then lead to an estimate of the total error $e_{u}^{(n+1)}$.

Theorem 17

Let the assumptions of Lemma 12 be satisfied. Let the solution of the semi-discrete equation (7) satisfy $u_{S}\in W^{5,\infty}\left(\left[0,T\right];L^{2}\left(\Omega\right)\right)$ and the right-hand side $f_{S}\in W^{3,\infty}\left(\left[0,T\right];L^{2}\left(\Omega\right)\right)$. Then the fully discrete solution $u_{S}^{\left(n+1\right)}$ of (12) satisfies the error estimate

[TABLE]

with

[TABLE]

and a constant $C$ which is independent of $n$, $\Delta t$, $T$, $h$, $p$, $f_{S}$, and $u_{S}$.

Proof. We apply the stability estimate to the second component of the error representation (61). From Theorem 16 and (49) we obtain222For a pair of functions $\mathbf{v}=\left(v_{1},v_{2}\right)^{\intercal}\in S^{2}$ we use the notation $\left\|\mathbf{v}\right\|_{\ell^{1}}:=\left\|v_{1}\right\|+\left\|v_{2}\right\|$.

[TABLE]

For the summands in the second term of the right-hand side in (66), we obtain by a Taylor argument and Corollary 15

[TABLE]

with

[TABLE]

and

[TABLE]

Now, let $\psi$ denote the second component of the first term in the right-hand side of (67),

[TABLE]

By using $\ddot{u}_{S}+A_{S}u_{S}=f_{S}$ (cf. (7a) and (10)) we obtain

[TABLE]

We employ (27) and argue as in the proof of Corollary 15 to obtain

[TABLE]

This yields

[TABLE]

In summary we have proved

[TABLE]

Next, we estimate the remaining terms in (66). We employ the discrete wave equation and a Taylor argument to obtain

[TABLE]

The estimate of the last term in (66) follows by setting $n=1$ in (68)

[TABLE]

Inserting these estimates into (66) leads to

[TABLE]

It remains to estimate the initial error $\mathbf{e}_{S,\Delta t}^{\left(1\right)}$. Let $u_{S}^{\left(0\right)}:=u_{S}\left(0\right)$ and $v_{S}^{\left(0\right)}:=\dot{u}_{S}\left(0\right)\in S$ be as in (7b). A Taylor argument for some $0\leq\theta\leq\tau\leq\Delta t$ and the definition of $u_{S}^{\left(0\right)}$, $u_{S}^{\left(1\right)}$ as in (12) lead to

[TABLE]

For the initial error in $v_{S}$ we obtain by a similar Taylor argument

[TABLE]

In summary, we have estimated the initial error by

[TABLE]

The combination of (71) and (76) leads to the assertion.

Theorem 17 can be combined with known error estimates for the semi-discrete error $\mathbf{e}_{S}^{\left(n+1\right)}$ to obtain an error estimate of the total error.

Theorem 18

Let the bilinear form $a\left(\cdot,\cdot\right)$ satisfy (1) and let the CFL condition (19) hold. Assume that the exact solution satisfies $u\in W^{1,\infty}\left(\left[0,T\right];H^{m+1}\left(\Omega\right)\right)\cap W^{5,\infty}\left(\left[0,T\right];L^{2}\left(\Omega\right)\right)$. Then, the corresponding fully discrete Galerkin FE formulation with local time-stepping (12) has a unique solution $u_{S}^{\left(n+1\right)}$ which satisfies the error estimate

[TABLE]

with

[TABLE]

and a constant $C$ which is independent of $n$, $\Delta t$, $h$, $p$, $f_{S}$, $u_{S}$, and the final time $T$.

Proof. The existence of the semi-discrete solution $u_{S}$ follows from [3, Theorem 3.1], which directly implies the existence of our fully discrete LTS-Galerkin FE solution.

Next, we split the total error

[TABLE]

according to (36). Following [40], we note that the semi-discrete solution $u_{S}$ inherits the same regularity from $u\in W^{5,\infty}\left(\left[0,T\right];L^{2}\left(\Omega\right)\right)$; thus, we can apply Theorem 17.

To estimate the remaining error from the semi-discretization,

[TABLE]

we use [3, Theorem 3.1] to obtain

[TABLE]

Inspection of the proof in [3, Theorem 3.1] shows that the constant in (77) can be estimated by $C\left(1+\sqrt{T}\right)$. Using a Hölder inequality in the second summand of the right-hand side in (77) thus results in

[TABLE]

from which we conclude that

[TABLE]

with a constant $C^{\prime}$ which is independent of the final time $T$. Finally, the triangle inequality leads to the assertion.

4 Numerical Experiments

Numerical experiments that corroborate the convergence rates and illustrate the stability properties of the LTS-LF scheme when combined with continuous or discontinuous Galerkin FEM [28] were presented in [18]. Together with its higher order versions, the LTS-LF method was also successfully applied to other (vector-valued) second-order wave equations from electromagnetics [26] and elasticity [36, 42] . Here we demonstrate the versatility of the LTS approach in the presence of adaptive mesh refinement near a re-entrant corner.

To illustrate the usefulness of the LTS approach, we consider the classical scalar wave equation (Example 1) in the L-shaped domain $\Omega$ shown in Fig. 1. The re-entrant corner is located at $(0.5,0.5)$ and we set $c=1$, $f=0$ and the final time $T=2$. Next, we impose homogeneous Neumann boundary conditions on all boundaries and choose as initial conditions the vertical Gaussian plane wave

[TABLE]

of width $\delta=10^{-5}$ centered about $x_{0}=0.25$ . For the spatial discretization we opt for ${\mathcal{P}}^{2}$ continuous finite elements with mass lumping [10].

First, we partition $\Omega$ into equal triangles of size $h_{\mbox{\scriptsize init}}$ – see Fig. 1 (a). Then we bisect the six elements nearest to the corner and subsequently bisect in the resulting mesh all elements with a vertex at $(0.5,0.5)$. Starting from that intermediate mesh, shown in Fig. 1 (b), we repeat this procedure again with the six elements adjacent to the corner, which finally yields the mesh shown in Fig. 1 (c). Hence the mesh refinement ratio, that is the ratio between smallest elements in the ”coarse” and the ”fine” regions, in the resulting mesh is 4:1. We therefore choose a four times smaller time-step $\Delta\tau=\Delta t/p$ with $p=4$ inside the fine region.

Clearly, this refinement strategy is heuristic, as optimal mesh refinement in the presence of corner singularities generally requires hierarchical mesh refinement [39]. However, when the region of local mesh refinement itself contains a sub-region of even smaller elements, and so forth, any local time-step will again be overly restricted due to even smaller elements inside the ”fine” region. To remedy the repeated bottleneck caused by hierarchical mesh refinement, multi-level local time-stepping methods were proposed in [19, 42], which permit the use of the appropriate time-step at every level of mesh refinement. For simplicity, we restrict ourselves here to the standard (two-level) LTS-LF scheme.

In Fig. 2 we display snapshots of the numerical solution at different times: the plane wave splits into two wave fronts travelling in opposite directions. The lower half of the right propagating wave is reflected while the upper half proceeds into the upper left quadrant. To avoid any loss in the global CFL condition and reach the optimal global time-step, we always include an overlap by one element, that is, we also advance the numerical solution inside those elements immediately next to the ”fine” region with the fine time-step.

In Fig. 3 we compare the runtime of the LTS-LF($p$) on a sequence of meshes using the refinement strategy depicted in Fig. 1, with the runtime of a standard LF scheme with a time-step $\Delta t/4$ on the entire domain. As expected, the LTS-LF method is faster than the standard LF scheme, in fact increasingly so, as the number of refinements increases. Indeed, as the number of degrees of freedom in the ”coarse” region grows much faster than in the ”fine” region, where it remains essentially constant, the use of local time-stepping becomes increasingly beneficial on finer meshes.

Acknowledgements

We thank Loredana Gaudio for useful comments and suggestions during the initial stages of this work and Maximillian Matthäus for his Matlab program.

Appendix A Some Auxiliary Estimates

Lemma 19

For $p\geq 2$ let $\alpha_{j}^{p}$, $j=1,\dots,p-1$, be recursively defined as in (11). Then, the constants $\alpha_{j}^{p}$ are given by

[TABLE]

Moreover, for $\kappa\in\left[0,4p^{2}\right]$ it holds

[TABLE]

Proof. To show that the constants $\alpha_{j}^{p}$ are in fact given by (78), we first use the identity

[TABLE]

to rewrite (78) as

[TABLE]

By using (80) it is then straightforward to verify that $\alpha_{j}^{p}$ satisfies the recursive definition in (11).

Next, one proves by induction that

[TABLE]

with the Čebyšev polynomials $T_{p}$ of the first kind. We recall that

[TABLE]

where the first relation follows from [43, (1.97)] and the second one from [43, Theorem 2.24], see also [44, Corollary 7.3.1].

Now, let $x=\kappa/p^{2}$. The condition $\kappa\in\left[0,4p^{2}\right]$ implies $\left[1-\frac{x}{2},1\right]\subset\left[-1,1\right]$. Hence, a Taylor argument shows that there exists $\xi\in\left[-1,1\right]$ such that

[TABLE]

where we have also used (81). Similarly, we get

[TABLE]