Contrast-Independent Partially Explicit Time Discretizations for Quasi Gas Dynamics

Abstract

In the paper, we study a design and stability of contrast-independent partially explicit time discretizations for Quasi-Gas-Dynamics (QGD) Equations in multiscale high-contrast media. In our previous works, we have introduced contrast-independent partially explicit time discretizations. In this paper, we extend these ideas to multiscale QGD problems. Because of high contrast, explicit methods require a very small time stepping. By designing appropriate spatial splitting and temporal splitting, partially explicit methods remove this constraint. The proposed partially explicit time discretization consists of two steps. First, we split the space into contrast dependent (fast) and contrast independent (slow) components on a coarse grid that is much larger compared to spatial heterogeneities. Secondly, we design a temporal splitting algorithm in a such way that it is stable and the time step is independent of the contrast and only depends on the coarse mesh size. Using proposed method, a few degrees of freedom are treated implicitly and the approach is mostly explicit. We prove that the proposed splitting is unconditionally stable under some suitable conditions formulated for the second space (slow). We present numerical results and show that the proposed methods provide results similar to implicit methods with the time step that is independent of the contrast.

1. Introduction

Complex flows that combine different processes occur in many applications. Some important classes of flow models include kinetic models and continuum models. There are several intermediate-scale and intermediate-physics models that are used in the literature. One of such models is the quasi-gas dynamic (QGD) system of equations, which has a form

$$\frac{\partial u}{\partial t}+\alpha \frac{{\partial }^{2}u}{\partial t^2}+A\left(u\right)=f,$$

which combines both parabolic and wave phenomena. The QGD model is multiphysics as it includes both diffusive and wave phenomena, QGD model equations can be derived from kinetic equations. This derivation assumes that the distribution function is similar to a locally Maxwellian representation. The QGD model has a smoothing effect for the solution at the free path distance. One can find more detailed information about QGD model in the literature [1]. In this paper, a simplified QGD system is considered. The model involves second derivatives with respect to the time. The resulting model regularizes parabolic equations by adding a hyperbolicity. The latter can be used in designing efficient time stepping algorithms [1].

In the paper, we consider the QGD model in a multiscale and high-contrast environment (see Equation (2)). The equation parameters represent the media properties and spatially vary. The applications of these equations can be considered in porous media for compressible flows. The heterogeneities of the coefficients represent the media properties, which can have large variations. Because of large variations, the resulting system is stiff and requires very small time steps for explicit methods. In this paper, we use our previous works on multiscale methods [2] to design efficient temporal splitting (cf. [3,4]). The resulting scheme is a partially explicit method and has many advantages in simulations.

Explicit methods are commonly used for QGD models [1] due to the finite speed of propagation. Explicit methods have advantages: (1) they provide a fast marching in time; (2) they can easily preserve physical quantities; and (3) they have small communication between degrees of freedom (DOF). However, explicit methods require a small time step which is due to the mesh size and the contrast. Implicit methods typically give unconditionally stable schemes, but with a higher computational cost. It is therefore desirable to develop a practical compromise that takes advantage of both explicit and implicit methods. We propose to solve the problem on a coarse grid, where the mesh size is chosen to be much larger compared to heterogeneities. It is important to remove the contrast dependency in the time stepping. For this, we extend recently developed novel splitting algorithm [3,4] and design a partial explicit approach, which gives unconditional stability under some suitable conditions. Our partial explicit discretization (under some conditions, see (15)) splits the solution into two parts $u=u_1+u_2$ (formally, understood as projections to corresponding spaces) and can be formally written as

$$\begin{array}{c}\hfill \frac{u_1^{n+1}-u_1^{n-1}}{\Delta t}+\alpha \frac{u_1^{n+1}-2u_1^n+u_1^{n-1}}{\Delta t^2}+\frac{1}{2}A(u_1^{n+1}+u_1^{n-1}+2u_2^n)=f_1\\ \hfill \frac{u_2^{n+1}-u_2^{n-1}}{\Delta t}+\alpha \frac{u_2^{n+1}-2u_2^n+u_2^{n-1}}{\Delta t^2}+A(u_1^n+u_2^n)=f_2.\end{array}$$

(1)

We see that this discretization has a special form, which is needed for the stability. Because the problems are solved on a coarse grid that is much larger compared to heterogeneities, we use multiscale methods.

The proposed approach extends our previous work on partial explicit methods for parabolic and wave Equations [3,4] to deal with QGD models. These approaches use multiscale space decomposition coupled with temporal splitting to perform partial explicit discretization. In particular, it uses earlier developed spatial multiscale methods for the spatial discretization to identify the degrees of freedom that require implicit treatment. These degrees of freedom typically correspond to fast time scales and their identification requires a careful spatial decomposition. Furthermore, we propose a splitting algorithm. We note that proposed partial explicit approach for parabolic and wave equations differ in their form of discretization. For this reason, a different discretization is required for splitting that is stable. In the paper, we discuss such approach. Next, we briefly discuss some ingredients of our method, which include multiscale methods and temporal splitting algorithms.

