Generalized Multiscale Finite Element Method and Balanced Truncation for Parameter-Dependent Parabolic Problems

: We propose a generalized multiscale finite element method combined with a balanced truncation to solve a parameter-dependent parabolic problem. As an updated version of the standard multiscale method, the generalized multiscale method contains the necessary eigenvalue computation, in which the enriched multiscale basis functions are picked up from a snapshot space on users’ demand. Based upon the generalized multiscale simulation on the coarse scale, the balanced truncation is applied to solve its Lyapunov equations on the reduced scale for further savings while ensuring high accuracy. A θ -implicit scheme is utilized for the fully discretization process. Finally, numerical results validate the uniform stability and robustness of our proposed method.


Introduction
Many scientific and engineering problems in porous media flow and composite material involve inherently multiscale model equations.Basically, the model coefficient often exhibits high-contrast and parameter-dependent behaviors in a sense.Conventional methods such as the finite difference method, or finite element method, which are applied to solve these problems on a fully-resolved scale would become unbearably expensive, especially when a large amount of computations may be consumed for uncertainty parameter quantification, model reduction and optimization.So an urgent demand for developing multiscale model reductions with properties of both accuracy and efficiency is of particular interest.
In the research history, the multiscale finite element method (MsFEM) was first proposed by Hou and Chen in [1,2].Then it evolved to a generalized multiscale version and had many applications in scientific fields, see [3][4][5][6][7].For example, in [3] the generalized multiscale finite element (GMsFEM) approximated the solution space locally by employing the enriched multiscale basis functions.This approximation constructed a proper snapshot space and a local spectral decomposition to realize an efficient model reduction.Chung, Efendiev and Hou proposed an adaptive multiscale model reduction with GMsFEM in [5], an adaptively local approach was utilized to update the global snapshot when online computations proceeded.The local-global couple is shown to be effective in many applications.As for the parameter-dependent model and the thermoelastic model with phase flow, the strategy of multiscale model reduction was addressed in [8][9][10][11][12].Besides, there were other techniques such as multiscale reduced-basis method [8,10,13,14], localized orthogonal decomposition [15], weak Galerkin GMsFEM [16], mixed GMsFEM [17], multiscale-spectral GFEM [18], meshfree GMsFEM [19], higher-order three-scale method [20] and local discontinuous Galerkin method [21], were applied in a large number of interested fields.
In literature [22][23][24], a balanced truncation (BT) is directly addressed to their problems.The aim of BT is to study the principal component analyses in linear systems and in the time or frequency domains, and it identifies the weakly controllable and observable components which are used for model reductions.Moreover, the eigenvalues of the Lyapunov solutions in frequency-limited BT often decay faster than those in the standard BT.However, as for the parameter-dependent parabolic problem, it would be annoyingly costing or infeasible since Lyapunov equations are hard to solve.In particular, solving these equations from the high-contrast permeability, boundary condition, and forcing term on the fine scale would quickly become prohibitive.Therefore, in this work we propose the joint usages of a generalized multiscale and balanced truncation to obtain accurate and economical results.
The paper is characterized by the following highlights: As for the parameter-dependent parabolic problem, we present a combination strategy of the generalized multiscale finite element method with the balanced truncation, namely GMsFEM-BT.The GMsFEM is skilled in the random-parameterized model and different coefficients and boundary conditions, it constructs an enriched projection matrix in the global reduction from the snapshot space.And it offers a good balance for the accuracy and the efficiency, which is compared with the fully-resolved FEM and the coarse version of MsFEM.Based on this, the BT is addressed on the correspondingly reduced scale for further reduction.Then moving from a current time step to the next time step, an effective implicit full-discretization process is utilized for the spatial-temporal problem.The uniform stability and convergent results are to be validated in our work.
The paper proceeds as follows.In Section 2, a preliminary description of the parameterdependent parabolic problem is provided and we brief the idea of the local-global multiscale reduction.In Section 3, the GMsFEM is applied to the spatial variables through an eigenvalue computation in the spectral sub-problems, and an offline-online coupling is accomplished by a projection matrix.In Section 4, the detail of BT is offered for further reduction for the model.We present a fully discretization process both for time and for space in Section 5.In Section 6 numerical experiments validate the ability of the combined GMsFEM-BT strategy.Finally, the conclusion is presented in Section 7.

