constraint energy minimization, multiscale finite element methods, inhomogeneous boundary value problem, convection-diffusion equation

1 Introduction↩︎

Convection diffusion equation is involved in many physical applications of partial differential equations. Computational difficulty may arise in two-fold: (1) coefficients in high contrast and multiple scales and (2) demanding discretization for a high Péclet number. A lot of multiscale effort has contributed to the problems such as multiscale finite element method [1], variational multiscale method [2]–[5], multiscale discontinuous Galerkin method [6], [7], multiscale stablization [8], [9]. In particular, the Generalized Multiscale Finite Element Method (GMsFEM) aims to create multiscale basis to apply Galerkin approximation and it has been applied to an array of partial differential equations [10]–[16]. Similar to the finite element method, it consists of two stages: first the offline stage where basis functions are generated and used to span the approximation solution manifold; and second the online stage where the actual approximation is found in the generated space. However, due to the complexity of the problem, it is necessary to further reduce the computational cost of the offline stage.

A spectral decomposition method, Constraint Energy Minimization Generalized Multiscale Finite Element Method (CEM-GMsFEM), is then applied to the scheme to generate a sufficiently large solution manifold by using a few basis functions [17]–[24]. Each basis function is formed by the eigenfunction in a local spectral problem and captures some information about the medium and velocity, and thereby dependent on the Péclet number. One main important property of such eigenfunction is exponential decay. The span of such basis therefore can capture the cell decaying part of the solution while the non-decaying part is left to be dealt with in the online stage. However, it is only common to see the above results conducted in a homogeneous setting whilst, in practice, inhomogeneous boundary conditions are necessary. We propose a numerical scheme to apply CEM-GMsFEM to solve the convection-diffusion equation with inhomogeneous boundary conditions. There are two versions of CEM-GMsFEM: constrained and relaxed. The relaxed-CEM-GMsFEM is used to develop the error analysis, which takes advantage of the elliptic projection of the solution as a bridge between our multiscale space and the continuous space. An account of relaxed-CEM-GMsFEM applied to some inhomogeneous Boundary Value Problems was given already in [25]. Our goal is to extend this idea to the convection-diffusion equation. We will study three main cases, Dirichlet, Neumann, and Robin conditions. We also give a theoretical account of the convergence analysis and prove the scheme has first-order convergence in energy norm and second-order in \(L^2-\)norm, as verified by numerical examples.

For non-time independent problems, some work of applying CEM-GMsFEM has been done for parabolic equations [23] on homogeneous conditions. Our goal is to extend our method to the convection-diffusion equation with inhomogeneous boundary conditions. Since we assume the medium and velocity are independent of time, we can reuse the multiscale space in the time-independent case. For time-invariant boundary conditions, the corrector can be pre-computed once. However, for time-variant boundary conditions, the time derivative of the corrector and the boundary conditions need to be taken into account. Therefore, we give a new formulation to update the corrector at each time step in a relaxed CEM-GMsFEM fashion. The error of such approximation is also proved to be exponentially decaying in space. In addition, there are two versions of Backward Euler schemes to compute the next steps of the solution: the diffusion approach (D-approach) and the convection-diffusion approach (CD-approach). The former is less expensive but can be shown to bear a lower accuracy. We compare them via convection diffusion IBVPs and verify that our proposal gives a more accurate result at a higher computational expense.

Moreover, we introduce a nonlinearity into the IBVPs, which greatly increases the difficulty. A classical approach is to apply the Strang splitting method which considers convection, diffusion, and the nonlinearity terms in separate intermediate steps. This method utilizes the symmetric property of the algorithm to split operators in intermediate timesteps. This has been applied to a variety of problems including parabolic equations [26], diffusion-reaction equations [27], [28] and diffusion-reaction-advection equations with a homogeneous boundary conditions [29]. However, it was numerically tested that the inhomogeneous boundary conditions would drag the overall accuracy of the algorithm [29]–[31]. To resolve this, some classical approaches can be found in [31], [32]. More recent approaches involve designing a time-dependent boundary corrector in the intermediate steps [33]–[36], which is aligned with our previous sections. In particular, in the classical Strang splitting method applied to convection diffusion equations, the convection and diffusion operators are split and considered separately with their corresponding nonlinear terms [37]. Now with our scheme, we can consider them in the same step. We tested our choice of boundary corrector, the same as the previous part, to attain both spatial and temporal convergence.

The paper is organized as follows: we first give the problem setting and some preliminaries in section 2. The convergence analysis of its application to the time-independent convection diffusion equation is given in section 3, along with numerical results on the Dirichlet, Neumann, and Robin conditions. In section 4, both time-variant and time-invariant IBVPs are presented with analysis and numerical results, along with the comparison of different finite difference schemes for temporal discretization. A demonstration of applying this to nonlinear problems with Strang splitting is presented in this section as well.

2 Preliminaries↩︎

2.1 Problem setting↩︎

We consider the following convection diffusion initial boundary value problem: \[\begin{cases} \partial_t u + \boldsymbol{\beta}(x) \cdot \nabla u = \nabla\cdot(\boldsymbol{A}(x)\nabla u)+f, &\text{ in }\Omega\times (0,T],\\ u = g(x,t), &\text{ on }\Gamma_D\times (0,T],\\ b(x) u + \boldsymbol{\nu}\cdot(\boldsymbol{A}(x)\nabla u -\boldsymbol{\beta}(x) u) = q(x,t), &\text{ on }\Gamma_N\times (0,T],\\ u(\cdot,0)=u_\mathup{init}, &\text{ in }\Omega, \end{cases}\] where \(\Omega\subset \mathbb{R}^d\) is the computational domain, \(\boldsymbol{\beta}\in L^{\infty}(\Omega)^2\), and \(0<T<\infty\). The medium \(\boldsymbol{A}\in L^{\infty}(\Omega;\mathbb{R}^{d\times d})\) and the velocity \(\boldsymbol{\beta}\) are heterogeneous coefficients with multiple scales and potentially high contrast. \(\boldsymbol{A}\) is a positive definite matrix. There exist \(\kappa_0\) and \(\kappa_1\) such that \(0<\kappa_0\leq\lambda_{\min} (\boldsymbol{A}(x))\leq\lambda_{\max}(\boldsymbol{A}(x))\leq\kappa_1\) and \(\kappa_1/\kappa_0\) can be large. We denote \(\boldsymbol{\nu}\) the outward unit normal vectors to \(\partial\Omega\) and \(\Gamma_D\) and \(\Gamma_N\) two nonempty disjoint part of \(\partial\Omega\). We assume the velocity flows inward on the Neumann boundary \(\Gamma_N\) and \(\boldsymbol{\beta}\) is incompressible, i.e., \(\nabla\cdot\boldsymbol{\beta}=0\). i.e., \(\boldsymbol{\beta}\cdot\boldsymbol{\nu}\leq 0\) on \(\Gamma_N\). Denote \(\beta_0\geq 1\) and \(\beta_1\) as the infimum and supremum of \(|\boldsymbol{\beta}|\) respectively. The function \(b(x)\geq 0\) a.e.\(x\in\partial\Omega\) and there exists a positive constant \(b_0>0\) and a subset \(\Gamma \subset \Gamma_N\) with positive measure such that \(b(x)\geq b_0\) for a.e.\(x\in \Gamma\). The Dirichlet boundary value term \(g\in\mathcal{H}^{1/2}(\Gamma_D\times[0,T])\) and the Neumann boundary value term \(q\in L^2(\Gamma_N\times[0,T])\).

