Implicit-Explicit schemes for decoupling multicontinuum problems in porous media
Maria Vasilyeva



SFMViT : SlowFast Meet ViT in Chaotic World
    [
    April 29, 2024
===========================================




In this work, we present an efficient way to decouple the multicontinuum problems. To construct decoupled schemes, we propose Implicit-Explicit time approximation in general form and study them for the fine-scale and coarse-scale space approximations. 
We use a finite-volume method for fine-scale approximation, and the nonlocal multicontinuum (NLMC) method is used to construct an accurate and physically meaningful coarse-scale approximation. The NLMC method is an accurate technique to develop a physically meaningful coarse scale model based on defining the macroscale variables. The multiscale basis functions are constructed in local domains by solving constraint energy minimization problems and projecting the system to the coarse grid. The resulting basis functions have exponential decay properties and lead to the accurate approximation on a coarse grid. We construct a fully Implicit time approximation for semi-discrete systems arising after fine-scale and coarse-scale space approximations. 
We investigate the stability of the two and three-level schemes for fully Implicit and Implicit-Explicit time approximations schemes for multicontinuum problems in fractured porous media. We show that combining the decoupling technique with multiscale approximation leads to developing an accurate and efficient solver for multicontinuum problems. 




§ INTRODUCTION



Multicontinuum models occur in many real-world applications, such as geothermal energy production, unconventional oil and gas production, carbon capture and storage (CCS) in geological formations, deposition of nuclear waste, wastewater treatment, and many more <cit.>. 
In reservoir simulation, the multicontimuum models are used to describe complex interaction between multiple scales of heterogeneity such as hydraulic and natural fractures, vugs and cavities, and porous matrix <cit.>. In <cit.>, the general form of the dual-continuum model is derived from homogenization theory. 
Multicontinuum systems are characterized by high contrast properties between continua, and each continua may itself have a complex heterogeneous structure. Moreover, fractured porous media usually have complex fracture geometries and minimal thicknesses compared to typical reservoir sizes. A typical model approach for fractured media is based on the lower-dimensional representation of the fracture objects   <cit.>. In such mathematical models, we have a coupled mixed-dimensional multicontinuum system of equations.



Numerical simulation of flow in multicontinuum fracture porous media is a challenging task. Due to the multiple scales and complex geometries of the fractures, the simulation of the processes in fractured porous media requires fine grids for accurate approximation, which is computationally expensive. Moreover, the significant contrast in the continuum's permeability makes numerical simulations even more computationally expensive.
To construct efficient solvers, the multiscale methods are developed to reduce the dimension of the discrete system and perform simulations with fewer degrees of freedom than classical approximation methods such as finite volume, finite element, or finite difference methods <cit.>.
Several multiscale methods are developed to solve problems in fractured porous media. In <cit.>, an iterative multiscale finite volume (i-MSFV) method is developed for flow in fractured porous media. The process is based on the Multiscale Finite Volume Method <cit.>. 
A multiscale restriction smoothed basis method (F-MsRSB) for multiphase flow in heterogeneous fractured porous media is devised in <cit.>. 
In our previous works, we presented the multiscale method for flow in fractured porous media <cit.>. The proposed approach is based on the Generalized Multiscale Finite Element Method (GMsFEM) and utilizes local spectral problems in multiscale basis construction <cit.>. The GMsFEM is effectively extended to solve problems in multicontinuum fractured porous media  <cit.>. 
Recently, the Constraint Energy Minimization Generalized Multiscale Finite Element Method (CEM-GMsFEM) was proposed in <cit.>. In CEM-GMsFEM, constructing multiscale basis functions starts with an auxiliary multiscale space that is defined by solving local spectral problems. Then, a constraint energy minimization is used to build multiscale basis functions in the oversampling domain. 
In CEM-GMsFEM, the choice of the auxiliary space can be made using careful design of the macroscale parameters. Moreover, the resulting construction can lead to a physically meaningful coarse scale model. In <cit.>, we recently presented a nonlocal multicontinuum (NLMC) method for problems in fractured media. In the NLMC approximation, the resulting system is very similar to the traditional finite volume method but provides an accurate approximation due to the nonlocal coupling in and between each continuum. 
The NLMC method effectively solves different applied problems in multicontinuum media <cit.>. 



To develop an efficient numerical implementation for multicontinuum problems, we construct decoupled (splitted) schemes to separate equations for each continuum. Decoupled discrete schemes can significantly reduce the computational work and allow using different methods and software for the subproblem solution.  
The main idea of the splitting technique is to separate the original problem into several smaller subproblems, which can be designed for efficient computational implementation. In time approximation, additive schemes are a helpful tool for solving unsteady equations <cit.>.
The application of the additive schemes for multiscale approximation has been considered in <cit.>. The partially-explicit time discretization for nonlinear multiscale problems is presented in <cit.>. The method was combined with the machine learning techniques applied for the implicit part of the operator <cit.>. 
The next extension of the partially explicit scheme is presented in <cit.>, where the mutirate partially explicit scheme for multiscale flow problems is developed and analyzed. 
In <cit.>, we proposed an efficient decoupled scheme for multicontinuum flow problems in fractured porous media. The presented approach is based on the additive representation of the operator with implicit-explicit approximation by time to decoupled equations for each continuum. We developed, analyzed, and investigated three first-order decoupled schemes for solving the classical multicontinuum problems in fractured porous media on fine grids with finite volume approximation by space. We extend the continuum decoupling approach for multiscale multicontinuum problems with nonlocal multicontinuum approximation on the coarse grid.



In this work, we develop decoupled multiscale schemes for solving multiscale multicontinuum problems in general form. With the help of an implicit-explicit time approximation, we decouple the fine-scale and coarse-scale discrete systems to construct an efficient and robust numerical algorithm. This work is motivated by <cit.> where implicit-explicit schemes are utilized to the splitting of time-dependent problems by an additive representation of the problem operator.  

The Implicit-Explicit (ImEx) schemes are widely used for diffusion-reaction or diffusion-convection problems. For the convection-diffusion problems, an explicit scheme is applied for the convection part, and an implicit scheme is utilized for the diffusion part because it gives a more substantial restriction on the time step size in the explicit approximation. In the convection-diffusion operator, the reason to approximate explicitly can be delivered to the non-symmetric part of the matrix or Jacobian for nonlinear systems for a general case.
In reaction-diffusion problems, the reaction part can be approximated explicitly <cit.>.  
In <cit.> an  implicit-explicit scheme is constructed for radiation hydrodynamics, where coupled radiation equations are treated using implicit integration and hydrodynamics is treated by explicit integration.
Moreover, in such applications, ImEx schemes are used to decouple the slow and fast parts of the operator. The stiff part uses the implicit approximation to avoid small-time stepping. Furthermore, the implicit part is typically associated with a  diffusion operator and can be solved effectively by iterative methods. 
In <cit.>, linear multistep ImEx schemes are developed for general equations. Such schemes can be constructed for spatially discretized partial differential equations. In <cit.>, the fundamental aspects of the splitting schemes are discussed for the problems described by time-dependent partial differential equations. The construction of the splitting schemes is based on the additive representation of the operator of the time-dependent problem, where various schemes are developed based on the implicit-explicit approximation. 
Due to the additive representation of an operator of a time-dependent problem as a sum of operators with a simpler structure, the transition to a new time level can be performed as a sequence of subproblems. 



Motivated by Implicit-Explicit schemes in <cit.> and additive schemes for time-dependent problems in  <cit.>, we develop decoupling schemes in general form for multicontinuum problems with fine-scale and coarse-scale approximations by space. The operator is split into the block diagonal and off-block-diagonal parts to apply it to decouple the high-contrast continua.  
In this work, we continue the development of the efficient decoupling schemes proposed in our previous work <cit.> and propose a general form of the Implicit-Explicit schemes. 
We note that the first-order schemes are similar to the scheme proposed in <cit.>. The paper's main novelty is developing a general form of the decoupling scheme for fine-scale and coarse-scale problems. We systematically analyze the stability and computational performance of the proposed scheme for two and three-level schemes. Finally, numerical experiments demonstrate that the proposed implicit-explicit two and three-level schemes give an efficient decoupling approach for multicontinuum problems, reduce the number of iterations and size of the discrete problems, and lead to a highly accurate and computationally efficient algorithm. 




The paper is organized as follows. Section 2 describes problem formulation for flow in multicontinuum media in a general form and presents a semi-discretization based on finite-volume method (FVM) and nonlocal multicontinua method (NLMC). The Implicit time approximation that leads to the coupled system of equations is given in Section 3. Section 4 presents the Implicit-Explicit schemes to decouple the multicontinuum system. Stability analysis is provided for two- and three-level Implicit and Implicit-Explicit schemes. Numerical results are given in Section 5. Finally, the conclusion is presented in Section 6.





§ PROBLEM FORMULATION

 
The multicontinuum models are widely used in reservoir simulations. For example,  in gas production from shale formation, we have a highly heterogeneous and complex mixture of organic matter, inorganic matter, and multiscale fractures <cit.>. Another example is the fractured vuggy reservoirs, where multicontinuum models are used to characterize the complex interaction between vugs, fractures, and porous matrix <cit.>.
The mathematical model of flow in multicontinuum media is described by a coupled system of equations

    c_α∂  u_α/∂ t
    - ∇·  (k_α∇ u_α)
    +  ∑_β≠ασ_αβ (u_α - u_β)  = f_α,     0 < t  ≤ T,    x ∈Ω_α,

where Ω_α is a computational domain, α, β = 1,2,...,L and L is the number of continuum. 
Here 
u_α = u_α(x, t) is the pressure of α continuum, 
c_α and k_α are the problem coefficients, 
σ_αβ = σ_βα  is the coupling coefficient between continuum α and β that characterize a flow between them,  
f_α is the source terms.  

The system of equations (<ref>) is considered with some given initial conditions and homogeneous Neumann boundary conditions for each continuum

    - k_α∇ u_α· n = 0,    0 < t  ≤ T,    x ∈∂Ω_α,

where n is the outer normal vector to the domain boundary ∂Ω_α.


Let u ∈ V and V = V_1 ⊕ V_2 ⊕ ... ⊕ V_L  be a direct sum of spaces V_α, where u_α∈ V_α and  V_α is a Hilbert space.   
Then for u(t) = (u_1, u_2, ...u_L)^T ∈ V, we have the following system of equations

    𝒞du/dt + 𝒜 u = f(t),      0 < t  ≤ T,

with 

    𝒜 = 
    𝒟 + 𝒬, 
      𝒟 =
     [ 𝒟_1   0 ...   0;   0 𝒟_2 ...   0; ... ... ... ...;   0   0 ... 𝒟_L ], 
      𝒞 = 
    [ c_1   0 ...   0;   0 c_2 ...   0; ... ... ... ...;   0   0 ... c_L ],

and 

    𝒬 =  
     [ ∑_β≠ 1σ _1β       -σ_12         ...       -σ_1L;       -σ_21  ∑_β≠ 2σ_2β         ...       -σ_2L;         ...         ...         ...         ...;       -σ_L1       -σ_L2         ...  ∑_β≠ Lσ_Lβ;             ], 
       
    f = [ f_1; f_2; ...; f_L ],

where f_α∈ L^2(0,T;L^2(Ω_α)) and 𝒟_α u_α = - ∇· (k_α∇ u_α) is the diffusion operator for component α.  

Let (u,v) and (u,v)_𝒜 = (𝒜u, v) be the scalar products for u,v ∈ V,  ||u|| = √((u,u)) and ||u||_𝒜 = √( (𝒜u, v)) be the norms in V with (u,v) = ∑_α=1^L (u_α, v_α), u_α∈ L^2(0,T;L^2(Ω_α)). 
For the system of equations (<ref>) that describe the flow problems in multicontinuum media,  we have the following physical properties

    c_α≥ c_0,α > 0,    
    k_α≥ k_0,α > 0,   σ_αβ = σ_βα≥ 0,

with c_α, σ_αβ∈ L^∞(Ω_α) and k_α∈ W^1, ∞(Ω_α). 
Therefore  𝒜 and 𝒞  are self-adjoint and positive definite operators.  





 §.§ Fine-scale approximation by space (FVM)

 
In domain Ω, we construct a structured grid 𝒯_h(Ω)  = ∪_i ς_i, where ς_i is the square cell with length h in each direction. 
For the problem (<ref>), we use a finite volume approximation and obtain the following semi-discrete form

    c_α,i∂ u^h_α, i/∂ t |ς^α_i|
     + ∑_j  T_α,ij  (u^h_α, i - u^h_α, j)
     + ∑_α≠β∑_jσ_αβ,ij (u^h_α, i - u^h_β, j )
     =  f_α,i   |ς^α_i|,   ∀ i = 1, N^α_h,

where 
T_α,ij = k_α |E_ij|/d_ij (|E_ij| is the length of facet between cells ς_i and ς_j, d_ij is the distance between midpoint of cells ς_i and ς_j),
σ_il = σ |E^i_l|/d^i_l if ι_l ∩ς_i ≠ 0 and zero else (|E^i_l| is the length of continuum interface and d^i_l is the distance between continuum), 
N^α_h is the number of cells and α = 1,2,...,L. 

Let u^h = (u^h_1, u^h_2, ...,u^h_L) and 
u^h_α = (u^h_α, 1, u^h_α, 2,..., u^h_α, N^α_H)^T. Then we have the following system of coupled equations  in the matrix form

    M^h ∂ u^h/∂ t + A^h u^h = F^h,

with

    M^h = [ M^h_1     0   ...     0;     0 M^h_2   ...     0;   ...   ...   ...   ...;     0     0   ... M^h_L ],    
    D^h = [ D^h_1     0   ...     0;     0 D^h_2   ...     0;   ...   ...   ...   ...;     0     0   ... D^h_L ],  
    F^h = [ F^h_1; F^h_2;   ...; F^h_L ],

    Q^h = 
    [  ∑_β≠ 1 Q^h_1β        -Q^h_12            ...        -Q^h_1L;        -Q^h_21  ∑_β≠ 2 Q^h_2β            ...        -Q^h_2L;            ...            ...            ...            ...               ;        -Q^h_L1        -Q^h_L2            ... ∑_β≠ L Q^h_L β ],

and 

    A^h = [ A^h_11 A^h_12    ... A^h_1L; A^h_21 A^h_22    ... A^h_2L;    ...    ...    ...    ...; A^h_L1 A^h_L2    ... A^h_LL ]
    = [  D^h_1 + ∑_β≠ 1 Q^h_1β                -Q^h_12                    ...                -Q^h_1L;                -Q^h_21  D^h_2 + ∑_β≠ 2 Q^h_2β                    ...                -Q^h_2L;                    ...                    ...                    ...                    ...                       ;                -Q^h_L1                -Q^h_L2                    ... D^h_L + ∑_β≠ L Q^h_L β ]

where 
M^h_α = {m_α,ij}, 
D^h_α = {a_α,ij}, 
Q^h_αβ = {q_αβ, ij}, 
F^h_α = {f_α,j |ς^α_j|}

    m_α,ij =
    {
    c_α,i |ς^α_i|      i = j, 
    
    0     otherwise
    . ,  
      
    a_α,ij =
    {∑_n ≠ i T_α, in    i = j, 
    
    -T_α, ij    otherwise
    . ,  
       
    q_αβ,ij =
    {σ_αβ,ij   ς^α_i ∩ς^β_j ≠ 0, 
    
    0     otherwise
    . .

We note that under assumptions (<ref>), we see that matrices D^h_α and Q^h_αβ are symmetric and diagonally dominant with non-negative diagonal entries and therefore are positive semidefinite. Then, the block matrices A^h and Q^h are symmetric and positive semidefinite <cit.>. Moreover, it is well-known that the given finite volume method with two-point flux approximation provides a solution with second-order accuracy by space.





 §.§ Multiscale approximation by space (NLMC)


A general multiscale finite element approach for multicontinuum upscaling was presented in <cit.>. The proposed method is based on the Constrained Energy Minimization Generalized Finite Element Method (CEM-GMsFEM) <cit.> and is called the nonlocal multicontinua method (NLMC). The MLMC method preserves the physical meaning of the coarse grid approximation by a particular way of defining constraints in the CEM-GMsFEM.
The system obtained using the NLMC method has the same size as a regular finite volume approximation on a coarse grid but provides an accurate approximation with a significant reduction of the discrete system size.  











Let 𝒯_H(Ω) be the coarse grid with cells K_i, 
K^+_i be an oversampled region for the coarse cell K_i obtained by enlarging K_i by several coarse cell layers,  
K_i^α =  K_i ∩Ω_α and 
K_i^α, + =  K^+_i ∩Ω_α.
We  construct a set of basis functions ψ^i,α = (ψ^i,α _1, ψ^i,α _2,...,ψ^i,α _L) for α-continuum in local domain K_i^+ with ψ^i,α_β∈ K_i^β,+ using the following constrains

    1/|K^β_j|∫_K^β_jψ^i,α_β dx = δ_ijδ_αβ,  
      ∀ K^β_j ∈ K^+_i,

with β = 1,...,L.
Given constraints provide meaning to the coarse scale solution: the local solution has zero means in another continuum except for the one for which it is formulated.

Then, we solve the following constrained energy minimizing problem in the oversampled local domain (K_i^+) using a fine-grid approximation for the coupled system

    [  D_1^K^+_i+ ∑_γ≠ 1 Q^K^+_i_1γ                           ...                   -Q^K^+_i_1L                         C^T_1                           ...                             0;                           ...                           ...                           ...                           ...                           ...                           ...;                   -Q^K^+_i_L1                           ... D^K^+_i_L + ∑_γ≠ L Q^K^+_i_Lγ                             0                           ...                         C^T_L;                           C_1                           ...                             0                             0                           ...                             0;                           ...                           ...                           ...                           ...                           ...                           ...;                             0                           ...                           C_L                             0                           ...                             0;                               ][ ψ^i,α_1;     ...; ψ^i,α_L;     μ_1;     ...;     μ_L;         ] =
    [       0;     ...;       0; F^α,i_1;     ...; F^α,i_L;         ]

with  the zero Dirichlet boundary conditions on ∂ K^+_i for ψ^i,α. 
Here a Lagrange multipliers μ_β = {μ^β_j} are used to impose the constraints with  F_β^α i = {F^α,i_β,j},  where F^α,i_β,j is  related to the K^β_j and  F^α,i_β,j =  δ_ijδ_αβ. 

By combining multiscale basis functions, we obtain the following multiscale space 

    V_H = span{ψ^i,α=(ψ^i,α_1, ψ^i,α_2,...,ψ^i,α_L),   α = 1,L,   
     i = 1,N_c},

and form the projection matrix

    R = [ R_11 R_12  ... R_1L; R_21 R_22  ... R_2L;  ...  ...  ...  ...; R_L1 R_L2  ... R_LL ],
       
    R_αβ = [   ψ^α, 1_β;   ψ^α, 2_β;        ...; ψ^α, N_c_β ].

 
Then, we have the following coarse scale coupled system for u^H = (u^H_1, u^H_2, ...,u^H_L)

    M^H ∂ u^H/∂ t + A^H u^H  = F^H,

where 

    D^H = R D^h R^T = [ D^H_11 D^H_12    ... D^H_1L; D^H_21 D^H_22    ... D^H_2L;    ...    ...    ...    ...; D^H_L1 D^H_L2    ... D^H_LL ],    
    F^H = [ F^H_1; F^H_2;   ...; F^H_L ],

with D^H_αβ = ∑_γ R_αγ D_γ R_βγ^T and F^H_α = R_αα F^h_α.   

Based on the properties of the multiscale basis functions for mass matrices, we have 

    M^H = [ M^H_1     0   ...     0;     0 M^H_2   ...     0;   ...   ...   ...   ...;     0     0   ... M^H_L ],    
    Q^H = [  ∑_γ≠ 1 Q^H_1γ        -Q^H_12            ...        -Q^H_1L;        -Q^H_21  ∑_γ≠ 2 Q^H_2γ            ...        -Q^H_2L;            ...            ...            ...            ...               ;        -Q^H_L1        -Q^H_L2            ... ∑_γ≠ L Q^H_L γ ],

with M^H_α = R_αα M^h_α R_αα^T,  Q^H_αβ = R_αα Q^h_αβ R_ββ^T
and 

    A^H = [ A^H_11 A^H_12    ... A^H_1L; A^H_21 A^H_22    ... A^H_2L;    ...    ...    ...    ...; A^H_L1 A^H_L2    ... A^H_LL ]
    = [  D^H_11 + ∑_γ≠ 1 Q^H_1γ         D^H_12 - Q^H_12                     ...         D^H_1L - Q^H_1L;         D^H_21 - Q^H_21  D^H_22 + ∑_γ≠ 2 Q^H_2γ                     ...         D^H_2L - Q^H_2L;                     ...                     ...                     ...                     ...                        ;         D^H_L1 - Q^H_L1         D^H_L2 - Q^H_L2                     ... D^H_LL + ∑_γ≠ L Q^H_L γ ].


Here for c_α, σ_αβ, f_α = in each coarse cell K_i^α,  the mass matrix, continuum coupling matrix, and right-hand side vector can be directly calculated on the coarse grid and similar to the finite volume approximation on the coarse grid

    M^H_α = {m_α,ij},    
    Q^H_αβ = {q_αβ, ij},    
    F^H_α = {f_α,j |K^α_j|},

    m_α,ij =
    {c_α,i |K^α_i|     i = j, 
    
    0     otherwise
    . ,  
      
    q_αβ,ij =
    {σ_αβ,ij     K^α_i ∩ K^β_j ≠ 0, 
    
    0     otherwise
    .

Therefore, we have M^H=(M^H)^T ≥ 0 and A^H = (A^H)^T ≥ 0 in the NLMC method. 





§ IMPLICIT APPROXIMATION BY TIME (COUPLED SCHEME)


We set u^n = u(t_n), where t_n = n τ,  n=1,2, ... and τ > 0 be the uniform time step size. 
To construct time approximations for (<ref>) and (<ref>) in a general form, we can use the following form of linear multistep scheme

    1/τ M  (∑_j=-1^s-1 a_j u^n-j) 
    + A (∑_j=-1^s-1 c_j u^n-j)  = F, 
       s ≥ 1
       n = 1,2,...

where M and A are linear operators. 
This system is coupled, and the size of the system is N_h = ∑_α N^α_h on the fine grid and N_H = ∑_α N^α_H on the coarse grid.
By Taylor expansion about t_n, we can obtain constraints for time approximation schemes of different order <cit.>.
Note that the time approximation techniques work with the discrete by-space systems and can be applied to fine and coarse-scale systems. 





 §.§ Two-level scheme
 

For two-level  scheme (s=1), we have the following constrains a_-1 + a_0 = 0 and a_-1 = c_-1 + c_0. 
With a_-1 = -a_0 = 1, c_-1 = θ, c_0 = (1-θ), we obtain the weighted scheme (θ-scheme)  <cit.>

    M u^n+1 - u^n/τ 
    + A (θ u^n+1 + (1-θ) u^n)  = F,

where we have the forward and backward Euler schemes for θ=0 and 1.




The solution of the discrete problem (<ref>)  is stable with θ≥ 1/2 and satisfies the following estimate

    ||u^n+1||_A^2 ≤ ||u^n||_A^2 
    + τ/2 ||F||^2_M^-1.



To find a stability estimate for the two-level  scheme, we write the equation (<ref>) as follows

    (M + τθ A )u^n+1 - u^n/τ + A u^n = F.

By scalar multiply to (u^n+1 - u^n)/τ with u^n = (u^n+1 + u^n)/2 - (u^n+1 - u^n)/2, we obtain 

    ( (M + τθ A) u^n+1 - u^n/τ, u^n+1 - u^n/τ) 
     + ( A u^n,  u^n+1 - u^n/τ)
       = ( (M + τθ A) u^n+1 - u^n/τ, u^n+1 - u^n/τ)  
    + ( A u^n+1 + u^n/2,   u^n+1 - u^n/τ) 
    - ( A u^n+1 - u^n/2,   u^n+1 - u^n/τ) 
       = ( (M + τ(θ - 1/2) A ) u^n+1 - u^n/τ, u^n+1 - u^n/τ)  
    + ( A u^n+1 + u^n/2,   u^n+1 - u^n/τ) 
    = (F,  u^n+1 - u^n/τ)

Here for A = A^T ≥ 0, we have ( A (u^n+1 + u^n),  u^n+1 - u^n) = (A u^n+1, u^n+1) - (A u^n, u^n).

By Cauchy inequality, 
we have

    (F,  u^n+1 - u^n/τ) 
    ≤ε u^n+1 - u^n/τ_M^2 
    + 1/4 ε F _M^-1^2.

Therefore the two-level  scheme  is stable with (θ - 1/2) ≥ 0 and stability estimate holds.






 §.§ Three-level scheme
 

For the three-level schemes (s=2), we have the following constrains a_-1 + a_0 + a_1 = 0, a_-1 - a_1 = c_-1 + c_0 + c_1 and (a_-1 + a_1)/2 = c_-1 - c_1. 
With 
a_-1 = μ, a_0 = 1-2 μ, a_1 = μ-1, 
c_-1=μ-1/2+σ, c_0 = 3/2-μ-2σ and c_1 = σ, we obtain <cit.>

    M    ( μu^n+1 - u^n/τ  + (1-μ) u^n - u^n-1/τ)
       + A ( 
    (σ + μ - 1/2) u^n+1 - (2 σ + μ - 3/2) u^n + σ u^n-1)  
    = F.




The solution of the discrete problem (<ref>)  is stable with μ≥ 1/2 and σ≥ (1-μ)/2 and satisfies the following estimate

    1/4||u^n+1 + u^n||_A^2 + ||u^n+1 - u^n||_S^Im^2 
    ≤1/4||u^n + u^n-1||_A^2 + ||u^n - u^n-1||_S^Im^2 
    + τ/2 ||F||^2_M^-1

with S^Im = 1/τ(μ - 1/2) M  
+ ( σ + μ - 1/2) A. 

    
Equation (<ref>) can be written in the following way

    M    ( μu^n+1 - u^n/τ  + (1-μ) u^n - u^n-1/τ)
    + A ( 
    (σ + μ - 1/2) (u^n+1 - u^n) - σ (u^n - u^n-1) + u^n
    )  
        =
    M ( μu^n+1 - 2 u^n + u^n-1/τ  +  u^n - u^n-1/τ)
    + A (  
    σ (u^n+1 - 2u^n + u^n-1)
    + (μ - 1/2) (u^n+1 - u^n)  + u^n
    )  
        =
     1/τ(μ M  + τσ A ) (u^n+1 - 2 u^n + u^n-1)_I
    + 1/τM (u^n - u^n-1)_II
    + (μ - 1/2) A (u^n+1 - u^n)_III  
    + A u^n_IV 
    = F.

Let w^n = u^n - u^n-1,   2 y^n= u^n + u^n-1, then

    w^n+1 - w^n = u^n+1 - 2 u^n + u^n-1,    
    w^n+1 + w^n = u^n+1 - u^n-1 = 2 (y^n+1 - y^n),

    w^n = (w^n+1 + w^n)/2 - (w^n+1 - w^n)/2,    
    w^n+1 = (w^n+1 + w^n)/2 + (w^n+1 - w^n)/2,

and 

    u^n = (y^n+1 + y^n)/2 - (w^n+1 - w^n)/4.

Next, we multiply to w^n+1 + w^n (in a scalar way) and obtain

    I:   ( ( μ M  + τσ A ) (w^n+1 - w^n), w^n+1 + w^n) 
    = ( ( μ M  + τσ A )w^n+1, w^n+1 )  
    - ( ( μ M  + τσ A ) w^n, w^n),

    II:   (M w^n, w^n+1 + w^n ) 
        = 1/2 ( M (w^n+1 + w^n), w^n+1 + w^n)  
    - 1/2 ( M (w^n+1 - w^n), w^n+1 + w^n) 
        = 1/2 ( M (w^n+1 + w^n), w^n+1 + w^n)  
    - 1/2(M w^n+1, w^n+1) + 1/2(M w^n, w^n),

    III:   (A w^n+1, w^n+1 + w^n ) 
        = 1/2 ( A (w^n+1 + w^n), w^n+1 + w^n)  
    + 1/2 ( A (w^n+1 - w^n), w^n+1 + w^n) 
        = 1/2 ( A (w^n+1 + w^n), w^n+1 + w^n)  
    + 1/2 (A w^n+1, w^n+1) - 1/2(A w^n, w^n),

    IV:   (A u^n, w^n+1 + w^n ) 
        = ( A (y^n+1 + y^n), (y^n+1 - y^n))  
    - 1/4 ( A (w^n+1 - w^n), w^n+1 + w^n) 
        = (A y^n+1, y^n+1) - (A y^n, y^n)
    - 1/4(A w^n+1, w^n+1) + 1/4 (A w^n, w^n),

for A=A^T and M=M^T.

Then, we have

    1/2   ( 
    ( 1/τ M + (μ - 1/2) A
    ) (w^n+1 + w^n), w^n+1 + w^n )
        + ( 
    ( μ/τ M  + σ A - 1/2 τ M  + ( μ/2 - 1/4) A - 1/4 A 
    ) w^n+1, w^n+1)  + (A y^n+1, y^n+1)
        -
    ( 
    ( μ/τ M  + σ A - 1/2 τ M  + ( μ/2 - 1/4) A - 1/4 A 
    ) w^n, w^n)  - (A y^n, y^n)= (F, w^n+1 + w^n).

Then using Cauchy inequality (<ref>), we have 

    1/2   (μ - 1/2) ( A (w^n+1 + w^n), w^n+1 + w^n )
       + 
    ( 
    ( 
    1/τ(μ - 1/2) M  
    + ( σ + μ - 1/2) A ) w^n+1, w^n+1)  
    + (A y^n+1, y^n+1)
        -
    ( 
    (1/τ(μ - 1/2) M  
    + ( σ + μ - 1/2) A ) w^n, w^n)  
    - (A y^n, y^n) 
    ≤τ/2 F _M^-1^2.

Let S^Im = 1/τ(μ - 1/2) M  
+ ( σ + μ - 1/2) A and for μ≥ 1/2, we have

    (S^Im w^n+1, w^n+1) - (A y^n+1, y^n+1)
    ≤ 
    (S^Im w^n, w^n) + (A y^n, y^n) 
    +
    τ/2 F _M^-1^2.

Then

    (S^Im (u^n+1 - u^n), u^n+1 - u^n) 
       + 1/4(A (u^n+1 + u^n), u^n+1 + u^n)
    ≤   (S^Im (u^n - u^n-1), u^n - u^n-1) 
    + 1/4(A (u^n + u^n-1), u^n + u^n-1) 
    + τ/2 F _M^-1^2.

Therefore, the three-level  scheme is stable with μ≥ 1/2 and σ≥ (1-μ)/2 and stability estimate holds.






For example, with μ = 1/2, we obtain

    M u^n+1 - u^n-1/2 τ
    + A ( 
    σ u^n+1 + (1 - 2 σ) u^n + σ u^n-1)  = F,

where with σ=0 we have LF (Leap Frog). 
A similar three-level scheme with weights was presented in <cit.> with σ = σ_1 = σ_2

    M u^n+1 - u^n-1/2 τ
    + A ( 
    σ_1 u^n+1 + (1 - σ_1 - σ_2) u^n + σ_2 u^n-1)  = F,

with stability for σ_1 ≥σ_2 and σ_1 + σ_2 > 1/2.   

For μ=1, we have

    M u^n+1 - u^n/τ 
    + A ( 
    ( σ + 1/2) u^n+1 - ( 2 σ - 1/2) u^n + σ u^n-1)  = F,

where with σ=0 we have CN (Crank-Nicolson).

For μ=3/2, we have

    M  3 u^n+1 - 4 u^n + u^n-1/2 τ 
    + A ( 
    (σ + 1) u^n+1 - 2 σ u^n + σ u^n-1)  = F,

where with σ=0 we have BDF2 (Backward Differentiation Formula).

Presented implicit approximation leads to the large coupled system of equations with high contrast properties. To separate equations by contrasts, we decouple the system of equations by constructing an Implicit-Explicit time approximation. The special additive representation of the operator with an explicit time approximation of the coupling term will lead us to the solution of the smaller problems associated with each continuum.  




§ IMPLICIT-EXPLICIT SCHEME (DECOUPLED SCHEME)


Let M be a diagonal matrix and A = A^(1) + A^(2) with A^(1) = diag(A_11, A_22, ..., A_LL) and A^(2) = A - A^(1). 
Such additive representation helps us to separate continuum coupling terms and, by approximating it explicitly, construct decoupled schemes

    1/τ M  (∑_j=-1^s-1 a_j u^n-j)
    +A^(1) (∑_j=-1^s-1 c_j u^n-j)
    _Im 
    + 
    A^(2) (∑_j=0^s-1 d_j u^n-j)
    _Ex 
    = F.
 
Note that we will have additional coupling for the non-diagonal matrix M and an additive representation should also be utilized to matrix M. The properties of the matrix M depend on the approximation method; for example, in the finite element method, we will have a non-diagonal matrix. This work considers finite volume approximation and the NLMC method, which lead to the diagonal matrix M. 

In considered multicontinuum problem  (<ref>), we have  (D-scheme)

    A^(1) = [ A_11    0  ...    0;    0 A_22  ...    0;  ...  ...  ...  ...;    0    0  ... A_LL ],
       
    A^(2) = A - A^(1) = 
     [    0 A_12  ... A_1L; A_21    0  ... A_2L;  ...  ...  ...  ...; A_L1 A_L2  ...    0 ].

We have the following properties for operators 

    A = A^(1) + A^(2)≥ 0,   
    A^(1) - A^(2)≥ 0,

where the property of operator A^(1)-A^(2) follows from the diagonally dominance with non-negative diagonal entries <cit.>. 

Furthermore, the additive representation of the operator A with diagonal operator A^(1) is similar to the Jacobi iterations <cit.>. Therefore, the most straightforward improvement of the method can be made by incorporating the available information of the solution of continuum α in computing solution for continuum α+1 like done in Gauss-Seidel iterations <cit.>. 
We sort the continuum based on their permeability,  k_1 < k_2 <...<k_L, and use the calculated continuum solution in the solution of the following equation. 
Therefore, we can construct the following modified schemes.


  * (L-scheme): We first calculate continuum with smaller permeability

    A^(1) = 
    [ A_11    0  ...    0; A_21 A_22  ...    0;  ...  ...  ...  ...; A_L1 A_L2  ... A_LL ],    
    A^(2) = A - A^(1) = 
    [    0 A_12  ... A_1L;    0    0  ... A_2L;  ...  ...  ...  ...;    0    0  ...    0 ],

where A̅_0 is the lower-triangular matrix and we have a forward substitution. 


  * (U-scheme) We first calculate continuum with larger permeability

    A^(1) =
    [ A_11 A_12  ... A_1L;    0 A_22  ... A_2L;  ...  ...  ...  ...;    0    0  ... A_LL ],    
    A^(2) =  A - A^(1) = 
    [    0    0  ...    0; A_21    0  ...    0;  ...  ...  ...  ...; A_L1 A_L2  ...    0 ],

where A̅_0 is the upper-triangular matrix and we have a backward substitution. 

The resulting system is decoupled, but we utilize information about the calculated solution for the previous continuum without affecting the solution time. 





 §.§ Two-level scheme
 

For two-level  scheme (s=1), we have the following constrains a_-1 + a_0 = 0 and a_-1 = c_-1 + c_0 = d_0. 
With 
a_-1 = -a_0 = 1, 
c_-1 = θ, c_0 = (1-θ) and d_0 = 1, we obtain the weighted  Implicit-Explicit  scheme

    M u^n+1 - u^n/τ 
    + A^(1) (θ u^n+1 + (1-θ) u^n)  
    + A^(2) u^n = F.




The solution of the discrete problem (<ref>) is stable with 

    (θ - 1/2) A^(1)≥1/2 A^(2)

and satisfies the following estimate

    ||u^n+1||_A^2 ≤ ||u^n||_A^2 
    +τ/2 ||F||^2_ M^-1.




 To find stability estimate for the two-level  scheme, we write the equation (<ref>) as follows

    (M + τθ A^(1))u^n+1 - u^n/τ + A u^n = F.

Similarly to the Theorem <ref> for A = A^T ≥ 0, we have

    ( (M + τθ A^(1)) u^n+1 - u^n/τ, u^n+1 - u^n/τ) 
     + ( A u^n,  u^n+1 - u^n/τ)
       = ( (M + τ(θ A^(1) - 1/2 A ) ) u^n+1 - u^n/τ, u^n+1 - u^n/τ)  
    + 1/2 τ ||u^n+1||_A^2
    - 1/2 τ ||u^n||_A^2
    = (F,  u^n+1 - u^n/τ)

Finally, by Cauchy inequality (<ref>) and for ( θ A^(1) - 1/2 A ) ≥ 0, the two-level  scheme  is stable and stability estimate holds.


We note that the scheme is similar to the one presented in <cit.>. However, in this paper we present a general way of the implicit-explicit scheme construction as a multistep implicit-explicit scheme presented in <cit.>. 





 §.§ Three-level scheme
 

For the three-level  schemes (s=2), we have 
a_-1 + a_0 + a_1 = 0, 
a_-1 - a_1 = c_-1 + c_0 + c_1 = = d_0 + d_1 and 
(a_-1 + a_1)/2 =  c_-1 - c_1 = -d_1. 
With 
a_-1 = μ, a_0 = 1-2 μ, a_1 = μ-1, 
c_-1=μ-1/2+σ, c_0 = 3/2-μ-2σ, c_1 = σ, 
d_0 = 1/2+μ, d_1 = 1/2-μ, 
we obtain

    M    ( μu^n+1 - u^n/τ  + (1-μ) u^n - u^n-1/τ)
       + A^(1)( 
    (σ + μ - 1/2) u^n+1 - (2 σ + μ - 3/2) u^n + σ u^n-1) 
       + A^(2)( 
    (μ + 1/2) u^n -  (μ - 1/2) u^n-1)
    = F.





The solution of the discrete problem (<ref>)  is stable with 

    μ≥ 1/2,   ( σ + μ - 1/2) A^(1)≥μ/2 A^(2)

and satisfies the following estimate

    1/4||u^n+1 + u^n||_A^2 + ||u^n+1 - u^n||_S^ImEx^2 
    ≤1/4||u^n + u^n-1||_A^2 + ||u^n - u^n-1||_S^ImEx^2 
    + τ/2 ||F||^2_M^-1

with S^ImEx = 1/τ(μ - 1/2) M  
+ ( σ + μ - 1/2) A^(1) - μ/2 A^(2). 



Equation (<ref>) can be written in the following way

    1/τA_I (u^n+1 - 2 u^n + u^n-1)_I
    + 1/τA_II (u^n - u^n-1)_II
    + (μ - 1/2) 
     A_III(u^n+1 - u^n)_III  
    +  A_IV u^n_IV 
    = F.

with 

    A_I =  ( μ M  + τσ A^(1)),     
    A_II = ( M + τ(μ -1/2) A^(2)),    
    A_III = A^(1),    
    A_IV = A^(1) + A^(2) = A.

Next, we do a scalar multiply to w^n+1 + w^n (w^n = u^n - u^n-1,   2 y^n= u^n + u^n-1) and obtain

    I:   ( A_I (w^n+1 - w^n), w^n+1 + w^n) 
    = (A_I w^n+1, w^n+1 )  
    - (A_I w^n, w^n),  
    
       II:   (A_II w^n, w^n+1 + w^n ) 
     = 1/2 ( A_II (w^n+1 + w^n), w^n+1 + w^n)  
    - 1/2(A_II w^n+1, w^n+1) 
    + 1/2(A_II w^n, w^n),
    
       III:   (A_III w^n+1, w^n+1 + w^n ) 
     = 1/2 (A_III (w^n+1 + w^n), w^n+1 + w^n)  
    + 1/2 (A_III w^n+1, w^n+1)
    - 1/2(A_III w^n, w^n),
    
       IV:   (A_IV  u^n, w^n+1 + w^n ) 
     = (A_IV  y^n+1, y^n+1) 
     - (A_IV y^n, y^n)
    - 1/4(A_IV w^n+1, w^n+1)
    + 1/4 (A_IV w^n, w^n).

Then, we have

    1/2   ( 
    ( 
    1/τ A_II + (μ - 1/2) A_III) (w^n+1 + w^n), w^n+1 + w^n )
        + ( 
    ( 1/τ  A_I - 1/2 τ A_II  
    + ( μ/2 - 1/4) A_III 
    - 1/4 A_IV) w^n+1, w^n+1)  
    + (A_IV y^n+1, y^n+1)
        -
    ( 
    ( 1/τ  A_I  - 1/2 τ A_II  
    + ( μ/2 - 1/4) A_III 
    - 1/4 A_IV) w^n, w^n)  - (A_IV y^n, y^n)
    = (F, w^n+1 + w^n)

Here

    1/τ A_II + (μ - 1/2) A_III  
    = 1/τ( M + τ(μ - 1/2) A^(2))
    + (μ - 1/2)  A^(1)
    = 1/τ M + (μ - 1/2) A,

and 

    S^ImEx   = 
    1/τ  A_I  - 1/2 τ A_II  
    + ( μ/2 - 1/4) A_III 
    - 1/4 A_IV
       = 1/τ( μ M  + τσ A^(1))
     - 1/2 τ( M + τ(μ - 1/2) A^(2))
    + ( μ/2 - 1/4) A^(1) 
    - 1/4 A
    
       = 1/τ( μ - 1/2) M
    + ( σ + μ-1/2) A^(1) 
    - μ/2 A^(2),

then using inequality (<ref>) we obtain a stable scheme for μ≥ 1/2 and 
( σ + μ - 1/2) A^(1)≥μ/2 A^(2).
   

For the particular choice of the parameter μ≥ 1/2, we can obtain the following schemes:


  * μ = 1/2:

    M u^n+1 - u^n-1/2 τ
    + A^(1)( 
    σ u^n+1 + (1 - 2 σ) u^n + σ u^n-1) + A^(2) u^n = F,

where with σ=0 we have LF (Leap Frog), and with σ = 1/2 we have CN (Crank-Nicolson) for the implicit part and LF for the explicit part (CNLF) <cit.>. 
A similar three-level scheme with weights was presented in <cit.> that is stable with (4 σ -1) A^(1)≥ A^(2).  



  * μ=1:

    M u^n+1 - u^n/τ 
    + A^(1)( 
    ( σ + 1/2) u^n+1 - ( 2 σ - 1/2) u^n + σ u^n-1) 
    + A^(2)( 3/2 u^n - 1/2 u^n-1)= F,

where with σ=0, we have CN for the implicit part and AB (Adams-Bashforth) for the explicit part (CNAB) <cit.>.


  * μ=3/2:

    M  3 u^n+1 - 4 u^n + u^n-1/2 τ 
    + A^(1)( 
    (σ + 1) u^n+1 - 2 σ u^n + σ u^n-1) 
    + A^(2) (2 u^n - u^n-1)= F,

where with σ=0 we have BDF2 for implicit part (semi-implicit BDF, SBDF) <cit.>.






§ NUMERICAL RESULTS


We investigate two test problems with two and three-continuum in fractured porous media, where one of the continuums represents hydraulic fracture. Let γ⊂ℝ^d-1 be the lower dimensional domain for fractures and Ω∈ℝ^d be the domain for a porous matrix with d = 2. 
To model fractured porous media, we use a mixed dimensional formulation <cit.>. Let u_α be the α-continuum pressure defined in domain Ω and u_f be the fracture continuum pressure in lower-dimensional domain γ.

We study two cases:


  * Two-continuum media (2C). 

We consider porous matrix continuum (u_m ∈Ω) and fracture continuum (u_f ∈γ)  

    c_m ∂ u_m /∂ t 
    - ∇· (k_m ∇ u_m) +  σ_mf (u_m - u_f) =   f_m, 
        x ∈Ω, 
        c_f ∂ u_f /∂ t 
    - ∇· (k_f ∇ u_f) +  σ_mf (u_f - u_m)=   f_f.
        x ∈γ.

where we set c_f = 1 and k_f = 10^6 for fracture continuum,  c_m = 0.1 and k_m = 1  for porous matrix continuum.  


  * Three-continuum media (3C). 

We have porous matrix continuum (u_1 ∈Ω), natural fracture continuum (u_2 ∈Ω), and hydraulic fracture continuum  (u_f ∈γ)  

    c_1∂ u_1/∂ t
    - ∇· ( k_1∇ u_1)
    +   σ_1 2 (u_1 - u_2)
    +   σ_1 f (u_1 - u_f) = f_1,   
    x ∈Ω, 
        c_2∂ u_2/∂ t
    - ∇· (k_2∇ u_2)
    +   σ_21 (u_2 - u_1)
    +   σ_2f (u_2 - u_f)  = f_2,   
    x ∈Ω, 
        c_f∂  u_f/∂ t
    - ∇· (k_f∇ u_f)
    +   σ_f1 (u_f - u_1)
    +   σ_f2 (u_f - u_2)  = f_f,   
    x ∈γ,

where we set c_f = 1 and k_f = 10^6 for hydraulic fracture continuum,  c_2 = 0.1 and k_2 = 1 for natural fracture continuum,  and  c_1 = 0.05 and k_1 = 10^-3 for porous matrix continuum.  


The problem is considered in domain Ω = [0,2] × [0,1] with 10 fracture lines γ_l, γ = ∪_l=1^10γ_l (Figure <ref>). 
For domain Ω, the fine grid is 400 × 200 grid with N_h = 80,000 square cells (h = 1/200) and the coarse grid is 40 × 20 grid with N_H = 800 square cells (H = 1/20). 
For the lower-dimensional fracture domain γ, the fine grid contains 1,474 cells, and the coarse grid contains 156 cells. 
Lower-dimensional fractures are considered an overlaying continuum and approximated using embedded fracture model (EFM) <cit.>. The fine-scale approach is accurate for a sufficiently fine grid. We construct accurate multiscale basis functions in the NLMC method for coarse-scale approximation to incorporate complex flow behavior on a coarse grid model. 





We set source term in fractures located in the coarse fracture cell related to the domain [1.85, 1.9] × [0.35, 0.4]).  We set f_m=0, f_1=f_2=0 and f_f(x) = q_w(x) (u_f - u_w) for x ∈γ with u_w = 1.2 and q_w = 10^5 and zero otherwise.  
As initial condition, we set u_α,0 = 1 (α = m, 1,2,f) and perform simulations with T_max = 0.005. 





To perform a numerical investigation of the proposed Implicit and Implicit-Explicit schemes for finite-volume and multiscale approximations, we use the solution on the fine grid with N_t = 1024 using coupled scheme (implicit scheme (<ref>) with θ=1) as a reference solution. 
In Figure <ref> and <ref>,  we present a solution for two- and three–continuum media, respectively. The solution u^n is shown at three time layers n=256, 512 and 1024 (N_t = 1024).  
In the coupled scheme, we solve a large coupled system of equations that have DOF_h = 81,474 for a two-continuum problem (2C) and DOF_h = 161,474 for a three-continuum problem (3C).  

Next, we numerically study schemes for various time step sizes N_t = 4, 8, 16, 32 and 64.128 with τ = T_max/N_t. We investigate the influence of the time step size on the method accuracy for coupled and decoupled schemes. 
We consider the following schemes:


  * Coupled schemes:


  * The two-level implicit scheme (<ref>) with θ=1 (Im1).

  * The three-level  implicit schemes (<ref>) with (μ, σ) = (1,0) (Im2-CN,  Crank-Nicholson scheme) and (μ, σ) = (1.5,0) (Im2-BDF, second order backward differentiation formula).



  * Decoupled schemes:


  * The two-level  implicit-explicit scheme (<ref>) with θ=1 (ImEx1).

  * The three-level implicit schemes (<ref>) with (μ, σ) = (1,0) (ImEx2-CNAB,  Crank-Nicholson-Adams-Bashforth scheme) and (μ, σ) = (1.5,0) (ImEx2-SBDF, second order semi-implicit backward differentiation formula).


We also consider three variations of decoupled schemes: L, D, and U-scheme. The schemes are based on the different additive representations of the operator A (see (<ref>), (<ref>) and (<ref>)) and related to the order of computations.  
Implementation is performed using the PETSc library <cit.>, using a conjugate gradient (CG) iterative solver with ILU preconditioner. Simulations are performed on MacBook Pro (2.3 GHz Quad-Core Intel Core i7 with 32 GB 3733 MHz LPDDR4X).





 §.§ Implicit-Explicit for fine-scale system



In this section, we consider a fine grid problem, where a finite-volume approximation is used to construct space approximation. We numerically investigate coupled and decoupled schemes. 








To compare schemes, we calculate L_2 and energy relative errors in percentage  at time t 

    e_h,1(t) = ||u_h(t) - ũ_h(t)||/||u_h(t)||× 100 %, 
      
    e_h,2(t) = ||u_h(t) - ũ_h(t)||_A/||u_h(t)||_A× 100 %,

with 

    ||u|| =  √((u, u)),    
    ||u||_A =  √((Au, u))

where u_h is the reference solution (coupled Im1-scheme with N_t=1024), and ũ is the approximate solution using different schemes (Im1, Im2-CN, Im2-BDF, ImEx1, ImEx2-CNAB and ImEx2-SBDF) and N_t = 4,8,16,32,64 and 128.














In Tables <ref> and <ref>, we present errors for two and three-continuum media, respectively. We show an error in percentage at the final time. We observe the following convergence behavior of the proposed schemes from the result. 

For the three-level  schemes, we have a smaller error. 
For example, in Table <ref> for a two-continuum test problem and in Table <ref> for a three-continuum test problem, we have less than 1% of energy error for N_t = 128 time steps in the two-level scheme, Im1. For three-level  implicit schemes, we have less than 1%  energy error for N_t = 8 and N_t = 16 time steps in Im2-CN and Im2-BDF, respectively. 
In the uncoupled case (D-schemes) for two-continuum case in Table <ref>, we have near 2% of energy error in two-level  decoupled scheme (ImEx1) for N_t = 128 and N_t = 32 for three-level  decoupled scheme (ImEx2-SBDF). 
For the three-continuum case in Table <ref>, we have nearly 3% of energy error in the first-order decoupled scheme (ImEx1) for N_t = 128 and N_t = 32 for second-order decoupled scheme (ImEx2-SBDF). 
Therefore, we can use a 4-times bigger time step size in the second-order schemes.


Next, we compare the coupled and decoupled schemes.
For the two-level  schemes in Table <ref>, we have nearly 1%  of energy error for N_t = 64 time steps in both coupled and decoupled schemes (Im1 and ImEx1-U). Moreover, we observe that the U-decoupling scheme gives better results compared with L and D-schemes in ImEx1. We see that continuum order matters in calculations: (1) It is better to use previous continuum solution information in calculations, and (2) First, calculate the continuum with higher permeability or faster flow. 
For the second-order decoupled scheme, we observe that the errors are smaller for the L-scheme in ImEx2-CNAB and ImEx2-SBDF. However, it is better to use a U-scheme for a smaller time step size. Moreover, we also observe that the D-scheme gives more considerable errors for ImEx2-CNAB with fewer time steps and can lead to a big energy error.
By comparing the coupled and uncoupled schemes of the three-level, we observe a more significant influence on the method accuracy. For example, we can use N_t = 16 in the Im2-BDF coupled scheme to obtain near 1%  of energy error, but in the ImEx2-SBDF decoupled scheme, we should use N_t = 32. 


In Figure <ref>, we depicted the dynamics of the error by time. We show the errors in L_2 and energy norms for N_t = 16 on the first row and for N_t = 64 on the second row. The figures illustrate a good performance of the second-order schemes compared with a first-order for a coupled case of calculations. We observe a well-known zigzag behavior of the Crank-Nicholson scheme at the beginning of simulations. When we use a more significant time step size (N_t=16, first row in Figure <ref>), we observe more minor errors for the Im2-CN scheme compared with the Im2-BDF for both two-continuum and three-continuum test problems. However, for N_t=64, we follow a more minor error for the Im2-BDF scheme. 
Finally, we discuss the computational performance of the presented schemes for the fine grid problem. We show the total time of computations in seconds with an average number of iterations for the preconditioned conjugate gradient method in Tables <ref> and <ref>.  In both tables, time_tot denotes the total time of computations for coupled schemes for given N_t and time_tot^α  is the total time of calculations related to the α-continuum in decoupled schemes, α=m, f for two-continuum test and α=1,2,f for three-continuum test. N_AvIt denotes an average number of iterations in one time layer, and N_AvIt^α is related to the α-continuum in the decoupled schemes.  
For the decoupled schemes, we have the same computational time and number of iterations in L, D, and U-schemes. Therefore, the results in Tables <ref> and <ref> are presented for U-scheme. 
To compare the computational efficiency, we present the reference solution's time and number of iterations.


  * Two-continuum media (2C): N_h = 81474, time_tot = 103.76 sec. and N_AverIt = 73.98  using coupled (Im1) scheme for N_t = 1024.

  * Three-continuum media (3C): N_h = 161474, time_tot = 234.37 sec. and N_AverIt = 73.90  using coupled (Im1) scheme for N_t = 1024.

We observe the same computation time for first- and second-order schemes from the presented results. 
For the coupled schemes for the two-continuum test, we have 12-15 sec for N_t=64 and 20-25 sec for N_t=128. 
In the coupled case on each time iteration, we solve the linear system with N_h = 81474 and N_h=161474 unknowns for two and three-continuum tests, respectively. In the decoupled case, we decoupled the system into systems with smaller sizes that are faster to solve. We have 3-4 times faster calculations for the two-continuum test with nearly 4 sec for N_t=64 and 6-7 sec for N_t=128. In the three-continuum cases, we obtain better results with 46-53 sec for the coupled schemes and 6-10 sec for decoupled implicit-explicit schemes for N_t=128. Then, we have 7.8 times faster calculations for the two-level  scheme and 4.8 times faster calculations for the second order scheme. The computational time is reduced due to the small size of the system in the decoupled case that requires a smaller number of iterations. The implicit-explicit schemes also decouple continuums and significantly affect the number of iterations. Furthermore, we observe that the most time is taken by a first continuum defined in Ω in 2C test and by second continuum p_2 ∈Ω in 3C test. The fracture continuum requires a more significant number of iterations in both test cases. However, the time of calculations is very short due to the size of the fracture continuum system N_h^f = 1474 compared with the continuum defined in Ω with N_h^m = 80,000 for 2C and N_h^1=N_h^2 = 80,000 for 3C.







 §.§ Implicit-Explicit for coarse-scale system


This section considers a coarse grid problem constructed with NLMC approximation. 
Similar to the previous subsection, we numerically investigate coupled and decoupled schemes. In this case, we are more interested in the coarse grid system errors; however, we can reconstruct the fine-scale solution using the constructed projection operator. 
We use a L_2 relative error in percentage at time t  on fine and coarse grids to compare the proposed schemes.

    e_ms,1(t) = ||u_h(t) - u_ms(t)||/||u_h(t)||× 100 %, 
      
    e_H,1(t) = ||u̅_h(t) - u_H(t)||/||u̅_h(t)||× 100 %,

with 

    u̅_h(t) = 1/|K_i|∫_K_i u_h(t) dx

where u_h is the reference solution,  u̅_h is the coarse grid reference solution (average over coarse cell K_i), u_ms is the multiscale solution on the fine grid, and u_H is the multiscale solution on the coarse grid. 







In NLMC approximation, we significantly influence the local domain size used in multiscale basis calculations to the method's accuracy. This study considers a two-continuum test problem to choose an optimal number of oversampling layers in basis calculations. In Table <ref>, we present the influence of the number of oversampling layers K_i^+,l (l=2,3,4 and 5) on the method accuracy for the fine and coarse grid multiscale solutions, u_ms and u_H at the final time. 
In Figure <ref>,  we present fine grid and coarse grid multiscale solutions at the final time on the first and second rows, respectively. Combining the visual representation of the solution from Figure <ref> and relative errors from Table <ref>, we can observe how several oversamplings affect the accuracy. Here, we observe that one cannot obtain a solution for the NLMC method by using one oversampling layer in local domain construction. When we use l=2 oversampling layers, we have a good approximation of the coarse grid solution. However, two oversampling layers are not enough for good fine-grid resolution of the solution. We should use five oversampling layers to obtain smooth results on the fine grid with a small e_ms,1 error. However, because we are more interested in the coarse grid model, we can say that the solution is good with less than one % of error if we use three oversampling layers. We should mention that the number of oversampling layers directly affects the time of multiscale basis calculations. However, the calculations are performed in a parallel way, i.e., for each local domain independently. 






Next, in Figures <ref> and <ref>, we depict the solution of the coarse grid system for two and three-continuum test problems, respectively. We use an NLMC approximation with three oversampling layers on a 40× 20 coarse grid. The resulting unknowns on the coarse grid are N_H = 956 and N_H = 1756 for 2C and 3C tests in the coupled scheme. The results are shown for the coupled scheme using the two-level  coupled scheme (Ms-Im1) with N_t=128 and represented for n=32, 64, and 128.
For the 2C test in Figure <ref>, we have e_H,1 = 0.12, 0.8 and 0.04 for t = T_max/4, T_max/2 and T_max, respectively. 
For the 3C test in Figure <ref>, we have e_H,1 = 0.09, 0.08 and 0.05 for t = T_max/4, T_max/2 and T_max, respectively. 
Solution time is 0.9 sec for 2C and 2.8 for 3C using N_t = 128 on coarse grid with N_H = 956 and N_H = 1756 for 2C and 3C tests, respectively. 
We have a speedy solution compared to the fine grid system using the same scheme (Im1). On the fine grid, we have 25.3 sec for  2C and 53.6 for 3C for N_t = 128. 








 
Tables <ref> and <ref> show multiscale method errors at final time for two and three-continuum test problems, respectively.
We observe the convergence behavior of the proposed Implicit and Implicit-Explicit schemes for multiscale approximation from the result. 

We observe excellent performance on the coarse grid approximation using the NLMC method for the first-order coupled and decoupled schemes. We can use fewer time steps and a coarse grid formulation to produce an accurate solution. For example, we have near 0.5 % of e_H,1 error using coupled Ms-Im1 scheme and decoupled version Ms-imEx1 using only eight time steps for 2C test (see Table <ref>).

However, for both two and three-continuum tests, we observe a lousy performance of the second-order coupled scheme based on the Crank-Nicolson approximation (Ms-Im2-CN), where we have fundamental error e_H,1 = 1.3 % for 2C and  e_H,1 = 3.5 % for 3C for N_t=128. However, we can obtain good results using decoupled implicit-explicit variation (Ms-Im2-CNAB-D) with  e_H,1 = 0.04 % for 2C and  e_H,1 = 0.06 % for 3C for N_t=128. Furthermore, we can obtain a perfect solution using only eight-time steps with less than one percent of e_H,1 error using the Ms-Im2-CNAB-L scheme.

For the coupled Ms-Im2-BDF scheme and decoupled Ms-ImEx2-SBDF schemes, we obtain a little bit better results in decoupled case with near 0.5 % of e_H,1 error using N_t=8 in   Ms-ImEx2-SBDF and using N_t=16 in Ms-Im2-BDF in 2C test. In 3C test, we observe a similar behavior with  near 0.8 % of e_H,1 error using N_t=8 in   Ms-ImEx2-SBDF and using N_t=16 in Ms-Im2-BDF. 
By comparing the L, D, and U versions of schemes, we see that the U-scheme works better for  Ms-ImEx1 and Ms-ImEx2-SBDF schemes. 









Finally, we discuss the computational performance of the presented combination of the time and space approximation schemes, where we have a substantial computational time reduction based on the coarse grid calculations using the NLMC method. Moreover, we can decouple a system and use a more significant time step size using presented Implicit-Explicit time approximation schemes. 
Similarly to the fine grid calculations, we present the total time of computations in seconds with an average number of iterations for the preconditioned conjugate gradient method in Tables <ref> and <ref>.  For the decoupled schemes, the results are presented for U-scheme.
For the coupled first-order scheme (Im1), we have the following:


  * Two-continuum media (2C): 


  * Reference solution with N_h = 81474 (fine grid) with N_t = 1024 (Im1): time_tot = 103.76 sec. and N_AverIt = 73.98. 

  * Fine-scale solution with N_h = 81474 (fine grid) with N_t = 128 (Im1): time_tot = 25.31 sec. and N_AverIt = 146.5. 

  * Multiscale solution with N_h = 956 (coarse grid) with N_t = 128 (Im1): 
time_tot = 0.90 sec. and N_AverIt = 13.4. 


  * Three-continuum media (3C): 


  * Reference solution with N_h = 161474 (fine grid) with N_t = 1024 (Im1): 
time_tot = 234.37 sec. and N_AverIt = 73.90

  * Fine-scale solution with N_h = 161474 (fine grid) with N_t = 128 (Im1): time_tot = 53.68 sec. and N_AverIt = 140.1. 

  * Multiscale solution with N_h = 1756 (coarse grid) with N_t = 128 (Im1): 
time_tot = 2.8 sec. and N_AverIt = 13.1. 


By constructing the accurate coarse-scale approximation, we reduce the size of the system significantly with 28.1 and 19.2 times faster calculations compared with a fine-scale solution with the same time step size (N_t = 128). 


Moreover, from the Tables <ref> and <ref>, we observe a significant influence of the decoupling to the time of computations. 
For N_t = 128, we have  for the two-continuum test (2C): 


  * Fine-scale solution with N_h = 81474 (fine grid): 


  * Coupled scheme (Im2-BDF): time_tot = 20.02 sec. with N_AverIt = 129.7. 

  * Decoupled scheme (ImEx2-SBDF): time_tot = 7.52 sec. with N_AverIt^m = 29.0 and N_AverIt^f = 53.0. 


  * Multiscale solution with N_h = 956 (coarse grid):


  * Coupled scheme (Ms-Im2-BDF): time_tot = 0.75 sec. with N_AverIt = 13.0. 

  * Decoupled scheme (Ms-ImEx2-SBDF): time_tot = 0.15 sec. with N_AverIt^m = 2.0 and N_AverIt^f = 13.2.


In a second-order scheme with (μ, σ) = (1.5, 0), we obtain a 2.7 times faster solution on the fine grid for the decoupled system and five times faster solution on the coarse grid with 0.08 % of e_H,1 error compared with a reference solution.
We also observe a significant reduction in the number of iterations for the continuum in Ω, significantly influencing the solution time. 

For the three-continuum media (3C): 


  * Fine-scale solution with N_h = 161474 (fine grid):


  * Coupled scheme (Im2-BDF): time_tot = 49.03 sec. with N_AverIt = 125.1. 

  * Decoupled scheme (ImEx2-SBDF): time_tot = 10.39 sec. with N_AverIt^1 = 3.0, N_AverIt^1 = 29.0 and N_AverIt^f = 52.0. 


  * Multiscale solution with N_h = 1756 (coarse grid):


  * Coupled scheme (Ms-Im2-BDF): time_tot = 2.60 sec. with N_AverIt = 13.0. 

  * Decoupled scheme (Ms-ImEx2-SBDF): time_tot = 0.27 sec. with N_AverIt^1 = 1.0, N_AverIt^1 = 2.0 and N_AverIt^f = 13.3. 


In the second-order scheme for 3C, we obtain a 4.7 times faster solution on the fine grid for the decoupled system and a 9.6 times faster solution on the coarse grid with 0.1 % of e_H,1 error compared with a reference solution. 

We observe a substantial computational time reduction by combining Implicit-Explicit schemes and accurate multiscale approximation by the NLMC method. We can solve a system in less than a second compared with 103.76 sec for two-continuum and  234.37 sec for three-continuum problems. 




§ CONCLUSION


We presented efficient decoupled schemes for multicontinuum problems in porous media. The decoupled schemes are constructed in a general way using Implicit-explicit time approximation schemes. An additive representation of the operator is used to decouple equations for each continuum. We used a finite-volume method for fine-scale approximation, and the nonlocal multicontinuum (NLMC) method was used to construct an accurate and physically meaningful coarse-scale approximation. The NLMC method is based on defining the macroscale variables on the coarse grid for the general multicontinua problem in porous media. We investigated the stability of the two and three-level time approximation schemes. 
We observed that combining decoupling techniques with multiscale approximation leads to developing an efficient solver for multicontinuum problems. An extensive numerical investigation was given for two and three continuum problems that describe flow processes in fractured porous media. By decoupled calculations, we obtain 3-4 times faster calculations for two-continuum media and 5-8 times faster calculations for three-continuum media than a regular implicit approximation leading to the coupled system on the fine grid. Furthermore, by combining two techniques (nonlocal multicontimuum method and continuum decoupling), the simulation time becomes 140-210 times faster for two-continuum and 180-280 times faster for three continuum media compared with a fine-scale coupled schemes of the same order of time approximation with similar time step size.




plain



