# Contrast-independent partially explicit time discretizations for nonlinear
multiscale problems

Eric T. Chung 111Department of Mathematics, The Chinese University of Hong
Kong (CUHK), Hong Kong SAR, Yalchin Efendiev222Department of Mathematics,
Texas A&M University, College Station, TX 77843, USA & North-Eastern Federal
University, Yakutsk, Russia, Wing Tat Leung333Department of Mathematics,
University of California, Irvine, USA, Wenyuan Li444Department of Mathematics,
Texas A&M University, College Station, TX 77843, USA

## Abstract

This work continues a line of works on developing partially explicit methods
for multiscale problems. In our previous works, we have considered linear
multiscale problems, where the spatial heterogeneities are at subgrid level
and are not resolved. In these works, we have introduced contrast-independent
partially explicit time discretizations for linear equations. The contrast-
independent partially explicit time discretization divides the spatial space
into two components: contrast dependent (fast) and contrast independent (slow)
spaces defined via multiscale space decomposition. Following this
decomposition, temporal splitting is proposed that treats fast components
implicitly and slow components explicitly. The space decomposition and
temporal splitting are chosen such that it guarantees a stability and
formulate a condition for the time stepping. This condition is formulated as a
condition on slow spaces. In this paper, we extend this approach to nonlinear
problems. We propose a splitting approach and derive a condition that
guarantees stability. This condition requires some type of contrast-
independent spaces for slow components of the solution. We present numerical
results and show that the proposed methods provide results similar to implicit
methods with the time step that is independent of the contrast.

##  1 Introduction

Nonlinear problems arise in many applications and they are typically described
by some nonlinear partial differential equations. In many applications, these
problems have multiscale nature and contain multiple scales and high contrast.
Examples include nonlinear porous media flows (e.g., Richards’ equations,
Forchheimer flow, and so on, e.g., [21, 4]), where the media properties
contain many spatial scales and high contrast. Because of high contrast in the
media properties, these processes also occur on multiple times scales. E.g.,
for nonlinear diffusion, the flow can be fast in high conductivity regions,
while the flow is slow in low conductivity regions. Because of disparity of
time scales, special temporal discretizations are often sought, which is the
main goal of the paper in the context of multiscale problems.

When the media properties are high, the flow and transport become fast and
requires small time step to resolve the dynamics. Implicit discretization can
be used to handle fast dynamics; however, this requires solving large-scale
nonlinear systems. For nonlinear problems, explicit methods are used when
possible to avoid solving nonlinear systems. The main drawback of explicit
methods is that they require small time steps that scale as the fine mesh and
depend on physical parameters, e.g., the contrast. To alleviate this issue, we
propose a novel nonlinear splitting algorithm following our earlier works [23,
24] for linear equations. The main idea of our approaches is to use multiscale
methods on a coarse spatial grid such that the time step scales as the coarse
mesh size.

Next, we give a brief overview of multiscale methods for spatial
discretizations that are used in our paper. Multiscale spatial algorithms have
been extensively studied for linear and nonlinear problems. For linear
problems, many multiscale methods have been developed. These include
homogenization-based approaches [18, 34], multiscale finite element methods
[18, 27, 33], generalized multiscale finite element methods (GMsFEM) [6, 7, 8,
11, 16], constraint energy minimizing GMsFEM (CEM-GMsFEM) [9, 10], nonlocal
multi-continua (NLMC) approaches [12], metric-based upscaling [38],
heterogeneous multiscale method [15], localized orthogonal decomposition (LOD)
[26], equation-free approaches [39, 41], multiscale stochastic approaches [29,
30, 28], and hierarchical multiscale method [5]. For high-contrast problems,
approaches such as GMsFEM and NLMC are developed. For example, in the GMsFEM
[9], multiple basis functions or continua are designed to capture the
multiscale features due to high contrast [10, 12]. These approaches require a
careful design of multiscale dominant modes. For nonlinear problems, linear
multiscale basis functions can be replaced by nonlinear maps [19, 20, 17].

Our proposed approaches follow concepts developed in [23, 24] for linear
equations. In these works, we design splitting algorithms for solving flow and
wave equations. In both cases, the solution space is divided into two parts,
coarse-grid part, and the correction part. Coarse-grid solution is computed
using multiscale basis functions with CEM-GMsFEM. The correction part uses
special spaces in the complement space (complement to the coarse space). A
careful choice of these spaces guarantees that the method is stable. Our
analysis in [23, 24] shows that for the stability, the correction space should
be free of contrast and thus, this requires a special multiscale space
construction. These splitting algorithms take their origin from earlier works
[36, 43]. In this paper, we extend the linear concepts to nonlinear problems.

Splitting algorithms for nonlinear problems have often been used in the
literature. Earlier approaches include implicit-explicit approaches and other
techniques [3, 35, 1, 22, 2, 37, 42, 14, 25, 32, 40, 31, 13]. In many
approaches, nonlinear contributions are roughly divided into two parts
depending on whether it is easy to implicitly solve discretized system. For
easy to solve part, implicit discretization is used, while for the rest,
explicit discretization is used. However, in general, one can not separate
these parts for the problems under the consideration. Our goal is to use
splitting concepts and treat implicitly and explicitly some parts of the
solution. As a result, we can use larger time steps that scale as the coarse
mesh size.

Our approach starts a nonlinear dynamical system

$u_{t}+f(u)+g(u)=0,$

where $f(u)$ represents diffusion-like operator, while $g(u)$ represents
reaction-like terms. In linear problems, for the stability, we formulate a
condition that involves the time step, the energy and $L^{2}$ norm of the
solution in complement space. This is a constraint for the time step. With an
appropriate choice of the complement space, this condition guarantees the
stability for the time steps that scale as the coarse mesh size. To obtain
similar conditions for nonlinear problems, we carry out the analysis for
nonlinear $f(u)$ and $g(u)$ functions. The analysis reveals conditions that
are required for stability. The conditions of multiscale spaces turn out to
share some similarities to those for nonlinear multiscale methods [17]. We
remark several observations.

  * •

Additional degrees of freedom are needed for dynamic problems, in general, to
handle missing information.

  * •

We note that restrictive time step scales as the coarse mesh size and, thus,
much coarser.

We present several numerical results. In our numerical results, we consider
two cases. In the first case, we take the reaction term, $g(u)$, to be
nonlinear, while the diffusion term, $f(u)$, to be linear. In the second case,
we consider both to be nonlinear and use for $f(u)$ the form
$f(u)=-div(\kappa(x,u)\nabla u)$. The media properties for the diffusion is
taken to be heterogeneous and we choose smooth and singular source terms. We
compare the proposed approach to the approach, where all degrees of freedom
are handled implicitly. We show that the proposed methods provide an
approximation similar to fully implicit methods.

The paper is organized as follows. In the next section, we provide some
preliminaries. In Section 3, we present our main assumptions and stability
estimates for fine-grid problem. In Section 4, we describe the partially
explicit method and its stability, following some discussions in Section 5.
Numerical results are presented in Section 6.

##  2 Problem Setting

We consider the following equation

$\displaystyle u_{t}=-f(u)-g(u),$ (2.1)

where $f=\cfrac{\delta F}{\delta u}$ and $g=\cfrac{\delta G}{\delta u}$ which
are the variational derivative of energies $F(u):=\int_{\Omega}E_{1}(u)$ and
$G(u):=\int_{\Omega}E_{2}(u)$. Here, $f$ is assumed to be contrast dependent
nonlinear (or linear) (i.e., $f$ introduces stiffness in the system) while $g$
is contrast independent (i.e., $g$ does not introduce stiffness).

We assume $f(u)\in V^{*}$ and $g(u)\in L^{2}(\Omega)$ for all $u\in V$. We can
then consider the weak formulation of the problem, namely, finding $u\in V$
such that

$(u_{t},v)=-(f(u),v)_{V^{*},V}\;\forall v\in V.$

To simplify the notation, we will simply write $(\cdot,\cdot)$ instead of
$(\cdot,\cdot)_{V^{*},V}$ in the following discussion.

Example 1: For $F(u)=\cfrac{1}{2}\int_{\Omega}\kappa|\nabla u|^{2}$ and
$G(u)=0$, we will have $f(u^{\prime})=\cfrac{\delta F}{\delta u}\in
H^{1}(\Omega)=(H^{1}(\Omega))^{*}$ for $V=H^{1}(\Omega)$.

$(\cfrac{\delta F}{\delta u}(u^{\prime}),v)=\int_{\Omega}\kappa\nabla
u^{\prime}\cdot\nabla v,$

and thus, we have

$-f(u^{\prime})=\nabla\cdot(\kappa\nabla u^{\prime}).$

In strong form, we have

$u_{t}=\nabla\cdot(\kappa\nabla u).$

Example 2: For $F=\cfrac{1}{p}\int\kappa|\nabla u|^{p}$ and $G(u)=0$, we will
have $f(u^{\prime})=\cfrac{\delta F}{\delta u}\in
W^{1,\frac{p}{p-1}}(\Omega)=(W^{1,p}(\Omega))^{*}$ for $V=W^{1,p}(\Omega)$