From now on, we denote \[\begin{align} a(w,v)&=\int_\Omega \boldsymbol{A}\nabla w\cdot \nabla v + \int_{\Gamma_N}(b-\boldsymbol{\beta}\cdot\boldsymbol{\nu}) wv \di \sigma,\\ \mathcal{A}(w,v)&=\int_\Omega \boldsymbol{A}\nabla w \cdot \nabla v +\int_{\Gamma_N}(b-\boldsymbol{\beta}\cdot\boldsymbol{\nu}) wv \di \sigma + \int_{\Omega}(\boldsymbol{\beta}\cdot\nabla w) v . \end{align}\]

So the variational form of the problem is: for \(t\in (0,T]\), find \(u_0(\cdot,t)\in V\mathrel{\vcenter{:}}=\{v\in\mathcal{H}^1(\Omega)\colon v=0 \text{ on }\Gamma_D\}\) such that at \(t\), \[\begin{align} (\partial_t u_0,v)+\mathcal{A}(u_0,v) &= (f,v) - \mathcal{A}(\widetilde{g},v) + (q,v)_{\Gamma_N} &\text{ for }v\in V\\ (u_0(\cdot,0),v)&=(u_\mathup{init}-\widetilde{g}(\cdot,0),v) &\text{ for }v\in V \end{align}\] where \(\widetilde{g}\in\mathcal{H}^1(\Omega\times[0,T])\) with \(\widetilde{g}=g\) on \(\Gamma_D\). The solution \(u\) to equation 4 is then \(u=u_0+\widetilde{g}\). Denote the energy norms on \(V\) \[\|v\|_{a(\omega)}=\left(\int_\omega \boldsymbol{A}\nabla v\cdot \nabla v + \int_{\Gamma_N\cap \partial\omega}(b-\boldsymbol{\beta}\cdot\boldsymbol{\nu}) v^2 \di \sigma \right)^{1/2}\text{ and } \|v\|_a = \|v\|_{a(\Omega)}.\] Also, notice that under the assumptions on \(\boldsymbol{\beta}\) and \(b\), for \(v\in V\), \[\begin{align} \mathcal{A}(v,v) &= \int_\Omega \boldsymbol{A}\nabla v\cdot \nabla v + \int_\Omega (\boldsymbol{\beta}\cdot \nabla v) v + \int_{\Gamma_N} (b-\boldsymbol{\beta}\cdot \boldsymbol{\nu}) v^2 \di \sigma\\ &= \int_\Omega \boldsymbol{A}\nabla v\cdot\nabla v+ \int_{\Gamma_N} \left(b-\frac{1}{2}\boldsymbol{\beta}\cdot\boldsymbol{\nu}\right) v^2 \di \sigma \geq 0. \end{align}\] So, the following quasi-norms are also well defined on \(V\): \[\|v\|^2_{\mathcal{A}(\omega)} = \int_\omega \boldsymbol{A} \nabla v \cdot \nabla v + \int_{\partial\omega\cap\Gamma_N}(b - \frac{1}{2}\boldsymbol{\beta}\cdot\boldsymbol{\nu}) v^2\di \sigma,\] for \(\omega \subset \Omega\) and \(\|v\|_\mathcal{A}=\|v\|_{\mathcal{A}(\Omega)}\). It is also easy to show that \(\|\cdot\|_\mathcal{A}\) and \(\|\cdot\|_a\) are equivalent on \(V\).

2.2 Constraint energy minimizing generalized multiscale finite element method↩︎

2.2.1 Auxiliary basis functions↩︎

We first introduce \(\mathcal{T}^H\) to be a conforming partition of \(\Omega\) into rectangular elements. The resulting elements are called the coarse grids. Let \(H\) be the meshsize of the coarse grid and \(N\) be the number of elements. For each \(K_i\in \mathcal{T}^H\) with \(1\leq i \leq N\), we define an oversampled domain \(K^m_i\) [25] by \[K^m_i\mathrel{\vcenter{:}}=\mathop{\mathrm{int}}\left\{\left(\bigcup_{K\in\mathcal{T}^H, \mathop{\mathrm{cl}}(K)\cap \mathop{\mathrm{cl}}(K^{m-1}_i)\neq \emptyset} \mathop{\mathrm{cl}}(K)\right)\cup \mathop{\mathrm{cl}}(K^{m-1}_i)\right\}\] where \(\mathop{\mathrm{int}}(S)\) and \(\mathop{\mathrm{cl}}(S)\) are the interior and the closure of a set \(S\). We also set \(K^0_i\mathrel{\vcenter{:}}= K_i\) for the consistency of notations. Let \(N_\nu\) be the number of vertices contained in an element. We can construct a set of Lagrange bases \(\{\eta^1_i,\dots, \eta^{N_\nu}_i\}\) of the element \(K_i\in\mathcal{T}^H\).

On each \(K_i\), we solve a local spectral problem: find \(\lambda^{j}_i\geq 0\), \(\phi^{j}_i\in \mathcal{H}^1(K_i)\) such that for all \(v\in\mathcal{H}^1(K_i)\), \[\mathcal{A}_{(K_i)}(\phi^{j}_i,v)=\lambda^{j}_i s_i(\phi^{j}_i,v), \label{eq:eigenvalueProblem}\tag{1}\] where \[s_i(w,v) = CH^{-2}\int_{K_i} \kappa_1 |\boldsymbol{\beta}|^2 wv \mathrel{\vcenter{:}}= \int_{K_i}\widetilde{\kappa}wv\] and \(C\) is a constant that depends on the choice of basis \(\{\eta_i\}\) and for the local auxiliary space on a structured mesh such that \(|\nabla \eta_i|^2 \leq CH^{-2}\). We take \(C = 24\) since we will use the Lagrange polynomials as the basis. For simplicity, we will omit the constant \(C\) in the analysis. Arrange the eigenvalues \(\{\lambda^j_i\}^{\infty}_{j=0}\) in an ascending order. Define the local auxiliary space \(V^\mathup{aux}_i\mathrel{\vcenter{:}}=\mathop{\mathrm{span}}\{\phi^{1}_i,\dots,\phi^{l_i}_i\}\) for some \(l_i\). Let \(s(w,v) = \sum^N_{i=1}s_i(w,v)\) on \(V\times V\). Denote \(\|w\|_{s(\omega)}=\sqrt{\int_{\omega} \widetilde{\kappa}w^2\di x }\) and \(\|w\|_{s}=\sqrt{s(w,w)}\).

Now, define the orthogonal projection \(\pi_i:L^2(K_i)\rightarrow V^\mathup{aux}_i\) by \[\pi_i(v)\mathrel{\vcenter{:}}=\sum_{j=0}^{l_i}\frac{s(\phi^{j}_i,v)}{s(\phi^{j}_i,\phi^{j}_i)}\phi^{j}_i.\]

Let \(V^\mathup{aux}\mathrel{\vcenter{:}}=\bigoplus^N_{i=1}V^\mathup{aux}_i\). Then \(s(\cdot,\cdot)\) and \(\|\cdot\|_s\) are an inner product and a norm on \(V^\mathup{aux}\) respectively. Also \(\pi \mathrel{\vcenter{:}}=\sum^N_{i=1}\pi_i\) maps \(L^2(\Omega)\) to \(V^\mathup{aux}\). Now, we can derive the following lemma [38]:

2.2.2 Multiscale basis functions↩︎

