Adaptive Discontinuous Galerkin Finite Elements for Advective Allen-Cahn Equation

TL;DR

This paper develops an adaptive interior penalty discontinuous Galerkin method to solve the advective Allen-Cahn equation with complex velocity fields, improving accuracy in capturing sharp layers.

Contribution

It introduces a residual-based a posteriori error estimator tailored for non-divergence-free velocity fields in the adaptive DG framework.

Findings

Effective resolution of sharp layers in numerical examples

High accuracy achieved with adaptive method

Error estimator accounts for non-divergence-free velocities

Abstract

We apply a space adaptive interior penalty discontinuous Galerkin method for solving advective Allen-Cahn equation with expanding and contracting velocity fields. The advective Allen-Cahn equation is first discretized in time and the resulting semi-linear elliptic PDE is solved by an adaptive algorithm using a residual-based a posteriori error estimator. The a posteriori error estimator contains additional terms due to the non-divergence-free velocity field. Numerical examples demonstrate the effectiveness and accuracy of the adaptive approach by resolving the sharp layers accurately.

Figures

Equations

$$\frac{\partial u}{\partial t}+\nabla \cdot (\mathbf{V}u)$$

$$u(\mathbf{x},0)$$

$$\frac{\partial u}{\partial n}$$

$$\begin{array}{rl}\frac{u^k(\mathbf{x})-u^{k-1}(\mathbf{x})}{\tau }& -ϵ\Delta u^k(\mathbf{x})+\mathbf{V}\cdot \nabla u^k(\mathbf{x})\\ & +(\nabla \cdot \mathbf{V})u^k(\mathbf{x})+\frac{1}{ϵ}f(u^k(\mathbf{x}))=0.\end{array}$$

$$\alpha u^k-ϵ\Delta u^k+\mathbf{V}\cdot \nabla u^k+r(u^k)=l(u^{k-1}),$$

$$\alpha :=\alpha (\mathbf{x})$$

$$r(u^k)$$

$$l(u^{k-1})$$

$$|r(s)|$$

$$\Vert r(s_1)-r(s_2){\Vert }_{L^2(\Omega )}$$

$$\alpha (\mathbf{x})-\frac{1}{2}\nabla \cdot \mathbf{V}(\mathbf{x})$$

$$\Vert -\nabla \cdot \mathbf{V}(\mathbf{x})+\alpha (\mathbf{x}){\Vert }_{L^{\infty }(\Omega )}$$

$$\frac{1}{\tau }+\frac{1}{2}\nabla \cdot \mathbf{V}(\mathbf{x})$$

$$\Vert {\tau }^{-1}{\Vert }_{L^{\infty }(\Omega )}$$

$$\underset{E\in {\mathcal{T}}_h^k}{\max }\frac{h_E^2}{\left|E\right|}\le c_0,$$

$$V_h^k=\left\{\upsilon \in L^2(\Omega ):\upsilon {|}_E\in {\mathbb{P}}^q(E),\forall E\in {\mathcal{T}}_h^k\right\},$$

$$[\upsilon {]}_e={\upsilon }_{{|}_{E_i}}{\mathbf{n}}_e-{\upsilon }_{{|}_{E_j}}{\mathbf{n}}_e,\{\upsilon {\}}_e=\frac{1}{2}({\upsilon }_{{|}_{E_i}}+{\upsilon }_{{|}_{E_j}}),$$

$$\begin{array}{rl}[\nabla \upsilon {]}_e& =\nabla {\upsilon }_{{|}_{E_i}}\cdot {\mathbf{n}}_e-\nabla {\upsilon }_{{|}_{E_j}}\cdot {\mathbf{n}}_e,\\ \{\nabla \upsilon {\}}_e& =\frac{1}{2}(\nabla {\upsilon }_{{|}_{E_i}}+\nabla {\upsilon }_{{|}_{E_j}}).\end{array}$$

$$\begin{array}{rl}[\upsilon {]}_e={\upsilon }_{{|}_{E_i}}\mathbf{n},& \{\upsilon {\}}_e={\upsilon }_{{|}_{E_i}},\\ [\nabla \upsilon {]}_e=\nabla {\upsilon }_{{|}_{E_i}}\cdot \mathbf{n},& \{\nabla \upsilon {\}}_e=\nabla {\upsilon }_{{|}_{E_i}},\end{array}$$

