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

: This work continues a line of work on developing partially explicit methods for multiscale problems. In our previous works, we considered linear multiscale problems where the spatial heterogeneities are at the 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 deﬁned via multiscale space decomposition. Following this decomposition, temporal splitting was proposed, which treats fast components implicitly and slow components explicitly. The space decomposition and temporal splitting are chosen such that they guarantees stability, and we formulated a condition for the time stepping. This condition was 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 a time step that is independent of the contrast.


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.
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 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(κ(x, u)∇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.

Problem Setting
We consider the following equation where f = δF δu and g = δG δu which are the variational derivative of energies F (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) ∈ V * and g(u) ∈ L 2 (Ω) for all u ∈ V .We can then consider the weak formulation of the problem, namely, finding u ∈ V such that To simplify the notation, we will simply write (•, •) instead of (•, •) V * ,V in the following discussion.
• The second variational derivative δ 2 F and δ 2 G satisfy • The second variational derivative δ 2 F and δ 2 G are bounded.That is, where 0 < C(u) < ∞ and 0 < B < ∞ are independent on v, w.

Discretization
To solve the problem, a standard method is finite element approach.We can consider the numerical solution 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 k H } N k=0 ⊂ V H such that For the backward Euler method, we consider Next, we would like to derive stability conditions for backward and forward Euler methods.Since and for some We have for any ∆t and thus, backward Euler method is stable if ∆tb ≤ 1.
Similarly, for forward Euler method, we can use and obtain We can see that although forward Euler method is easier for implementation, we require a small time is large.We remark that in typical cases we will consider b and B are not too large.Therefore, we have ∆tb and ∆tB is small and the energy G will not affect the stability too much.

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, The finite element solution is then satisfying where u H = u H,1 + u H,2 .We can use a partially explicit time discretization.For example, we can consider Energy stability where Proof.By substitute Summing up the above two equations, we have and obtain To prove the stability of the method, we can consider for some and Thus, we obtain then we have where c = inf u∈V H c(u), C2 = sup ξ∈V H C 2 (ξ) and we have 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 . Similar to the proof in previous lemma, we consider ) ) with λ ∈ (0, 1).We have we have we have the condition formulated in Lemma 4.1 By Lemma 4.1, we get the result.

G = 0 case
First, we present some discussions for G = 0 case.In this case, the first stability condition for partial explicit scheme is This condition can be understood as a nonlinear constraint on the "second space" that represents u n H,2 and in order to have a small bound, one needs to guarantee that u n H,1 captures important degrees of freeom.Indeed, the smallness of is a condition on u n H,1 (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 5.2 G = 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.

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 and nonlinear case In all examples, we will use two heterogeneous high contrast κ(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 1 10 and the fine mesh size is 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.

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 κ(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 Ṽ , 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 1 0 (S) for a proper subset S ⊂ Ω.

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 T H be a coarse grid partition of Ω.For K i ∈ T H , we first need to define a set of auxiliary basis functions in V (K i ).We solve Ki κ∇ψ (i) where with {χ 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 where s(u, v) := Ne i=1 s i (u| Ki , v| Ki ) 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 ψ (i) j , we can find a local basis function φ for some µ We then define the space V cem as Let Ṽ := {v ∈ V : Π(v) = 0} and we can now construct V H,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 (ξ and rearranging the eigenvalues by γ For each K i , we choose the first few J i eigenfunctions corresponding to the smallest J i eigenvalues.We define V aux,2 := span{ξ where we use the notation V aux,1 to denote the space V aux defined in Section 6.1.1.We define

Linear f (u)
In this subsection, we discuss the numerical results for Equation (2.1) becomes For the time discretization, we consider the time step ∆t = T 500 = 10 −4 .Let u h be the fine mesh solution for Equation (6.6).We use Newton's method to solve the following implicit equation.
• ∇v and (•, •) is the L 2 inner product.In finite element methods, let {ϕ i } i be fine mesh basis functions.Let m be the step number in Newton's method.We have u n+1,m+1 Le M and A be the mass and stiffness matrices, respectively.Let U n+1,m+1 h ) and U n h = (U n h,i ).We define Then where JG = ((JG) ij ) .
Then we have Newton's method for coarse mesh is similar.The partially explicit scheme is: In our first example, we consider In Figure 6.1, the permeability field κ 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.In this case, The difference is that we use a smooth source term.In Figure 6.4, the permeability field κ 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.In our third test, and we use a more complicated permeability field with more high contrast channels.Figure 6.7 shows the permeability field κ 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.In this example, g(u) = −(10 • u • (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.In this case, we use a new reaction term g = −(1 + cos(a 1 • u) + g 0 ), a 1 (x, y) = 2 cos(20πx) cos(20πy).
In numerical experiments, we set a 1 to be constant inside every fine element.Figure 6.13 shows the permeability field κ, 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.The permeability field κ, 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.

Nonlinear f (u)
We want to explore the case in which the diffusion operator is nonlinear.So in Equation (2.1), we set Let u h be the fine mesh solution for Equation (6.7).We use the Picard-Newton iteration to solve the implicit equation.
where m is the Picard-Newton step number and (•, •) is the L 2 inner product.In finite element methods, let {ϕ i } i be fine mesh basis functions.We have u n+1,m+1 ) and 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.
In the following three examples, we use g(u n+1 H,1 + u n H,2 ) 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 α(u) = 1 + u 2 and time step ∆t = T 500 = 10 −4 .The permeability field κ 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.In this second example, we consider α(u) = 2 + cos(u) and time step ∆t = 0.05 1500 .Figure 6.22 shows the permeability field κ 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.In this case, we have α(u) = 2 + cos(u) and time step ∆t = 0.05 1500 .The permeability field κ 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.

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.

Figure 6 . 2 :
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 .
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 .
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 .
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 .
Figure 6.20: Left: Reference solution at t = T .Middle: Implicit CEM solution (with additional basis) at = T .Right: Partially explicit solution at t = T .

Figure 6 .
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 .
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 .