Although lemma 1 has shown that the elliptic projection can approximate any vector \(v\) close enough given sufficient number of eigenfunctions used, this function in \(V^\mathup{aux}\) may not be continuous in \(\Omega\). For this, We now construct the multiscale local basis functions using the auxiliary space \(V^\mathup{aux}\). Define \[\psi^{j}_i=\mathop{\mathrm{argmin}}\left\{\|\psi\|_\mathcal{A}^2+\|\pi\psi-\phi^{j}_i\|_s^2 \colon \psi\in V\right\},\] \[\psi^{j,m}_i=\mathop{\mathrm{argmin}}\left\{\|\psi\|_\mathcal{A}^2+\|\pi\psi-\phi^{j}_i\|_s^2 \colon \psi\in V^m_i\right\}.\] where \[V^m_i\mathrel{\vcenter{:}}=\left\{v\in\mathcal{H}^1(K^m_i)\colon v=0 \text{ on }\Gamma_D\cap\partial K^m_i\text{ or }\Omega\cap\partial K^m_i\right\}.\]

With these definitions, we can obtain \[\mathcal{A}(\psi^j_i,v)+s(\pi\psi^j_i,\pi v)=s(\phi^j_i,\pi v),\text{ } \forall v\in V, \label{eq:psi}\tag{2}\] \[\mathcal{A}(\psi^{j,m}_i,v)+s(\pi\psi^{j,m}_i,\pi v)=s(\phi^j_i,\pi v), \text{ }\forall v\in V^m_i. \label{eq:psi95m}\tag{3}\]

Denote \(V^\mathup{glo}_\mathup{ms}=\mathop{\mathrm{span}}\{\psi^j_i\colon 0\leq j\leq l_i,1\leq i\leq N\}\) and \(V^{m}_\mathup{ms}\mathrel{\vcenter{:}}=\mathop{\mathrm{span}}\{\psi^{j,m}_i\colon 0\leq j\leq l_i,1\leq i\leq N\}\).

It is worth-mentioning that the bilinear form \(\mathcal{A}\) is not symmetric, and therefore is not an inner product on \(V\). However, the following lemma still holds, without stating the proof here for simplicity [38].

3 Time-independent convection diffusion boundary value problems↩︎

As the multiscale space \(V^{m}_{ms}\) is independent of \(t\), the method is applicable to time independent problems. To illustrate the idea, consider Find \(u\in \mathcal{H}^1(\Omega)\) such that \[\begin{cases} -\nabla \cdot(\boldsymbol{A} \nabla u)+ \boldsymbol{\beta}\cdot\nabla u = f, &\text{ in }\Omega,\\ u=g, &\text{ on }\Gamma_D,\\ b u + \boldsymbol{\nu}\cdot (\boldsymbol{A} \nabla u- \boldsymbol{\beta}u)=q, &\text{ on }\Gamma_N, \end{cases} \label{eqn:BVP}\tag{4}\] where \(f\in L^2(\Omega)\) is the source term independent of \(u\).

3.1 Derivation of the Method↩︎

The methods are the following steps:

1. Find \(\mathcal{D}^m_i \widetilde{g}\in V^m_i\) and \(\mathcal{N}^m_i q\in V^m_i\) such that for all \(v\in V^m_i\), \[\mathcal{A}(\mathcal{D}^m_i\widetilde{g},v)+s(\pi \mathcal{D}^m_i\widetilde{g},\pi v)=\mathcal{A}_{(K_i)}(\widetilde{g},v),\] \[\mathcal{A}(\mathcal{N}_i^m q,v)+s(\pi\mathcal{N}_i^m q,\pi v)=\int_{\partial K_i\cap \Gamma_N}qv \di \sigma.\] Further denote \(\mathcal{D}^m\widetilde{g}=\sum^N_{i=1}\mathcal{D}^m_i\widetilde{g}\) and \(\mathcal{N}^m q=\sum^N_{i=1}\mathcal{N}_i^m q\).

2. Construct the multiscale function space \(V^m_\mathup{ms}\) according to equation 3 .

3. Solve \(w^m\in V^m_\mathup{ms }\) such that for all \(v\in V^m_\mathup{ms}\), \[\mathcal{A}(w^m,v)=(f,v)-\mathcal{A}(\widetilde{g},v)+\int_{\Gamma_N}qv \di \sigma +\mathcal{A}(\mathcal{D}^m\widetilde{g},v)-\mathcal{A}(\mathcal{N}^m q,v). \label{eq:w94m}\tag{5}\]

4. Construct the numerical solution \(u^\mathup{ms}\) to approximate the actual solution \(u\) of equation 4 by \[u^\mathup{ms}_0 = w^m-\mathcal{D}^m\widetilde{g}+\mathcal{N}^m q\text{ and } u^\mathup{ms}\approx u^\mathup{ms}_0+\widetilde{g}.\]

3.2 Analysis↩︎

We will use the following notation: for \(w,v\in V\) and \(\omega \subset \Omega\), \[\mathcal{B}_{(\omega)}(w,v) = \mathcal{A}_{(\omega)}(w,v)+s_{(\omega)}(\pi w,\pi v)\] and \(\|v\|_{\mathcal{B}(\omega)} \mathrel{\vcenter{:}}= \sqrt{\mathcal{B}_{(\omega)}(v,v)}\). Before we give an account of the analysis, we will summarize all the quasi-norms in the following lemma,

With the above lemmas, we can now start our error analysis.

3.2.1 Global Approximations↩︎

We will approximate the error using the global approximation \(w^\mathup{glo}\) of the real solution via \(V^\mathup{glo}_\mathup{ms}\). Define \(\mathcal{D}^\mathup{glo}\widetilde{g}=\sum^N_{i=1}\mathcal{D}^\mathup{glo}_i\widetilde{g}\) and \(\mathcal{N}^\mathup{glo} q=\sum^N_{i=1}\mathcal{N}^\mathup{glo}_i q\) where \(\mathcal{D}^\mathup{glo}_i\widetilde{g}\in V^\mathup{glo}_\mathup{ms}\) satisfies that for all \(v\in V\), \[\mathcal{B}(\mathcal{D}^\mathup{glo}_i\widetilde{g},v) = \mathcal{A}(\widetilde{g},v), \label{eq:Dgtg}\tag{6}\] and \(\mathcal{N}^\mathup{glo}_i q\in V^\mathup{glo}_\mathup{ms}\) satisfies that for all \(v\in V\), \[\mathcal{B}(\mathcal{N}^\mathup{glo}_i q,v)=\int_{\partial K_i\cap \Gamma_N} q v \di \sigma. \label{eq:Ngq}\tag{7}\]

Also define \(w^\mathup{glo}\in V^\mathup{glo}_\mathup{ms}\) such that, \[\mathcal{A}(w^\mathup{glo},v) - \mathcal{A}(\mathcal{D}^\mathup{glo}\widetilde{g}, v) + \mathcal{A}(\mathcal{N}^\mathup{glo} q, v) = \mathcal{A}(u_0,v) \text{ for } v\in V^\mathup{glo}_\mathup{ms}; \label{eq:w94glo95u0}\tag{8}\]

We first show the convergence of the elliptic projection defined using \(V^\mathup{glo}_\mathup{ms}\).

3.2.2 Abstract Problem↩︎

Following the approach in [25], we summarise the analysis of CEM-GMsFEM by considering the following abstract problem:

To solve this problem, we prepare ourselves with the following lemmas 4, 5 and 6. To prove them, we need to define cutoff functions \(\{\chi_i^{n,m}\}\). Let \(V^H\) be the Lagrange basis function space of \(\mathcal{T}^H\). For \(K_i\in\mathcal{T}^H\), a cutoff function \(\chi^{n,m}_i\in V^H\) with \(n<m\) satisfies that: \[\chi^{n,m}_i(x)=1\text{ in }K^n_i;\] \[\chi^{n,m}_i(x)=0\text{ in }\Omega\backslash K^m_i;\] \[0\leq\chi^{n,m}_i\leq 1\text{ in }K^m_i\backslash K^n_i.\]