$${\Gamma }_h^{-}=\left\{\mathbf{x}\in \partial \Omega :\mathbf{V}(\mathbf{x})\cdot \mathbf{n}(\mathbf{x})<0\right\},{\Gamma }_h^{+}=\partial \Omega \setminus {\Gamma }_h^{-}.$$

$$\partial E_h^{-}=\left\{\mathbf{x}\in \partial E:\mathbf{V}(\mathbf{x})\cdot {\mathbf{n}}_e(\mathbf{x})<0\right\},\partial E_h^{+}=\partial E\setminus \partial E_h^{-},$$

$$a_h(u_h^k,{\upsilon }_h^k)+r_h(u_h^k)=l_h^{k-1}({\upsilon }_h^k),$$

$$a_h(u,\upsilon ):=D_h(u,\upsilon )+O_h(u,\upsilon )+K_h(u,\upsilon )+J_h(u,\upsilon ),$$

$$\begin{array}{rl}{D}_h(u,\upsilon )=& \sum _{E\in {\mathcal{T}}_h^k}{\int }_E\left(ϵ\nabla u\cdot \nabla \upsilon +(\alpha -\nabla \cdot \mathbf{V})u\upsilon \right)d\mathbf{x},\\ O_h(u,\upsilon )=& -\sum _{E\in {\mathcal{T}}_h^k}{\int }_E\mathbf{V}u\cdot \nabla \upsilon d\mathbf{x}+\sum _{E\in {\mathcal{T}}_h^k}{\int }_{\partial E_h^{+}\cap {\Gamma }_h^{+}}\mathbf{V}\cdot {\mathbf{n}}_{E}u\upsilon ds\\ & +\sum _{E\in {\mathcal{T}}_h^k}{\int }_{\partial E_h^{+}\setminus \partial \Omega }\mathbf{V}\cdot {\mathbf{n}}_{E}u(\upsilon -{\upsilon }^{out})ds,\\ K_h(u,\upsilon )=& -\sum _{e\in {\mathcal{F}}_{h,0}^k}{\int }_e\left(\{ϵ\nabla u\}[\upsilon ]+\{ϵ\nabla \upsilon \}[u]\right)ds,\\ J_h(u,\upsilon )=& \sum _{e\in {\mathcal{F}}_{h,0}^k}\frac{\sigma ϵ}{h_e}{\int }_e[u][\upsilon ]ds,\end{array}$$

$$\begin{array}{rl}{r}_h(u)& =\sum _{E\in {\mathcal{T}}_h^k}{\int }_{E}r(u)\upsilon d\mathbf{x},\\ l_h^{k-1}(\upsilon )& =\sum _{E\in {\mathcal{T}}_h^k}{\int }_E\frac{1}{\tau }{u}_h^{k-1}\upsilon d\mathbf{x}.\end{array}$$

$$\Vert \upsilon {\Vert }_{dG}^2:=|||\upsilon ||{|}^2+|\upsilon {|}_C^2,$$

$$|||\upsilon ||{|}^2=$$

$$+\sum _{e\in {\mathcal{F}}_{h,0}^k}\frac{\sigma ϵ}{h_e}\Vert [\upsilon ]{\Vert }_{L^2(e)}^2,$$

$$|\upsilon {|}_C^2=|\mathbf{V}\upsilon {|}_{\ast }^2+\sum _{e\in {\mathcal{F}}_{h,0}^k}({\kappa }_0{h}_e+\frac{h_e}{ϵ})\Vert [\upsilon ]{\Vert }_{L^2(e)}^2,$$

$$|\mathbf{V}\upsilon {|}_{\ast }={\sup }_{w\in H_0^1(\Omega )\setminus \{0\}}\frac{{\int }_{\Omega }\mathbf{V}\upsilon \cdot \nabla wdx}{|||w|||},$$

$$({\eta }_E^k{)}^2$$

$$\begin{array}{rl}({\eta }_{E,R}^k{)}^2=& {\rho }_E^2\Vert l(u_h^{k-1})-{\alpha }_h{u}_h^k+ϵ{\Delta }_h{u}_h^k\\ & -{\mathbf{V}}_h\cdot {\nabla }_h{u}_h^k-r(u_h^k){\Vert }_{L^2(E)}^2,\end{array}$$