$(\cfrac{\delta F}{\delta u}(u^{\prime}),v)=\int\kappa|\nabla
u^{\prime}|^{p-2}\nabla u^{\prime}\cdot\nabla v,$

and thus, we have

$-f(u^{\prime})=\nabla\cdot(\kappa|\nabla u^{\prime}|^{p-2}\nabla
u^{\prime}).$

In strong form, we have

$u_{t}=\nabla\cdot(\kappa|\nabla u|^{p-2}\nabla u).$

In the following discussion, we will make the following assumptions on the
second variational derivative of $F$ and $G$.

  * •

The second variational derivative $\delta^{2}F$ and $\delta^{2}G$ satisfy

$\delta^{2}F(u)(v,v)\geq c(u)\|v\|_{V}^{2}\;\forall u,v\in V$
$\delta^{2}G(u)(v,v)\geq b(u)\|v\|^{2}\;\forall u,v\in V,$

where $0\leq c(u)<\infty$ and $-\underline{b}\leq b(u)<\infty$ are independent
on $v$.

  * •

The second variational derivative $\delta^{2}F$ and $\delta^{2}G$ are bounded.
That is,

$|\delta^{2}F(u)(w,v)|\leq C(u)\|v\|_{V}\|w\|_{V}\;\forall u,v,w\in V$
$|\delta^{2}G(u)(w,v)|\leq B\|v\|_{L_{2}}\|w\|_{L_{2}}\;\forall u,v,w\in V,$

where $0<C(u)<\infty$ and $0<B<\infty$ are independent on $v,w$.

##  3 Discretization

To solve the problem, a standard method is finite element approach. We can
consider the numerical solution $u_{H}(t,\cdot)\in V_{H}$ satisfies

$\displaystyle(u_{H,t},v)=-(f(u)+g(u),v)\;\forall v\in V_{H},$ (3.1)

where $V_{H}$ is a finite element space in $V$.

For the time discretization, we can consider two simplest discretizations
which are forward Euler and backward Euler methods. For the forward Euler
method, we consider $\\{u_{H}^{k}\\}_{k=0}^{N}\subset V_{H}$ such that

$(\cfrac{u_{H}^{n+1}-u_{H}^{n}}{\Delta
t},v)+(f(u_{H}^{n})+g(u_{H}^{n}),v)=0\;\forall v\in V_{H}.$

For the backward Euler method, we consider $\\{u_{H}^{k}\\}_{k=0}^{N}\subset
V_{H}$ such that

$(\cfrac{u_{H}^{n+1}-u_{H}^{n}}{\Delta
t},v)+(f(u_{H}^{n+1})+g(u_{H}^{n+1}),v)=0\;\forall v\in V_{H}.$

Next, we would like to derive stability conditions for backward and forward
Euler methods.

Since

$F(u_{H}^{n})=F(u_{H}^{n+1})-(f(u_{H}^{n+1}),u_{H}^{n+1}-u_{H}^{n})+\cfrac{1}{2}\delta^{2}F(\xi_{1}^{n})(u_{H}^{n+1}-u_{H}^{n},u_{H}^{n+1}-u_{H}^{n})$

and

$G(u_{H}^{n})=G(u_{H}^{n+1})-(g(u_{H}^{n+1}),u_{H}^{n+1}-u_{H}^{n})+\cfrac{1}{2}\delta^{2}G(\xi_{2}^{n})(u_{H}^{n+1}-u_{H}^{n},u_{H}^{n+1}-u_{H}^{n})$

for some $\xi_{i}^{n}=(1-\lambda_{i})u_{H}^{n+1}+\lambda_{i}u_{H}^{n}$ with
$\lambda_{i}\in(0,1)$ and $i=1,2$, we have

$\displaystyle 0=$ $\displaystyle(\cfrac{u_{H}^{n+1}-u_{H}^{n}}{\Delta
t},u_{H}^{n+1}-u_{H}^{n})+(f(u_{H}^{n+1})+g(u_{H}^{n+1}),u_{H}^{n+1}-u_{H}^{n})$
$\displaystyle=$ $\displaystyle\cfrac{1}{\Delta
t}\|u_{H}^{n+1}-u_{H}^{n}\|^{2}+F(u_{H}^{n+1})-F(u_{H}^{n})+G(u_{H}^{n+1})-G(u_{H}^{n})$
$\displaystyle+\cfrac{1}{2}\delta^{2}F(\xi_{1}^{n})(u_{H}^{n+1}-u_{H}^{n},u_{H}^{n+1}-u_{H}^{n})+\cfrac{1}{2}\delta^{2}G(\xi_{2}^{n})(u_{H}^{n+1}-u_{H}^{n},u_{H}^{n+1}-u_{H}^{n})$
$\displaystyle\geq$ $\displaystyle\Big{(}\cfrac{1}{\Delta
t}+b(u)\Big{)}\|u_{H}^{n+1}-u_{H}^{n}\|^{2}+F(u_{H}^{n+1})-F(u_{H}^{n})+G(u_{H}^{n+1})-G(u_{H}^{n})+\cfrac{c(u)}{2}\|u_{H}^{n+1}-u_{H}^{n}\|_{V}^{2}.$

We have

$\displaystyle F(u_{H}^{n+1})+G(u_{H}^{n+1})$ $\displaystyle\leq
F(u_{H}^{n+1})+G(u_{H}^{n+1})+\cfrac{c(u)}{2}\|u_{H}^{n+1}-u_{H}^{n}\|_{V}^{2}+\Big{(}\cfrac{1}{\Delta
t}-\underline{b}\Big{)}\|u_{H}^{n+1}-u_{H}^{n}\|^{2}$ $\displaystyle\leq
F(u_{H}^{n})+G(u_{H}^{n})$

for any $\Delta t$ and thus, backward Euler method is stable if $\Delta
t\underline{b}\leq 1$.

Similarly, for forward Euler method, we can use

$F(u_{H}^{n+1})=F(u_{H}^{n})+(f(u_{H}^{n}),u_{H}^{n+1}-u_{H}^{n})+\cfrac{1}{2}\delta^{2}F(\xi_{1}^{n})(u_{H}^{n+1}-u_{H}^{n},u_{H}^{n+1}-u_{H}^{n})$
$G(u_{H}^{n+1})=G(u_{H}^{n})+(g(u_{H}^{n}),u_{H}^{n+1}-u_{H}^{n})+\cfrac{1}{2}\delta^{2}G(\xi_{2}^{n})(u_{H}^{n+1}-u_{H}^{n},u_{H}^{n+1}-u_{H}^{n})$

and obtain

$\displaystyle 0=$ $\displaystyle(\cfrac{u_{H}^{n+1}-u_{H}^{n}}{\Delta
t},u_{H}^{n+1}-u_{H}^{n})+(f(u^{n})+g(u_{H}^{n}),u_{H}^{n+1}-u_{H}^{n})$
$\displaystyle=$ $\displaystyle\cfrac{1}{\Delta
t}\|u_{H}^{n+1}-u_{H}^{n}\|^{2}+F(u_{H}^{n+1})-F(u_{H}^{n})+G(u_{H}^{n+1})-G(u_{H}^{n})$
$\displaystyle-\cfrac{1}{2}\delta^{2}F(\xi_{1}^{n})(u_{H}^{n+1}-u_{H}^{n},u_{H}^{n+1}-u_{H}^{n})-\cfrac{1}{2}\delta^{2}G(\xi_{2}^{n})(u_{H}^{n+1}-u_{H}^{n},u_{H}^{n+1}-u_{H}^{n})$
$\displaystyle\geq$ $\displaystyle(\cfrac{1}{\Delta
t}-B)\|u_{H}^{n+1}-u_{H}^{n}\|^{2}+F(u_{H}^{n+1})-F(u_{H}^{n})+G(u_{H}^{n+1})-G(u_{H}^{n})-\cfrac{C(\xi^{n})}{2}\|u_{H}^{n+1}-u_{H}^{n}\|_{V}^{2}.$

Therefore, if $\Delta
t\Big{(}\cfrac{C(\xi)}{2}\cfrac{\|u_{H}^{n+1}-u_{H}^{n}\|_{V}^{2}}{\|u_{H}^{n+1}-u_{H}^{n}\|^{2}}+B\Big{)}\leq
1$ for any $\xi=(1-\lambda)u_{H}^{k+1}+\lambda u_{H}^{k}$ with $0\leq k\leq
N-1$, we have

$\displaystyle F(u_{H}^{n+1})+G(u_{H}^{n+1})$ $\displaystyle\leq
F(u_{H}^{n})+G(u_{H}^{n}).$

We can see that although forward Euler method is easier for implementation, we
require a small time step for stability if $\sup_{v\in
V_{H}}\cfrac{\|v\|_{V}^{2}}{\|v\|^{2}}$ or $C(\xi)$ is large.

We remark that in typical cases we will consider $\underline{b}$ and $B$ are
not too large. Therefore, we have $\Delta t\underline{b}$ and $\Delta tB$ is
small and the energy $G$ will not affect the stability too much.