3.2.3 Error Estimates for Boundary Correctors↩︎

We now directly apply these results to estimate \(\mathcal{D}^\mathup{glo}\widetilde{g}\) \(-\mathcal{D}^m\widetilde{g}\) and \(\mathcal{N}^\mathup{glo} q\) \(-\mathcal{N}^m q\).

To further analyse function spaces \(V^\mathup{glo}_\mathup{ms}\) and \(V^m_\mathup{ms}\), we define operators \(\mathcal{R}^\mathup{glo}\mathrel{\vcenter{:}}=\sum^N_{i=1}\mathcal{R}^\mathup{glo}_i\) \(\colon L^2(\Omega)\rightarrow V^\mathup{glo}_\mathup{ms}\) and \(\mathcal{R}^{m}\colon= \sum^N_{i=1}\mathcal{R}^{m}_i\colon L^2(\Omega)\rightarrow V^m_\mathup{ms}\) where \[\mathcal{A}(\mathcal{R}^\mathup{glo}_i\varphi,v)+s(\pi\mathcal{R}^\mathup{glo}_i\varphi,\pi v)= s(\pi_i\varphi,\pi v),\text{ for } v\in V; \label{eq:Rgloi}\tag{9}\] \[\mathcal{A}(\mathcal{R}^{m}_i\varphi,v)+s(\pi\mathcal{R}^{m}_i\varphi,\pi v)=s(\pi_i\varphi,\pi v),\text{ for } v\in V^m_i. \label{eq:Rmmi}\tag{10}\]

We also remark another lemma [38].

3.2.4 Error Estimate for the Multiscale Solution↩︎

We now state the main result.

3.3 Numerical Experiments↩︎

In this section, we will demonstrate the method via several numerical examples in a high-contrast setting and verify the significance of the inflow condition. For simplicity, we take point-wise isotropic coefficients, \(\boldsymbol{A}(x)=\kappa(x)\boldsymbol{I}\), the domain \(\Omega = [0,1]\times[0,1]\). We will calculate the reference solutions on a \(200\times200\) mesh with the bilinear Lagrange finite element method. The medium \(\kappa\) is presented in Figure 2 (a) and the source term in Figure 2 (b). The experiments are each tested for coarse mesh \(H = \frac{1}{10}, \frac{1}{20}\) and \(\frac{1}{40}\) with the fixed number \(l_m\) of eigenfunctions to generate the auxiliary space \(V^\mathup{aux}\). By our experiments, we tested \(l_m=3\) to be sufficient to verify our results.

3.3.1 Example 1↩︎

We first look at the Dirichlet condition, via considering the following problem: \[\begin{cases} -\nabla\cdot\left(\kappa(x_1,x_2) \nabla u\right)+\boldsymbol{\beta}(x_1,x_2)\cdot\nabla u=f &\text{ for } (x_1,x_2)\in\Omega\\ u(x_1,x_2)=\widetilde{g}(x_1,x_2)=x_1^2 + e^{x_1 x_2}&\text{ for }(x_1,x_2)\in \{0,1\}\times[0,1] \end{cases}\] where \[\boldsymbol{\beta}(x_1,x_2)=\left[\cos(18\pi x_2)\sin(18\pi x_1), -\cos(18\pi x_1)\sin(18\pi x_2)\right]^T.\]

To simplify notations, we denote the relative errors for the Dirichlet corrector \[D^m_a \mathrel{\vcenter{:}}=\frac{\|\mathcal{D}^m\widetilde{g}-\mathcal{D}^\mathup{glo}\widetilde{g}\|_\mathcal{A}}{\|\mathcal{D}^\mathup{glo}\widetilde{g}\|_\mathcal{A}} \text{ and } D^m_L \mathrel{\vcenter{:}}= \frac{\|\mathcal{D}^m\widetilde{g}-\mathcal{D}^\mathup{glo}\widetilde{g}\|_{L^2(\Omega)}}{\|\mathcal{D}^\mathup{glo}\widetilde{g}\|_{L^2}};\]

For the Neumann condition corrector, \[N^m_a \mathrel{\vcenter{:}}=\frac{\|\mathcal{N}^m q-\mathcal{N}^\mathup{glo} q\|_\mathcal{A}}{\|\mathcal{N}^\mathup{glo} q\|_\mathcal{A}} \text{ and } N^m_L \mathrel{\vcenter{:}}= \frac{\|\mathcal{N}^m q-\mathcal{N}^\mathup{glo} q\|_{L^2(\Omega)}}{\|\mathcal{N}^\mathup{glo} q\|_{L^2}};\] to measure errors and \(\Lambda^\prime = \max_i\lambda^{l_m}_i\). As for the error estimate, \[E^m_a \mathrel{\vcenter{:}}= \frac{\|u^\mathup{ms}-u\|_\mathcal{A}}{\|u\|_\mathcal{A}} \text{ and } E^m_L \mathrel{\vcenter{:}}= \frac{\|u^\mathup{ms}-u\|_{L^2}}{\|u\|_{L^2}}.\]

As can be seen in Table 1, the error of the Dirichlet corrector decays exponentially as we increase the number of oversampling layers. From Table 2, the number of eigenfunctions does improve the results for both the overall multiscale solutions and the boundary correctors. Most importantly, from Table 3, we can see a second-order convergence when increasing the oversampling layers and decreasing the coarse mesh. This also echoes the idea that the convergence of \(u^\mathup{ms}\) depends on the oversampling layers \(N_{ov}\).

3.3.2 Example 2↩︎

Another demonstration comes from the Neumann and Robin conditions. In this problem, we first consider the same velocity field \(\boldsymbol{\beta}\): \[\begin{cases} -\nabla\cdot\left(\kappa(x_1,x_2) \nabla u\right)+\boldsymbol{\beta}(x_1,x_2)\cdot\nabla u=f &\text{ for } (x_1,x_2)\in\Omega\\ u(x_1,x_2)=\widetilde{g}(x_1,x_2)=x_1^2 + e^{x_1 x_2}&\text{ for }(x_1,x_2)\in \{0,1\}\times[0,1]\\ bu + \boldsymbol{\nu}\cdot(A\nabla u -\boldsymbol{\beta}u) = -1 &\text{ for } x_1 = 0 \text{ and } x_2 \in [0,1]\\ bu + \boldsymbol{\nu}\cdot(A\nabla u -\boldsymbol{\beta}u) = 1 &\text{ for } x_1 = 1 \text{ and } x_2 \in [0,1]\\ bu + \boldsymbol{\nu}\cdot(A\nabla u -\boldsymbol{\beta}u) = 1 &\text{ for } x_1 = (0.5,1] \text{ and } x_2 =0\\ bu + \boldsymbol{\nu}\cdot(A\nabla u -\boldsymbol{\beta}u) = 0 &\text{ for } x_1 = [0, 0.5] \text{ and } x_2 = 0 \end{cases}\] where \(b(x_1,x_2)= \kappa(x_1,x_2)\) is the Robin coefficient.

For the following results, the numbers in brackets show the relative error for varying numbers of oversampling layers.

From Table 4, we also can observe the second-order convergence with respect to \(H\). Moreover, we compare the effect of the inflow conditions using the following two examples.

3.3.3 Example 3↩︎