$$\begin{array}{rl}({\eta }_{E,0}^k{)}^2=& \frac{1}{2}\sum _{e\in \partial E\cap {\mathcal{F}}_{h,0}^{k-1}}\{{ϵ}^{-\frac{1}{2}}{\rho }_e\Vert [ϵ{\nabla }_h{u}_h^k]{\Vert }_{L^2(e)}^2\\ & +\left(\frac{ϵ\sigma }{h_e}+{\kappa }_0{h}_e+\frac{h_e}{ϵ}\right)\Vert [u_h^k]{\Vert }_{L^2(e)}^2\}.\end{array}$$

$${\rho }_E=\min \{h_E{ϵ}^{-\frac{1}{2}},{\kappa }_0^{-\frac{1}{2}}\},{\rho }_e=\min \{h_e{ϵ}^{-\frac{1}{2}},{\kappa }_0^{-\frac{1}{2}}\},$$

Full text

MSC 2010: 65M20; 65M22; 65M50; 65M60

††journal: Numerical Algebra, Control & Optimization (Accepted)

1 Introduction

Interfacial dynamics has great importance in the modeling of multi-phase flow and it plays an important role in different scientific and industrial applications such as micro-structure evolution and grain growth in material science [1], binary fluids flow movement [2], and complex interfacial dynamics [3].

There have been various diffuse interface models for multi-phase flow [4, 5]. In this study, we consider a specific model of diffuse interface for two phase flow; Allen-Cahn (AC) equation with advection. AC equation without advection is the most known dynamical model for diffuse interface dynamics. In the past it was investigated extensively with many numerical and analytical methods. But few studies are dealing with the numerical solution of the advective AC equation [6, 7, 8, 9]. The AC model with advection is given by

[TABLE]

In (1), the bi-stable cubic nonlinearity $f(u)$ is such that $f(u)=F^{\prime}(u)=2u(1-u)(1-2u)$ for a double-well potential $F(u)$, and $\bm{V}=(V_{1},V_{2})^{T}$ is the prescribed velocity field. The velocity field is either given or computed from the Navier-Stokes equation [4, 5]. In many cases the velocity field $\bm{V}$ is divergence-free [6, 7] and satisfies incompressible Navier-Stokes equations. In this study we consider non-divergence-free velocity fields $\nabla\cdot\bm{V}\neq 0$, either expanding $\nabla\cdot\bm{V}>0$ or contracting $\nabla\cdot\bm{V}<0$ [5]. In other words we consider advective AC equation in compressible fluids.

Problems with surface tension in two-phase fluids are known as multi-scale problems with two different time scales, the small surface tension, and the convection time scale, which results in computational stiffness. Actually, there exists three main algorithms: the sharp interface algorithm method, the level-set algorithm method and the diffuse interface method [5]. The numerical simulations are illustrated using finite elements method in space and semi-implicit schemes or semi-implicit schemes with splitting in time [10, 11].

Numerical solutions of advective AC equation may exhibit unphysical oscillations at the interior layers due to convection and sharp fronts may occur due to the non-linear reaction term. Since the standard FEMs are known to produce strong oscillations around layers, adaptive algorithms are developed to tackle the unphysical oscillations and shocks. By refining the mesh locally at layers and sharp fronts, accurate solutions are obtained with less degrees of freedom (DoFs) and computational time. The major part of an adaptive algorithm is the estimation of the local errors and refine the elements with large estimated errors. A posteriori error estimation is the main tool to estimate the local errors which uses the approximate solution and the given problem data [12, 13, 14]. Since the discontinuous Galerkin (dG) methods have the flexibility on adaptive meshes, a posteriori error estimators are developed using dG discretization [15, 16, 17, 18, 19, 20].

We develop an adaptive strategy for the numerical solution of advective AC equation (1) for contracting and expanding flows. Because the solutions do not show strong variations with respect to time, only space adaptivity is applied. Usually the time dependent PDEs are discretized first in space and integrated in time. Here we discretized first in time by Rothe’s method [21]. Then, the resulting semi-linear elliptic equations are solved with the adaptive version of symmetric interior penalty Galerkin (SIPG) method using upwinding for the convective term. The adaptive strategy is based on a residual based a posteriori error estimation [15, 22]. In [22] we developed a posteriori error bounds with respect to the energy norm induced by the SIPG formulation of the system given for semi-linear diffusion-convection-reaction equations with divergence-free velocity field. In this study we develop a posteriori error estimates for the advective AC equation with non-divergence-free convective terms. The extra terms in the error estimates coming from the non-divergence-free vector field are added to the reaction terms. Numerical tests show that the proposed adaptive algorithm can efficiently detect the layers and accurate solutions around these layers can be obtained.