##  4 Partially explicit scheme with space splitting

To obtain an efficient method, one can consider partially explicit scheme by
splitting finite element space. We consider $V_{H}$ is a direct sum of two
subspace $V_{H,1}$ and $V_{H,2}$, namely, $V_{H}=V_{H,1}\oplus V_{H,2}$. The
finite element solution is then satisfying

$\displaystyle(u_{H,1,t}+u_{H,2,t},v_{1})+(f(u_{H,1}+u_{H,2})+g(u_{H,1}+u_{H,2}),v_{1})$
$\displaystyle=0\;\forall v_{1}\in V_{H,1},$
$\displaystyle(u_{H,1,t}+u_{H,2,t},v_{2})+(f(u_{H,1}+u_{H,2})+g(u_{H,1}+u_{H,2}),v_{2})$
$\displaystyle=0\;\forall v_{2}\in V_{H,2},$

where $u_{H}=u_{H,1}+u_{H,2}$. We can use a partially explicit time
discretization. For example, we can consider

$\displaystyle(\cfrac{u_{H,1}^{n+1}-u_{H,1}^{n}}{\Delta
t}+\cfrac{u_{H,2}^{n}-u_{H,2}^{n-1}}{\Delta
t},v_{1})+(f(u_{H,1}^{n+1}+u_{H,2}^{n})+g(u_{H,1}^{n}+u_{H,2}^{n}),v_{1})$
$\displaystyle=0\;\forall v_{1}\in V_{H,1},$
$\displaystyle(\cfrac{u_{H,1}^{n}-u_{H,1}^{n-1}}{\Delta
t}+\cfrac{u_{H,2}^{n+1}-u_{H,2}^{n}}{\Delta
t},v_{2})+(f(u_{H,1}^{n+1}+u_{H,2}^{n})+g(u_{H,1}^{n}+u_{H,2}^{n}),v_{2})$
$\displaystyle=0\;\forall v_{2}\in V_{H,2}.$

### Energy stability

######  Lemma 4.1.

If

$\begin{split}(f(u_{H}^{n+1})-f(u_{H,1}^{n+1}+u_{H,2}^{n}),u_{H}^{n+1}-u_{H}^{n})\leq\\\
\cfrac{\bar{c}}{2}\|u_{H}^{n+1}-u_{H}^{n}\|_{V}^{2}+\Big{(}\cfrac{(1-\gamma)}{\Delta
t}-(1+\gamma)\cfrac{B}{2}\Big{)}\sum_{i}\|u_{H,i}^{n+1}-u_{H,i}^{n}\|^{2},\end{split}$
(4.1)

where $\bar{c}=\inf_{u\in V_{H}}c(u)$ and $\gamma=\sup_{v_{1}\in
V_{H,1},v_{2}\in V_{H,2}}\cfrac{(v_{1},v_{2})}{\|v_{1}\|\|v_{2}\|}<1$, we have

$\cfrac{\gamma}{2\Delta
t}\sum_{i}\|u_{H,i}^{n+1}-u_{H,i}^{n}\|^{2}+F(u_{H}^{n+1})+G(u_{H}^{n+1})\leq\cfrac{\gamma}{2\Delta
t}\sum_{i}\|u_{H,i}^{n}-u_{H,i}^{n-1}\|^{2}+F(u_{H}^{n})+G(u_{H}^{n}).$

###### Proof.

By substitute $v_{1}=u_{H,1}^{n+1}-u_{H,1}^{n}$ and
$v_{2}=u_{H,2}^{n+1}-u_{H,2}^{n}$, we have

$\displaystyle\cfrac{1}{\Delta
t}\|u_{H,1}^{n+1}-u_{H,1}^{n}\|^{2}+\cfrac{1}{\Delta
t}(u_{H,2}^{n}-u_{H,2}^{n-1},u_{H,1}^{n+1}-u_{H,1}^{n})$
$\displaystyle+(f(u_{H,1}^{n+1}+u_{H,2}^{n})+g(u_{H}^{n}),u_{H,1}^{n+1}-u_{H,1}^{n})$
$\displaystyle=0,$

and

$\displaystyle\cfrac{1}{\Delta
t}\|u_{H,2}^{n+1}-u_{H,2}^{n}\|^{2}+\cfrac{1}{\Delta
t}(u_{H,1}^{n}-u_{H,1}^{n-1},u_{H,2}^{n+1}-u_{H,2}^{n})$
$\displaystyle+(f(u_{H,1}^{n+1}+u_{H,2}^{n})+g(u_{H}^{n}),u_{H,2}^{n+1}-u_{H,2}^{n})$
$\displaystyle=0.$

Summing up the above two equations, we have

$\displaystyle\cfrac{1}{\Delta
t}\sum_{i}\|u_{H,i}^{n+1}-u_{H,i}^{n}\|^{2}+\cfrac{1}{\Delta t}\sum_{i\neq
j}(u_{H,i}^{n}-u_{H,i}^{n-1},u_{H,j}^{n+1}-u_{H,j}^{n})$
$\displaystyle+(f(u_{H,1}^{n+1}+u_{H,2}^{n})+g(u_{H}^{n}),u_{H}^{n+1}-u_{H}^{n})$
$\displaystyle=0.$

We first use

$\displaystyle\cfrac{1}{\Delta t}|\sum_{i\neq
j}(u_{H,i}^{n}-u_{H,i}^{n-1},u_{H,j}^{n+1}-u_{H,j}^{n})|$
$\displaystyle\leq\cfrac{\gamma}{\Delta t}\sum_{i\neq
j}\|u_{H,i}^{n}-u_{H,i}^{n-1}\|\|u_{H,j}^{n+1}-u_{H,j}^{n}\|$
$\displaystyle\leq\cfrac{\gamma}{2\Delta
t}\sum_{i}\Big{(}\|u_{H,i}^{n+1}-u_{H,i}^{n}\|^{2}+\|u_{H,i}^{n}-u_{H,i}^{n-1}\|^{2}\Big{)}$

and obtain

$\displaystyle\cfrac{1}{\Delta
t}\sum_{i}\|u_{H,i}^{n+1}-u_{H,i}^{n}\|^{2}+\cfrac{1}{\Delta t}\sum_{i\neq
j}(u_{H,i}^{n}-u_{H,i}^{n-1},u_{H,j}^{n+1}-u_{H,j}^{n})$ $\displaystyle\geq$
$\displaystyle\cfrac{2-\gamma}{2\Delta
t}\sum_{i}\|u_{H,i}^{n+1}-u_{H,i}^{n}\|^{2}-\cfrac{\gamma}{2\Delta
t}\sum_{i}\|u_{H,i}^{n}-u_{H,i}^{n-1}\|^{2}.$

To prove the stability of the method, we can consider

$F(u_{H}^{n})=F(u_{H}^{n+1})-(f(u_{H}^{n+1}),u_{H}^{n+1}-u_{H}^{n})+\cfrac{1}{2}\delta^{2}F(\xi_{1}^{n})(u_{H}^{n+1}-u_{H}^{n},u_{H}^{n+1}-u_{H}^{n})$
$G(u_{H}^{n+1})=G(u_{H}^{n})+(g(u_{H}^{n}),u_{H}^{n+1}-u_{H}^{n})+\cfrac{1}{2}\delta^{2}G(\xi_{2}^{n})(u_{H}^{n+1}-u_{H}^{n},u_{H}^{n+1}-u_{H}^{n})$

for some $\xi_{i}^{n}=(1-\lambda_{i})u_{H}^{n+1}+\lambda_{i}u_{H}^{n}$ with
$\lambda_{i}\in(0,1)$ and $i=1,2$.

Therefore, we have

$\begin{split}&(f(u_{H,1}^{n+1}+u_{H,2}^{n}),u_{H}^{n+1}-u_{H}^{n})\\\
=&(f(u_{H,1}^{n+1}+u_{H,2}^{n})-f(u_{H}^{n+1}),u_{H}^{n+1}-u_{H}^{n})\\\
&+F(u_{H}^{n+1})-F(u_{H}^{n})+\cfrac{1}{2}\delta^{2}F(\xi_{1}^{n})(u_{H}^{n+1}-u_{H}^{n},u_{H}^{n+1}-u_{H}^{n})\end{split}$
(4.2)

and

$\displaystyle(g(u_{H}^{n}),u_{H}^{n+1}-u_{H}^{n})$ $\displaystyle=$
$\displaystyle
G(u_{H}^{n+1})-G(u_{H}^{n})-\cfrac{1}{2}\delta^{2}G(\xi_{2}^{n})(u_{H}^{n+1}-u_{H}^{n},u_{H}^{n+1}-u_{H}^{n}).$

Thus, we obtain