Parameter-Dependent Parabolic Model
We investigate a parameter-dependent parabolic problem in the two-dimensional heterogeneous media, where κ(x; µ) is a high-contrast coefficient with spatial variables x in the two-dimensional square domain Ω, µ is used to represent the parameter dependence and may be random.
T is the time termination, f (x, t) and g(x, t) are given functions, and u 0 (x) is an initial function at t = 0.The numerical solution with the most accuracy, efficiency and stability to u(x, t) is our target in this work.
In the finite-dimensional discretization, the spatial domain Ω is divided into many quadrilateral finite elements on coarse grids and then divided on fine grids, which is denoted by T H and T h , respectively.For a clearer illustration see Figure 1, firstly we divide the domain on coarse grids, in this way coarse elements K are formed.Denote the coarse nodes by {y i } N c i (N c is the number of coarse nodes), and denote the coarse neighborhood of y i by In the following, it will be applied in our generalized multiscale finite element method.Secondly, on the coarse grids T H we refine every coarse element K into more fine elements k, which is denoted by T h on those fine grids.In this paper, we denote the coarse and the fine (of every coarse element) partition number in each spatial direction by N and M, thus H(= 1 N ) and h(= H M ) are used to denote the coarse mesh size and the fine one, respectively.It is known that the classical finite element method is manipulated on very fine grids T h to achieve its maximum accuracy.The functional space is spanned by the piecewise bilinear polynomials Q 1 (k), which are subordinated to T h as In Galerkin weak form, it is applied to find u h ∈ V h such that where As for the initial-boundary value parabolic problem (1), the above (4) produces a discrete system to solve where is a right-side vector, and ψ i is the bilinear basis function in the 2D problem.As a consequence, in the finite element scheme with a very large partition number, a huge-size discrete system is obliged to solve where The subscript f means that the matrices and vectors are processed from all of the fine elements, which is different from the subscript c/r from the coarse/reduced elements in the following situations.For instance, the stiffness matrix A f (µ) is assembled from all of the fine elements, whose size is N f × N f , where N f = (NM + 1) 2 for the 2D spatial discretization.In the paper, we denote the nodal matrices by B, C (which are transposed to each other B = C T , and the size of B is N f × N 2 ), which are allocated to the domain boundary vertices.