The paper is organized as follows. We first derive the space/time discretization of the advective AC equation (1) in Section 2. A detailed explanation for space adaptivity algorithm is presented in Section 3, where we derive a posteriori error bounds for a stationary problem. Finally, we present two numerical examples with expanding and contracting flows in Section 4 to demonstrate the performance of the adaptive approach.

2 Discretization (Rothe’s method)

Since the solutions of advective AC equation (1) do not change much with the evolution of time, we apply only space adaptivity ($h$-adaptivity), which utilizes a posteriori error estimation or stationary (elliptic) problems at each time step. We first discretize (1) in time for obtaining by Rothe’s method [21]. Then, on each time interval, space discretization and adaptivity is applied to the elliptic stationary problem.

2.1 Time discretization

For the semi-discrete scheme, we consider the uniform partition $0=t_{0}<t_{1}<\ldots<t_{J}=T$ of the time interval $[0,T]$ with the uniform time step-size $\tau=T/J$, and with the time intervals $I_{k}=(t_{k-1},t_{k}]$. For $k=1,\ldots,J$, we denote by $u^{k}(\bm{x})$ the approximate solution at the time instance $t=t_{k}$, i.e. $u^{k}(\bm{x})\approx u(\bm{x},t_{k})$, and we set $u^{0}(\bm{x})=g(\bm{x})$. We use backward Euler as the time integrator to discretize the advective AC (1), for $k=1,\ldots,J$, given $u^{k-1}(\bm{x})$ find $u^{k}(\bm{x})$ satisfying

[TABLE]

For each $k=1,2,\ldots,J$, the equation (2) can be written in the form of a semi-linear elliptic problem as

[TABLE]

where

[TABLE]

We assume that the non-linear reaction term is bounded and locally Lipschitz continuous, i.e., satisfy for any $s,s_{1},s_{2}\geq 0$, $s,s_{1},s_{2}\in\mathbb{R}$ the following conditions

[TABLE]

Moreover, we assume that there is a non-negative constant $\kappa_{0}$ satisfying

[TABLE]

for a positive constant $c^{*}$. The identity (6a) guarantees the coercivity of the bilinear form $a_{h}(\cdot,\cdot)$ in (8), and (6b) is needed to prove the reliability of a posteriori error estimator in Section 3. According to our setting for $\alpha$ in (4), the assumptions in (6) are equivalent that

[TABLE]

2.2 Space discretization

For the space discretization, we use symmetric interior penalty Galerkin (SIPG) method [23, 24], which is a member of the family of discontinuous Galerkin (dG) finite elements methods, and we apply upwinding [25, 26] for the convective term. The dG methods exhibit attractive properties of both classical finite elements method (FEM) and finite volume method (FVM). The functions in dG spaces are discontinuous along the inter-element boundaries, which makes dG methods flexible. In addition, different order of basis functions on each element can be used within dG schemes ($p$-adaptivity). Hence, it allows to use hp-adaptivity methods [27] which arranges the mesh elements and also the order of polynomials on each element adaptively. Further, the dG methods locally conserve several physical quantities such as mass and energy, which plays an important role in the flow and transport problems. Moreover, the sharp gradients or the singularities in the mesh can be locally detected by the fully discontinuous polynomial representation of the solution in dG schemes. A sample archive of MATLAB codes together with a dGFEM tutorial can be found at https://github.com/muzunca/dGFEM-Tutorial.git.

On the $k$-th time interval, let $\{\mathcal{T}_{h}^{k}\}_{h}$ be a family of shape regular meshes with triangular elements $\{E_{i}\}$ such that $\overline{\Omega}=\cup\overline{E}_{i}$ and $E_{i}\cap E_{j}=\emptyset$ for $E_{i}$, $E_{j}$ $\in\mathcal{T}_{h}^{k}$. By shape regularity, we mean that there exists a constant $c_{0}$ such that

[TABLE]

where $h_{E}$ is the diameter and ${\left|E\right|}$ is the area of the triangular element $E$. We split the set of all edges $\mathcal{F}_{h}^{k}$ into the set of interior edges $\mathcal{F}^{k}_{h,0}$ and the set of boundary edges $\mathcal{F}^{k}_{h,\partial}$ so that $\mathcal{F}_{h}^{k}=\mathcal{F}^{k}_{h,0}\cup\mathcal{F}^{k}_{h,\partial}$. We set the finite dimensional solution and test function space by