In this model, we consider the inflow condition on \(\Gamma_N\). \[\boldsymbol{\beta}_\mathup{in}(x_1,x_2) = \boldsymbol{\beta}+ c_\mathup{flow}\left[\frac{1}{2}-x_1,x_1\right]^T.\] Notice that \(\boldsymbol{\beta}_\mathup{in}\cdot\boldsymbol{\nu}\leq 0\) on \(\Gamma_N\) as shown in Figure 3 (b). The constant \(c_\mathup{flow}\) is proportional to the magnitude of the velocity on the boundary \(\Gamma_N\),

As can be seen in Tables 5 and 6, similar observations of the results are shown, resembling our theoretical analysis and the numerical results in the examples.

3.3.4 Example 4↩︎

We consider the same setting as Example 3, but with the following velocity field \(\boldsymbol{\beta}_\mathup{out} = -\boldsymbol{\beta}_\mathup{in}\) with \(c_\mathup{flow}>0\). Note that \(\nabla\cdot\boldsymbol{\beta}_\mathup{out}=0\) on \(\Omega\) but \(\boldsymbol{\beta}_\mathup{out}\cdot\boldsymbol{\nu}> 0\) on \(x_1=0\) and \(x_1=1\), shown in Figure 3 (c).

In Table 7, second-order convergence in \(H\) is still observed when given a large enough number of oversampling layers and a fixed \(c_\mathup{flow}=3\). However, as we increase the velocity field \(c_\mathup{flow}=4\), we see that the problem is more demanding. Without a larger number of oversampling layers, reducing the coarse mesh \(H\) alone may not be enough to improve the approximation. However, when given a sufficient number of layers, the job can still be done. The comparison is more apparent in Figure 4 and Figure 5, where the outflow case converges much slower in both \(\|\cdot\|_{L^2}\) and \(\|\cdot\|_\mathcal{A}\) compared to the inflow case.

4 Time Dependent Convection-Diffusion Initial Boundary Value Problems↩︎

4.1 Derivation of the Method↩︎

We present the convection diffusion initial boundary value problem as follows: \[\begin{cases} \partial_t u(x,t) -\nabla \cdot(\boldsymbol{A}(x) \nabla u(x,t))+ \boldsymbol{\beta}(x) \cdot\nabla u(x,t) = f(x,t) &\text{ on }\Omega\times (0,T]\\ u(x,t) = g(x,t) &\text{ on }\Gamma_D\times (0,T]\\ b(x) u(x,t) + \boldsymbol{\nu}\cdot (\boldsymbol{A}(x) \nabla u(x.t)- \boldsymbol{\beta}(x) u(x,t))= q(x,t) &\text{ on }\Gamma_N\times (0,T]\\ u(x,0) = u_\mathup{init}(x) &\text{ on } \Omega. \end{cases} \label{eqn:unsteadyBVP}\tag{11}\] Notice \(f\) is still independent of the solution \(u\). The variational formulation becomes: find \(u_0(\cdot,t)\in V\) such that for \(v\in V\), \[\begin{align} (\partial_t u_0,v) + \mathcal{A}(u_0,v) &= (f,v) + (q,v)_{\Gamma_N} - \mathcal{A}(\widetilde{g},v) - (\widetilde{g}_t,v)\tag{12}\\ (u_0(\cdot, 0),v) &= (u_\mathup{init}(x)-\widetilde{g}(\cdot,0),v)\tag{13} \end{align}\] and the actual solution would be \(u(\cdot,t) = u_{0}(\cdot,t)+\widetilde{g}(\cdot,t)\). We will use the same auxiliary spaces \(V^\mathup{aux}\) and multiscale space \(V^m_\mathup{ms}\) as in the previous section, and thereby the same set of correctors \(\mathcal{D}^m \widetilde{g}\) and \(\mathcal{N}^m q\). In particular, for \(t\in (0,T]\),

1. Find \(\mathcal{D}^m_i \widetilde{g}(\cdot, t)\in V^m_i\) such that \(v\in V^m_i\), \[\begin{cases} (\mathcal{D}^m_i \widetilde{g}_t, v) + \mathcal{B}(\mathcal{D}^m_i \widetilde{g},v) = (\widetilde{g}_t, v)_{(K_i)} + \mathcal{A}_{(K_i)}(\widetilde{g},v)\\ \mathcal{B}(\mathcal{D}^m_i \widetilde{g}(\cdot, 0),v) = \mathcal{A}_{(K_i)}(\widetilde{g}(\cdot,0),v). \end{cases} \label{eq:dg95time95variant95ms}\tag{14}\] Denote \(\mathcal{D}^m \widetilde{g} = \sum_{i=1}^{N}\mathcal{D}^m_i \widetilde{g}\).

2. Find \(\mathcal{N}^m_i q(\cdot, t) \in V^m_i\) such that \(v\in V^m_i\), \[\begin{cases} (\mathcal{N}^m_i q_t, v) + \mathcal{B}(\mathcal{N}^m_i q,v)= \int_{\Gamma_N\cap \partial K_i}qv \di \sigma\\ \mathcal{B}(\mathcal{N}^m_i q(\cdot, 0),v) = \int_{\Gamma_N\cap\partial K_i }q(\cdot, 0)v \di \sigma. \end{cases} \label{eq:nq95time95variant95ms}\tag{15}\] Denote \(\mathcal{N}^m q = \sum^N_{i=1} \mathcal{N}^m_i q.\)

3. find \(w^m(\cdot, t) \in V^m_\mathup{ms}\) such that for \(v\in\) \(V^m_\mathup{ms}\) \[\begin{align} (w^m_t,v) + \mathcal{A}(w^m,v) &= (f,v) + (q,v)_{\Gamma_N} -\mathcal{A}(\widetilde{g},v)-(\widetilde{g}_t,v) \\ & + (\mathcal{D}^m \widetilde{g}_t, v) + \mathcal{A}(\mathcal{D}^m \widetilde{g}, v) -(\mathcal{N}^m q_t,v) -\mathcal{A}(\mathcal{N}^m q, v) , \end{align} \label{eqn:ms95t95weakBVP}\tag{16}\]

\[(w^m(\cdot, 0), v) = \left(u_\mathup{init}-\widetilde{g}(\cdot,0)+\mathcal{D}^m \widetilde{g}(\cdot,0)-\mathcal{N}^m q(\cdot,0),v\right) .\]

Then the multiscale approximation becomes \[u^\mathup{ms} = u^\mathup{ms}_0 + \widetilde{g} = w^m- \mathcal{D}^m \widetilde{g} + \mathcal{N}^m q+\widetilde{g}. \label{eq:time95dpdt95solution}\tag{17}\]

4.2 Analysis↩︎

To test the performance of the oversampling layers, we define \(\mathcal{D}^\mathup{glo}\widetilde{g}=\sum^N_{i=1}\mathcal{D}^\mathup{glo}_i\widetilde{g}\) and \(\mathcal{N}^\mathup{glo} q=\sum^N_{i=1}\mathcal{N}^\mathup{glo}_i q\) where \(\mathcal{D}^\mathup{glo}_i\widetilde{g},\mathcal{N}^\mathup{glo}_i q\in V^\mathup{glo}_\mathup{ms}\) satisfies that for all \(v\in V\), \[\begin{cases} (\mathcal{D}^\mathup{glo}_i \widetilde{g}_t, v) + \mathcal{B}(\mathcal{D}^\mathup{glo}_i \widetilde{g},v) = (\widetilde{g}_t, v)_{(K_i)} + \mathcal{A}_{(K_i)}(\widetilde{g},v)\\ \mathcal{B}(\mathcal{D}^\mathup{glo}_i \widetilde{g}(\cdot, 0),v) = \mathcal{A}_{(K_i)}(\widetilde{g}(\cdot,0),v), \end{cases} \label{eq:dg95time95variant}\tag{18}\] , and \[\begin{cases} (\mathcal{N}^\mathup{glo}_i q_t, v) + \mathcal{B}(\mathcal{N}^\mathup{glo}_i q,v)= \int_{\Gamma_N\cap \partial K_i}qv \di \sigma\\ \mathcal{B}(\mathcal{N}^\mathup{glo}_i q(\cdot, 0),v) = \int_{\Gamma_N\cap\partial K_i }q(\cdot, 0)v \di \sigma. \end{cases} \label{eq:nq95time95variant}\tag{19}\]