Local-Global Multiscale Reduction
In this subsection, we begin with the process of semi-discretization in the space.Different from the above-mentioned finite element method (5) on very fine grids, we present two types of local-global multiscale reductions on coarse grids by contrast.
First, we recall the standard multiscale finite element method (see [1,2,[25][26][27]).A homogeneous sub-problem is proposed for solving the multiscale basis functions ϕ i in every quadrilateral coarse element K as where the boundary condition l i satisfies that δ j (y It is known that useful local information is delivered through the multiscale basis functions into a global system.However, we also know that this type of standard multiscale scheme has its deficiencies, such as the relatively low accuracy, slow convergence even divergence in former literature.Second, we track the generalized multiscale finite element method (see [3][4][5]).Another version of the non-homogeneous local problem in every coarse neighborhood is solved for the (updated) multiscale basis functions, via the necessary eigenvalue computation.In this way, plenty of parameter-dependent details are captured through the over-sampling technique and this would be solved via a singular value decomposition.
Thus in contrast to the finite element method in the fine version (6), the generalized multiscale method in the coarse version is In this multiscale scheme, the size of the stiffness matrix A c (µ) is just N c × N c , where N c = (N + 1) 2 to realize a preliminary computation reduction.It also builds the counterpart of nodal matrices, which are denoted by B c , C c with B c = RB, C c = CR T .The projection matrix R with sufficient microscopic information is constructed from the multiscale basis functions, which would be provided in the following Formula (21).
If it is necessary for a further decline in computational costs, we could carry out a further matrix reduction with a strategy of balanced truncation in Section 4, which is in the reduced-version like The size of A r is N r × N r , and the number N r (<< N c << N f ) can be artificially chosen by users later on.Singular value decompositions would be employed to form the matrices U, V in (29) with UV = I N r .Set B r = UB c , C r = C c V, thus the matrices are further reduced to save consumptions.At the same time, accuracy and efficiency are granted in the novel method.
In a word, we are constantly motivated to chase a better cost-effective balance for the complicated multiscale parameter-dependent parabolic model.

Mechanism of GMsFEM
The GMsFEM is an extended framework that generalizes the MsFEM.Its key mechanism relies on that the multiple basis functions in every coarse neighborhood could be systematically updated, by flexibly enriching the functional space through necessary eigenvalue computations.
The GMsFEM has two stages: offline and online.At the offline stage, firstly a snapshot space is created then based on it a small dimensional construction is required for solving local spectral decompositions.One of the advantages of the GMsFEM is that the created snapshot space can be re-used to capture the necessary fine-scale features.After that, a space reduction is performed and the dominant modes (including sufficient eigenvalues and eigenvectors) are selected as the multiscale basis functions, for any input parameters which are fitting our parameter-dependent problem ideally.At the online stage, when the corresponding boundary condition and force term are defined, the available offline basis functions are utilized to pursue a high-precision numerical solution with the potential least costs.Users can also adaptively determine the various number of eigenfunctions needed in the coarse neighborhoods, to fulfill a proper balance of efficiency and accuracy.
The mathematical theories behind the offline and online spaces are the construction of snapshot space and the design of spectral sub-problems, we will present the details in the following subsections in turn.Since plenty of parameter-dependent multiscale problems may have small eigenvalues, we would like to comprise enough of the eigenvectors associated with the small eigenvalues to be enriched into the coarse functional space surely.As a consequence, the multiscale basis functions could be achieved by multiplying the dominant eigenmodes by the partition of unity functions.

Offline Stage
At this stage, we first create a snapshot space subordinated to each coarse neighborhood ω i in (2).The construction of snapshot space involves solving the local problems from a set of parameter-dependent inputs.Being different from the standard multiscale finite element method solving the local problem (7), an eigenvalue computation is utilized in this generalized multiscale finite element method.
First, we solve for the snapshot functions from the spectral sub-problems: seeking (λ where l = 1, • • • , L i (L i is the chosen number of eigenfunctions from small to large that correspond to the dominant eigenvalues).Note that we define zero Neumann boundary condition for (10), except in the discontinuous case the Dirichlet condition is defined instead.It is known that solving all the eigenfunctions of a large matrix is exhausting, so we choose the L i eigenfunctions to construct a snapshot space.The snapshot space V ω i ,snap is associated to the neighborhood ω i , which is produced from the multiscale eigenfunctions as Next, the snapshot functions are reordered by using a single index as the columns to create a projection matrix in a snapshot where L snap denotes the number of functions in the snapshot construction.It should be pointed out that L snap is not only dependent on the chosen number L i , but also it is related to the number of parameter µ.For brevity, we omit the superscript ω i to brief as V snap , R snap .However, keep in mind that throughout the paper, the offline and online eigenvalue decompositions are localized for the corresponding coarse neighborhoods.
Then, we perform an elementary dimension reduction of the snapshot space by applying an auxiliary spectral decomposition to build an offline space.We seek a subspace of the original snapshot space so that it may properly approximate any neighborhood.
The key is to efficiently build a set of multiscale basis functions for each parameter µ.We reiterate the offline bilinear forms to be parameter-independent, then there is no need to rebuild the offline space for each value µ.Then the eigenvalue problem can be described as where where κ(x; µ) and κ(x; µ) are averaged parameter-dependent coefficients.A and M denote similar fine scale matrices as in (13), except that averaged coefficients of them are used.The smallest number of eigenvalues L off are picked up in (14) and form the corresponding eigenvectors as , where Φ off l,j are the coordinates of the vector Φ off l .In this way, the offline projection matrix is constructed as

Online Stage
With the offline space available, next we build an online space on every coarse neighborhood for a specified parameter µ.For the computational efficiency, it is assumed that the online space is a subspace with small dimensions.At the online stage, the reducedorder bilinear forms are chosen to be parameter-dependent.The eigenvalue problem is ready to solve where where κ(x; µ) and κ(x; µ) now rely on the specified value µ, A(µ) and M(µ) are specified fine-scale matrices in (13).
To construct the online space we pick the smallest number of eigenvalues L on from ( 16) and finalize the corresponding eigenvectors as , where Φ on l,j are the coordinates of the vector Φ on l .It should be noted that the matrices size in ( 16) solely relies on the dimension of the offline subspace.This results in the declining consumption of the online space, which in contrast to the original snapshot space would cost more expenses.

Offline-Online Coupling
In order to enrich the online basis functions into a reduced formula of (1), we denote χ i as the multiscale partition of unity functions which are from another local problem where g i is a continuous and bilinear boundary condition.A summed pointwise energy is defined as As above, H is the coarse mesh size and N c is the number of coarse nodes.To construct the global coarse-grid solution space, the partition of unity functions are multiplied by the online eigenfunctions to obtain the final multiscale bases Then it follows the online space which is quite different from the large-scale finite dimensional space in (3).Using a single index notation we denote V on = span{ϕ k } N c k=1 .Recalling in (8) now we finalize the projection matrix R to be utilized in the global reduction where ) represents the nodal vector of the generalized multiscale bases.The matrix R (whose size is N c × N f ) is now plugged in (8), which could be enriched with the dominant eigenfunctions and meaningful microscopic behaviors.This would lead to strengthening the numerical abilities later on.Now it completes the coarse model reduction in the GMsFEM.In the next section, a further model reduction will be carried on a balanced truncation (BT).

Related Lyapunov Equations
First recalling that for the standard Lyapunov equation where the matrix W is dense even though the system is sparse, and generally the eigenvalues of W decay quickly when the right-side Z has a low rank.This decay is typically associated with the approximation error from the balanced truncation.Next, a joint measure of controllability and observability is offered for solving Lyapunov equations.It should be pointed out that these equations are not based on the above-mentioned fine matrices, but instead, they are based on the coarse matrices A c (µ) defined in (8) and B c , C c below that.Then the related Lyapunov equations are in the form of where W con is a controllability matrix, W obs is an observability one, and Z is dependent on the B c , C c .The solutions of these Lyapunov equations play a key role in the model reduction and in the balanced truncation particularly.Generally, the computations of dense matrices are of order O(n 3 ).Some techniques have been evolved to efficiently compute the solutions by reducing the complexity.However, most of these techniques are often iterative and low-rank approximations, while in our high-contrast situation, it is hard to provide enough satisfactory approximations.In Matlab code lyap, the matrices are transformed into the complex Schur forms, and it computes the resulting triangular systems.In this way, the high-rank approximations are available and the solutions of Lyapunov equations have a close relationship to the dominant eigenvalues.Then are used to realize the numerical advantages of the coupled Lyapunov equations.Since the matrices A c (µ), B c , C c are based on the coarse-scale, then the sizes of W con , W obs are just N c × N c .
Once we obtain the two matrices W con , W obs , which may be expressed in the form of Cholesky singular value decompositions to get L con , L obs as and where U, V are to be specified soon in (29).Set , and the subscript r means the reduced case.Note that our BT is based on the GMsFEM.An independent set of online basis functions are constructed in which the coarse-scale solution may be sought.As a consequence, the global system would be processed on a coarse scale to avoid applying the BT on a fine scale directly.After the offline and online stages in the GMsFEM, it produces a subspace with small dimensions for each fixed parameter µ.It serves to further diminish the computational cost, which is associated with the balanced truncation in the model reductions.

Further Reduction
Based on coarse matrices (compared with the fine ones), now we have the platform to fulfill the potential further reduction.In our work, we use the first N r ordered eigenvalues and the normalized eigenvectors of the matrix L T con W obs L con , by solving The sets of eigenvalues and eigenvectors are given by {λ i }, and {ζ where is the singular values.The balanced truncation is validated to be more efficient in the numerical simulation.Since their a-priori error estimates rely on the quick decay of eigenvalues in (27) and the magnitude of all truncated eigenvalues summation.On the other hand, the balanced truncation has its weaknesses too.The Lyapunov Equations ( 23) and ( 24) require expensive computations for unknown W con and W obs , and the singular value decompositions are costing too, especially for the very fine case.
With the eigenvalues and corresponding eigenvectors in place, the above-mentioned matrices U and V are constructed from the reduced matrices in (9) and B r , C r below that as whose size is N r × N c and N c × N r , respectively.And the property of UV = I N r is ensured.As a consequence, the reduced version is A r (µ) = U A c (µ)V, where r depends on the eigenvalue decay and is determined by users in the numerical experiment.Now we achieve the reduced scale since A r (µ) has a size of N r ×N r , compared with the coarse scale A c (µ) with a size of N c × N c , while the fine-scale A f (µ) with a size of N f × N f , whose N r << N c << N f .

Fully Discretization for Time and Space
As for the parabolic problem, we are now in the procedure of fully discretization.A θ-implicit scheme is presented to the original Formula (5), we have We know the influence of the temporal parameter θ on the fully discretization.In the following numerical experiments, we set θ = 1, namely the backward Euler difference is applied.
Solving the above ordinary differential equations, the corresponding numerical results are obtained from different schemes.Thus in the paper, a proposed combination of GMsFEM-BT approach for the parameterdependent parabolic problem is evolved for both the spatial and the temporal discretizations.It avoids a huge system computation, while maintaining an ideal balance between accuracy and efficiency at the same time.The generalized multiscale method is addressed in an offline-online coarse scale system, which is further reduced via using the balanced truncation, whose computational cost would be saved greatly.

Numerical Experiment
Now numerical results are ready to reveal the strength of our proposed method.Matlab codes are executed on a Windows desktop with Intel(R) Core(TM) i9-10900K CPU 3.70 GHz and 32 GB RAM.We set the right side f = 1, so there is no exact solution or analytical solution for this random-parameterized parabolic problem.To solve the problem (1), first we divide the square domain Ω = [0, 1]×[0, 1] into N M × N M = 200 × 200 fine elements and use the classical FEM to solve on space and the FDM to solve on time for (6).This approach is the most expensive, and it serves as a reference solution.
A random-parameterized coefficient is combined like where the parameter µ = (µ 1 , µ 2 , µ 3 , µ 4 ) may be randomly picked from [0, 1] 4 , and r = (x − 1 2 ) 2 + (y − 1 2 ) 2 , and set ε 1 = 0.2, ε 2 = 0.08, ε 3 = 0.125, ε 4 = 0.0078125, ε 5 = 0.012, ε 6 = 0.004, ε 7 = 0.8.A resulting permeability field expresses the structure with high contrast and quasi-periodicity, see Figure 2. The homogeneous boundary condition g = 0 and the initial function u 0 = 0 are defined for simplicity.In the coarse-scale schemes, the MsFEM and the updated GMsFEM for (8) are applied to obtain the coarse-scale outputs.In both schemes, we set the double refinement of the coarse partition N = 5, 10, 20 in turn and fix the fine partition M = 10 in each coarse element.However the latter GMsFEM, with the contributions of eigenvalue computation and enriched multiscale basis functions, would offer more accurate simulations.Finally, we utilize the BT on the reduced scale for (9) to test the behaviors of our GMsFEM-BT approach.
To be more specifical, in the GMsFEM to create the snapshot space first for the offlineonline stages, three equally spaced points in each dimension are set and the eigenfunctions number L i = 10 is kept in (11), then the number of snapshot functions is L snap = 3 j × L i = 810 in (12), where j = 4 is the number of random parameter µ j .In each coarse neighborhood, the online-stage eigenfunctions number L on in (20) is smaller than the original L snap .Through a further reduction in the BT, after solving the Lyapunov Equations in ( 23) and ( 24) for a controllability W con and observability W obs on the coarse scale, we use the singular value decompositions to build U, V in (29), whose size is dependent on a smaller number N r on the reduced-scale.
We are quite interested in comparing the simulation behaviors from the fully-resolved fine system, the coarse system and the reduced system.We set the refinement of N = 5, 10, 20, and the fine partition is kept as M = 10 in each coarse element, then on both x, y directions in the 2D problem the number of coarse vertices is N c = (N + 1) 2 = 36, 121, 441.Then the matrix size of the fine system is N f = (NM + 1) 2 = 2601, 10,201, 40,401, while the corresponding matrix size of the generalized coarse system is 207, 988, 4293, which is much smaller than N f and is a bit larger than N c .In contrast, the matrix size of a reduced system is N r , which would be chosen as smaller integers in the following.
Figure 3 shows the discrete errors of MsFEM and GMsFEM, both are compared to the reference solution.Note that the reference is solved from a fully-resolved fine system on N M × N M = 200 × 200, while the solutions of MsFEM and GMsFEM are obtained from their coarse system.It is certified that with the contributions of eigenvalue computation for the dominant eigenfunctions and multiscale enrichments, the GMsFEM performs better than the MsFEM.We list the norm errors from the MsFEM and the GMsFEM with the refinement of the partition number in Table 1 and Figure    In Figure 5, a series of outputs are provided for the reference on the very fine partition N M = 200, the GMsFEM on the quite coarse partition N = 10, and the BT on the several reductions N r .In Table 2 we observe that on small reductions N r such as 20, the GMsFEM-BT does not behave well.However, with the increase of N r it shows the ability to deliver more necessary details of the reference, and its error percentage is less than 0.04% in the end.Note that the reduced size of GMsFEM-BT on N r = 260, which is compared to the generalized coarse size of GMsFEM with N c = 988 and the fine size of FEM with N f = 40,401.In the meantime, in fact the GMsFEM-BT works based on the GMsFEM's coarse scale, and it would take some time to solve the Lyapunov equations and the singular value decompositions.However on the contrary, if the GMsFEM-BT works based on the FEM's fully-resolved fine scale, the computational cost would increase inevitably.
Figure 6 shows the error percentages between the fine and coarse, and between the coarse and reduced, respectively.We find that with time goes by, it is reasonable that the error percentages between the reference and the GMsFEM with a coarse partition number N = 20 (which is fixed during the process) would increase with the time step, but the error remains at a relatively low level overall (less than 0.16%).However, it is exciting to observe that the error percentages between the generalized coarse and the reduced N r = 340 would decrease and remain stable at 0.2% during the time elapses, which validates the strength of GMsFEM-BT to simulate the real application with costing just an amount of computational resource.
Finally, we provide the CPU time among different methods based on their corresponding matrix sizes.In Table 3 it is apparent that the FEM consumes the maximum time on a fully fine matrix size of N f = 40,401 to solve for the reference solution.The MsFEM (on a coarse size of N c = 441) and the GMsFEM (on a coarse size of N c = 4293) spend less time to reach their numerical results.However, note that the GMsFEM has to construct the snapshot space (which could be reused) and couple the offline-online stages with the singular value decompositions.This process consumes some time, so three parts of time (for offline, online and solving stages, respectively) are shown in the table for GMsFEM, but it would bring good accuracy.A Lyapunov solver based on the GMsFEM with the generalized coarse size of N c = 4393 offers its strength.After that, the reduced matrix size of BT is chosen as the user's demand, and the satisfactory balance both for accuracy and efficiency of the GMsFEM-BT is shown in the numerical results.

Conclusions
In the paper, a parameter-dependent parabolic model is addressed by a combination strategy of GMsFEM-BT.Neither the FEM works on the fully-resolved fine scale, nor the MsFEM works on the coarse scale, the GMsFEM acts on the updated coarse scale by introducing the necessary eigenvalue computation to create the snapshot space and couple the offline-online stages.In this way, the selected dominant eigenfunctions are plugged into the projection matrix.The balanced truncation is then presented to solve the Lyapunov equations, whose results are served in the reduced matrices for further reduction.The fully discretization for time and space is applied to the parabolic problem.In the numerical experiments, the capabilities of the GMsFEM-BT are demonstrated convincingly.Its simulation is quite close to the reference result, and the error percentages and CPU time of the GMsFEM-BT are encouraging.Our method could be extended to other random parameter cases in the multiscale simulations.

Figure 1 .
Figure 1.Illustration of coarse grids, elements and neighborhoods.
4. The error percentages indicatethe fast convergence and stability of GMsFEM, which is much superior to the behavior of MsFEM.

Figure 6 .
Figure 6.With time elapsing, error percentages between fine and coarse N = 20 (left), and between coarse and reduced number N r = 340 (right).

Table 1 .
Percentages of L 2 and H 1 norm errors from MsFEM and GMsFEM on log-log scale.Percentages of norm error from MsFEM and GMsFEM, on different partition numbers N.

Table 3 .
Comparisons of CPU time.