$\displaystyle\cfrac{\gamma}{2\Delta
t}\sum_{i}\|u_{H,i}^{n+1}-u_{H,i}^{n}\|^{2}+\cfrac{(1-\gamma)}{\Delta
t}\sum_{i}\|u_{H,i}^{n+1}-u_{H,i}^{n}\|^{2}+F(u_{H}^{n+1})+G(u_{H}^{n+1})+\cfrac{c(\xi^{n})}{2}\|u_{H}^{n+1}-u_{H}^{n}\|_{V}^{2}$
$\displaystyle\leq$ $\displaystyle\cfrac{\gamma}{2\Delta
t}\sum_{i}\|u_{H,i}^{n}-u_{H,i}^{n-1}\|^{2}+F(u_{H}^{n})+G(u_{H}^{n})+\cfrac{B}{2}\|u_{H}^{n+1}-u_{H}^{n}\|^{2}$
$\displaystyle+(f(u_{H}^{n+1})-f(u_{H,1}^{n+1}+u_{H,2}^{n}),u_{H}^{n+1}-u_{H}^{n})$

and

$\cfrac{B}{2}\|u_{H}^{n+1}-u_{H}^{n}\|^{2}\leq(1+\gamma)\cfrac{B}{2}\sum_{i}\|u_{H,i}^{n+1}-u_{H,i}^{n}\|^{2}.$

If

$(f(u_{H}^{n+1})-f(u_{H,1}^{n+1}+u_{H,2}^{n}),u_{H}^{n+1}-u_{H}^{n})\leq\cfrac{c(\xi^{n})}{2}\|u_{H}^{n+1}-u_{H}^{n}\|_{V}^{2}+\Big{(}\cfrac{(1-\gamma)}{\Delta
t}-(1+\gamma)\cfrac{B}{2}\Big{)}\sum_{i}\|u_{H,i}^{n+1}-u_{H,i}^{n}\|^{2}$

then we have

$\cfrac{\gamma}{2\Delta
t}\sum_{i}\|u_{H,i}^{n+1}-u_{H,i}^{n}\|^{2}+F(u_{H}^{n+1})+G(u_{H}^{n+1})\leq\cfrac{\gamma}{2\Delta
t}\sum_{i}\|u_{H,i}^{n}-u_{H,i}^{n-1}\|^{2}+F(u_{H}^{n})+G(u_{H}^{n}).$

∎

######  Lemma 4.2.

If

$\displaystyle\cfrac{\bar{C}_{2}^{2}}{2\bar{c}}\sup_{v_{2}\in
V_{H,2}}\cfrac{\|v_{2}\|_{V}^{2}}{\|v_{2}\|^{2}}+(1+\gamma)\cfrac{B}{2}\leq\cfrac{(1-\gamma)}{\Delta
t},$ (4.3)

where $\bar{c}=\inf_{u\in V_{H}}c(u)$, $\bar{C}_{2}=\sup_{\xi\in
V_{H}}C_{2}(\xi)$ and

$C_{2}(\xi)=\sup_{v\in V_{H},w\in
V_{H,2}}\cfrac{1}{\|v\|_{V}\|w\|_{V}}\delta^{2}F(\xi)(w,v)\leq C(\xi),$

we have

$\cfrac{\gamma}{2\Delta
t}\sum_{i}\|u_{H,i}^{n+1}-u_{H,i}^{n}\|^{2}+F(u_{H}^{n+1})+G(u_{H}^{n+1})\leq\cfrac{\gamma}{2\Delta
t}\sum_{i}\|u_{H,i}^{n}-u_{H,i}^{n-1}\|^{2}+F(u_{H}^{n})+G(u_{H}^{n}).$

###### Proof.