We will give an overview of the analysis. Define for \(v\in V\), \(\|v\|_\mathcal{E}^2 = \|v(\cdot,T)\|_{L^2}^2 + \int^T_0 \|v\|_\mathcal{B}^2\). Note that since \(\|\cdot\|_\mathcal{A}\) is a quasi-norm, so are \(\|\cdot\|_\mathcal{B}\) and \(\|\cdot\|_\mathcal{E}\). Suppose \(\widetilde{u}(\cdot,t)\) is the elliptic projection of the solution \(u\), i.e. \[\mathcal{A}(u-\widetilde{u},v) =0 \text{ for } v\in V^\mathup{glo}_\mathup{ms}.\] Akin to the previous treatment, our strategy is to decompose \(u-v\) into two parts: \[u-u^\mathup{ms} = (u-\widetilde{u})+ (\widetilde{u}-u^\mathup{ms}).\] The error of the former term is computable while that of the latter term can be bounded by a specific choice \(v\in V^m_\mathup{ms}\). In particular, for any \(v\in V^m_\mathup{ms}\), \[\begin{align} &\int^T_0\left[((u-u^\mathup{ms})_t,u-u^\mathup{ms}) + \mathcal{A}(u-u^\mathup{ms},u-u^\mathup{ms})\right]\\ &=\int^T_0\left[((u-u^\mathup{ms})_t,u-v) + \mathcal{A}(u-u^\mathup{ms},u-v)\right]\\ &\leq \left.(u-u^\mathup{ms},u-v)\right|^T_0 -\int^T_0 ((u-v)_t,u-u^\mathup{ms}) + \int^T_0 \mathcal{A}(u-u^\mathup{ms},u-v)\\ &\leq \|(u-u^\mathup{ms})(\cdot,T)\|_{L^2}\|(u-v)(\cdot,T)\|_{L^2} +\|(u-u^\mathup{ms})(\cdot,0)\|_{L^2}\|(u-v)(\cdot,0)\|_{L^2}\\ &+\sqrt{\int^T_0\|(u-v)_t\|_{L^2}^2}\sqrt{\int^T_0\|u-u^\mathup{ms}\|_{L^2}^2} +\frac{3}{4}{\int^T_0\|u-u^\mathup{ms}\|_\mathcal{A}^2} +\frac{1}{3}\overline{C}^2{\int^T_0\|u-v\|_\mathcal{B}^2}. \end{align}\]

Hence, by the repeated use of Young’s inequality and the Cauchy Schwartz inequality, \[\begin{align} \|(u-u^\mathup{ms})(\cdot,T)\|_{L^2}^2+\int^T_0\|u-u^\mathup{ms}\|_\mathcal{A}^2 &\leq 4\|(u-u^\mathup{ms})(\cdot,0)\|_{L^2}^2 + 2\|(u-v)(\cdot,0)\|_{L^2}^2+ 4\overline{C}^2\|u-v\|_{\mathcal{E}}^2\\ &+4\sqrt{\int^T_0\|(u-v)_t\|_{L^2}^2}\sqrt{\int^T_0\|u-u^\mathup{ms}\|_{L^2}^2}\\ &=: (i) + (ii) + (iii) + (iv). \end{align} \label{eq:main95t}\tag{20}\]

Note that \((i)\) has been investigated in the time-independent case. Now, define \(w^\mathup{glo}\in V^\mathup{glo}_\mathup{ms}\) such that for any \(v\in V\), \[\mathcal{A}(w^\mathup{glo},v) = \mathcal{A}(\widetilde{u} - \mathcal{D}^\mathup{glo}\widetilde{g} + \mathcal{N}^\mathup{glo} q,v). \label{eq:elliptic95projection95t}\tag{21}\]

Then, with the same set of \(\mathcal{R}^\mathup{glo}\mathrel{\vcenter{:}}=\sum^N_{i=1}\mathcal{R}^\mathup{glo}_i\) and \(\mathcal{R}^{m}\mathrel{\vcenter{:}}=\sum^N_{i=1}\mathcal{R}^{m}_i\) operators in the last section, there is \(\varphi(\cdot,t) \in L^2(\Omega)\) such that \[\begin{align} (\mathcal{R}^\mathup{glo}_i\varphi_t,v)+\mathcal{B}(\mathcal{R}^\mathup{glo}_i\varphi,v)=s(\pi_i \varphi, \pi v)&\text{ for }v\in V,\\ (\mathcal{R}^{m}_i\varphi_t,v)+\mathcal{B}(\mathcal{R}^{m}_i\varphi,v)= s(\pi_i \varphi,\pi v)&\text{ for }v\in V^m_\mathup{ms}. \end{align}\] with initial conditions \(\mathcal{B}(\mathcal{R}^\mathup{glo}_i\varphi(\cdot,0),v)=s(\pi_i\varphi(\cdot,0),\pi v)\) for \(v\in V\) and \(\mathcal{B}(\mathcal{R}^{m}_i\varphi(\cdot,0),v)=s(\pi_i\varphi(\cdot,0),\pi v)\) for \(v\in V^m_\mathup{ms}\). Then, by the surjectivity of \(\mathcal{R}^\mathup{glo}\), we can find \(\varphi_*(\cdot,t)\in L^2(\Omega)\) such that \[\mathcal{R}^\mathup{glo}\varphi_* = w^\mathup{glo} \text{ and define } w^m_*\mathrel{\vcenter{:}}= \mathcal{R}^{m}\varphi_*.\]

By putting \(v=w^m_* - \mathcal{D}^m \widetilde{g} + \mathcal{N}^m q\), we can further decompose \[u - v = (u-\widetilde{u}) + (w^\mathup{glo}-w^m_*)+ (\mathcal{D}^\mathup{glo}\widetilde{g}-\mathcal{D}^m \widetilde{g}) + (\mathcal{N}^\mathup{glo} q-\mathcal{N}^m q).\]

The error analysis of these terms will suffice that of the remaining terms \((ii), (iii)\), and \((iv)\) in equation 20 . \(u-\widetilde{u}\) is first dealt with and the rest in another abstract problem.

4.2.1 Elliptic Projection↩︎

We first give an error bound for the elliptic projection.

We also can obtain the \(L^2\)-norm of the global estimate \(\widetilde{u}\).

4.2.2 Abstract Problem↩︎

We now move on to the analysis of the corrector. The main idea is that the error propagation has an exponential decay with respect to the oversampling layers, similar to the time-independent case.

We aim to estimate \[\left\|\sum^N_{i=1}{\mathcal{P}_i t_i}-{\mathcal{P}^m_i t_i}\right\|_\mathcal{E}^2 =\left\|\sum^N_{i=1}({\mathcal{P}_i t_i}-{\mathcal{P}^m_i t_i})(\cdot,T)\right\|_{L^2}^2+ \int^T_0\left\|\sum^N_{i=1}{\mathcal{P}_i t_i}-{\mathcal{P}^m_i t_i}\right\|_\mathcal{B}^2 \di t.\]