In the paper, we use Constraint Energy Minimizing Generalized Multiscale Finite Element Method (CEM-GMsFEM) [5,6], though other approaches can possibly be also used. We note that many multiscale methods have been developed and analyzed. Some of them formulate coarse-grid problems using effective media properties and based on homogenization [7,8,9]. Other approaches are based on constructing multiscale basis functions and formulating coarse-grid equations. These include multiscale finite element methods [9,10,11], generalized multiscale finite element methods (GMsFEM) [12,13,14,15], constraint energy minimizing GMsFEM (CEM-GMsFEM) [5,6], nonlocal multi-continua (NLMC) approaches [16], metric-based upscaling [17], heterogeneous multiscale method [18], localized orthogonal decomposition (LOD) [19], equation-free approaches [20,21], and so on. Because proposed problems are high contrast, GMsFEM and NLMC are proposed to extract macroscopic quantities associated with the degrees of the freedom that the operator “cannot see”. For this reason, for GMsFEM and related approaches [5], multiple basis functions or continua are constructed to capture the multiscale features due to high contrast [6,16]. As we mentioned, the contrast introduces a stiffness in the forward problems. When treating explicitly, one needs to take very small time steps when the contrast is high. Our approach allows taking the time step to be independent of the contrast by handling some degrees of freedom implicitly and some explicitly.

Our approaches use some fundamental splitting algorithms [22,23]. Unlike algorithms that split various physics, our algorithms identify spatial components of the solution that corresponds to fast and slow time scales and treat them implicitly and explicitly, correspondingly. Though our proposed approaches share some common concepts with implicit-explicit approaches (e.g., [24,25,26,27,28,29,30,31,32,33,34,35,36]), they differ from these approaches. Previous works on explicit-implicit methods treat nonlinearities in an explicit fashion, other methods treat globally defined stiff terms implicitly. Our approaches identify local features that need implicit treatment and design splitting algorithms for implicit-explicit discretization on a coarse grid. Another important aspect is that the explicit time step scales with the coarse-mesh size.

We summarize our contributions.

• We design a partial explicit scheme for QGD model, which differs from previously schemes in [3,4] and use both centered difference and one-sided differences.
• We use multiscale spatial decomposition that identifies fast and slow components of the solution.
• Proposed approach uses implicit treatment for a few degrees of freedom and the rest is treated explicitly. The resulting time constraint is independent of the contrast and scales as the coarse-mesh size.

In the paper, we present some representative numerical results. We consider a heterogeneous conductivity field and a source term that has a small spatial support, which introduces additional scales. In our numerical results, we consider various values of $\alpha$ that represents the interaction between wave and parabolic phenomena. In all cases, our numerical results show that the proposed partial explicit approach performs very similar to fully implicit approach when handing additional degrees of freedom. In particular, in a strong mixed cases, our approach also performs well.

The paper is organized as follows. In next section, we present preliminaries. Section 3 is devoted to partial explicit method description and analysis. In Section 4, we present numerical results.

2. Preliminaries

We consider the QGD model in heterogeneous domain. The problem consists of finding u such that

$$\frac{\partial }{\partial t}u+\alpha \frac{{\partial }^2}{\partial t^2}u=\nabla ·(\kappa \nabla u)\mathrm{in}\Omega ,$$

(2)

where $\kappa \in L^{\infty }(\Omega )$ is a high contrast heterogeneous field and $\alpha$ is constant (all our findings apply for heterogeneous spatial field $\alpha \left(x\right)$). We assume both $\kappa \left(x\right)$ and $\alpha \left(x\right)$ are positive. Equation (2) is equipped with initial and boundary conditions, $u(0,·)=u_0\left(x\right)$, $u_t(0,·)=u_{00}\left(x\right)$, and $u(t,x)=g\left(x\right)$ on $\partial \Omega$.

The weak formulation of the problem is to find $u(t,·)\in V:=H^1(\Omega )$ such that

$$(\frac{\partial }{\partial t}u,v)+\alpha (\frac{{\partial }^2}{\partial t^2}u,v)=-a(u,v)\forall v\in V,$$

(3)

where

$$a(u,v)={\int }_{\Omega }\kappa \nabla u·\nabla v.$$

We take homogeneous boundary conditions.

We consider the energy

$$E\left(u\right)=\frac{1}{2}\left({\int }_{\Omega }\kappa {|\nabla u|}^2+\alpha {\int }_{\Omega }{\left(\frac{\partial u}{\partial t}\right)}^2\right).$$