[TABLE]

where $\mathbb{P}^{q}(E)$ denotes the set of all polynomials of degree at most $q$ on the element $E\in\mathcal{T}_{h}^{k}$. Note that the functions in $V_{h}^{k}\not\subset H^{1}_{0}(\Omega)$ are allowed to be discontinuous along the inter-element boundaries, thus, in contrast to continuous FEM, the dG methods are suitable to use a non-conforming space.

The discontinuities of the functions in $V_{h}^{k}$ along the inter-element boundaries lead to different traces from the neighboring elements sharing an edge $e$. Let the edge $e$ be a common edge for two elements $E_{i}$ and $E_{j}$ (w.l.o.g. assume $i<j$). Then for a scalar function $\upsilon\in V_{h}^{k}$, there are two traces of $\upsilon$ along $e$, denoted by $\upsilon_{|_{E_{i}}}$ from inside $E_{i}$ and $\upsilon_{|_{E_{j}}}$ from inside $E_{j}$. Then, the jump and average of $\upsilon$ across the edge $e$ are defined as

[TABLE]

where $\bm{n}_{e}$ is the unit normal to the edge $e$ oriented exterior to $E_{i}$. Similarly, we set the jump and average values of a vector field $\nabla\upsilon$ on e as

[TABLE]

For a boundary edge $e\subset E_{i}\cap\partial\Omega$, we set

[TABLE]

where $\bm{n}$ is the unit outward normal to the boundary at $e$. We further define the sets of inflow and outflow boundary parts, respectively, by

[TABLE]

The set of inflow and outflow boundary edge parts for an element $E$ is defined in a similar way by

[TABLE]

where $\bm{n}_{e}(\bm{x})$ is the unit outward normal vector to the element boundary $\partial E$ at $\bm{x}$. Moreover, for an interior edge $e$, we denote the trace of a scalar function $\upsilon$ from inside the element $E$ by $\upsilon^{in}$ and from outside the element $E$ by $\upsilon^{out}$.

We multiply the continuous (in space) equation (3) by a test function $\upsilon_{h}^{k}\in V_{h}^{k}$, integrate over $\Omega$, apply Green’s identity for diffusive term together with the upwinding for the convective term, and we obtain the following weak problem: set $u_{h}^{0}$ be the $L^{2}$-projection of $u^{0}(\bm{x})$ onto the dG space $V_{h}^{k}$, for $k=1,2,\ldots,J$, given $u_{h}^{k-1}$, $\forall\upsilon_{h}^{k}\in V_{h}^{k}$ find $u_{h}^{k}\in V_{h}^{k}$ satisfying

[TABLE]

where the dG bilinear form is such that

[TABLE]

where the bilinear terms are given by

[TABLE]

with $\sigma$ denoting the penalty parameter which should be sufficiently large for the stability of the dG scheme [24]. The nonlinear form $r_{h}(u)$ and the linear right hand side $l_{h}^{k-1}(\upsilon)$ in (8) are given by

[TABLE]

3 Space adaptivity

In this section, a space-adaptive procedure is constructed for the SIPG discretized system (8) of the stationary problem (3) which mimics the $k$-th time step of the semi-discretized advective AC equation (1). We use a residual-based a posteriori error approach to apply the adaptivity in which not only we employ the refinement but also we consider the coarsening phenomena. In [28], an adaptive scheme using a posteriori error estimates is constructed for stationary linear diffusion-convection-reaction equation using dG. It is extended to stationary diffusion-convection equations with non-linear reaction using dG in [22], and to parabolic diffusion-convection equations with non-linear reaction in [15, 29].

The adaptive algorithm starts with a sufficiently coarse initial mesh $\mathcal{T}_{h}^{0}$ together with an initial vector $u_{h}^{0}$ which is the $L^{2}$-projection of the initial condition $u^{0}(\bm{x}):=g(\bm{x})$ onto the initial solution space $V_{h}^{0}$. Then, on each time interval $I_{k}=(t_{k-1},t_{k}]$, we consider a single stationary problem (3) with its SIPG discretized system (8). Firstly, we solve the discrete system (8) for the solution $u_{h}^{k}$ on the space $V_{h}^{k-1}$ given that $u_{h}^{k-1}\in V_{h}^{k-1}$ is the known solution vector from the previous time interval. Then, it follows the estimation step providing information about the elements for refinement/coarsening. In the estimation procedure, we utilize a posteriori error estimation. An adaptive scheme with the use of an a posteriori error estimation requires two crucial ingredients. One is an error indicator to compute the local errors on each element, and the other is a compatible norm to measure the error. Here, in order to measure the error, we use the following so-called dG norm

