Generalized multiscale finite element method for the steady state linear Boltzmann equation

TL;DR

This paper develops and analyzes a generalized multiscale finite element method for the steady state linear Boltzmann equation, effectively handling highly oscillatory media and small Knudsen numbers with spectral convergence.

Contribution

The paper introduces a GMsFEM approach with proven well-posedness and spectral convergence, capable of homogenizing and preserving asymptotic limits across scales.

Findings

Spectral convergence of the basis functions independent of media oscillation

Method achieves numerical homogenization and asymptotic preservation

Proven well-posedness of the GMsFEM for the Boltzmann equation

Abstract

The Boltzmann equation, as a model equation in statistical mechanics, is used to describe the statistical behavior of a large number of particles driven by the same physics laws. Depending on the media and the particles to be modeled, the equation has slightly different forms. In this article, we investigate a model Boltzmann equation with highly oscillatory media in the small Knudsen number regime, and study the numerical behavior of the Generalized Multi-scale Finite Element Method (GMsFEM) in the fluid regime when high oscillation in the media presents. The Generalized Multi-scale Finite Element Method (GMsFEM) is a general approach to numerically treat equations with multi-scale structures. The method is divided into the offline and online steps. In the offline step, basis functions are prepared from a snapshot space via a well-designed generalized eigenvalue problem (GEP), and these basis functions are then utilized to patch up for a solution through DG formulation in the online step to incorporate specific boundary and source information. We prove the wellposedness of the method on the Boltzmann equation, and show that the GEP formulation provides a set of optimal basis functions that achieve spectral convergence. Such convergence is independent of the oscillation in the media, or the smallness of the Knudsen number, making it one of the few methods that simultaneously achieve numerical homogenization and asymptotic preserving properties across all scales of oscillations and the Knudsen number.

Figures

Tables

Table 1. Table 1: Errors for Example 1 with ϵ = 10 − 1 italic-ϵ superscript 10 1 \epsilon=10^{-1} and a δ = κ 4 superscript 𝑎 𝛿 superscript 𝜅 4 a^{\delta}=\kappa^{4} .

$L$	snapshot ratio	$e_1$	$e_2$
1	0.79%	19.64%	9.41%
2	1.59%	17.68%	8.53%
3	2.38%	14.41%	7.40%
5	3.97%	8.11%	4.92%
7	5.56%	6.16%	3.62%
10	7.94%	3.44%	1.62%
15	11.90%	2.24%	1.04%
20	15.87%	1.64%	0.68%

Table 2. Table 3: Errors for Example 1 with ϵ = 10 − 3 italic-ϵ superscript 10 3 \epsilon=10^{-3} and a δ = κ 4 superscript 𝑎 𝛿 superscript 𝜅 4 a^{\delta}=\kappa^{4} .

$L$	snapshot ratio	$e_1$	$e_2$
1	0.79%	12.80%	12.80%
2	1.59%	26.43%	26.42%
3	2.38%	17.86%	17.85%
5	3.97%	4.45%	4.43%
7	5.56%	3.60%	3.59%
10	7.94%	3.55%	3.53%
15	11.90%	3.20%	3.18%
20	15.87%	3.19%	3.17%

Table 3. Table 5: Errors for Example 2 with ϵ = 5 × 10 − 2 italic-ϵ 5 superscript 10 2 \epsilon=5\times 10^{-2} .

$L$	snapshot ratio	$e_1$	$e_2$
1	0.79%	22.70%	9.73%
2	1.59%	20.36%	8.43%
3	2.38%	16.97%	8.13%
5	3.97%	11.94%	6.86%
7	5.56%	8.09%	4.64%
10	7.94%	4.70%	1.99%
15	11.90%	2.48%	1.22%
20	15.87%	1.86%	0.91%

Table 4. Table 7: Errors for Example 2 with ϵ = 5 × 10 − 4 italic-ϵ 5 superscript 10 4 \epsilon=5\times 10^{-4} .

$L$	snapshot ratio	$e_1$	$e_2$
1	0.79%	14.12%	14.11%
2	1.59%	20.87%	20.86%
3	2.38%	11.69%	11.69%
5	3.97%	2.95%	2.95%
7	5.56%	2.71%	2.71%
10	7.94%	2.71%	2.71%
15	11.90%	2.88%	2.88%
20	15.87%	2.93%	2.92%

Equations

$${\partial }_{t}u+v\cdot \nabla u=\sigma \mathcal{S}u(x,v)-\eta u,(x,v)\in \Omega \times V.$$

$$\mathcal{R}u(x,v)={\int }_{V}u(x,v^{\prime })d{v}^{\prime }-u(x,v),$$

$$v\cdot \nabla u+ϵu=\frac{1}{ϵ}\sigma \mathcal{R}u(x,v).$$

$$u(x,v)\sim \rho (x),$$

$$C\nabla \cdot (\frac{1}{\sigma }\nabla \rho )=\rho ,$$

$$v\cdot \nabla u+ϵu=\frac{1}{ϵ}{\sigma }^{\delta }\mathcal{R}u(x,v),$$

$$\begin{array}{rl}\mathbf{v}\cdot \nabla u(x,\mathbf{v})+ϵu(x,\mathbf{v})& =\frac{1}{ϵa^{\delta }}\mathcal{R}u(x,\mathbf{v})\text{ in }\Omega \times S^1\\ u(x,\mathbf{v})& =g(x,\mathbf{v})\text{ on }{\Gamma }^{-}\end{array},$$

$${\Gamma }^{-}=\left\{(x,v)\in \partial \Omega \times S^1|\mathbf{v}\cdot {\mathbf{n}}_x<0\right\}.$$

$$\mathcal{R}u(x,\mathbf{v})=\overline{u}(x)-u(x,\mathbf{v})=\frac{1}{2\pi }{\int }_{S^1}u(x,\mathbf{v})d\mathbf{v}-u(x,\mathbf{v}).$$

$$\sum _{i=1}^m{\alpha }_i=1,\text{ and }\frac{1}{2\pi }{\int }_{S^1}u(x,\mathbf{v})d\mathbf{v}\approx \sum _{i=1}^m{\alpha }_{i}u(x,{\mathbf{v}}_i).$$

$$\begin{array}{rl}{\mathbf{v}}_i\cdot \nabla u_i+ϵu_i+\frac{1}{ϵa^{\delta }}\left(u_i-\sum _{j=1}^m{\alpha }_j{u}_j\right)& =0\text{in}\Omega ,\\ u_i& =g_i\text{on}{\Gamma }^{-},\end{array}$$