Then, if we assume $f=0$, we have (after multiplying by $\frac{\partial u}{\partial t}$.

$${\int }_{\Omega }\frac{\partial u}{\partial t}\frac{\partial u}{\partial t}+\alpha {\int }_{\Omega }\frac{{\partial }^{2}u}{\partial t^2}\frac{\partial u}{\partial t}+{\int }_{\Omega }\kappa \nabla u·\nabla \frac{\partial u}{\partial t}=0.$$

From here, we have

$$E\left(u\right)\left(t\right)=E\left(u\right)\left(0\right)-{\int }_{\Omega }{u}_t^2\le E\left(u\right)\left(0\right).$$

(4)

Our goal is to preserve this energy inequality in the full temporal discretization.

3. Partially Explicit Temporal Splitting Scheme

In this section, we will discuss a temporal Splitting scheme for solving problem (2). We consider a numerical solution $u_H$ is defined as

$$(\frac{\partial u_H}{\partial t},v)+\alpha (\frac{{\partial }^2{u}_H}{\partial t^2},v)+a(u_H,v)=(f,v)\forall v\in V_H$$

(5)

where $V_H$ is a coarse grid finite element space. We consider $V_H$ can be decomposed into two subspaces $V_{H,1}$, $V_{H,2}$ namely,

$$V_H=V_{H,1}\oplus V_{H,2}.$$

We will use a time discretization scheme: finding ${\left\{u_{H,1}^n\right\}}_{n=1}^N\in V_{H,1},{\left\{u_{H,2}^n\right\}}_{n=1}^N\in V_{H,2}$

$$\begin{array}{c}\hfill (\frac{u_H^{n+1}-u_H^{n-1}}{\Delta t},v)+\alpha (\frac{u_H^{n+1}-2u_H^n+u_H^{n-1}}{\Delta t^2},v)+\frac{1}{2}a(u_{H,1}^{n+1}+u_{H,1}^{n-1}+2u_{H,2}^n,v)\\ \hfill =(f,v),v\in V_{H,1}\\ \hfill (\frac{u_H^{n+1}-u_H^{n-1}}{\Delta t},v)+\alpha (\frac{u_H^{n+1}-2u_H^n+u_H^{n-1}}{\Delta t^2},v)+a(u_{H,1}^n+u_{H,2}^n,v)\\ \hfill =(f,v),v\in V_{H,2}\end{array}$$

(6)

where the numerical solution $u_H^n\in V_H$ is the sum of $u_{H,1}^n$ and $u_{H,2}^n$, namely $u_H^n=u_{H,1}^n+u_{H,2}^n$.

3.1. Stability

We define

$$E^{n+1}=\alpha \parallel \frac{u_H^{n+1}-u_H^n}{\Delta t}{\parallel }^2+\frac{1}{2}\left(\parallel u_{H,1}^{n+1}+u_{H,2}^n{\parallel }_a^2+{\parallel u_{H,2}^{n+1}+u_{H,1}^n\parallel }_a^2\right)-\frac{{(\Delta t)}^2}{2}{\parallel \frac{u_{H,2}^{n+1}-u_{H,2}^n}{\Delta t}\parallel }_a^2$$

Lemma 1.

For $f=0$, we have

$$E^{n+1}\le E^{n+1}+\Delta t^{-1}{\parallel u_H^{n+1}-u_H^{n-1}\parallel }^2=E^n.$$

Proof. 

We consider the test function $u_{H,1}^{n+1}-u_{H,1}^{n-1}$ and $u_{H,2}^{n+1}-u_{H,2}^{n-1}$. We have

$$\begin{array}{c}\hfill (\frac{u_H^{n+1}-u_H^{n-1}}{\Delta t}+\alpha \frac{u_H^{n+1}-2u_H^n+u_H^{n-1}}{\Delta t^2},u_{H,1}^{n+1}-u_{H,1}^{n-1})+\\ \hfill \frac{1}{2}a(u_{H,1}^{n+1}+u_{H,1}^{n-1}+2u_{H,2}^n,u_{H,1}^{n+1}-u_{H,1}^{n-1})=0,v\in V_{H,1}\\ \hfill (\frac{u_H^{n+1}-u_H^{n-1}}{\Delta t}+\alpha \frac{u_H^{n+1}-2u_H^n+u_H^{n-1}}{\Delta t^2},u_{H,2}^{n+1}-u_{H,2}^{n-1})+\\ \hfill a(u_{H,1}^n+u_{H,2}^n,u_{H,2}^{n+1}-u_{H,2}^{n-1})=0,v\in V_{H,2}\end{array}$$

(7)

We have

$$\begin{array}{cc}\hfill & (\frac{u_H^{n+1}-u_H^{n-1}}{\Delta t}+\alpha \frac{u_H^{n+1}-2u_H^n+u_H^{n-1}}{\Delta t^2},u_H^{n+1}-u_H^{n-1})+\frac{1}{2}a(u_{H,1}^{n+1}+u_{H,1}^{n-1},u_{H,1}^{n+1}-u_{H,1}^{n-1})+\hfill \\ \hfill & a(u_{H,2}^n,u_{H,1}^{n+1}-u_{H,1}^{n-1})+a(u_{H,1}^n+u_{H,2}^n,u_{H,2}^{n+1}-u_{H,2}^{n-1})=0.\hfill \end{array}$$

(8)

We first have

$$\begin{array}{c}\hfill (\frac{u_H^{n+1}-u_H^{n-1}}{\Delta t}+\alpha \frac{u_H^{n+1}-2u_H^n+u_H^{n-1}}{\Delta t^2},u_H^{n+1}-u_H^{n-1})=\Delta t^{-1}{\parallel u_H^{n+1}-u_H^{n-1}\parallel }^2+\\ \hfill \alpha \Delta t^{-2}\left(\parallel u_H^{n+1}-u_H^n{\parallel }^2-{\parallel u_H^n-u_H^{n-1}\parallel }^2\right)\end{array}$$

(9)

We next estimate the term $\frac{1}{2}a(u_{H,1}^{n+1}+u_{H,1}^{n-1},u_{H,1}^{n+1}-u_{H,1}^{n-1})$ and have

$$\begin{array}{c}\hfill \frac{1}{2}a(u_{H,1}^{n+1}+u_{H,1}^{n-1},u_{H,1}^{n+1}-u_{H,1}^{n-1})=\frac{1}{2}\left(\parallel u_{H,1}^{n+1}{\parallel }_a^2-{\parallel u_{H,1}^{n-1}\parallel }_a^2\right)\\ \hfill =\frac{1}{2}\left(\parallel u_{H,1}^{n+1}{\parallel }_a^2+{\parallel u_{H,1}^n\parallel }_a^2\right)-\frac{1}{2}\left(\parallel u_{H,1}^n{\parallel }_a^2+{\parallel u_{H,1}^{n-1}\parallel }_a^2\right).\end{array}$$

(10)

We next consider the term $a(u_{H,2}^n,u_{H,1}^{n+1}-u_{H,1}^{n-1})+a(u_{H,1}^n,u_{H,2}^{n+1}-u_{H,2}^{n-1})+a(u_{H,2}^n,u_{H,2}^{n+1}-u_{H,2}^{n-1})$ and obtain

$$\begin{array}{c}\hfill a(u_{H,2}^n,u_{H,1}^{n+1}-u_{H,1}^{n-1})+a(u_{H,1}^n,u_{H,2}^{n+1}-u_{H,2}^{n-1})+a(u_{H,2}^n,u_{H,2}^{n+1}-u_{H,2}^{n-1})\\ \hfill =\left(a(u_{H,2}^{n+1},u_{H,1}^n)+a(u_{H,1}^{n+1},u_{H,2}^n)+a(u_{H,1}^{n+1},u_{H,2}^n)\right)-(a(u_{H,2}^n,u_{H,1}^{n-1})+\\ \hfill a(u_{H,1}^n,u_{H,2}^{n-1})+a(u_{H,2}^n,u_{H,2}^{n-1})).\end{array}$$

(11)

Thus, we have

$$\begin{array}{cc}\hfill & \Delta t^{-1}{\parallel u_H^{n+1}-u_H^{n-1}\parallel }^2+\alpha \Delta t^{-2}\left(\parallel u_H^{n+1}-u_H^n{\parallel }^2\right)+\frac{1}{2}\left(\parallel u_{H,1}^{n+1}{\parallel }_a^2+{\parallel u_{H,1}^n\parallel }_a^2\right)\hfill \\ \hfill & +\left(a(u_{H,2}^{n+1},u_{H,1}^n)+a(u_{H,1}^{n+1},u_{H,2}^n)+a(u_{H,1}^{n+1},u_{H,2}^n)\right)\hfill \\ \hfill =& \alpha \Delta t^{-2}{\parallel u_H^n-u_H^{n-1}\parallel }^2+\frac{1}{2}\left(\parallel u_{H,1}^n{\parallel }_a^2+{\parallel u_{H,1}^{n-1}\parallel }_a^2\right)+\hfill \\ \hfill & \left(a(u_{H,2}^n,u_{H,1}^{n-1})+a(u_{H,1}^n,u_{H,2}^{n-1})+a(u_{H,2}^n,u_{H,2}^{n-1})\right).\hfill \end{array}$$

(12)

We observe that

$$\begin{array}{cc}\hfill & \frac{1}{2}\left(\parallel u_{H,1}^{k+1}{\parallel }_a^2+{\parallel u_{H,1}^k\parallel }_a^2\right)+\left(a(u_{H,2}^{k+1},u_{H,1}^k)+a(u_{H,1}^{k+1},u_{H,2}^k)+a(u_{H,2}^{k+1},u_{H,2}^k)\right)\hfill \\ \hfill =& \frac{1}{2}\sum _{m=1,2}\left(\parallel u_{H,m}^{k+1}{\parallel }_a^2+{\parallel u_{H,m}^k\parallel }_a^2\right)+\left(a(u_{H,2}^{k+1},u_{H,1}^k)+a(u_{H,1}^{k+1},u_{H,2}^k)\right)-\frac{1}{2}{\parallel u_{H,2}^{k+1}-u_{H,2}^k\parallel }_a^2\hfill \\ \hfill =& \frac{1}{2}\left(\parallel u_{H,1}^{k+1}+u_{H,2}^k{\parallel }_a^2+{\parallel u_{H,2}^{k+1}+u_{H,1}^k\parallel }_a^2\right)-\frac{1}{2}{\parallel u_{H,2}^{k+1}-u_{H,2}^k\parallel }_a^2\hfill \end{array}$$

(13)

for $k=n,n-1$. Therefore, we have

$$\begin{array}{cc}\hfill E^{n+1}& \le E^{n+1}+\Delta t^{-1}{\parallel u_H^{n+1}-u_H^{n-1}\parallel }^2\hfill \\ \hfill & =E^n.\hfill \end{array}$$

(14)

□

We will discuss two cases. In the first case, we assume $V_{H,1}$ and $V_{H,2}$ are orthogonal and in the second case, we will consider the case when they are not orthogonal.

3.2. Case $V_{H,1}$ and $V_{H,2}$ Are Orthogonal

When $V_{H,1}$ and $V_{H,2}$ are orthogonal, the partial explicit scheme has a form

$$\begin{array}{c}\hfill (\frac{u_{H,1}^{n+1}-u_{H,1}^{n-1}}{\Delta t},v)+\alpha (\frac{u_{H,1}^{n+1}-2u_{H,1}^n+u_{H,1}^{n-1}}{\Delta t^2},v)+\frac{1}{2}a(u_{H,1}^{n+1}+u_{H,1}^{n-1}+2u_{H,2}^n,v)\\ \hfill =(f,v),v\in V_{H,1}\\ \hfill (\frac{u_{H,2}^{n+1}-u_{H,2}^{n-1}}{\Delta t},v)+\alpha (\frac{u_{H,2}^{n+1}-2u_{H,2}^n+u_{H,2}^{n-1}}{\Delta t^2},v)+a(u_{H,1}^n+u_{H,2}^n,v)\\ \hfill =(f,v),v\in V_{H,2}.\end{array}$$

(15)

Lemma 2.

The partially explicit scheme (6) is stable if

$$\alpha \parallel v_2{\parallel }^2\ge \frac{{(\Delta t)}^2}{2}{\parallel v_2\parallel }_a^2\forall v_2\in V_{H,2}.$$

(16)

Proof. 

We have

$$\begin{array}{c}\hfill E^{n+1}=\alpha {\parallel \frac{u_H^{n+1}-u_H^n}{\Delta t}\parallel }^2+\frac{1}{2}\left(\parallel u_{H,1}^{n+1}+u_{H,2}^n{\parallel }_a^2+{\parallel u_{H,2}^{n+1}+u_{H,1}^n\parallel }_a^2\right)\\ \hfill -\frac{{(\Delta t)}^2}{2}\parallel \frac{u_{H,2}^{n+1}-u_{H,2}^n}{\Delta t}{\parallel }_a^2\ge \alpha {\parallel \frac{u_{H,1}^{n+1}-u_{H,1}^n}{\Delta t}\parallel }^2+\frac{1}{2}\left(\parallel u_{H,1}^{n+1}+u_{H,2}^n{\parallel }_a^2+{\parallel u_{H,2}^{n+1}+u_{H,1}^n\parallel }_a^2\right).\end{array}$$

(17)

The latter defines a norm and in this norm, we have a stability. □

3.3. Case $V_{H,1}$ and $V_{H,2}$ Are Non-Orthogonal

In this case, we note that (see (6)) a coupling in mass matrix term remains. This can be removed through some mass lumping (see [4]). We define $\gamma ,{\gamma }_a<1$ and $\alpha \in \mathbb{R}$ as

$$\gamma :=\underset{v_1\in V_{H,1},v_2\in V_{H,2}}{sup}\frac{(v_1,v_2)}{\parallel v_1{\parallel }_{}{\parallel v_2\parallel }_{}},{\gamma }_a:=\underset{v_1\in V_{H,1},v_2\in V_{H,2}}{sup}\frac{a(v_1,v_2)}{\parallel v_1{\parallel }_a{\parallel v_2\parallel }_a}$$

and

$$\beta =\underset{v_2\in V_{H,2}}{sup}\frac{\parallel v_2{\parallel }_a}{\parallel v_2{\parallel }_{}}.$$

Lemma 3.

If

$$2\alpha (1-{\gamma }^2){\beta }^{-2}\ge {(\Delta t)}^2,$$

(18)

we have

$$\begin{array}{c}\hfill (\alpha (1-{\gamma }^2)-\frac{{\beta }^2{(\Delta t)}^2}{2}){\parallel \frac{u_{H,2}^{n+1}-u_{H,2}^n}{\Delta t}\parallel }_{}^2+\frac{(1-{\gamma }_a)}{2}\sum _{i=1,2}\left(\parallel u_{H,i}^{n+1}{\parallel }_a^2+{\parallel u_{H,i}^n\parallel }_a^2\right)\le E^1\end{array}$$

(19)

Proof. 

Since we have

$$E^1\ge E^{n+1}\mathrm{for}\mathrm{any}n\ge 0,$$

we have

$$\begin{array}{cc}\hfill E^{n+1}& =\alpha \parallel \frac{u_H^{n+1}-u_H^n}{\Delta t}{\parallel }^2+\frac{1}{2}\left(\parallel u_{H,1}^{n+1}+u_{H,2}^n{\parallel }_a^2+{\parallel u_{H,2}^{n+1}+u_{H,1}^n\parallel }_a^2\right)-\frac{{(\Delta t)}^2}{2}{\parallel \frac{u_{H,2}^{n+1}-u_{H,2}^n}{\Delta t}\parallel }_a^2\hfill \\ \hfill & \le E^1.\hfill \end{array}$$

We have

$$\begin{array}{cc}\hfill \parallel u_H^{n+1}-u_H^n{\parallel }^2& =\sum _{i=1,2}{\parallel u_{H,i}^{n+1}-u_{H,i}^n\parallel }^2+2(u_{H,1}^{n+1}-u_{H,1}^n,u_{H,2}^{n+1}-u_{H,2}^n)\hfill \\ \hfill & \ge \sum _{i=1,2}\parallel u_{H,i}^{n+1}-u_{H,i}^n{\parallel }^2-2\gamma \parallel u_{H,1}^{n+1}-u_{H,1}^n\parallel \parallel u_{H,2}^{n+1}-u_{H,2}^n\parallel \hfill \\ \hfill & \ge (1-{\gamma }^2){\parallel u_{H,2}^{n+1}-u_{H,2}^n\parallel }^2\hfill \end{array}$$

If $2\alpha (1-{\gamma }^2){\beta }^{-2}\ge {(\Delta t)}^2$, we have

$$\begin{array}{cc}\hfill & \alpha (1-{\gamma }^2)\parallel u_{H,2}^{n+1}-u_{H,2}^n{\parallel }^2-\frac{{(\Delta t)}^2}{2}{\parallel u_{H,2}^{n+1}-u_{H,2}^n\parallel }_a^2\hfill \\ \hfill \ge & \left(\alpha (1-{\gamma }^2)-\frac{{\beta }^2{(\Delta t)}^2}{2}\right){\parallel u_{H,2}^{n+1}-u_{H,2}^n\parallel }^2.\hfill \end{array}$$

Thus, we have

$$\begin{array}{cc}\hfill E^1& \ge \alpha \parallel \frac{u_H^{n+1}-u_H^n}{\Delta t}{\parallel }^2+\frac{1}{2}\left(\parallel u_{H,1}^{n+1}+u_{H,2}^n{\parallel }_a^2+{\parallel u_{H,2}^{n+1}+u_{H,1}^n\parallel }_a^2\right)-\frac{{(\Delta t)}^2}{2}{\parallel \frac{u_{H,2}^{n+1}-u_{H,2}^n}{\Delta t}\parallel }_a^2\hfill \\ \hfill & \ge \left(\alpha (1-{\gamma }^2)-\frac{{\beta }^2{(\Delta t)}^2}{2}\right){\parallel \frac{u_{H,2}^{n+1}-u_{H,2}^n}{\Delta t}\parallel }^2+\frac{1}{2}\left(\parallel u_{H,1}^{n+1}+u_{H,2}^n{\parallel }_a^2+{\parallel u_{H,2}^{n+1}+u_{H,1}^n\parallel }_a^2\right).\hfill \end{array}$$

We also obtain that

$$\begin{array}{cc}\hfill & \frac{1}{2}\left(\parallel u_{H,1}^{n+1}+u_{H,2}^n{\parallel }_a^2+{\parallel u_{H,2}^{n+1}+u_{H,1}^n\parallel }_a^2\right)\hfill \\ \hfill =& \frac{1}{2}\sum _{i=1,2}\left(\parallel u_{H,i}^{n+1}{\parallel }_a^2+{\parallel u_{H,i}^n\parallel }_a^2\right)+a(u_{H,2}^n,u_{H,1}^{n+1})+a(u_{H,1}^n,u_{H,2}^{n+1})\hfill \\ \hfill \ge & \frac{1}{2}\sum _{i=1,2}\left(\parallel u_{H,i}^{n+1}{\parallel }_a^2+{\parallel u_{H,i}^n\parallel }_a^2\right)-{\gamma }_a\left(\parallel u_{H,2}^n{\parallel }_a\parallel u_{H,1}^{n+1}{\parallel }_a+\parallel u_{H,1}^n{\parallel }_a{\parallel u_{H,2}^{n+1}\parallel }_a\right)\hfill \end{array}$$

and

$${\gamma }_a\left(\parallel u_{H,2}^n{\parallel }_a\parallel u_{H,1}^{n+1}{\parallel }_a+\parallel u_{H,1}^n{\parallel }_a{\parallel u_{H,2}^{n+1}\parallel }_a\right)\le \frac{{\gamma }_a}{2}\sum _{i=1,2}\left(\parallel u_{H,i}^{n+1}{\parallel }_a^2+{\parallel u_{H,i}^n\parallel }_a^2\right).$$

Therefore, we have

$$\frac{1}{2}\left(\parallel u_{H,1}^{n+1}+u_{H,2}^n{\parallel }_a^2+{\parallel u_{H,2}^{n+1}+u_{H,1}^n\parallel }_a^2\right)\ge \frac{1-{\gamma }_a}{2}\sum _{i=1,2}\left(\parallel u_{H,i}^{n+1}{\parallel }_a^2+{\parallel u_{H,i}^n\parallel }_a^2\right)$$

and obtain

$$\begin{array}{cc}\hfill E^1& \ge \left(\alpha (1-{\gamma }^2)-\frac{{\beta }^2{(\Delta t)}^2}{2}\right)\parallel \frac{u_{H,2}^{n+1}-u_{H,2}^n}{\Delta t}\parallel +\frac{1}{2}\left(\parallel u_{H,1}^{n+1}+u_{H,2}^n{\parallel }_a^2+{\parallel u_{H,2}^{n+1}+u_{H,1}^n\parallel }_a^2\right)\hfill \\ \hfill & \ge \left(\alpha (1-{\gamma }^2)-\frac{{\beta }^2{(\Delta t)}^2}{2}\right){\parallel \frac{u_{H,2}^{n+1}-u_{H,2}^n}{\Delta t}\parallel }^2+\frac{1-{\gamma }_a}{2}\sum _{i=1,2}\left(\parallel u_{H,i}^{n+1}{\parallel }_a^2+{\parallel u_{H,i}^n\parallel }_a^2\right).\hfill \end{array}$$

This completes the proof. □

3.4. Remarks

We make several remarks. First, we note that the above studies and analysis apply to the case when $f\left(x\right)\ne 0$. Secondly, we note that our studies and analysis also apply to the case—$\alpha$ is a spatial functions, $\alpha =\alpha \left(x\right)$. In this case, one needs to assume that $\alpha \left(x\right)\ge {\alpha }_0>0$ and defines the norm

$$E^{n+1}=\parallel \sqrt{\alpha }\frac{u_H^{n+1}-u_H^n}{\Delta t}{\parallel }^2+\frac{1}{2}\left(\parallel u_{H,1}^{n+1}+u_{H,2}^n{\parallel }_a^2+{\parallel u_{H,2}^{n+1}+u_{H,1}^n\parallel }_a^2\right)-\frac{{(\Delta t)}^2}{2}{\parallel \frac{u_{H,2}^{n+1}-u_{H,2}^n}{\Delta t}\parallel }_a^2.$$

Finally, we would like to remark that the value of $\gamma$ is difficult to determine in practical examples. Our numerical results show the condition for the time step is close to that in standard discretizations (though scales as the coarse mesh size) and we choose it two times larger.

We would like to remark that the computational advantage of our approach is due to the fact that a small system is solved implicitly. Given that the implicit solution procedure is more expensive, the computational gain is in solving a smaller system that corresponds to the degrees of freedom in $V_{H,1}$. The proposed method brings additional computational gains for nonlinear problems. We also would like to point out that the proposed splitting requires the inequality (16) for $V_{H,2}$. In this regard, some other multiscale space constructions (that are not designed for high contrast problems) may not suitable.

The analysis of our approach can be carried out. Here, we give a brief overview of the analysis. We can define $P:V\to V_H$ as the a-projection operator such that

$$a(Pv,w)=a(v,w)\forall w\in V_H.$$

Considering u to be the solution of the QGD equation, we have

$$2({\partial }_{t}u\left(t_n\right),v)+\alpha ({\partial }_{tt}u\left(t_n\right),v)+a(u\left(t_n\right),v)=(f\left(t_n\right),v)\forall v\in V.$$

We define $u^n=u\left(t_n\right)$ and obtain

$$(\frac{u^{n+1}-u^{n-1}}{\Delta t},v)+\alpha (\frac{u^{n+1}-2u^n+u^{n-1}}{\Delta t^2},v)+a(u^n,v)=(f^n,v)+(r^n,v)\forall v\in V_H,$$

where

$$r^n=\left(\frac{u^{n+1}-u^{n-1}}{\Delta t}-2{\partial }_{t}u\right)+\alpha \left(\frac{u^{n+1}-2u^n+u^{n-1}}{\Delta t^2}-{\partial }_{tt}u\right).$$

Assuming u is smooth enough in time, we have

$$\parallel r^n\parallel \le C\Delta t^2\left(\parallel u_{ttt}\parallel +\alpha \parallel u_{tttt}\parallel \right).$$

These arguments can be made rigorous.

3.5. $V_{H,1}$ and $V_{H,2}$ Constructions

In this section, we introduce a possible way to construct the spaces satisfying (18). Here, we follow our previous work [3]. We will show that the constrained energy minimization finite element space is a good choice of $V_{H,1}$ since the CEM basis functions are constructed such that they are almost orthogonal to a space $\tilde{V}$ which can be easily defined. To obtain a $V_{H,2}$ satisfying the condition (18), one of the possible way is using an eigenvalue problem to construct the local basis function. Before, discussing the construction of $V_{H,2}$, we will first introduce the CEM finite element space. In the following, we let $V\left(S\right)=H_0^1\left(S\right)$ for a proper subset $S\subset \Omega$.

3.5.1. CEM Method

In this section, we will discuss the CEM method for solving the problem (2). We will construct the finite element space by solving a constrained energy minimization problem. We let ${\mathcal{T}}_H$ be a coarse grid partition of $\Omega$. For each element $K_i\in {\mathcal{T}}_H$, we consider a set of auxiliary basis functions ${\left\{{\psi }_j^{\left(i\right)}\right\}}_{j=1}^{L_i}\in V\left(K_j\right)$. We then can define a projection operator ${\Pi }_{K_i}:L^2\left(K_i\right)\mapsto V_{aux}^{\left(i\right)}\subset L^2\left(K_i\right)$ such that

$$s_i({\Pi }_{i}u,v)=s_i(u,v)\forall v\in V_{aux}^{\left(i\right)}:=\mathrm{span}\{{\psi }_j^{\left(i\right)}:1\le j\le L_i\},$$

where

$$s_i(u,v)={\int }_{K_i}\tilde{\kappa }uv.$$

(20)

and $\tilde{\kappa }=\kappa H^{-2}$ or $\tilde{\kappa }=\kappa {\sum }_i{|\nabla {\chi }_i|}^2$ with some partition of unity ${\chi }_i$.

We next define a projection operator by $\Pi :L^2(\Omega )\mapsto V_{aux}\subset L^2(\Omega )$

$$s(\Pi u,v)=s(u,v)\forall v\in V_{aux}:=\sum _{i=1}^{N_e}{V}^{\left(i\right)},$$

where $s(u,v):={\sum }_{i=1}^{N_e}{s}_i{\left(u\right|}_{K_i},v{|}_{K_i})$. For each auxiliary basis functions ${\psi }_j^{\left(i\right)}$, we can define a local basis function ${\varphi }_j^{\left(i\right)}\in V\left(K_i^{+}\right)$ such that

$$\begin{array}{cc}\hfill a({\varphi }_j^{\left(i\right)},v)+s({\mu }_j^{\left(i\right)},v)& =0\forall v\in V\left(K_i^{+}\right)\hfill \\ \hfill s({\varphi }_j^{\left(i\right)},\nu )& =s({\psi }_j^{\left(i\right)},\nu )\forall \nu \in V_{aux}\left(K_i^{+}\right)\hfill \end{array}$$

where $K_i^{+}$ is an oversampling domain of $K_i$, which is a few coarse blocks larger than $K_i$ [5]. We then define the space $V_{cem}$ as

$$\begin{array}{cc}\hfill V_{cem}& :=\mathrm{span}\{{\varphi }_j^{\left(i\right)}:1\le i\le N_e,1\le j\le L_i\},\hfill \end{array}$$

where $N_e$ is the number of coarse elements. The CEM solution $u_{cem}$ is given by

$$\begin{array}{cc}\hfill (\frac{\partial }{\partial t}\left(u_{cem}\right),v)+\alpha (\frac{{\partial }^2}{\partial t^2}\left(u_{cem}\right),v)& =-a(u_{cem},v)\forall v\in V_{cem}.\hfill \end{array}$$

We remark that the $V_{glo}$ is $a-$orthogonal to a space $\tilde{V}:=\{v\in V:\Pi \left(v\right)=0\}$. We also know that $V_{cem}$ is closed to $V_{glo}$ and therefore it is almost orthogonal to $\tilde{V}$. Thus, we can choice $V_{cem}$ to be $V_{H,1}$ and construct a space $V_{H,2}$ in $\tilde{V}$.

3.5.2. Construction of $V_{H,2}$

We discuss a choice for the space $V_{H,2}\subset \tilde{V}$. We will investigate the stability properties numerically in Section 4.

The choice of $V_{H,2}$ is based on the CEM type finite element space. For each coarse element $K_i$, an eigenvalue problem is solved to obtain the auxiliary basis. We find eigenpairs $({\eta }_j^{\left(i\right)},{\gamma }_j^{\left(i\right)})\in (V\left(K_i\right)\cap \tilde{V})\times \mathbb{R}$ via

$$\begin{array}{cc}\hfill {\int }_{K_i}\kappa \nabla {\eta }_j^{\left(i\right)}·\nabla vs.& ={\gamma }_j^{\left(i\right)}{\int }_{K_i}{\eta }_j^{\left(i\right)}v,\forall v\in V\left(K_i\right)\cap \tilde{V}.\hfill \end{array}$$

(21)

For each $K_i$, we choose the first few $J_i$ eigenfunctions corresponding to the smallest $J_i$ eigenvalues. The span of these functions form a space which is called $V_{aux,2}$ and $V_{aux,1}$ is defined as the auxiliary space of the CEM space, namely,

$$\begin{array}{cc}\hfill V_{aux,1}& :=\mathrm{span}\left\{{\psi }_j^{\left(i\right)}\right|\forall K_i,\forall 1\le j\le L_i\},\hfill \\ \hfill V_{aux,2}& :=\mathrm{span}\left\{{\eta }_j^{\left(i\right)}\right|\forall K_i,\forall 1\le j\le J_i\}.\hfill \end{array}$$

For each auxiliary basis function ${\eta }_j^{\left(i\right)}\in V_{aux,2}$, we define a basis function ${\zeta }_j^{\left(i\right)}\in V\left(K_i^{+}\right)$ such that ${\mu }_j^{\left(i\right)}\in V_{aux,1}$, ${\mu }_j^{\left(i\right),2}\in V_{aux,2}$ and

$$\begin{array}{cc}\hfill a({\zeta }_j^{\left(i\right)},v)+s({\mu }_j^{\left(i\right),1},v)+({\mu }_j^{\left(i\right),2},v)& =0,\forall v\in V\left(K_i^{+}\right),\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill s({\zeta }_j^{\left(i\right)},\nu )& =0,\forall \nu \in V_{aux,1},\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill ({\zeta }_j^{\left(i\right)},\nu )& =({\xi }_j^{\left(i\right)},\nu ),\forall \nu \in V_{aux,2}.\hfill \end{array}$$

(24)

where we use the notation $V_{aux,1}$ to denote the space $V_{aux}$ defined in Section 3.5.1, and $K_i^{+}$ is an oversampling domain a few coarse blocks larger than $K_i$ (see [5]). We define

$$V_{H,2}=\mathrm{span}\left\{{\zeta }_j^{\left(i\right)}\right|\forall K_i,\forall 1\le j\le J_i\}.$$

4. Numerical Result

We present representative numerical results. We consider the following mesh and time step parameters in all examples.

$$H=1/10,h=1/100.$$

We consider the medium parameter $\kappa \left(x\right)$ depicted in Figure other (see Figure 1). Medium properties are high contrast and multiscale. We choose the source term as a source distrubuted in a small region as shown in Figure 1). It, f, is given by