[TABLE]

which is composed of the energy-like norm

[TABLE]

and the semi-norm

[TABLE]

where

[TABLE]

and the parameter $\kappa_{0}$ is the lower bound in (7). As the indicator on the mesh $\mathcal{T}_{h}^{k-1}$, we use the following a posteriori error indicator [29, 28, 22]:

[TABLE]

where $\eta_{E,R}^{k}$ stands for the volume residual on the element $E$, given by

[TABLE]

while $\eta_{E,0}^{k}$ denotes the edge residuals coming from the jump of the numerical solution on the interior edges, given by

[TABLE]

In (12) and (13), $\Delta_{h}$ and $\nabla_{h}$ denote the discrete versions of the Laplace and gradient operators, respectively, and $\alpha_{h}$ and $\bm{V}_{h}$ are the $L^{2}$-projections of the data $\alpha$ and $\bm{V}$ onto $V_{h}^{k-1}$, respectively. The positive weights $\rho_{E}$ and $\rho_{e}$ are defined by

[TABLE]

for $\kappa_{0}\neq 0$, or $\rho_{E}=h_{E}\epsilon^{-\frac{1}{2}}$ and $\rho_{e}=h_{e}\epsilon^{-\frac{1}{2}}$ when $\kappa_{0}=0$. Then, the global error indicator is set as

[TABLE]

We also introduce the data approximation terms,

[TABLE]

and the data approximation error

[TABLE]

Then, using the definitions in (14) and (15), we can obtain the a posteriori error bounds

[TABLE]

where $\lesssim$ mimics a bound up to a constant. The proof of the a posteriori error bounds in (16) and (17) proceeds as the following: first the dG solution $u_{h}^{k}$ is rewritten as the direct sum of its conforming part $u_{h,c}^{k}\in H_{0}^{1}(\Omega)\cap V_{h}^{k-1}$ and the non-conforming remainder term $u_{h,r}^{k}\in V_{h}^{k-1}$:

[TABLE]

and from triangle inequality,

[TABLE]

where the term $\|u^{k}-u_{h,c}^{k}\|_{dG}$ is now well-defined. Then, using coercivity of the bilinear form $a_{h}$, continuity of the bilinear forms $D_{h}$, $O_{h}$, $K_{h}$ and $J_{h}$, inf-sup condition, error bounds for approximation and interpolation operator, the assumptions (5) on the nonlinearity, we derive the following bounds

[TABLE]

which finishes the proof of reliability. For further details on the proof of reliability and efficiency of the a posteriori error estimation we refer to [15, 29, 28, 22].

After collecting information about the elements in the estimation step, we mark the elements for refinement/coarsening. We prescribe two tolerances $\bf{stol^{r}}$ and $\bf{stol^{c}}$ related to the refinement and coarsening, respectively. We refine an element $E$ for which the corresponding error indicator $\eta_{E}^{k}$ is greater than $\bf{stol^{r}}$, and we coarsen an element $E$ if it is less than $\bf{stol^{c}}$, provided that $E$ is not an element of the initial mesh. Accordingly, we form the sets $M_{R}$ and $M_{C}$ of the elements in $\mathcal{T}_{h}^{k-1}$ to be refined and coarsened defined, respectively, by

[TABLE]

Then, we create the new mesh $\mathcal{T}_{h}^{k}$ by refining the elements $E\in M_{R}$ using the newest vertex bisection method [30], and by coarsening the elements $E\in M_{C}$. As the final stage, we resolve the discrete system (8) for $u_{h}^{k}\in V_{h}^{k}$ on the new mesh $\mathcal{T}_{h}^{k}$. However, the known solution vector $u_{h}^{k-1}\in V_{h}^{k-1}$ does not belong to the new solution space $V_{h}^{k}$. Before resolving the system on the new mesh, we simply recover (by interpolation or projection) the known solution vector $u_{h}^{k-1}$ to be used on the new solution space $V_{h}^{k}$.