$$a_{ij}=\{\begin{array}{cc}{\alpha }_i-{\alpha }_i^2,& i=j,\\ -{\alpha }_i{\alpha }_j,& i\ne j.\end{array}.$$

$${\mathbf{v}}_i\cdot \nabla u_i+ϵu_i+\frac{1}{ϵa^{\delta }{\alpha }_i}\sum _j{a}_{ij}{u}_j=0.$$

$${\mathsf{u}}^{\top }\mathsf{Au}=\sum _{i,j=1}^m{a}_{ij}{u}_j{u}_i=$$

$$=$$

$${\mathsf{v}}^{\top }\mathsf{Au}=\sum _{i,j=1}^m{a}_{ij}{u}_j{v}_i=$$

$$=$$

$$\mathcal{L}({\eta }_l^{(i)})=0\text{in}{K}_i$$

$$\begin{array}{r}{a}_{\Omega }(u-u_H,u-u_H)⪯\sum _K{a}_K(u^K-u_H^K,u^K-u_H^K),\end{array}$$

$$a_{K_i}(\eta -{\eta }_0,\eta -{\eta }_0)⪯\beta s_{K_i}(\eta -{\eta }_0,\eta -{\eta }_0),$$

$$V_h=$$

$$u_h=(u_{h,1},u_{h,2},...,u_{h,m})\in (V_h{)}^m.$$

$$\begin{array}{rl}& -{\int }_{\Omega }{u}_{h,i}\nabla w_i\cdot {\mathbf{v}}_i+\sum _{e\in {\mathcal{E}}_H^0}{\int }_e{u}_{h,i}^{+}[w_i]\cdot {\mathbf{v}}_i+\sum _{e\in {\Gamma }^{+}}{\int }_e{u}_{h,i}{w}_i{\mathbf{v}}_i\cdot \mathbf{n}+{\int }_{\Omega }ϵu_{h,i}{w}_i\\ +& {\int }_{\Omega }\frac{1}{ϵa^{\delta }}\left(u_{h,i}-\sum _{j=1}^m{\alpha }_j{u}_{h,j}\right)w_i=-\sum _{e\in {\Gamma }^{-}}{\int }_e{g}_i{w}_i{\mathbf{v}}_i\cdot \mathbf{n},\forall w_i\in V_h,\end{array}$$

$$\begin{array}{rl}& \sum _{i=1}^m{\alpha }_i\left(-{\int }_{\Omega }{u}_{h,i}\nabla w_i\cdot {\mathbf{v}}_i+\sum _{e\in {\mathcal{E}}_H^0}{\int }_e{u}_{h,i}^{+}[w_i]\cdot {\mathbf{v}}_i+\sum _{e\in {\Gamma }^{+}}{\int }_e{u}_{h,i}{w}_i{\mathbf{v}}_i\cdot \mathbf{n}+{\int }_{\Omega }ϵu_{h,i}{w}_i\right)\\ & +{\int }_{\Omega }\frac{1}{ϵa^{\delta }}\sum _{i,j=1}^m{a}_{ij}{u}_{h,j}{w}_i=-\sum _{i=1}^m{\alpha }_i\sum _{e\in {\Gamma }^{-}}{\int }_e{g}_i{w}_i{\mathbf{v}}_i\cdot \mathbf{n},\forall w\in (V_h{)}^m,\end{array}$$

$$[w]=\{\begin{array}{cc}{w}^{-}{\mathbf{n}}^{-}+w^{+}{\mathbf{n}}^{+}& \text{on }{\mathcal{E}}_H^0\\ w^{-}{\mathbf{n}}^{-}& \text{on }{\Gamma }^{-}\\ w^{+}{\mathbf{n}}^{+}& \text{on }{\Gamma }^{+}\end{array}.$$

$$\begin{array}{rl}a(u,w)=& \sum _{i=1}^m{\alpha }_i{a}_i(u_i,w_i),\text{with}{a}_i(u_i,w_i)=-{\int }_{\Omega }{u}_i\nabla w_i\cdot {\mathbf{v}}_i+\sum _{e\in {\mathcal{E}}_H^0}{\int }_e{u}_i^{+}[w_i]\cdot {\mathbf{v}}_i+\sum _{e\in {\Gamma }^{+}}{\int }_e{u}_i{w}_i{\mathbf{v}}_i\cdot \mathbf{n},\\ l(u,w)=& \sum _{i=1}^m{\alpha }_i{\int }_{\Omega }\frac{1}{ϵa^{\delta }}\left(u_i-\sum _{j=1}^m{\alpha }_j{u}_j\right)w_i+\sum _{i=1}^m{\alpha }_i{\int }_{\Omega }ϵu_i{w}_i={\int }_{\Omega }\frac{1}{ϵa^{\delta }}\sum _{i,j=1}^m{a}_{ij}{u}_j{w}_i+\sum _{i=1}^m{\alpha }_i{\int }_{\Omega }ϵu_i{w}_i,\end{array}$$

$$F(w)=\sum _{i=1}^m{\alpha }_i{F}_i(w_i),\text{with}{F}_i(w_i)=-\sum _{e\in {\Gamma }^{-}}{\int }_e{g}_i{w}_i{\mathbf{v}}_i\cdot \mathbf{n}.$$

$$a(u_h,w)+l(u_h,w)=F(w),\forall w\in (V_h{)}^m.$$

$$u_H=(u_{H,1},u_{H,2},...,u_{H,m})\in V_H$$

$$\begin{array}{rl}& \sum _{i=1}^m{\alpha }_i\left(-{\int }_{\Omega }{u}_{H,i}\nabla w_i\cdot {\mathbf{v}}_i+\sum _{e\in {\mathcal{E}}_H^0}{\int }_e{u}_{H,i}^{+}[w_i]\cdot {\mathbf{v}}_i+\sum _{e\in {\Gamma }^{+}}{\int }_e{u}_{H,i}{w}_i{\mathbf{v}}_i\cdot \mathbf{n}+{\int }_{\Omega }ϵu_{H,i}{w}_i\right)\\ & +{\int }_{\Omega }\frac{1}{ϵa^{\delta }}\sum _{i,j=1}^m{a}_{ij}{u}_{H,j}{w}_i=-\sum _{i=1}^m{\alpha }_i\sum _{e\in {\Gamma }^{-}}{\int }_e{g}_i{w}_i{\mathbf{v}}_i\cdot \mathbf{n},\forall w\in V_H.\end{array}$$

Full text

1 Introduction

The Boltzmann equation is a fundamental model in statistical mechanics. It traces the evolution of the distribution function on the phase space, and describes the dynamics of a large number of particles that follow the same physics rules via a statistical manner. The equation encodes the particles’ free transport and their interactions with the media and each other. Depending on the physics the particles follow, the interaction term may differ, but to a large extent, many particles, including neutrons, photons and phonons interact mainly with the media, making the collision term linear. The dynamics then can be described by the linear Boltzmann equation:

[TABLE]

In the equation, $u$ is a function on the phase space $(x,v)\in\Omega\times V$. The evolution is governed by $v\cdot\nabla u$, a free transport term, and the terms on the right side of the equation that represent the “collision” and quantify the particles’ interactions. These interactions include a pure absorption term $\eta u$ where $\eta$ is the absorption coefficient, and a scattering term $\sigma\mathcal{S}u$. The specific form of the operator $\mathcal{S}$ varies from particle to particle, and it is typically a functional independent of $x$. The strength of the interactions are governed by the size of $\sigma$ and $\eta$. For photons specifically, the radiative transfer equation is used, these coefficients are termed the optical thickness. In this paper, for simplicity, we take $\mathcal{R}$ to be:

[TABLE]

where $dv$ is a normalized measure: $\int dv=1$, and we set $\eta=1$.

The equation demonstrates different behavior in different regimes. One particularly interesting regime is called the diffusion regime, in which the scattering coefficient is extremely strong and the pure absorption term is weak. Mathematically, consider the steady state case, we model the equation to:

[TABLE]

In this equation, $\epsilon$ is termed the Knudsen number, and it characterizes the ratio of the mean free path and the typical domain length. Physically it reflects the number of collision a normal particle experiences inside $\Omega$ before emitting. When $\epsilon$ is small, the number of collision per particle is large, meaning the particle gets scattered many times before emitting, and thus some kind of averaging effects take place, and the local equilibrium is achieved. In the case of (2), the equilibrium reads:

[TABLE]

and through asymptotic analysis, one could mathematically derive that $\rho$ satisfies the diffusion equation:

[TABLE]

where $C$ depends on the dimension.

The convergence from (3) to the asymptotic limit (5) was conjectured in [4] and was made rigorous in [3] for periodic boundary condition. In [29] the authors studied the boundary layer effect with geometric corrections, and the asymptotic convergence rate was shown to degrade [26, 27].

However, all the rigorous proofs are done assuming certain smoothness of $\sigma$. In particular, it is assumed that $\sigma$ is sufficiently smooth. At the current stage, very limited works have been done when oscillations present in the media. Denote $\delta$ the small scale in the media, we rewrite our equation as:

[TABLE]

where $\sigma^{\delta}(x)=\sigma(x,\frac{x}{\sigma})$ to explicit reflects the fast variable $\frac{x}{\sigma}$ dependence. On the theoretical level, to our best knowledge, except a few cases [22], the theory is largely in lack, except a few cases [22], and to a large extend, we do not yet know the resonance of the two parameters, and how they contribute in the asymptotic limits of the equation. And on the computation level, the only numerical study aware to the authors is presented in [25] where the limits are taken in order: $\delta\ll\epsilon\ll 1$.

The problem is very challenging on the numerical level. The small $\epsilon$ makes the collision term extremely stiff, bringing ill-conditioning to the associated discrete system, and thus severe stability issue; and the small $\delta$ brings wild oscillations to the media and the solution, and for accuracy of the numerical solution, high resolution is needed and small discretization is necessary, driving up the numerical cost.

This is certainly unaffordable, especially in the zero limit of $\epsilon$ and $\delta$. The main goal of this paper is to develop a general numerical treatment that could deal with the equation with a wide range of $\epsilon$ and $\delta$, and perform uniformly well, with the error term independent on the small parameters.

The approach we take is in line of Generalized Multi-scale Finite Element Method (GMsFEM). This is an offline-online framework that builds a good set of local basis functions during the offline step and patches local solutions up in the online step, similar to the original Multi-scale Finite Element Method (MsFEM). One main feature of GMsFEM is its basis selection procedure in the offline step where a special generalized eigenvalue problem is designed. This special generalized eigenvalue problem encodes the oscillations and the ill-conditioning of the problem.

More explicitly, like many other multi-scale methods, we build nested grids with coarse grid $H$ and fine grid $h$ satisfying $H\gg\epsilon\gg h$. In the offline step, local basis functions are constructed within coarse mesh $H$ on fine mesh $h$ that capture fine scale structure and preserve the heterogeneities in the media; and in the online step, the basis functions are patched up through Galerkin framework [2, 6, 17, 14, 15, 16, 20, 10, 11, 9, 19, 1, 24]. Online step is rather standard and different methods give various algorithms in the offline step. What makes GMsFEM favorable is indeed its offline step, in which the full list of $a$-harmonic functions are collected, and then the most “representative” modes are selected through a specially designed generalized eigenvalue problem (GEP). The definition of the matrices in the GEP is associated with the final error term, which permits certain spectral error decay. We should mention GMsFEM was initially used for elliptic equations containing strong heterogeneous media, a topic about which the literature is extremely rich. For this particular problem, there is another category of method: upscaling-type methods. In upscaling methods, either locally or globally an effective media is numerically prepared so that equations could be computed on coarser grids with the effective media replacing the heterogeneous one [12, 30, 18, 21]. But this approach is not going to be pursued in this paper.

As a framework, the GMsFEM approach is rather easy to use, and the main mathematical challenge, when utilized for tackle different equations, is to develop the right GEP. For the linear Boltzmann equation with heterogeneous media, we frame the problem in the discontinuous Galerkin setting, and are able to find two matrices that resemble the mass and stiffness matrices in the GEP of the elliptic equations, which allows us to show the optimality of the basis functions with physically meaning definition of the norm for the error. As the standard approach, these basis functions are then used in the online computation.

The paper will be organized in the following way. In Section 2, we introduce some preliminaries. Both discrete ordinates, the standard kinetic solver, and GMsFEM for elliptic equation will be presented. Some properties will be presented along. We present the algorithm in Section 3, which is further divided into two subsection introducing offline and online procedures. Section 4 contains the analysis where we present the wellposedness, and convergence results. The small $\epsilon$ limit of the method will also be discussed. Numerical results are to be shown in Section 5.

To end the introduction, we comment that the scaling problem studied in this paper is not mathematically artificial. In fact, as one redefines $x\to\frac{x}{\epsilon}$, $\sigma(x)$ should have been automatically changed to highly oscillatory media $\sigma(x/\epsilon)$. Another practical example is to inject light into crystals, where the radiative transfer equation (one particular linear Boltzmann equation) is utilized. In this case, the periodic crystal structure should be encoded in the media and the period that corresponds to $\delta$ in our math formulation, is expected to be small.

2 Preliminaries

In this section we prepare some basic important concepts. In particular, we will first present the discrete ordinate method for the linear Boltzmann equation, and then give a brief account of the Generalized Multiscale Finite Element Method (GMsFEM). They are the building blocks for the algorithm designed in this paper.

2.1 $S_{N}$ Boltzmann Equation

The Boltzmann equation gives a statistical description of particle dynamics. Its extensive use in all kinds of engineering problems brought its great popularity, and literature on both theory and numerics has been very rich. Among all numerical methods developed for the Boltzmann equation, the discrete ordinate method stands out for its simplicity and intuitiveness, and is the method we will use in our GMsFEM. Essentially it discretizes the velocity domain, and the semi-discrete system is a coupled PDE in the physical space.

We start the discussion with the following model equation:

[TABLE]

where $x\in\Omega\subset\mathbb{R}^{2}$ is a bounded domain with a Lipschitz boundary $\partial\Omega$. The velocity is $v\in S^{1}$, the unit circle. The media $a^{\delta}$ presents fine scale structure at $\delta$ order, and the stiffness of the collision operator $\mathcal{R}$ is determined by $\frac{1}{\epsilon}\gg 1$. We have the inflow boundary condition, with the inflow data $g(x,\mathbf{v})$ defined on $\Gamma^{-}$, a collection of coordinates on the boundary with velocity pointing into the domain:

[TABLE]

Here $\mathbf{n}_{x}$ is the unit outer normal direction at $x\in\partial\Omega$. For simplicity, we use the model collision operator with homogeneous scattering coefficient:

[TABLE]

The discrete ordinate method, denoted by $S_{N}$, is a standard method to discretize the velocity domain. One first sample $m$ quadrature points on $S^{1}$ and each sample point is associated with a weight, denoted by: $\{(\mathbf{v}_{i},\alpha_{i}),i=1,...,m\}$, where $\mathbf{v}_{i}$ are the quadrature points and $\alpha_{i}$ are the corresponding positive weights. These quadrature points and weights are chosen so that:

[TABLE]

The equation then will be discretized into semi-discrete system. Let $u_{i}(x)=u(x,\mathbf{v}_{i})$, the integro-differential (7) is then transformed into a system of $m$ coupled partial differential equations:

[TABLE]

where $g_{i}=g(x,\mathbf{v}_{i})$ is the inflow boundary data. Denote

[TABLE]

Then (9) is further simplified to:

[TABLE]

Since $\{a_{ij}\}$ is basically the discrete version of the collision operator $-\mathcal{R}$, it resembles the properties of $\mathcal{R}$. In particular, the matrix is positive semi-definite with a known kernel.

Proposition 1**.**

Define a matrix $\mathsf{A}$ so that $\mathsf{A}_{ij}=a_{ij}$, we claim:

• •

$\mathsf{A}$* is positive semi-definite.*

• •

$\mathsf{u}^{\top}\mathsf{A}\mathsf{u}=0$* if and only if $\mathsf{u}=(u_{1},u_{2},...,u_{m})$ is isotropic, i.e., $u_{1}=u_{2}=\cdots=u_{m}$.*

• •

$\mathsf{v}^{\top}\mathsf{A}\mathsf{u}=0$* if either $\mathsf{u}$ or $\mathsf{v}$ is isotropic.*

Proof.

The computation is straightforward:

[TABLE]

The equal sign is achieved only when $u_{i}=u_{j}$ for all $i\neq j$.

To show the third bullet point, we note that the matrix $\mathsf{A}$ is symmetric, it suffices to assume that $\mathsf{u}$ is isotropic: $u_{1}=u_{2}=\cdots=u_{m}=\overline{u}$. Then:

[TABLE]

where we have used the weight condition (8): $\sum_{i}\alpha_{i}=1$. ∎

2.2 Generalized Multi-scale Finite Element Method

The discrete ordinate method is used to discretize the velocity domain, and for the spatial domain, we follow the GMsFEM approach, which, by choosing “optimal basis functions” via a special design of a GEP, we can obtain a reduced model that is robust for all values of $\epsilon$ and $\delta$. For the completeness of the paper, we now present a general idea of GMsFEM, and its application to the heterogeneous Boltzmann equation will be discussed in details in Section 3.

The GMsFEM uses two stages: offline and online. In the offline stage, a small dimensional approximation space is constructed to solve the global problem for any external source on a coarse grid, which does not need to resolve any scales of the media and solution. The offline stage consists of two main concepts. The snapshot space, $V_{\text{snap}}^{i}$, is constructed for a generic coarse element $K_{i}$. The snapshot solutions are used to compute local multiscale basis functions. An appropriate snapshot space can

• •

provide a faster convergence,

• •

provide problem relevant restrictions on the coarse spaces (e.g., divergence free solutions),

• •

reduce the cost associated with constructing the offline spaces.

Standard choices of snapshot spaces (see [7]) are (1) all fine-grid functions; (2) snapshots of local solutions; (3) oversampling snapshots of local solutions; and (4) force-based snapshots. In this paper, we will use snapshots of local solutions.

More specifically, these are functions $\eta_{l}^{(i)}$ that satisfy

[TABLE]

subject to some boundary conditions, where $\mathcal{L}$ is the differential operator under consideration, and $l$ is the index for the boundary condition. One can use all fine grid delta functions as boundary conditions or randomized boundary conditions [7, 5].

The offline space, $V_{H}$, is computed for each $K_{i}$ (with elements of the space denoted $\psi_{l}^{(i)}$). We perform a spectral decomposition in the snapshot space and select the dominant eigenfunctions (corresponding to the smallest eigenvalues) to construct the offline (multiscale) space. The convergence rate of the resulting method is proportional to $1/\Lambda_{*}$, where $\Lambda_{*}$ is the smallest eigenvalue that the corresponding eigenvector is not included in the multiscale space. We would like to select local spectral problem such that we can remove many small eigenvalues with fewer multiscale basis functions. The choice of spectral problems is usually problem dependent and is based on convergence analysis. In general, the error is decomposed into coarse subdomains. The energy functional corresponding to the domain $\Omega$ is denoted by $a_{\Omega}(u,u)$, e.g., $a_{\Omega}(u,u)=\int_{\Omega}\kappa\nabla u\cdot\nabla u$. Then,

[TABLE]

where $K$ are coarse regions ($K_{i}$), $u^{K}$ is the localization of the solution. The local spectral problem is chosen to bound $a_{K}(u^{K}-u_{H}^{K},u^{K}-u_{H}^{K})$. We seek the subspace $V_{H}^{i}$ such that for any $\eta\in V_{\text{snap}}^{i}$, there exists $\eta_{0}\in V_{H}^{i}$ with,

[TABLE]

where $s_{K_{i}}(\cdot,\cdot)$ is an auxiliary bilinear form, and $\beta$ is an accuracy parameter. The auxiliary bilinear form needs to be chosen such that the solution is bounded in the corresponding norm.

Finally, in the online stage, the space $V_{H}$ is used together with a suitable coarse grid discretization to solve the problem. The same space $V_{H}$ is used for all input sources.

3 GMsFEM for heterogeneous Boltzmann equation

We now apply the GMsFEM approach to numerically study the heterogeneous Boltzmann equation, expressed in the discrete ordinate system (9).

The numerical difficulties in solving this equation are summarized as follows. First, the media $a^{\delta}$ is highly oscillatory. The fine structure oscillates at the scale of $\delta$ injects high heterogeneities to $u_{i}$. In order to capture these details, the mesh size $h$ has to be smaller than $\delta$, which in turn brings prohibitive numerical cost. Secondly, the operator $\mathcal{L}$ is scaled by $\frac{1}{\epsilon}$, and in the zero limit of $\epsilon$, the term is extremely stiff, and this brings concern in stability. It is our aim in this paper to develop a multiscale method that can address these issues. In particular, inspired by GMsFEM, we will design a numerical method that relies on offline basis construction and online basis patching procedure, and its numerical error has limited dependence on the two small parameters.

We will construct nested grids and call $\mathcal{T}^{h}$ the partition of $\Omega$ into fine finite elements, and $\mathcal{T}^{H}$ the partition into coarse elements, where $h$ and $H$ are the fine and coarse mesh sizes respectively. For simpler notation, we consider rectangular coarse elements as shown in Figure 1. The basis functions and discretization are based on the coarse grid, and the fine grid is used to numerically compute the basis functions. We also denote the collection of coarse edges $\mathcal{E}_{H}$, and $\mathcal{E}_{H}^{0}=\mathcal{E}_{H}\backslash\partial\Omega$ the collection of coarse edges in the interior of the domain.

The discontinuous Galerkin method allows one to pick different values of the solution on different sides of the edges. Suppose two adjacent coarse blocks $\tau_{i}$ and $\tau_{j}$ share an edge, and that $\tau_{i}$ is the upwind block, then we denote $w^{+}=w|_{\tau_{i}}$ and $w^{-}=w|_{\tau_{j}}$. Notice that depending on the direction of a specific $\mathbf{v}_{i}$, different block could be picked as the upwind block, as shown in Figure 2.

For the fine scale approximation, we choose the discrete function space to be:

[TABLE]

and we seek for numerical solution such that

[TABLE]

This means the numerical solution for each $u_{h,i}$, when confined in each fine grid, is a linear function, and is continuous function across coarse grids. In the variational formulation: for all $i=1,2,...,m$, we have:

[TABLE]

or with the definition of $\{a_{ij}\}$ in (10), they could be summed up to

[TABLE]

In the equation, we have used upwind approximation for $\mathbf{v}_{i}\cdot\nabla u_{i}$ and the jump operator $[\cdot]$ is defined by

[TABLE]

For notational simplicity, we define two bilinear operators

[TABLE]

and a linear operator:

[TABLE]

With these notations, equation (15) now writes

[TABLE]

With $h\ll\text{min}\{\epsilon\,,\delta\}$, it is a standard result that $u_{h}\approx u$, with an error term of size $\mathcal{O}(h^{2}/\text{min}(\epsilon,\delta))$. For significantly small $h$, the function $u_{h}$ is considered as a reference solution in accessing the performance of our method.

However, using small $h$ that resolves $\epsilon$ and $\delta$ leads to a very big system that is numerically very costly. We would like to develop an algorithm that seeks for solution only on the coarse grid $H$ and the corresponding solution $u_{H}\approx u_{h}\approx u$. To do that, an offline-online procedure developed in [5] for elliptic equation, termed GMsFEM (Generalized Multiscale Finite Element Method) will be pursued. In the offline step, an approximate space $V_{H}$ is constructed to replace $V_{h}$. This newly constructed space would have much less degrees of freedom but preserves $V_{h}$’s important factors. The final multiscale solution will be computed in the online step where the boundary condition $g(x,\mathbf{v})$ will be taken into account to determine the degrees of freedom in $V_{H}$.

We quickly review the online stage in Section 3.1, and the complicated offline step will be discussed in detail in Section 3.2.

3.1 Online computation

In online stage, we will use the multiscale basis functions together with a coarse grid discretization to solve the given problem. The coarse grid discretization we used is a discontinuous Galerkin method with upwind flux. Assume that a multiscale finite element space $V_{H}=\text{span}\{\phi_{p}\}$ is determined, and this space, in some sense, approximates $(V_{h})^{m}$. Then similar to the formulation as in (15), the solution will be sought in

[TABLE]

so that

[TABLE]

Similar to (16), we use a compact notation:

[TABLE]

To implement the scheme above, we define the following matrices

[TABLE]

Then the multiscale solution $u_{H}$ is formulated as

[TABLE]

where the coefficient vector $U$ solves $(A+L)U=b$.

3.2 Construction of $V_{H}$

The key to the success of our method is the construction of the space $V_{H}$ on the coarse mesh during the offline stage. We will give the details here.

As discussed in Section 2.2, the offline step is further decomposed into two substages: constructing the snapshot space, and selecting modes associated with small eigenvalues. These two stages will be presented in Section 3.2.1 and Section 3.2.2 respectively.

In the snapshot space construction stage, in each coarse region, the Boltzmann solution will be solved multiple times together with all possible boundary conditions resolved by the fine grid. This give a high dimensional space. However, some modes in the snapshot space are more important than the others, and they dominate the numerical solution. To identify these basis functions, a specially designed local spectral problem (generalized eigenvalue problem GEP) is formulated and solved. The modes that correspond to the smallest eigenvalues are selected for form $V_{H}$. The number of modes to be selected depends on the error tolerance and the eigenvalues of the GEP. The design of the local spectral problem is to encode the convergence error that is to be discussed in Section 4.

3.2.1 Snapshot Space

We present the construction of the snapshot spaces in this subsection. The procedure is the same in each coarse element, and we take the coarse element $K_{j}$ as an example. The snapshot space for this particular element is denoted by $V_{\text{snap}}^{j}$. We use the notation $J^{i}(D)$ to denote the set of all nodes of the fine mesh $\mathcal{T}^{h}$ lying in the upwind part of $\partial D$ associated with velocity $\mathbf{v}_{i}$. And we also use $J(D)=\bigoplus_{i=1}^{m}J^{i}(D)$ to denote the union. Then the snapshot space is simply the linear span of solutions to the local Boltzmann equation with delta function as boundary condition, namely:

[TABLE]

where $\leftidx{{}_{n}}\psi_{l}^{j}=(\leftidx{{}_{n}}\psi_{l,1}^{j},\leftidx{{}_{n}}\psi_{l,2}^{j},...,\leftidx{{}_{n}}\psi_{l,m}^{j})$ solves:

[TABLE]

Here we use multi-index Kronecker delta function $\delta_{l}\mathbf{e}_{n}$, where $\mathbf{e}_{n}$ is the standard basis in $\mathbb{R}^{m}$ and $\delta_{l}$ is the standard Kronecker delta function:

[TABLE]

This strategy is summarized in Algorithm DetLocal.

The full snapshot space is given by

[TABLE]

Remark 1**.**

Numerically to prepare all snapshot basis functions is hard. It requires the computation of local Boltzmann equation with a large number of possible incoming delta functions. To reduce the cost of computation, we use the idea of oversampling [15]. To do so, the local computational domain is slightly enlarged to $K_{j}^{+}$ (see Figure 1), and a collection of random boundary condition is imposed on $K_{j}^{+}$. The low rank structure of the solution space allows one to correctly capture the range, even with a limited sampling. In particular, we define the snapshot space

[TABLE]

where $k_{j}$ is the number of snapshot functions we could customize, and $\leftidx{{}_{n}}\psi_{l}^{j,+}=(\leftidx{{}_{n}}\psi_{l,1}^{j,+},\leftidx{{}_{n}}\psi_{l,2}^{j,+},...,\leftidx{{}_{n}}\psi_{l,m}^{j,+})$ solves:

[TABLE]

where $r_{l}$ are random i.i.d. Gaussian sampling on $J(K^{+}_{j})$. The solutions $\leftidx{{}_{n}}\psi_{l}^{j,+}$ confined on $K_{j}$ are then used to form the snapshot spaces. We remark that the use of randomized boundary conditions on oversampling domains is able to reduce the offline computational cost as there is no need to impose delta function boundary conditions as in (24).

This strategy is summarized in Algorithm RanLocal.

Similar to (16) and (18), we can solve the snapshot solution $u_{\text{snap}}\in V_{\text{snap}}$ by the following equation:

[TABLE]

We note that the snapshot solution can be considered as a reference solution. The error of the snapshot solution is related to the approximation property of the snapshot space in the fine scale space.

3.2.2 Offline Space

Now, we will present the construction of the solution space $V_{H}$, with the property we mentioned in (13). In the end $V_{H}$, when confined on each coarse element, say $K_{j}$, will be a subspace of $V^{j}_{\text{snap}}$, and the construction of $V_{H}$ is amount to finding the most appropriate basis functions in $V^{j}_{\text{snap}}$ to be included. The procedure is further divided into two sub-steps: the energy minimizing oversampling (EMO), and a design of a generalized eigenvalue problem, as used in [8].

We first denote the local oversampled snapshot space $\bigoplus_{K_{i}\subset K_{j}^{+}}V_{\text{snap}}^{i}$ by $V_{\text{snap}}^{j,+}$. Notice that, for a given coarse element $K_{j}$ and its corresponding oversampling region $K_{j}^{+}$, the space $V_{\text{snap}}^{j,+}$ is the union of all snapshot spaces $V_{\text{snap}}^{i}$ with the condition that $K_{i}\subset K_{j}^{+}$.

Then the energy minimizing snapshots are calculated. For any snapshot function $\psi\in V_{\text{snap}}^{j}$, its energy minimizing extension $\widetilde{\psi}$ has the smallest energy in some norm $a_{\text{Energy}}^{j}(\cdot,\cdot)$ and is sought in the local oversampled snapshot space $V_{\text{snap}}^{j,+}$ with the constraint $\widetilde{\psi}|_{K_{j}}=\psi|_{K_{j}}$. In mathematical expression, for any $\psi\in V_{\text{snap}}^{j}$, we seek for $\widetilde{\psi}\in V_{\text{snap}}^{j,+}$ so that

[TABLE]

in which

[TABLE]

We notice that this construction is well-defined and the strategy is summarized in Algorithm EnergyMin. As one can see, $\widetilde{\psi}$ is an extension of $\psi$ onto the oversampling domain that achieves the minimum energy, defined in (32) and this extension is crucial, as will be seen in the later analysis.

Remark 2**.**

This is about a stable decomposition property. It is important that the local basis functions satisfy a stable decomposition property. More precisely, the sum of local energies is bounded by the global energy.

Next, we define the two bilinear operators $a_{K_{j}}(\cdot,\cdot)$ and $s_{K_{j}}(\cdot,\cdot)$, mentioned in (13). For simplicity of notation, we use $a^{j}(\cdot,\cdot)$ and $s^{j}(\cdot,\cdot)$ instead. For the element $K_{j}$, define:

[TABLE]

Using the above bilinear forms, a spectral problem is defined. On $K_{j}$, we look for $\left(\phi_{k}^{j},\lambda_{k}^{j}\right)\in V_{\text{snap}}^{j}\times\mathbb{R}$ such that

[TABLE]

where the eigenvalues are ordered in the ascending way:

[TABLE]

For implementation, we define the following matrices

[TABLE]

Then the pair $\left(\phi_{k}^{j},\lambda_{k}^{j}\right)$ is computed by solving

[TABLE]

Suppose $L_{j}$ modes are used for each $K_{j}$. This strategy is summarized in Algorithm LocalGEP. The offline space $V_{H}$ is given by

[TABLE]

This will be the approximation space for solving the system (9) in the scheme (18).

3.3 Algorithm summary

We finally summarize the algorithm. Largely speaking, we prepare the basis functions in the offline step, and patch them up in the online step. The offline step is further divided into preparing a snapshot space in which either all Green’s functions are accumulated or a good random selection is obtained, and the basis selection step, in which the local generalized eigenvalue problem is computed and eigenfunctions with highest energies are chosen. These basis functions are ultimately used in the online step via the weak formulation (20). We summarize the procedure in Algorithm 1.

4 Analysis of the GMsFEM

In this section, we will present some analysis of our GMsFEM. In Section 4.1, we will prove the well-posedness of the discrete system resulting from the GMsFEM, and in Section 4.2, we will prove the convergence of the method. Finally, in Section 4.3, we will analyze the behavior of the method when $\epsilon$ is small.

4.1 Well-posedness

We first show the well-posedness of the GMsFEM (18).

Theorem 1**.**

Problem (18) has a unique solution, and the solution $u_{H}$ satisfies the following stability condition

[TABLE]

Proof.

Since the system (18) is a square linear system, showing the existence and uniqueness is amount to proving that $a(\widehat{u},w)+l(\widehat{u},w)=0$ for all $w\in V_{H}$ only for trivial solution $\widehat{u}=0$.

We will first prove the following inequalities

[TABLE]

First, $l$ is non-negative since the matrix $(a_{ij})$ is a positive semi-definite matrix, as discussed in Proposition 1. Next, the non-negativity of $a(\cdot,\cdot)$ is shown below:

[TABLE]

Assuming $a(\widehat{u},w)+l(\widehat{u},w)=0$ for any $w\in(V_{h})^{m}$, then setting $w=\widehat{u}$, we have

[TABLE]

According to Proposition 1, we have

[TABLE]

meaning $\widehat{u}=0$, and the solution to (18) is thus unique. For stability, we start with

[TABLE]

Considering $a(u_{H},u_{H})+l(u_{H},u_{H})=F(u_{H})$, we conclude with the stability inequality (37). ∎

We notice that the snapshot equation (30) has the same structure, and the wellposedness is proved in the same way.

4.2 Convergence Analysis

We now analyze the convergence of the proposed method. The goal of this section is to estimate the difference between the snapshot solution, $u_{\text{snap}}$, computed in (30), and the multiscale coarse solution, $u_{H}$, computed in (18). To do so, we first define the following norms. We define the $V$-norm as:

[TABLE]

and $W$-norm as:

[TABLE]

We also extend them by incorporating the collision term:

[TABLE]

The total energy is now defined by:

[TABLE]

Note that we have following propositions

Proposition 2**.**

$a(u,u)=\|u\|^{2}_{{V}},$* and  $a(u,u)+l(u,u)=\|u\|^{2}_{\widetilde{V}}.$*

Proof.

This proposition simply comes from the calculations in (38). ∎

Proposition 3**.**

If $u\in V_{\text{snap}}$, we have

[TABLE]

Here $a^{j}$ and $s^{j}$ are bilinear operator defined in (33). $M=\max_{K,E}\{M_{K},M_{E}\}$ where $M_{K}$ is the number of oversampled regions $K_{j}^{+}$’s which have nonempty intersection with coarse block $K$, and $M_{E}$ is the number of oversampled regions $K_{j}^{+}$’s whose interior coarse edges $\mathcal{E}_{H}^{0}(K_{j}^{+})$ contains coarse edge $E$. They are both small numbers.

Proof.

We denote $u|_{K_{j}}$ by $u^{j}$. So $u^{j}\in V_{\text{snap}}^{j}$. According to (31), $u^{j}$ has an energy minimizing extension $\widetilde{u}^{j}\in V_{\text{snap}}^{j,+}$ that satisfies $\widetilde{u}^{j}=u^{j}$ in $K_{j}$. Then we have

[TABLE]

Combining with the definition of $\left\|\cdot\right\|^{2}_{\widetilde{W}}$, we proved (44).

Next, we denote $u^{j,+}=u|_{K_{j}^{+}}\in V_{\text{snap}}^{j,+}$. By the definition of the energy minimizing extension in (31), we have

[TABLE]

Hence,

[TABLE]

and thus we have (45). ∎

For the convergence analysis, we first examine the best approximation property. For that, we have the following:

Lemma 1**.**

Let $u_{\text{snap}}$ be the snapshot solution to the equation (30) and let $u_{H}$ be the multiscale solution to the equation (18). Then:

[TABLE]

where $C$ is a constant independent of $\epsilon$, $a^{\delta}$ and the mesh size.

Proof.

Using (30) and (18), and the fact that $V_{H}\subset V_{\text{snap}}$, we have:

[TABLE]

Then for all $w\in V_{H}$:

[TABLE]

Using Proposition 2, we have:

[TABLE]

To obtain (46), noticing that $u_{\text{snap}}-u_{H}$ and $u_{\text{snap}}-w$ are both in $V_{\text{snap}}$, it amounts to show that:

[TABLE]

In fact it suffices to show that

[TABLE]

since it is obvious that $\|u\|_{V}\leq\|u\|_{\widetilde{V}}$.

To show (50), we first use integration by parts to obtain

[TABLE]

where we have used the continuity accross fine scales $\partial\tau$, and the assumption that, in each $K_{j}$, $u$ satisfies the following equation

[TABLE]

for all $i=1,2,\cdots,m$. Here $w$ could be any function in $(V_{h})^{m}$ restricted on $K_{j}$. In particular, (51) works for $w\in V_{\text{snap}}^{j}$.

Using the definition of $a(\cdot,\cdot)$ and direct calculations, we have

[TABLE]

Then, applying the Cauchy-Schwarz inequality, we have

[TABLE]

The two terms on the right hand side are taken care of separately. To handle the first term, recall the definition of $V$-norm in equation (41), we have

[TABLE]

And to compute the second term, we notice that

[TABLE]

which in turn gives

[TABLE]

Plug (53) and (4.2) into (4.2), we have proved the desired boundedness condition (50) which concludes the proof of (46). ∎

Now, we are ready to prove our main convergence result in this section.

Theorem 2**.**

Let $u_{\text{snap}}$ be the snapshot solution to problem (30) and let $u_{H}$ be the multiscale solution to problem (18). Then:

[TABLE]

where $\Lambda_{*}=\min_{j}\lambda_{j,L_{j}+1}$, $C$ is the same constant from Lemma 1, and $M$ is the same constant from Proposition 3.

Proof.

We first denote

[TABLE]

where $\phi_{l}^{j}$ is the $l$-th multiscale basis function for the coarse element $K_{j}$ (35). Note that $\text{span}\{\phi_{l}^{j}\}$ covers the entire snapshot space. We then define a projection of $u_{\text{snap}}^{j}$ into $V_{H}^{j}$, as well as a projection of $u_{\text{snap}}$ into $V_{H}$:

[TABLE]

It is easy to see that $P^{j}(u_{\text{snap}}^{j})=P(u_{\text{snap}})|_{K_{j}}$. Combining with Proposition 3, we have

[TABLE]

Combining with Lemma 1, we proved the theorem. ∎

In the above theorem, we estimate the error between the snapshot solution $u_{\text{snap}}$ and the multiscale solution $u_{H}$. We see that the error is inversely proportional to the eigenvalues. This shows that the multiscale space gives a good approximation property in the snapshot space. In our analysis, we assume that the snapshot functions satisfy the PDE in the strong sense, that is, (51). On the other hand, there is an error between the snapshot solution $u_{\text{snap}}$ and the fine scale solution $u_{h}$ if we use Algorithm RanLocal in Section 3.2. This amounts to an irreducible error, and the analysis of this is beyond the scope of this paper.

We should emphasize that the difficulty brought by small $\delta$ is encoded in the quality of $\Lambda_{\ast}$ and thus is not explicitly expressed in the error analysis.

4.3 Small $\epsilon$ regime

An important property the algorithm satisfies is that it is robust with respect to the parameters. In the limiting regime of $\epsilon\to 0$, $\Lambda_{*}$ has a positive lower bound, and this serves as the stability argument that allows the algorithm to be effective across regimes. In particular we will show:

Theorem 3**.**

Denote $\lambda_{j,k}$ the $k$-th eigenvalue of the GEP defined in (35) for coarse element $K_{j}$. It has a asymptotic limit in the zero limit of $\epsilon$, meaning there is a constant $\lambda^{0}_{j,k}$ so that:

[TABLE]

This theorem, when combined with our main Theorem 2, indicates that the error bound, which is controlled by $\Lambda_{\ast}=\frac{1}{\text{min}_{j}\{\lambda_{j,L_{j+1}}\}}$, will not grow in $\epsilon$ and thus the error is uniformly bounded.

To show the theorem, we first start with a lemma.

Lemma 2**.**

For every coarse element $K_{j}$, we have

[TABLE]

where entries in $A^{j,0}$ and $S^{j,0}$ are defined by:

[TABLE]

and

[TABLE]

where $V_{\text{snap}}^{j}=\text{span}\{\psi_{p}^{j}\}$, and $\widetilde{\psi}_{p}^{j}$ is the basis functions’ energy minimizing extension. We further denote $\widetilde{\psi}_{p}^{j,0}$ the leading order asymptotic expansion of $\widetilde{\psi}_{p}^{j}$.

Proof.

To proceed we notice that $\widetilde{\psi}_{p}^{j}\in V_{\text{snap}}^{j,+}$ can be written as the sum of some ${\psi}^{s}$’s, where ${\psi}^{s}\in V_{\text{snap}}^{s}$ and $K_{s}\subset K_{j}^{+}$. Recall the assumption of equation (51). Then in each $K_{s}$ and for all $i=1,2,...,m$, we have:

[TABLE]

Here $w$ could be any function in $(V_{h})^{m}$ restricted on $K_{j}$.

Therefore we have:

[TABLE]

Take the asymptotic expansion for $\widetilde{\psi}_{p}^{j}$:

[TABLE]

set $w=\widetilde{\psi}_{p}^{j,0}$, and plug them back into (58). We have, in the leading order of $\frac{1}{\epsilon}$:

[TABLE]

meaning $\widetilde{\psi}_{p}^{j,0}$ is isotropic in each $K_{s}$ due to Proposition 1. Therefore $\widetilde{\psi}_{p}^{j,0}$ is isotropic. The same analysis is applied to $\widetilde{\psi}_{q}^{j,0}$.

Recall the definition of $A^{j}$ and $a^{n}$ in (33), we have:

[TABLE]

Due to Proposition 1, we have:

[TABLE]

and thus

[TABLE]

The proof for $s^{j}$ is the same and is omitted here. ∎

Theorem 3 is straightforward consequence of the following perturbation theorem:

Proof for Theorem 3.

According to Lemma 2, $A^{j}$ and $S^{j}$ have expansions $A^{j}=A^{j,0}+\mathcal{O}(\epsilon)=:A^{j,0}+\epsilon A^{j,1}$ and $S^{j}=S^{j,0}+\mathcal{O}(\epsilon)=:S^{j,0}+\epsilon S^{j,1}$. We also define $\mathbf{x}_{k}^{j,0}$ as the $k$-th generalized eigenvector of the two matrices $A^{j,0}$ and $S^{j,0}$, i.e.

[TABLE]

Using Absolute Weyl theorem for generalized eigenvalue problems in [28], when $\epsilon$ is small enough such that $\epsilon\left\|S^{j,1}\right\|_{2}<\lambda_{\min}(S^{j,0})$, we have:

[TABLE]

where $\left\|\cdot\right\|_{2}$ is the spectral norm of a matrix. ∎

According to the formula for $A^{j,0}$ and $S^{j,0}$ in (55) and (56), the eigenvalues are positive except that the smallest one is 0 with constant as corresponding eigenvector. So $\Lambda_{*}$ has positive limit in the limiting regime of $\epsilon\to 0$.

5 Numerical results

We take boundary condition $g(x,\mathbf{v})=\cos(2\pi(x_{1}+x_{2}))+1$. And we set $m=6$, and use Gaussian quadrature rule to define $\{(\mathbf{v}_{i},\alpha_{i}),i=1,...,m\}$. As for $a^{\delta}$, we give two examples. In the first example, we will choose $a^{\delta}$ to be based on a high contrast media $\kappa$, shown in Figure 3(Left), and choose $a^{\delta}$ to be an oscillatory function for the second example used in [25, 13, 23], shown in Figure 3(Right), with expression

[TABLE]

The space domain $\Omega$ is taken as the unit square $[0,1]\times[0,1]$ and is divided into $10\times 10$ coarse blocks consisting of uniform squares. Each coarse element is then divided into $10\times 10$ fine elements consisting of uniform squares. That is, $\Omega$ is partitioned by $100\times 100$ square fine elements. And we use oversampling technique in equation (26)-(29) to obtain the snapshot space. We define an oversampling region $K_{j}^{+}$ by enlarging $K_{j}$ by one coarse grid layer.

To compare the accuracy, we will use the following error quantities:

[TABLE]

where $\overline{u}$ is defined as $\overline{u}=\sum_{i=1}^{m}\alpha_{i}u_{i}.$

For Example 1, we first fix $a^{\delta}=\kappa^{4}$ and give the error tables for Knudsen number $\epsilon=10^{-1},10^{-2},10^{-3}$, respectively. And $L$ is the number of multiscale basis chosen from each coarse element, and snapshot ratio is define by

[TABLE]

From Table 2, 2, 4, we can see this framework works for all Knudsen number $\epsilon$, which verifies our proved conclusion. In addition, we see clearly the reduction of error when more basis functions are used, and the reduction of error is more rapid when fewer basis functions are used. We also observe that the method gives reasonable error levels with small snapshot ratios. On the other hand, Figures 4 show the fine and multiscale solutions with $\epsilon=10^{-2}$ and $L=5$. From these figures, we observe very good agreements between the fine-scale and multiscale solutions

Next, we fix $\epsilon=10^{-2}$ and change the high contrast value of $a^{\delta}$. We set $a^{\delta}=\kappa^{2},\;\kappa^{4},\;\kappa^{6}$, respectively. From Table 4, we can see that contrast values do not affect the error.

For Example 2, we give the error tables for $\epsilon=5\times 10^{-2},5\times 10^{-3},5\times 10^{-4}$, respectively. We present the errors for using various choices of number of basis functions in Table 6, 6, 7. We clearly see that, with a very small snapshot ratio, our method is able to obtain solutions with very good accuracy. Furthermore, we observe a faster decay of the error when smaller number of basis functions are used. In Figures 5, we present the fine and multiscale solutions with $\epsilon=5\times 10^{-3}$ and $L=5$. We observe very good agreement of both solutions.

Acknowledgements

Q. L.’s research is supported by NSF-DMS 1619778, 1750488, and UW-Madison Data Science Initiative. E. C.’s research is partially supported by Hong Kong RGC General Research Fund (Projects 14304217 and 14302018) and CUHK Direct Grant for Research 2018-19.