One should note that the initial condition here is exactly the abstract problem in the previous section. The results are carried over here. In addition, we will further define two norms for our analysis, \[\|t_i\| = \max_{v\in V}\frac{\left<t_i,v\right>}{\|v\|_\mathcal{B}} \text{ and } \|t^0_i\| = \max_{v\in V}\frac{\left<t^0_i,v\right>}{\|v\|_\mathcal{B}}.\]

In parallel, we can have \[\|{\mathcal{P}^m_i t_i}\|^2_{\mathcal{E}^2(\Omega\backslash K^m_i)}\leq \theta^m \left((m+1)H^2c_\#^2\|t^0_i\|^2+\int^T_0 \|t_i\|^2\right). \label{eq:32Pmiti95E}\tag{24}\]

4.2.3 Multiscale Solutions↩︎

We will first bound \(\int^T_0\|u-u^\mathup{ms}\|_{L^2}^2\).

We can now combine all the results to bound the terms in equation 20 .

4.3 Temporal Discretization↩︎

We apply the Backward Euler method to the scheme. Explicitly, we let \(\tau\) be the time step and \(U^n_\mathup{ms}\mathrel{\vcenter{:}}= u^\mathup{ms}(t_n)\). Then, the variational formulation 16 becomes \[\begin{align} \left(\frac{U^n_\mathup{ms}-U^{n-1}_\mathup{ms}}{\tau},v\right) + \mathcal{A}(U^n_\mathup{ms},v) = (f(t_n),v) &\text{ for } v\in V^m_\mathup{ms},\\ U^0_\mathup{ms} = u^0_\mathup{ms}. & \end{align} \label{eqn:beuler95weakBVP}\tag{25}\]

4.3.1 Backward Euler Schemes↩︎

We compare two versions of the Backward Euler Scheme for convection diffusion, using the Dirichlet boundary condition as an illustration. The Neumann and Robin conditions follow similar treatments. Let \(\tau\) be the timestep.

1. \[\frac{u^{n+1}-u^n}{\tau}+ (\boldsymbol{\beta}\cdot\nabla u^{n+1})=\nabla \cdot (A \nabla u^{n+1})\]

2. \[\frac{u^{n+1}-u^n}{\tau}+ (\boldsymbol{\beta}\cdot\nabla u^{n})=\nabla \cdot (A \nabla u^{n+1}).\]

Respectively, the application of the method becomes: given the multiscale solution \(u^{n}\) at step \(n\),

1. Convection-Diffusion approach (CD-approach)

   1. Find \(\mathcal{D}^m_i \widetilde{g}^{n+1}\in V^{m}_\mathup{ms}\) such that for \(v\in V^{m}_\mathup{ms}\) such that \[\begin{align} \left(\frac{\mathcal{D}^m_i \widetilde{g}^{n+1}-\mathcal{D}^m_i \widetilde{g}^n}{\tau},v\right) + \mathcal{B}(\mathcal{D}^m_i \widetilde{g}^{n+1},v)=\mathcal{A}_{(K_i)}(\widetilde{g}^{n+1},v)+(\widetilde{g}_t^{n+1},v)_{(K_i)}. \end{align}\] Then, set \(\mathcal{D}^m \widetilde{g}^{n+1} = \sum^N_{i=1}\mathcal{D}^m_i \widetilde{g}^{n+1}\).

   2. Find \(w^{n+1}\in V^{m}_\mathup{ms}\) such that for \(v\in V^m_\mathup{ms}\), \[\begin{align} \left(\frac{w^{n+1}-w^{n}}{\tau},v\right)+\mathcal{A}(w^{n+1},v) &= (f^{n+1},v)-\mathcal{A}(\widetilde{g}^{n+1},v)-(\widetilde{g}^{n+1}_t,v)\\ &+\mathcal{A}(\mathcal{D}^m \widetilde{g}^{n+1},v)+(\mathcal{D}^m \widetilde{g}^{n+1}_t,v). \end{align}\]

   3. Set \(u_\mathup{ms}^{n+1}=w^{n+1}-\mathcal{D}^m \widetilde{g}^{n+1}+\widetilde{g}^{n+1}\).

2. Diffusion approach (D-approach)

   1. Find \(\mathcal{D}^m_i \widetilde{g}^{n+1}\in V^{m}_\mathup{ms}\) such that for \(v\in V^{m}_\mathup{ms}\) such that \[\begin{align} \left(\frac{\mathcal{D}^m_i \widetilde{g}^{n+1}-\mathcal{D}^m_i \widetilde{g}^n}{\tau},v\right) + a(\mathcal{D}^m_i \widetilde{g}^{n+1},v)+s(\pi \mathcal{D}^m_i \widetilde{g}^{n+1},\pi v)=a_{(K_i)}(\widetilde{g}^{n+1},v)+(\widetilde{g}_t^{n+1},v)_{(K_i)}. \end{align}\] Then, set \(\mathcal{D}^m \widetilde{g}^{n+1} = \sum^N_{i=1}\mathcal{D}^m_i \widetilde{g}^{n+1}\).

   2. Find \(w^{n+1}\in V^{m}_\mathup{ms}\) such that for \(v\in V^m_\mathup{ms}\), \[\begin{align} \left(\frac{w^{n+1}-w^{n}}{\tau},v\right)+a(w^{n+1},v) &= (f^{n+1},v)-a(\widetilde{g}^{n+1},v)-(\widetilde{g}^{n+1}_t,v)\\ &-(\boldsymbol{\beta}\cdot\nabla u^n,v)+a(\mathcal{D}^m \widetilde{g}^{n+1},v)+(\mathcal{D}^m \widetilde{g}^{n+1}_t,v). \end{align}\]

   3. Set \(u_\mathup{ms}^{n+1}=w^{n+1}-\mathcal{D}^m \widetilde{g}^{n+1}+\widetilde{g}^{n+1}\).

4.4 Nonlinear Convection Diffusion IBVPs↩︎

In this subsection, we are interested in the convection diffusion with a nonlinearity term \(f(u)\): \[\begin{cases} \partial_t u + \boldsymbol{\beta}\cdot\nabla u = f(u) + \nabla\cdot(\boldsymbol{A}u), &\text{ in } \Omega\times[0,T],\\ u = g, &\text{ on }\Gamma_D\times[0,T],\\ b u + \boldsymbol{\nu}\cdot(\boldsymbol{A}\nabla u -\boldsymbol{\beta}u) = q, &\text{ on }\Gamma_N\times[0,T],\\ u(\cdot,0) = u_\mathup{init}, &\text{ on }\Omega. \end{cases}\]

The traditional approach is to perform operator splitting. It is to decompose the convection-diffusion operator into two sub-problems, each targeting one operator [37]. However, with inhomogeneous boundary conditions, some operators are not left-invariant [33]. Therefore, a correction term \(\eta_n\) needs to be introduced. On the other hand, with the current scheme on convection diffusion equations, the convection and diffusion operators can be considered at once. In other words, our ultimate goal is to construct the solution at the next step via: \[u^{n+1} = (\mathcal{S}^{\tau/2}_{f-\eta_n}\circ \mathcal{S}^\tau_{CD+\eta_n}\circ \mathcal{S}^{\tau/2}_{f-\eta_n})(u^n)\] where \(\mathcal{S}^{\tau/2}_{f-\eta_n}\) maps to the solution considering the adjusted nonlinearity term \(f-\eta_n\) for a stepsize \(\tau/2\) and \(\mathcal{S}^\tau_{CD+\eta_n}\) maps to the solution considering the convection diffusion equation with the source term \(\eta_n\). Also, this corrector term \(\eta_n\) is to be carefully selected and get cancelled out in the process; otherwise can be accumulated over time and become a stiff term. The problem becomes even more apparent with time-variant boundary conditions. In light of this, we will decompose the nonlinearity into two parts, one dependent on the boundary and one independent. We will illustrate this idea using time-invariant Dirichlet boundary conditions. The time-variant version extends naturally by following a similar fashion. The proof of this well-known method can be found in [33], [34], [36]. To summarise: at each step \(n\),

1. Define \(z\in L^2(\Omega)\) such that \(\left.z\right|_{\partial \Omega}=g(\cdot,t_n)\).

2. Find \(w(\frac{\tau}{2})\) by \[\begin{cases} \partial_t w = f(w+z) - f(z)\\ w(0) = u^n-z. \end{cases}\]

3. Find \(v(\tau)\) by \[\begin{cases} \partial_t v +\mathcal{A}v = f(z) - \mathcal{A}z -\partial_t z\\ v(0) = w(\frac{\tau}{2}). \end{cases}\]

4. Find \(w(\frac{\tau}{2})\) by \[\begin{cases} \partial_t w = f(w+z) - f(z)\\ w(0) = v(\tau). \end{cases}\]

5. Define \(z\in L^2(\Omega)\) such that \(\left.z\right|_{\partial \Omega}=g(\cdot,t_{n+1})\).

6. Define \(u^{n+1} = w(\frac{\tau}{2}) + z\).

Above, the introduction of the corrector \(f(z)\) is independent of the solution in Step 3. Therefore, \(f(z)-\mathcal{A}z-\partial_t z\) can be treated as a source term of a convection diffusion equation, which has been discussed already. We note that for the first step, one has the flexibility to determine the interior value of \(z\). Aligned with our previous direction, we will use our choice of \(\mathcal{D}^m\widetilde{g}-\widetilde{g}\) and \(\mathcal{N}^m q\) for \(z\), catering to the corresponding boundary conditions.

4.5 Numerical Results↩︎

We demonstrate the numerical experiments for time-invariant and time-variant boundary conditions. \(\Omega\), \(\kappa\), and \(\boldsymbol{\beta}\) are taken as in section 3. For simplicity, we assume the source term vanishes \(f=0\). We compare the CD-approach and D-approach. Experiments will be run on the coarse mesh \(H=\frac{1}{10},\frac{1}{20}\) and \(\frac{1}{40}\) with the number of eigenfunctions \(l_m=3\) and a fixed timestep \(\tau = \frac{1}{10}\). Reference solutions are generated on a \(200\times200\) mesh with \(1000\) timesteps with the bilinear Lagrange finite element method. Without further mentioning, the error terms are recorded at the final time \(T=1\).

4.5.1 Time-invariant Boundary Conditions↩︎

We first consider time-invariant boundary conditions: \[\begin{cases} u_t - \nabla\cdot(\boldsymbol{A}\nabla u) + \boldsymbol{\beta}\cdot\nabla u = 0 & \text{ in } [0,1]^2 \times (0,1]\\ u(x,t) = x_1^2+e^{x_1x_2} &\text{ for } (\partial [0,1]^2 )\times (0,1]\\ u(x,0) = x_1^2 + e^{x_1 x_2} &\text{ on } [0,1]^2. \end{cases}\]

Since the boundary condition is independent of time, \(\mathcal{D}^m \widetilde{g}_t\) and \(\widetilde{g}_t\) vanish in the numerical method. The error analysis followed directly from the time-independent problem in section 3.

From Table 8, both second-order convergence in \(L^2\)-norm and at least first-order convergence in the energy norm can be observed for both D- and CD-approaches. Similar running times are recorded for both cases. However, since the convection term is dependent on the velocity field, the CD-approach outperforms D-approach in Table 9. Not only does it perform relatively poorly, the D-approach also fails to achieve the convergence in \(\|\cdot\|_\mathcal{A}\) with respect to \(H\).

4.5.2 Time-variant Dirichlet Boundary Conditions↩︎

Consider the following: \[\begin{cases} u_t - \nabla\cdot(\kappa\nabla u) + \boldsymbol{\beta}\cdot\nabla u = 0 & \text{ in } [0,1]^2\times [0,1]\\ u(x,t) = (x_1^2+e^{x_1 x_2})e^{-t} &\text{ for } (\partial [0,1]^2 )\times [0,1]\\ u(x,0) = x_1^2 + e^{x_1 x_2} &\text{ on } [0,1]^2. \end{cases}\]

As can be seen in Table 10, the D-approach is more efficient than the CD-approach, at the cost of a higher sensitivity to the Péclet number. CD-approach is more accurate at the cost of longer computational time. Even though both present spatial convergence but the CD-approach is more robust. The time-variant case also verifies our theoretical error estimates by showing second-order (or higher) convergence in \(L^2\)-norm and first-order convergence in the energy norm with respect to \(H\). As in Table 11, the increase in the number of oversampling layers does allow decay in the corrector error, and therefore verifies our theoretical convergence in \(H\) as shown in Table 12. The solution profiles can be found in Figure 6.

4.5.3 Nonlinear Example↩︎

We demonstrate combining CEM-GMsFEM with Strang Splitting to solve the following nonlinear convection-diffusion problem with a time-invariant Dirichlet boundary condition: \[\begin{cases} \partial_t u-\nabla \cdot(\boldsymbol{A}(x)\nabla u)+ \boldsymbol{\beta}(x) \cdot\nabla u = u-u^3 &\text{ on }[0,1]^2\times (0,1]\\ u = x_1^2+e^{x_1x_2} &\text{ on }(\partial [0,1]^2 )\times (0,1]\\ u(\cdot,0) = x_1^2 + e^{x_1x_2} &\text{ on } [0,1]^2 \end{cases} \label{eqn:splitting}\tag{26}\] with \(\boldsymbol{\beta}= \boldsymbol{\beta}_{in}\) and \(c_\mathup{flow}=\frac{1}{4}\).

The reference solutions are generated on a \(200\times 200\) mesh with \(1000\) steps using the bilinear Lagrange finite element method and Backward Euler for the time discretization. Running tests on a combinations of stepsize \(\tau\in \{\frac{1}{10},\frac{1}{20},\frac{1}{40}\}\), mesh size \(H\in\{\frac{1}{10},\frac{1}{20},\frac{1}{40}\}\) and oversampling layers \(N_\mathup{ov}\in \{7,8,9\}\).

As shown in Table 13, the first order convergence in energy norm with respect to space and time are guaranteed. However, to achieve second order convergence in \(L^2-\)norm, higher oversampling layer is needed. The results could potentially be improved via adapting other temporal discretization scheme such as exponential integration [26].

5 Conclusions↩︎

In this paper, we propose an application of CEM-GMsFEM to solve convection diffusion equations under various types of inhomogeneous boundary conditions along with high-contrast coefficients. The method begins with constructing an auxiliary space and builds a multiscale space upon it. Boundary correctors are built upon this multiscale space by solving local oversampled cell problems. For both time independent and dependent problems, we provide convergence anaylsis and show second order convergence in \(L^2\) and first order convergence in energy norm with respect to the coarse mesh, given sufficient oversampling, as agreed with numerical results. We also compare different time discretization strategies using the Backward Euler scheme. For nonlinear problems, we apply this modified method with Strang Splitting and demonstrate using a Dirichlet initial boundary value problem.

6 Acknowledgements↩︎

The research of Eric Chung is partially supported by the Hong Kong RGC General Research Fund (Projects: 14304021 and 14305423).

References↩︎