The accuracy of the solutions is through a better choice for the initial mesh $\mathcal{T}_{h}^{0}$ by pre-refining the uniform initial mesh so that the approximation error of the initial condition is well enough. Thus, in order to determine the initial mesh $\mathcal{T}_{h}^{0}$, we form a sequence of initial meshes $\{\mathcal{T}_{h}^{0,k}\}$ for $k=0,1,\ldots$, and we set the initial mesh $\mathcal{T}_{h}^{0}:=\mathcal{T}_{h}^{0,K}$ for some integer $K\geq 1$, where $\mathcal{T}_{h}^{0,0}$ is the initial uniform mesh and each $\mathcal{T}_{h}^{0,k}$ is obtained by refining $\mathcal{T}_{h}^{0,k-1}$. To select the elements to be refined in $\mathcal{T}_{h}^{0,k-1}$, we simply form the set

[TABLE]

where $\bf{stol^{0}}$ is a user defined tolerance, and $\rho_{E}^{k-1}$ is the approximation error given by

[TABLE]

We select the initial mesh when it is satisfied that $\max_{E\in\mathcal{T}_{h}^{0,k}}\leq\bf{stol^{0}}$, or a maximum number $K$ of pre-refinement level is reached. The adaptive procedure is given in Algorithm 1.

4 Numerical results

In this section, we present numerical examples for advective AC equation under homogenous Neumann boundary conditions to demonstrate the effectiveness of the adaptive SIPG method to recapture sharp layers. All simulations are performed on the spatial domain $\Omega=[-1,1]^{2}$ using linear polynomials in the dG space. The final time is $T=0.06$, and the step-size is taken as $\tau=0.001$. In the figures, we indicate the (spatial) dimension of dG space by DoFs which is the product of the number of triangular elements in the mesh and the local dimension on the elements, i.e. the number of dG basis functions defined on each element.

4.1 Expanding flow

We consider the advective AC equation with an expanding velocity field as [5]

[TABLE]

with $v_{0}=10$. The initial condition is taken as the symmetric data

[TABLE]

and the interface length is $\epsilon=0.001$. We first obtain the solutions on a uniform mesh with the mesh size $\Delta x=\Delta y=1/64$. Figure 1, top, shows the uniform solutions at different time instances, where there seem spurious oscillations. For the adaptive case, we prescribe the tolerances $\bf{stol^{r}}=1\times 10^{-2}$ and $\bf{stol^{c}}=1\times 10^{-5}$, and we take the initial uniform mesh with the mesh size $\Delta x=\Delta y=1/4$. The adaptive solutions and corresponding adaptive meshes are given in Figure 1, middle-bottom. We can see that with much less DoFs, the oscillations are almost disappeared, and that the internal layers are well-captured together with the expanding behavior. In Figure 2, we present the propagation of the maximum element error $\eta_{E}^{k}$ over the time, and also the evaluation of DoFs during the time progression. We see that maximum element error is bounded in a small band according to the prescribed tolerances, and both refinement and coarsening phenomena works well.

4.2 Sheer flow

Now, we test the advective AC equation with a sheering flow given by [5]

[TABLE]

with $v_{0}=100$, and the initial condition is the square data

[TABLE]

The interface length is taken as $\epsilon=0.01$. As the expanding case, the uniform solutions are obtained on the uniform mesh with the mesh size $\Delta x=\Delta y=1/32$, and for the adaptive case we use the same settings. In Figure 3, top, uniform solutions are shown with small oscillations, whereas they are dumped out in the adaptive solutions, Figure 3, middle. The layers are well-captured by the related adaptive meshes in Figure 3, bottom, in coherence with sheering behavior of the solutions. Moreover, the maximum element errors again lies in a small band, and the number of DoFs needed for the adaptive scheme is pretty less then it requires for the uniform case, Figure 4.

5 Conclusion

Most of the adaptive methods are developed for linear and nonlinear diffusion-convection-reaction equations with non-divergence-free velocity fields. In this paper we have developed an adaptive procedure for AC equation in compressible fluids using interior penalty discontinuous Galerkin method. Numerical results show that comparing with the uniform meshes, the sharp layers can be resolved accurately using less number of DoFs on adaptive meshes.

Acknowledgments

The authors would like to thank the reviewer for the comments and suggestions that helped to improve the manuscript.