Modeling of flow and transport in complicated porous media in various physical and engineering applications encounters problems with multiscale features. In particular, the properties of the underlying media, such as thermal diffusivity or hydraulic conductivity, have values across different magnitudes. This poses challenges in the numerical simulation since the high contrast feature of the heterogeneous media introduces stiffness for the system. In terms of temporal discretization, the time-stepping depending on the magnitude of the multiscale coefficient is needed for explicit schemes. For the spatial discretization, many multiscale methods [27, 25, 10, 22, 5, 20] are introduced to handle the issue. The multiscale model reduction methods include both local [22, 2, 1, 4, 21] and global [26, 7, 9, 8] approaches to reduce computational expenses. The idea is to construct reduced order models to approximate the full fine-scale model and achieve efficient computation. Among these methodologies, the family of generalized multiscale finite element methods (GMsFEM) [20, 13, 16, 17]
are proposed to effectively address multiscale problem with high-contrast parameters. It first formulates some local problems on coarse grid regions to get snapshot basis that can capture the heterogeneous properties, and then designs appropriate spectral problems to get important modes in the snapshot space. The GMsFEM approach share some similarities with multi-continuum methods. The basis functions can recognize the high-contrast features such as channels that need to be represented individually. The convergence of the GMsFEM depends on the eigenvalue decay, and the small eigenvalues correspond to the high permeable channels.
To construct multiscale basis such that the convergence of the method is independent of the contrast and linearly decreases with respect to mesh size under suitable assumptions, the constraint energy minimizing GMsFEM (CEM-GMsFEM) was initiated[14, 11]. This approach begins with a suitable choice of auxiliary space, where some local spectral problems in coarse blocks are solved. The auxiliary space includes the minimal number of basis functions to identify the essential information of the channelized media. Then it will be used to compute the solutions of constraint energy minimizing problem in some oversampling coarse regions to handle the non-decaying property. The resulting localized solutions form the multiscale space.
To adapt the CEM-GMsFEM for flow-based upscaling, the nonlocal multicontinuum upscaling method (NLMC)  is proposed by modifying the above framework. The idea is to use simplified auxiliary space by assuming that each separate fracture network within a coarse grid block is known. The auxiliary basis are piecewise constants corresponding to fracture networks and matrix, which are called continua. Then the local problems are formulated for each continuum by minimizing the local energy subject to appropriate constraints. This construction returns localized basis functions which can automatically identify each continuum. Further, due to the property of the NLMC basis, this approach will provide non-local transmissibilities which describe the transfer among coarse blocks in an oversampled region and among different continua.
Consider the time-dependent problem with high-contrast coefficients, there have been various approaches to handle multiscale stiff systems [3, 6, 23, 24]. Recently, a temporal splitting method is combined with the spatial multiscale method  to produce a contrast-independent partially explicit time discretization scheme. It splits the solution of the problem into two subspaces which can be computed using implicit and explicit methods, instead of splitting the operator of the equation directly based on physics [32, 33, 34, 28, 29]. The multiscale subspaces are carefully constructed. The dominant basis functions stem from CEM-GMsFEM which have very few degrees of freedom and are treated implicitly. The additional space as a complement will be treated explicitly. It was shown that with the designed spaces, the proposed implicit-explicit scheme is unconditionally stable in the sense that the time step size is independent of the contrast. Following a similar idea in , in this work, we will propose a multirate time-stepping method for the multiscale flow problem.
Multirate time integration method has been studied extensively for the flow problems [18, 19, 30, 31]. Based on different splittings of the target equation, multiple time stepping are utilized in different parts of the spatial domain according to computational cost or complexity of the physics.
In this paper, we integrate the multirate approach with multiscale space construction. Due to the high contrast property of the coefficients, the solutions pass through different regions of the porous medium with different speeds in the flow problem. Different from the previous approach , where the multiscale basis functions are formulated for dominant features and missing information, we propose to design multiscale spaces to handle the fast and slow components of the slow separately. We use the simplified auxiliary space containing piecewise constant functions as in the NLMC framework. We only keep the basis representing the high-diffusive component in the first space and adopt an implicit time discretization scheme. The second space consists of bases representing the remaining component, it will take care of the slow flow and the corresponding degrees of freedom are solved explicitly. Next, we introduce a multirate approach where different time step sizes are employed in the partially explicit splitting scheme, such that different parts of the solution are sought with time steps in line with the dynamics. We start with a coarse step size for both equations and refine local coarse time blocks based on some error estimators. With a finer discretization, the accuracy of the approximation can be improved. We analyze the stability of the multirate methods for all four cases when we use coarse or fine time step size alternatively for the implicit and explicit parts of the splitting scheme. It shows that the scheme is stable as long as the coarse time step size satisfies some suitable conditions independent of the contrast. Moreover, we propose an adaptive algorithm for the implicit-explicit scheme by deriving error estimators based on the residuals. The two error estimators corresponding to the two components of the solutions can provide an upper bound of the errors. Compared with uniform refinement, an adaptive refining algorithm can enhance the efficiency significantly. Several numerical examples are presented to demonstrate the effectiveness of the proposed adaptive method.
The paper is organized as follows. In Section 2, we describe the problem setup and the partially explicit scheme. The construction of the multiscale spaces is discussed in Section 3. In Section 4, the multirate method is presented, the subsection 4.1 is devoted to the stability analysis and the subsection 4.3 presents the adaptive algorithm. Numerical tests are shown in Section 5. A conclusion is drawn in Section 6.
2 Problem Setup
Consider the parabolic equation
where is a heterogeneous coefficient with high contrast, that is, the value of the conductivity/permeability in different regions of can differ in magnitudes.
For a finite element space , the discretization in space leads to seeking such that
Now consider a coarse spatial partition of the computational domain , we will construct suitable multiscale basis functions on and form a multiscale space which is an approximation of the finite element space. Let be the time step size. The discretization in the space with implicit backward Euler scheme in time reads
where is the number of time steps, and . It is well-known that this implicit scheme is unconditionally stable.
Suppose the multiscale space can be decomposed into two subspaces
then a partial explicit temporal splitting scheme  is to find and , for all satisfying
, where is a customized parameter. In the case , the two equations are decoupled, and can be solved simultaneously. In the case , the second equation depends on the solution , thus the two equations will be solved sequentially.
The solution at time step will be . It was shown in  that under appropriate choices of the multiscale spaces and , the above implicit-explicit scheme resulted from the temporal splitting method for multiscale problems are stable with time step independent of contrast. In , the dimension of is low and it contains some dominant multiscale basis functions, the second space includes additional bases representing the missing information. In this paper, we will construct multiscale spaces corresponding to different time scales, where the fast and slow parts of the solution are treated separately.
3 Construction of multiscale spaces
In this section, we will present the construction of multiscale spaces.
3.1 Multiscale basis for
To start with, we introduce some notations for the fine and coarse discretization of the computational domain . Let be a coarse partition with mesh size , and be a conforming refinement of with mesh size , where . Denote by () the set of coarse blocks in , and is an oversampled region with respect to each , where the oversampling part contains a few layers of coarse blocks neighboring . Let be the restriction of on .
Under the framework of CEM-GMsFEM, one first constructs an auxiliary space. Consider the spectral problem
where and are corresponding eigenpairs, and
with , and denotes the multiscale partition of unity function. Upon solving the spectral problem, we arrange the eigenvalues of (4) in an ascending order, and select the first eigenfunctions to form the auxiliary basis functions. Define , where and is the number of coarse elements. The the global auxiliary space . We note that the auxiliary space needs to be chosen appropriately in order to get good approximation results. That is, the first few basis functions corresponding to small eigenvalues (representing all the channels) have to be included in the space.
Define a projection operator satisfying
Then we let , and . Define the null space of to be :
Let the global basis be the solution of the optimization problem
Define . It can be seen that is -orthogonal to , that is
Then the CEM multiscale basis is a localization of , and are also computed using the auxiliary space . The idea is to solve the constraint energy minimization problem in a localized region
where is an auxiliary basis.
The multiscale space is then , it is an approximation to the global space .
Note that the construction of CEM basis which we have presented here is general and can handle complex heterogeneous permeability field with high contrast. In a fractured media where the media configurations are explicitly known, here the assumption is valid in many real applications, we can consider the simplified construction.
The domain for the media with fracture networks can be represented as follows
where the subscripts and denote the matrix and fractures correspondingly. In the fracture regions , the scalar denotes the aperture, and is the number of discrete fracture networks. The permeabilities of matrix and fractures usually have high contrast.
In this setting, some simplified auxiliary basis can be adopted, the constraint energy minimizing basis can be constructed via NLMC  and the resulting basis functions which can separate the continua such as matrix and fracture automatically. To be specific, for a given coarse block, we use constants for each separate fracture network, and then a constant for the matrix for the simplified auxiliary space.
Consider an oversampled region , for a coarse element , let be the set of discrete fractures, and be the number of elements in . The NLMC basis are obtained by solving the following localized constraint energy minimizing problem on the fine grid
We remark that the resulting basis separates the matrix and fractures automatically, and have spatial decay property. The NLMC basis for fractured media is then .
We denote the the average of the NLMC basis representing the fractures by
, where . Note that in this double summation, the subscript starts from , which indicates the basis functions corresponding to the matrix are not included here. Let . To simplify the notation, we omit the double scripts in and denote the set of basis by
Finally, we define the space as follows:
where we take away the last basis to remove linear dependency. We remark that, contains basis representing the high contrast regions only. The basis functions corresponding to the matrix and the basis will be included in the second subspace in the following section.
3.2 Multiscale basis for
In this section, we will present basis construction for . We first include which are computed from equations (6), and the average basis obtained from equation (7) in . Additionally, we can also conduct some spectral decomposition in the local coarse region to construct more basis corresponding to the information in the non-high-contrast regions. Here be a coarse neighborhood with respect to the -th coarse node. The basis in the spaces needs to satisfy the stability condition (16) and (20). The additional basis are obtained through finding the eigenvalues and corresponding eigenfunctions from the spectral problem
for all .
Again, rearranging eigenvalues in an ascending order, and select the first eigenfunctions correspondingly, we get
The space is then defined to be
We remark that, for simplicity, we only use for the second space in our numerical examples.
4 Multirate time stepping for partially explicit scheme
Based on the multiscale spaces constructed in Section 3, we introduce a multirate time stepping partially explicit temporal splitting scheme. Consider the coarse time step size and fine time step size , where . Denote by the fine partition of the time domain by
The coarse partition of the time domain is formed by
Further, we write each coarse time interval where .
The multirate scheme is then defined as follows. For each coarse interval , we are seeking for with a given , such that defined in one of the following four cases which are using coarse time step size in both (2) and (3) (coarse-coarse), using coarse time step size in (2) and using fine time step size in (3) (coarse-fine), using coarse time step size in (2) and using fine time step size in (3) (fine-coarse), using fine time step size in (2) and using fine time step size in both (2) and (3) (fine-fine).
We remark that
may not be defined in some of the fine time steps. In this case, we use the linear interpolation of the nearest two coarse time step solutionsto define intermediate time step solutions for .
4.1 Stability for different cases
Consider a coarse time block , the stability of the multirate method for the above mentioned four cases is proved in this subsection.
Let be a constant such that
and for case 1, for case 4.
We will show the stability for case 2 and 3 in the following.
4.1.1 Stability for case 2
Use the coarse time step for and use the fine time step for ,
The multirate partially explicit scheme in (12) satisfies the stability estimate
The equations in (12) can be written as
For the left hand side of (18), we have