$$f(t,x)=\frac{2-2/f_0}{4h^2}exp(-{\pi }^2{f}_0^2{(t-2/f_0)}^2)f_x\left(x\right),$$

(25)

where $f_0$ is a frequency, which will be taken to be $f_0=0.5$. We take the values for $\alpha$ to be $\alpha =5$, $\alpha =0.5$, $\alpha =0.05$, and $\alpha =0.005$, which balances wave terms and parabolic terms. For larger $\alpha$’s, we will observe slightly larger errors since the wave phenomena dominates. In all examples, we plot errors in time for three methods.

• Implicit CEM. In this, we use only implicit method with fewer degrees of freedom (without additional degrees of freedom) and compute the error associated with the multiscale approach.
• Implicit CEM with additional basis functions. In this approach, we take into account additional degrees of freedom and handle them in implicit manner. The method is more expensive and corresponds to full implicit approach.
• Partial Explicit Splitting CEM. In our approach, we take into account additional degrees of freedom and handle them in explicit manner. One approach for additional degrees is presented in Section 3.5.2. Partial Explicit Splitting CEM and Implicit CEM have the same degrees of freedom, while partial Explicit Splitting CEM handles additional degrees of freedom explicitly.

In Figure 2, we depict the results for the case $\alpha =5$ and time step $dt=0.006$. As we see from the results that additional degrees of freedom improves the results. We note that in all examples, we choose the source term such that the error without additional degrees of freedom is large and we can improve it. In general, this error can be improved and made smaller with additional degrees of freedom. This is partially because of the source term which brings additional multiscale nature. Moreover, we observe that the errors of proposed partial explicit approach and fully implicit approach are similar. This suggest that our proposed approach is stable and can be used for explicit time stepping. In Figure 3, we use smaller $\alpha$, $\alpha =0.5$. In this case, we again observe similar phenomena, where our proposed approach provides very similar results when additional degrees of freedom treated implicitly. In other cases, we consider other values of $\alpha$, $\alpha =0.05$ (see Figure 4) and $\alpha =0.005$ (see Figure 5). In all cases, the proposed partial explicit approach provides very similar results compared to fully implicit approach when treating additional degrees of freedom. We note that because of nonzero right hand side, the energy can increase.

We will next present a numerical example with the same medium as the previous examples. We consider the source term is smoother spacially and with the same frequency as the previous example. In this case, we will have a much smaller error in both $L_2$ and Energy norm since the $V_H$ space can approximte the source better. In Figure 6, we show the medium with marker indicating three locations for comparison and the source term.

In Figure 7, we compare the partially explicit solution with different grids in $L_2$ norm and energy norm. We see that as we decrease the mesh size, the error decreases. In Figure 8, we compare the partially explicit solution with at three different locations.

5. Conclusions

In this paper, we design and analyze contrast-independent partially explicit time discretization methods for quasi gas dynamics model. The proposed methods differ from our previous works [3,4] as parabolic and wave equations require different treatment for temporal and spatial discretization. Temporal splitting is based on spatial multiscale space splitting that we introduce in the paper. Using these multiscale spaces, time splitting identifies and treats fast components implicitly while slow component explicitly. Our proposed method is still implicit via mass matrix, which can be treated via mass lumping as in [4]. A stability of the proposed splitting is demonstrated under some conditions. The proposed method allows identifying the degrees of freedom that need explicit treatment and those that require implicit treatment. Numerical results show that the proposed methods provide very similar results as fully implicit methods using explicit methods with the time stepping that is independent of the contrast.