For the proof of this lemma, we will show that if the condition of Lemma 4.2
holds, then the condition of Lemma 4.1 holds. For this reason, we will need to
estimate
$(f(u_{H}^{n+1})-(f(u_{H,1}^{n+1}+u_{H,2}^{n}),u_{H}^{n+1}-u_{H}^{n})$.
Similar to the proof in previous lemma, we consider

$(f(u_{H}^{n+1}),u_{H}^{n+1}-u_{H}^{n})=(f(u_{H,1}^{n+1}+u_{H,2}^{n}),u_{H}^{n+1}-u_{H}^{n})+(\delta^{2}F(\tilde{\xi}^{n})(u_{H}^{n+1}-u_{H}^{n}),u_{H,2}^{n+1}-u_{H,2}^{n})$

for some
$\tilde{\xi}^{n}=(1-\tilde{\lambda})u_{H}^{n+1}+\tilde{\lambda}(u_{H,1}^{n+1}+u_{H,2}^{n})$
with $\tilde{\lambda}\in(0,1)$. We have

$\displaystyle(f(u_{H}^{n+1})-f(u_{H,1}^{n+1}+u_{H,2}^{n}),u_{H}^{n+1}-u_{H}^{n})$
$\displaystyle=(\delta^{2}F(\tilde{\xi}^{n})(u_{H}^{n+1}-u_{H}^{n}),u_{H,2}^{n+1}-u_{H,2}^{n})$
$\displaystyle\leq
C_{2}(\tilde{\xi}^{n})\|u_{H}^{n+1}-u_{H}^{n}\|_{V}\|u_{H,2}^{n+1}-u_{H,2}^{n}\|_{V}.$

Since

$C_{2}(\tilde{\xi}^{n})\|u_{H}^{n+1}-u_{H}^{n}\|_{V}\|u_{H,2}^{n+1}-u_{H,2}^{n}\|_{V}\leq\cfrac{\bar{c}}{2}\|u_{H}^{n+1}-u_{H}^{n}\|_{V}^{2}+\cfrac{C_{2}^{2}(\tilde{\xi}^{n})}{2\bar{c}}\|u_{H,2}^{n+1}-u_{H,2}^{n}\|_{V}^{2},$

we have

$(f(u_{H}^{n+1})-f(u_{H,1}^{n+1}+u_{H,2}^{n}),u_{H}^{n+1}-u_{H}^{n})\leq\cfrac{\bar{c}}{2}\|u_{H}^{n+1}-u_{H}^{n}\|_{V}^{2}+\cfrac{C_{2}^{2}(\tilde{\xi}^{n})}{2\bar{c}}\|u_{H,2}^{n+1}-u_{H,2}^{n}\|_{V}^{2}.$

If

$\cfrac{\bar{C}_{2}^{2}}{2\bar{c}}\cfrac{\|u_{H,2}^{n+1}-u_{H,2}^{n}\|_{V}^{2}}{\|u_{H,2}^{n+1}-u_{H,2}^{n}\|^{2}}+(1+\gamma)\cfrac{B}{2}\leq\cfrac{(1-\gamma)}{\Delta
t},$

we have the condition formulated in Lemma 4.1

$\begin{split}(f(u_{H}^{n+1})-f(u_{H,1}^{n+1}+u_{H,2}^{n}),u_{H}^{n+1}-u_{H}^{n})\leq\\\
\cfrac{\bar{c}}{2}\|u_{H}^{n+1}-u_{H}^{n}\|_{V}^{2}+\Big{(}\cfrac{(1-\gamma)}{\Delta
t}-(1+\gamma)\cfrac{B}{2}\Big{)}\sum_{i}\|u_{H,i}^{n+1}-u_{H,i}^{n}\|^{2}.\end{split}$
(4.4)

By Lemma 4.1, we get the result. ∎

Example 1. For $F=\cfrac{1}{2}\int_{\Omega}\kappa|\nabla u|^{2}$ and $G(u)=0$,
we have

$(\cfrac{\delta F}{\delta u},v)=\int_{\Omega}\kappa\nabla u\cdot\nabla v,$

and

$\delta^{2}F(u)(w,v)=\int_{\Omega}\kappa\nabla v\cdot\nabla w\;\forall u\in
V,$

and thus, we have

$\bar{c}=C_{2}^{2}=1,\;B=0$

and the partially explicit scheme is stable when

$\cfrac{\Delta t}{2}\sup_{v_{2}\in V_{H,2}}\cfrac{\|\kappa^{\frac{1}{2}}\nabla
v_{2}\|^{2}}{\|v_{2}\|^{2}}\leq(1-\gamma)\;\forall v_{2}\in V_{H,2}.$

##  5 Discussions

###  5.1 $G=0$ case

First, we present some discussions for $G=0$ case. In this case, the first
stability condition for partial explicit scheme is

$\begin{split}(f(u_{H,1}^{n+1}+u_{H,2}^{n+1})-f(u_{H,1}^{n+1}+u_{H,2}^{n}),u_{H}^{n+1}-u_{H}^{n})\leq\\\
\cfrac{\bar{c}}{2}\|u_{H}^{n+1}-u_{H}^{n}\|_{V}^{2}+\cfrac{(1-\gamma)}{\Delta
t}\sum_{i}\|u_{H,i}^{n+1}-u_{H,i}^{n}\|^{2}.\end{split}$ (5.1)

This condition can be understood as a nonlinear constraint on the “second
space” that represents $u_{H,2}^{n}$ and in order to have a small bound, one
needs to guarantee that $u_{H,1}^{n}$ captures important degrees of freeom.
Indeed, the smallness of

${(f(u_{H,1}^{n+1}+u_{H,2}^{n+1})-f(u_{H,1}^{n+1}+u_{H,2}^{n}),u_{H}^{n+1}-u_{H}^{n})\over\|u_{H}^{n+1}-u_{H}^{n}\|_{V}^{2}}$

is a condition on $u_{H,1}^{n}$ (on the coarse space) and requires that this
term is chosen such that the difference is independent of the contrast. This
condition is more evident in Lemma 4.2, where the condition on $V_{2}$ is

$\cfrac{\bar{C}_{2}^{2}}{2\bar{c}}\sup_{v_{2}\in
V_{H,2}}\cfrac{\|v_{2}\|_{V}^{2}}{\|v_{2}\|^{2}}\leq\cfrac{(1-\gamma)}{\Delta
t}.$

###  5.2 $G\not=0$ case.

In this case, we will first treat the nonlinear forcing explicitly and we will
also discuss the case when $g(u)$ is partially explicit.

##  6 Numerical Results

In this section, we will present numerical results for various cases. We will
consider several choices for $f(u)$ and $g(u)$. For $f(u)$, we will use
diffusion operator for linear case

$f(u)=-\nabla\cdot(\kappa\nabla u),$

and nonlinear case

$f(u)=-\nabla\cdot(\kappa\alpha(u)\nabla u).$ (6.1)

In all examples, we will use two heterogeneous high contrast $\kappa(x)$ that
represent the media, where one is more complex (more channels). As for $g(u)$,
we will consider several choices of nonlinear reaction terms as discussed
below. This term will contain nonlinear reaction and steady state spatial
source term. One source term will be more regular and the other more singular.
The singular source term is chosen so that CEM solution requires additional
basis functions as the source term contains subgrid features. In all numerical
examples, the coarse mesh size is $\frac{1}{10}$ and the fine mesh size is
$\frac{1}{100}$. For the time discretization, we will consider the final time
$T=0.05$. In our numerical tests, we will compare three methods.

  * •

First, we will use implicit CEM to compute the solution without additional
degrees of freedom (called “Implicit CEM” in our graphs).

  * •

Secondly, we will compute the solution with additional degrees of freedom
using implicit CEM (called “Implicit CEM with additional basis” in our
graphs).

  * •

Finally, we will compute the solution with additional degrees of fredoom using
our proposed partially explicit approach (called “Partially Explicit Splitting
CEM” in our graphs).

In all examples, we use Newton or Picard iterations to find the solution of
nonlinear equations. In all examples, we observe that our proposed partially
explicit method provides similar accuracy as the implicit CEM approach that
uses additional degrees of freedom.

###  6.1 $V_{H,1}$ and $V_{H,2}$ constructions

In this section, we present a way to construct the spaces satisfying (4.3)
based on linear problems. These spaces are constructed under the assumption
that linear multiscale structure can be used to accurately model coarse-grid
solution. It can be shown that the linear spaces satisfy (4.3) under some
assumptions on $\kappa(x,u)$ (see (6.1)). Here, we follow our previous work
[23]. As the constrained energy minimization basis functions are constructed
such that they are almost orthogonal to a space $\tilde{V}$, the CEM finite
element space is a good option for $V_{H,1}$. To find a $V_{H,2}$ satisfying
the condition (4.3), we can use an eigenvalue problem to construct the local
basis functions. We will first introduce the CEM finite element space,
followed by the discussion of constructing $V_{H,2}$. In the following, we let
$V(S)=H_{0}^{1}(S)$ for a proper subset $S\subset\Omega$.

####  6.1.1 CEM method

In this section, we introduce the CEM method for solving the problem (3.1). We
will construct the finite element space by solving a constrained energy
minimization problem. Let $\mathcal{T}_{H}$ be a coarse grid partition of
$\Omega$. For $K_{i}\in\mathcal{T}_{H}$, we first need to define a set of
auxiliary basis functions in $V(K_{i})$. We solve

$\displaystyle\int_{K_{i}}\kappa\nabla\psi_{j}^{(i)}\cdot\nabla
v=\lambda_{j}^{(i)}s_{i}(\psi_{j}^{(i)},v)\;\forall v\in V(K_{i}),$

where

$\displaystyle s_{i}(u,v)=\int_{K_{i}}\tilde{\kappa}uv,\;\tilde{\kappa}=\kappa
H^{-2}\;\text{or}\;\tilde{\kappa}=\kappa\sum_{i}\left|\nabla\chi_{i}\right|^{2}$

with $\\{\chi_{i}\\}$ being a partition of unity functions corresponding to an
overlapping partition of the domain. We then collect the first $L_{i}$
eigenfunctions corresponding to the first $L_{i}$ smallest eigenvalues. We
define

$V_{aux}^{(i)}:=\text{span}\\{\psi_{j}^{(i)}:\;1\leq j\leq L_{i}\\}.$

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

$s(\Pi u,v)=s(u,v)\;\forall v\in V_{aux}:=\sum_{i=1}^{N_{e}}V_{aux}^{(i)},$

where $s(u,v):=\sum_{i=1}^{N_{e}}s_{i}(u|_{K_{i}},v|_{K_{i}})$ and $N_{e}$ is
the number of coarse elements. We let $K_{i}^{+}$ be an oversampling domain of
$K_{i}$, which is a few coarse blocks larger than $K_{i}$ [9]. For each
auxiliary basis functions $\psi_{j}^{(i)}$, we can find a local basis function
$\phi_{j}^{(i)}\in V(K_{i}^{+})$ such that

$\displaystyle a(\phi_{j}^{(i)},v)+s(\mu_{j}^{(i)},v)$
$\displaystyle=0\;\forall v\in V(K_{i}^{+}),$ $\displaystyle
s(\phi_{j}^{(i)},\nu)$ $\displaystyle=s(\psi_{j}^{(i)},\nu)\;\forall\nu\in
V_{aux}(K_{i}^{+})$

for some $\mu_{j}^{(i)}\in V_{aux}$. We then define the space $V_{cem}$ as

$\displaystyle V_{cem}$ $\displaystyle:=\text{span}\\{\phi_{j}^{(i)}:\;1\leq
i\leq N_{e},1\leq j\leq L_{i}\\}.$

The CEM solution $u_{cem}$ is given by

$\displaystyle(\cfrac{\partial u_{cem}}{\partial
t},v)+(f(u_{cem})+g(u_{cem}),v)=0\;\forall v\in V_{cem}.$

Let $\tilde{V}:=\\{v\in V:\;\Pi(v)=0\\}$ and we can now construct $V_{H,2}$.

####  6.1.2 Construction of $V_{H,2}$

The construction of $V_{H,2}$ is based on the CEM type finite element space.
For each coarse element $K_{i}$, we will solve an eigenvalue problem to get
the second type of auxiliary basis. We obtain eigenpairs
$(\xi_{j}^{(i)},\gamma_{j}^{(i)})\in(V(K_{i})\cap\tilde{V})\times\mathbb{R}$
by solving

$\displaystyle\int_{K_{i}}\kappa\nabla\xi_{j}^{(i)}\cdot\nabla v$
$\displaystyle=\gamma_{j}^{(i)}\int_{K_{i}}\xi_{j}^{(i)}v,\;\ \forall v\in
V(K_{i})\cap\tilde{V}$ (6.2)

and rearranging the eigenvalues by
$\gamma_{1}^{(i)}\leq\gamma_{2}^{(i)}\leq\cdots$. For each $K_{i}$, we choose
the first few $J_{i}$ eigenfunctions corresponding to the smallest $J_{i}$
eigenvalues. We define $V_{aux,2}:=\text{span}\\{\xi_{j}^{(i)}:1\leq i\leq
N_{e},1\leq j\leq J_{i}\\}$. For each auxiliary basis function
$\xi_{j}^{(i)}\in V_{aux,2}$, we define a basis function $\zeta_{j}^{(i)}\in
V(K_{i}^{+})$ such that $\mu_{j}^{(i),1}\in V_{aux,1}$, $\mu_{j}^{(i),2}\in
V_{aux,2}$ and

$\displaystyle a(\zeta_{j}^{(i)},v)+s(\mu_{j}^{(i),1},v)+(\mu_{j}^{(i),2},v)$
$\displaystyle=0,\;\forall v\in V(K_{i}^{+}),$ (6.3) $\displaystyle
s(\zeta_{j}^{(i)},\nu)$ $\displaystyle=0,\;\forall\nu\in V_{aux,1},$ (6.4)
$\displaystyle(\zeta_{j}^{(i)},\nu)$
$\displaystyle=(\xi_{j}^{(i)},\nu),\;\forall\nu\in V_{aux,2},$ (6.5)

where we use the notation $V_{aux,1}$ to denote the space $V_{aux}$ defined in
Section 6.1.1. We define

$V_{H,2}=\text{span}\\{\zeta_{j}^{(i)}|\;1\leq i\leq N_{e},\;1\leq j\leq
J_{i}\\}.$

###  6.2 Linear $f(u)$

In this subsection, we discuss the numerical results for

$f(u)=-\nabla\cdot(\kappa\nabla u).$

Equation $(\ref{eq1})$ becomes

$\displaystyle u_{t}-\nabla\cdot(\kappa\nabla u)+g(u)=0.$ (6.6)

For the time discretization, we consider the time step $\Delta
t=\frac{T}{500}=10^{-4}$.

Let $u_{h}$ be the fine mesh solution for Equation $(\ref{eq71})$. We use
Newton’s method to solve the following implicit equation.

$(\frac{u^{n+1}_{h}-u^{n}_{h}}{\Delta
t},v)+a(u^{n+1}_{h},v)+(g(u^{n+1}_{h}),v)=0\quad\forall v\in V_{h},$

where $a(u^{n+1}_{h},v)=\int_{\Omega}\kappa\nabla u^{n+1}_{h}\cdot\nabla v$
and $(\cdot,\cdot)$ is the $L^{2}$ inner product. In finite element methods,
let $\\{\varphi_{i}\\}_{i}$ be fine mesh basis functions. Let $m$ be the step
number in Newton’s method. We have
$u^{n+1,m+1}_{h}=\sum\limits_{i}U^{n+1,m+1}_{h,i}\varphi_{i}$,
$u^{n+1,m}_{h}=\sum\limits_{i}U^{n+1,m}_{h,i}\varphi_{i}$ and
$u^{n}_{h}=\sum\limits_{i}U^{n}_{h,i}\varphi_{i}$. Le $M$ and $A$ be the mass
and stiffness matrices, respectively. Let
$U^{n+1,m+1}_{h}=(U^{n+1,m+1}_{h,i})$, $U^{n+1,m}_{h}=(U^{n+1,m}_{h,i})$ and
$U^{n}_{h}=(U^{n}_{h,i})$. We define

$\displaystyle P(U^{n+1,m}_{h})=MU^{n+1,m}_{h}+\Delta t\cdot
AU^{n+1,m}_{h}+\Delta t\cdot\mathcal{G}-MU^{n}_{h},$

where $\mathcal{G}=(\mathcal{G}_{i})$

$\mathcal{G}_{i}=(g(u^{n+1,m}_{h}),\varphi_{i}).$

Then

$\displaystyle(JP)(U^{n+1,m}_{h})=M+\Delta t\cdot A+\Delta
t\cdot(J\mathcal{G}),$

where $J\mathcal{G}=((J\mathcal{G})_{ij})$

$(J\mathcal{G})_{ij}=\frac{\partial(g(u^{n+1,m}_{h}),\varphi_{i})}{\partial
U^{n+1,m}_{h,j}}.$

Then we have

$U^{n+1,m+1}_{h}=U^{n+1,m}_{h}-(JP)^{-1}(U^{n+1,m}_{h})P(U^{n+1,m}_{h}).$

Newton’s method for coarse mesh is similar. The partially explicit scheme is:

$\displaystyle(\frac{u^{n+1}_{H,1}-u^{n}_{H,1}}{\Delta
t}+\frac{u^{n}_{H,2}-u^{n-1}_{H,2}}{\Delta
t},v_{1})+a((u^{n+1}_{H,1}+u^{n}_{H,2}),v_{1})+(g(u^{n}_{H,1}+u^{n}_{H,2}),v_{1})=0\quad\forall
v_{1}\in V_{H,1},$ $\displaystyle(\frac{u^{n+1}_{H,2}-u^{n}_{H,2}}{\Delta
t}+\frac{u^{n}_{H,1}-u^{n-1}_{H,1}}{\Delta
t},v_{2})+a((u^{n+1}_{H,1}+u^{n}_{H,2}),v_{2})+(g(u^{n}_{H,1}+u^{n}_{H,2}),v_{2})=0\quad\forall
v_{2}\in V_{H,2}.$

In our first example, we consider

$g(u)=-(10\cdot u\cdot(u^{2}-1)+g_{0}).$

In Figure 6.1, the permeability field $\kappa$ and $g_{0}$ are presented. As
is shown, this permeability field has heterogeneous high contrast channels and
$g_{0}$ is a singular source term. In Figure 6.2, we first present the
reference solution which is implicitly solved using fine grid basis functions.
The middle plot in Figure 6.2 is implicit CEM solution obtained with
additional basis functions and the solution in the right plot is obtained
using the partially explicit scheme presented above. These three plots all
show the solution at $t=T$. We present two relative error plots in Figure 6.3.
The first one is the relative $L^{2}$ error plot and the second one is the
relative energy error plot. The blue, red and black curves (in both plots)
stand for the relative error for implicit CEM solution, implicit CEM solution
(with additional basis) and partially explicit solution respectively. In each
of these two plots, there is a noticeable improvement for error when we use
additional basis functions. We find that the black curve coincides with the
red curve, which means that the partially explicit scheme can achieve similar
accuracy as the fully implicit scheme.

Figure 6.1: Left: $\kappa$. Right: $g_{0}$.

Figure 6.2: Left: Reference solution at $t=T$. Middle: Implicit CEM solution
(with additional basis) at $t=T$. Right: Partially explicit solution at $t=T$.

Figure 6.3: Left: Relative $L^{2}$ error. Right: Relative energy error.

In this case,

$g(u)=-(10\cdot u\cdot(u^{2}-1)+g_{0}).$

The difference is that we use a smooth source term. In Figure 6.4, the
permeability field $\kappa$ and source term $g_{0}$ are shown. The reference
solution at the final time, implicit CEM solution (with additional basis) at
the final time and partially explicit solution at the final time are presented
in Figure 6.5. We show the relative $L^{2}$ error plot and the relative energy
error plot in Figure 6.6. We see that the relative $L^{2}$ and energy error
curves for implicit CEM (with additional basis) and partially explicit scheme
almost coincide, which implies similar accuracy between them.

Figure 6.4: Left: $\kappa$. Right: $g_{0}$.

Figure 6.5: Left: Reference solution at $t=T$. Middle: Implicit CEM solution
(with additional basis) at $t=T$. Right: Partially explicit solution at $t=T$.

Figure 6.6: Left: Relative $L^{2}$ error. Right: Relative energy error.

In our third test,

$g(u)=-(10\cdot u\cdot(u^{2}-1)+g_{0}),$

and we use a more complicated permeability field with more high contrast
channels. Figure 6.7 shows the permeability field $\kappa$ and source term
$g_{0}$. The reference solution at $t=T$, implicit CEM solution (with
additional basis) at $t=T$ and partially explicit solution at $t=T$ are
presented in Figure 6.8. In Figure 6.9, we show the relative $L^{2}$ error
plot and the relative energy error plot. In this case, in both error plots,
the curves for implicit CEM (with additional basis) and partially explicit
scheme coincide.

Figure 6.7: Left: $\kappa$. Right: $g_{0}$.

Figure 6.8: Left: Reference solution at $t=T$. Middle: Implicit CEM solution
(with additional basis) at $t=T$. Right: Partially explicit solution at $t=T$.

Figure 6.9: Left: Relative $L^{2}$ error. Right: Relative energy error.

In this example,

$g(u)=-(10\cdot u\cdot(u^{2}-1)+g_{0}).$

We use the more complicated permeability field and the smooth source term
which are shown in Figure 6.10. In Figure 6.11, the reference solution at the
final time, implicit CEM solution (with additional basis) at the final time
and partially explicit solution at the final time are presented. Relative
$L^{2}$ error and energy error plots are shown in Figure 6.12. In this case,
the relative error for implicit CEM scheme is small and comparable to the two
schemes with additional basis. From Figure 6.12, we can find that the $L^{2}$
and energy error for implicit CEM (with additional basis) and partially
explicit scheme are nearly the same.

Figure 6.10: Left: $\kappa$. Right: $g_{0}$.

Figure 6.11: Left: Reference solution at $t=T$. Middle: Implicit CEM solution
(with additional basis) at $t=T$. Right: Partially explicit solution at $t=T$.

Figure 6.12: Left: Relative $L^{2}$ error. Right: Relative energy error.

In this case, we use a new reaction term

$g=-(1+\cos(a_{1}\cdot u)+g_{0}),\;a_{1}(x,y)=2\cos(20\pi x)\cos(20\pi y).$

In numerical experiments, we set $a_{1}$ to be constant inside every fine
element. Figure 6.13 shows the permeability field $\kappa$, the source term
$g_{0}$ and the function $a_{1}$. The reference solution, implicit CEM
solution (with additional basis) and partially explicit solution at the final
time are shown in Figure 6.14. The relative $L^{2}$ error plot and relative
energy error plot are presented in Figure 6.15. From the relative $L^{2}$
error plot, we find that the relative $L^{2}$ error for three schemes are
comparable. For the relative energy error, we can find a large improvement
when we introduce additional basis. The relative energy error for implicit CEM
solution (with additional basis) and partially explicit solution are nearly
the same.

Figure 6.13: Left: $\kappa$. Middle: $g_{0}$. Right: $a_{1}$.

Figure 6.14: Left: Reference solution at $t=T$. Middle: Implicit CEM solution
(with additional basis) at $t=T$. Right: Partially explicit solution at $t=T$.

Figure 6.15: Left: Relative $L^{2}$ error. Right: Relative energy error.

In this case,

$g=-(1+\cos(a_{1}\cdot u)+g_{0}),\;a_{1}(x,y)=2\cos(20\pi x)\cos(20\pi y).$

The permeability field $\kappa$, the source term $g_{0}$ and the funciotn
$a_{1}$ are presented in Figure 6.16. The reference solution, implicit CEM
solution (with additional basis) and partially explicit solution at the final
time are shown in Figure 6.17. We present the relative $L^{2}$ error plot and
the relative energy error plot in Figure 6.18. We observe that the $L^{2}$ and
energy error curves for implicit CEM (with additional basis) and partially
explicit scheme coincide.

Figure 6.16: Left: $\kappa$. Middle: $g_{0}$. Right: $a_{1}$.

Figure 6.17: Left: Reference solution at $t=T$. Middle: Implicit CEM solution
(with additional basis) at $t=T$. Right: Partially explicit solution at $t=T$.

Figure 6.18: Left: Relative $L^{2}$ error. Right: Relative energy error.

###  6.3 Nonlinear $f(u)$

We want to explore the case in which the diffusion operator is nonlinear. So
in Equation $(\ref{eq1})$, we set

$f(u)=-\nabla\cdot(\kappa\alpha(u)\nabla u).$

Equation $(\ref{eq1})$ becomes

$\displaystyle u_{t}-\nabla\cdot(\kappa\alpha(u)\nabla u)+g(u)=0.$ (6.7)

Let $u_{h}$ be the fine mesh solution for Equation $(\ref{eq72})$. We use the
Picard-Newton iteration to solve the implicit equation.

$\Big{(}\frac{u^{n+1,m+1}_{h}-u^{n}_{h}}{\Delta
t},v\Big{)}+\int_{\Omega}\kappa\alpha(u^{n+1,m}_{h})\nabla
u^{n+1,m+1}_{h}\cdot\nabla v+(g(u^{n+1,m+1}_{h}),v)=0\quad\forall v\in V_{h},$

where $m$ is the Picard-Newton step number and $(\cdot,\cdot)$ is the $L^{2}$
inner product. In finite element methods, let $\\{\varphi_{i}\\}_{i}$ be fine
mesh basis functions. We have
$u^{n+1,m+1}_{h}=\sum\limits_{i}U^{n+1,m+1}_{h,i}\varphi_{i}$,
$u^{n+1,m}_{h}=\sum\limits_{i}U^{n+1,m}_{h,i}\varphi_{i}$ and
$u^{n}_{h}=\sum\limits_{i}U^{n}_{h,i}\varphi_{i}$. Let $M$ be the mass matrix.
Let $U^{n+1,m+1}_{h}=(U^{n+1,m+1}_{h,i})$, $U^{n+1,m}_{h}=(U^{n+1,m}_{h,i})$
and $U^{n}_{h}=(U^{n}_{h,i})$. We define

$\displaystyle Q(U^{n+1,m}_{h})=MU^{n+1,m}_{h}+\Delta
t\cdot\mathcal{A}U^{n+1,m}_{h}-MU^{n}_{h}+\Delta t\cdot\mathcal{G},$

where $\mathcal{A}=(\mathcal{A}_{ij})$

$\mathcal{A}_{ij}=\int_{\Omega}\kappa\alpha(u^{n+1,m}_{h})\nabla\varphi_{j}\nabla\varphi_{i}$

and $\mathcal{G}=(\mathcal{G}_{i})$

$\mathcal{G}_{i}=(g(u^{n+1,m}_{h}),\varphi_{i}).$

Then

$(JQ)(U^{n+1,m}_{h})=M+\Delta t\mathcal{A}+\Delta t\cdot(J\mathcal{G}),$

where $J\mathcal{G}=((J\mathcal{G})_{ij})$

$(J\mathcal{G})_{ij}=\frac{\partial(g(u^{n+1,m}_{h}),\varphi_{i})}{\partial
U^{n+1,m}_{h,j}}.$

Thus,

$U^{n+1,m+1}_{h}=U^{n+1,m}_{h}-(JQ)^{-1}(U^{n+1,m}_{h})Q(U^{n+1,m}_{h}).$

Newton’s method for coarse mesh is similar. For the partially explicit scheme,
we use the following Picard-Newton iteration which can be solved similarly
using the method introduced above. Note that we change the reaction term to be
partially explicit.

$\displaystyle(\frac{u^{n+1,m+1}_{H,1}-u^{n}_{H,1}}{\Delta
t}+\frac{u^{n}_{H,2}-u^{n-1}_{H,2}}{\Delta
t},v_{1})+\int_{\Omega}\kappa\alpha(u^{n+1,m}_{H,1}+u^{n}_{H,2})\nabla(u^{n+1,m+1}_{H,1}+u^{n}_{H,2})\cdot\nabla
v_{1}$ $\displaystyle=(-g(u^{n+1,m+1}_{H,1}+u^{n}_{H,2}),v_{1})\quad\forall
v_{1}\in V_{H,1},$ $\displaystyle(\frac{u^{n+1}_{H,2}-u^{n}_{H,2}}{\Delta
t}+\frac{u^{n}_{H,1}-u^{n-1}_{H,1}}{\Delta
t},v_{2})+\int_{\Omega}\kappa\alpha(u^{n+1}_{H,1}+u^{n}_{H,2})\nabla(u^{n+1}_{H,1}+u^{n}_{H,2})\cdot\nabla
v_{2}$ $\displaystyle=(-g(u^{n+1}_{H,1}+u^{n}_{H,2}),v_{2})\quad\forall
v_{2}\in V_{H,2}.$

In the following three examples, we use $g(u_{H,1}^{n+1}+u_{H,2}^{n})$ in the
partially explicit algorithm. One can prove its stability as in Lemma 4.1 and
Lemma 4.2. The proof is similar and we omit it here.
In the following three cases,

$g(u)=-(10u(u^{2}-1)+g_{0}).$

In the first example, we consider $\alpha(u)=1+u^{2}$ and time step $\Delta
t=\frac{T}{500}=10^{-4}$. The permeability field $\kappa$ and the source term
$g_{0}$ are presented in Figure 6.19. The reference solution, implicit CEM
solution (with additional basis) and partially explicit solution at the final
time are shown in Figure 6.20. We present the relative $L^{2}$ error plot and
relative energy error plot in Figure 6.21. We can find that the curves for
implicit CEM solution (with additional basis) and partially explicit solution
coincide, which implies that our new partially explicit scheme is also
effective and has similar accuracy as the implicit CEM (with additional basis)
scheme.

Figure 6.19: Left: $\kappa$. Right: $g_{0}$.

Figure 6.20: Left: Reference solution at $t=T$. Middle: Implicit CEM solution
(with additional basis) at $t=T$. Right: Partially explicit solution at $t=T$.

Figure 6.21: Left: Relative $L^{2}$ error. Right: Relative energy error.

In this second example, we consider $\alpha(u)=2+\cos(u)$ and time step
$\Delta t=\frac{0.05}{1500}$. Figure 6.22 shows the permeability field
$\kappa$ and the source term $g_{0}$. The reference solution, implicit CEM
solution (with additional basis) and partially explicit solution at $t=T$ are
presented in Figure 6.23. The relative $L^{2}$ error plot and relative energy
error plot are shown in Figure 6.24. The partially explicit scheme also works
in this case and has similar accuracy as the implicit CEM (with additional
basis) scheme.

Figure 6.22: Left: $\kappa$. Right: $g_{0}$.

Figure 6.23: Left: Reference solution at $t=T$. Middle: Implicit CEM solution
(with additional basis) at $t=T$. Right: Partially explicit solution at $t=T$.

Figure 6.24: Left: Relative $L^{2}$ error. Right: Relative energy error.

In this case, we have $\alpha(u)=2+\cos(u)$ and time step $\Delta
t=\frac{0.05}{1500}$. The permeability field $\kappa$ and the source term
$g_{0}$ are presented in Figure 6.25. We show the reference solution, implicit
CEM solution (with additional basis) and partially explicit solution at $t=T$
in Figure 6.26. The relative $L^{2}$ error plot and the relative energy error
plot are presented in Figure 6.27. From Figure 6.27, we see that the curves
for implicit CEM (with additional basis) and partially explicit scheme
coincide, which implies that they have similar accuracy.

Figure 6.25: Left: $\kappa$. Right: $g_{0}$.

Figure 6.26: Left: Reference solution at $t=T$. Middle: Implicit CEM solution
(with additional basis) at $t=T$. Right: Partially explicit solution at $t=T$.

Figure 6.27: Left: Relative $L^{2}$ error. Right: Relative energy error.

##  7 Conclusions

In this work, we design and analyze contrast-independent time discretization
for nonlinear problems. The work continues our earlier works on linear
problems, where we propose temporal splitting and associated spatial
decomposition that guarantees a stability. We introduce two spatial spaces,
the first account for spatial features related to fast time scales and the
second for spatial features related to “slow” time scales. We propose time
splitting, where the first equation solves for fast components implicitly and
the second equation solves for slow components explicitly. We introduce a
condition for multiscale spaces that guarantees stability of the proposed
splitting algorithm. Our proposed method is still implicit via mass matrix;
however, it is explicit in terms of stiffness matrix for the slow component.
We present numerical results, which show that the proposed methods provide
very similar results as fully implicit methods using explicit methods with the
time stepping that is independent of the contrast.

## Acknowledgments

The research of Eric Chung is partially supported by the Hong Kong RGC General
Research Fund (Project numbers 14304719 and 14302018) and the CUHK Faculty of
Science Direct Grant 2020-21.

## References

  * [1] A. Abdulle.  Explicit methods for stiff stochastic differential equations.  In Numerical Analysis of Multiscale Computations, pages 1–22. Springer, 2012.
  * [2] G. Ariel, B. Engquist, and R. Tsai.  A multiscale method for highly oscillatory ordinary differential equations with resonance.  Mathematics of Computation, 78(266):929–956, 2009.
  * [3] U. M. Ascher, S. J. Ruuth, and R. J. Spiteri.  Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations.  Applied Numerical Mathematics, 25(2-3):151–167, 1997.
  * [4] J. Bear.  Dynamics of fluids in porous media.  Courier Corporation, 2013.
  * [5] D. L. Brown, Y. Efendiev, and V. H. Hoang.  An efficient hierarchical multiscale finite element method for stokes equations in slowly varying media.  Multiscale Modeling & Simulation, 11(1):30–58, 2013.
  * [6] E. T. Chung, Y. Efendiev, and T. Hou.  Adaptive multiscale model reduction with generalized multiscale finite element methods.  Journal of Computational Physics, 320:69–95, 2016.
  * [7] E. T. Chung, Y. Efendiev, and C. Lee.  Mixed generalized multiscale finite element methods and applications.  SIAM Multiscale Model. Simul., 13:338–366, 2014.
  * [8] E. T. Chung, Y. Efendiev, and W. T. Leung.  Generalized multiscale finite element methods for wave propagation in heterogeneous media.  Multiscale Modeling & Simulation, 12(4):1691–1721, 2014.
  * [9] E. T. Chung, Y. Efendiev, and W. T. Leung.  Constraint energy minimizing generalized multiscale finite element method.  Computer Methods in Applied Mechanics and Engineering, 339:298–319, 2018.
  * [10] E. T. Chung, Y. Efendiev, and W. T. Leung.  Constraint energy minimizing generalized multiscale finite element method in the mixed formulation.  Computational Geosciences, 22(3):677–693, 2018.
  * [11] E. T. Chung, Y. Efendiev, and W. T. Leung.  Fast online generalized multiscale finite element method using constraint energy minimization.  Journal of Computational Physics, 355:450–463, 2018.
  * [12] E. T. Chung, Y. Efendiev, W. T. Leung, M. Vasilyeva, and Y. Wang.  Non-local multi-continua upscaling for flows in heterogeneous fractured media.  Journal of Computational Physics, 372:22–34, 2018.
  * [13] J. Du and Y. Yang.  Third-order conservative sign-preserving and steady-state-preserving time integrations and applications in stiff multispecies and multireaction detonations.  Journal of Computational Physics, 395:489–510, 2019.
  * [14] L. Duchemin and J. Eggers.  The explicit–implicit–null method: Removing the numerical instability of pdes.  Journal of Computational Physics, 263:37–52, 2014.
  * [15] W. E and B. Engquist.  Heterogeneous multiscale methods.  Comm. Math. Sci., 1(1):87–132, 2003.
  * [16] Y. Efendiev, J. Galvis, and T. Hou.  Generalized multiscale finite element methods (GMsFEM).  Journal of Computational Physics, 251:116–135, 2013.
  * [17] Y. Efendiev, J. Galvis, G. Li, and M. Presho.  Generalized multiscale finite element methods. nonlinear elliptic equations.  Communications in Computational Physics, 15(3):733–755, 2014.
  * [18] Y. Efendiev and T. Hou.  Multiscale Finite Element Methods: Theory and Applications, volume 4 of Surveys and Tutorials in the Applied Mathematical Sciences.  Springer, New York, 2009.
  * [19] Y. Efendiev and A. Pankov.  Numerical homogenization of monotone elliptic operators.  SIAM J. Multiscale Modeling and Simulation, 2(1):62–79, 2003.
  * [20] Y. Efendiev and A. Pankov.  Homogenization of nonlinear random parabolic operators.  Advances in Differential Equations, 10(11):1235–1260, 2005.
  * [21] W. Ehlers.  Darcy, forchheimer, brinkman and richards: classical hydromechanical equations and their significance in the light of the tpm.  Archive of Applied Mechanics, pages 1–21, 2020.
  * [22] B. Engquist and Y.-H. Tsai.  Heterogeneous multiscale methods for stiff ordinary differential equations.  Mathematics of computation, 74(252):1707–1742, 2005.
  * [23] W. T. L. Eric Chung, Yalchin Efendiev and P. N. Vabishchevich.  Contrast-independent partially explicit time discretizations for multiscale flow problems.  arXiv:2101.04863.
  * [24] W. T. L. Eric Chung, Yalchin Efendiev and P. N. Vabishchevich.  Contrast-independent partially explicit time discretizations for multiscale wave problems.
  * [25] J. Frank, W. Hundsdorfer, and J. Verwer.  On the stability of implicit-explicit linear multistep methods.  Applied Numerical Mathematics, 25(2-3):193–205, 1997.
  * [26] P. Henning, A. Målqvist, and D. Peterseim.  A localized orthogonal decomposition method for semi-linear elliptic problems.  ESAIM: Mathematical Modelling and Numerical Analysis, 48(5):1331–1349, 2014.
  * [27] T. Hou and X. Wu.  A multiscale finite element method for elliptic problems in composite materials and porous media.  J. Comput. Phys., 134:169–189, 1997.
  * [28] T. Y. Hou, D. Huang, K. C. Lam, and P. Zhang.  An adaptive fast solver for a general class of positive definite matrices via energy decomposition.  Multiscale Modeling & Simulation, 16(2):615–678, 2018.
  * [29] T. Y. Hou, Q. Li, and P. Zhang.  Exploring the locally low dimensional structure in solving random elliptic pdes.  Multiscale Modeling & Simulation, 15(2):661–695, 2017.
  * [30] T. Y. Hou, D. Ma, and Z. Zhang.  A model reduction method for multiscale elliptic pdes with random coefficients using an optimization approach.  Multiscale Modeling & Simulation, 17(2):826–853, 2019.
  * [31] W. Hundsdorfer and S. J. Ruuth.  Imex extensions of linear multistep methods with general monotonicity and boundedness properties.  Journal of Computational Physics, 225(2):2016–2042, 2007.
  * [32] G. Izzo and Z. Jackiewicz.  Highly stable implicit–explicit runge–kutta methods.  Applied Numerical Mathematics, 113:71–92, 2017.
  * [33] P. Jenny, S. Lee, and H. Tchelepi.  Multi-scale finite volume method for elliptic problems in subsurface flow simulation.  J. Comput. Phys., 187:47–67, 2003.
  * [34] C. Le Bris, F. Legoll, and A. Lozinski.  An MsFEM type approach for perforated domains.  Multiscale Modeling & Simulation, 12(3):1046–1077, 2014.
  * [35] T. Li, A. Abdulle, et al.  Effectiveness of implicit methods for stiff stochastic differential equations.  In Commun. Comput. Phys. Citeseer, 2008.
  * [36] G. I. Marchuk.  Splitting and alternating direction methods.  Handbook of numerical analysis, 1:197–462, 1990.
  * [37] M. Narayanamurthi, P. Tranquilli, A. Sandu, and M. Tokman.  Epirk-w and epirk-k time discretization methods.  Journal of Scientific Computing, 78(1):167–201, 2019.
  * [38] H. Owhadi and L. Zhang.  Metric-based upscaling.  Comm. Pure. Appl. Math., 60:675–723, 2007.
  * [39] A. Roberts and I. Kevrekidis.  General tooth boundary conditions for equation free modeling.  SIAM J. Sci. Comput., 29(4):1495–1510, 2007.
  * [40] S. J. Ruuth.  Implicit-explicit methods for reaction-diffusion problems in pattern formation.  Journal of Mathematical Biology, 34(2):148–176, 1995.
  * [41] G. Samaey, I. Kevrekidis, and D. Roose.  Patch dynamics with buffers for homogenization problems.  J. Comput. Phys., 213(1):264–287, 2006.
  * [42] H. Shi and Y. Li.  Local discontinuous galerkin methods with implicit-explicit multistep time-marching for solving the nonlinear cahn-hilliard equation.  Journal of Computational Physics, 394:719–731, 2019.
  * [43] P. N. Vabishchevich.  Additive Operator-Difference Schemes: Splitting Schemes.  Walter de Gruyter GmbH, Berlin, Boston, 2013.