Generalized Multiscale Finite Element Method for discrete network (graph) models
Maria Vasilyeva








Generalized Multiscale Finite Element Method 

for discrete network (graph) models
    
Maria VasilyevaDepartment of Mathematics & Statistics, Texas A& M University - Corpus Christi, Corpus Christi, TX, USA.  
Email: maria.vasilyeva@tamucc.edu. 

    April 29, 2024
===================================================================================================================================================================




In this paper, we consider a time-dependent discrete network model with highly varying connectivity.  The approximation by time is performed using an implicit scheme.  We propose the coarse scale approximation construction of network models based on the Generalized Multiscale Finite Element Method.  An accurate coarse-scale approximation is generated by solving local spectral problems in sub-networks.  
Convergence analysis of the proposed method is presented for semi-discrete and discrete network models.  We establish the stability of the multiscale discrete network.  
Numerical results are presented for structured and random heterogeneous networks.   






§ INTRODUCTION



Multiphysics models on large networks are used in many applications, for example, pore network models in reservoir simulation, epidemiological models of disease spread, ecological models on multispecies interaction, medical applications such as multiscale multidimensional simulations of blood flow, fibrous materials, electric power systems, and many others <cit.>. 

In porous media flow simulation,  instead of direct numerical simulation of the flow in the exact pore geometry-based Navier-Stokes flow models, one can be approximated by a simplified network model of throats and pores <cit.>. This technique reduces the computational complexity and allows for simulation using larger computational domains.

In simulations of the class of insulation materials that are composed of a large number of fibers, the network models are used to represent the high-conductive fibrous <cit.>

The application of the discrete network model to river network simulations is presented in <cit.>. The implementation of the model is based on the PETSc library for high-performance computing systems (HPC).  
The application of the network models to transient hydraulic simulations is considered in <cit.> and performed for problems such as water distribution in urban distribution systems, oil distribution, and hydraulic generation. 

In <cit.>,  the authors present a graph-based computational framework that facilitates the construction and analysis of large-scale optimization and simulation applications of coupled infrastructure networks.

The dynamic optimal electric power flow is simulated using network models in <cit.>.  

In  <cit.>, traffic flows are considered in complex networks, where the model is based on a graph or network of streets in which vehicles can move.

The application of the spatial networks for disease transmission in epidemiological models is considered in <cit.>. 

In <cit.>, a cerebral blood flow is modeled as fluid flow driven through a network of resistors by pressure gradients. The authors introduce a vascular graph modeling framework based on these principles to compute blood pressure, flow, and scalar transport in realistic vascular networks. 
The simulation of blood flow in microvascular networks and the surrounding tissue is considered in <cit.>. To reduce the computational complexity of this issue, the network structures are modeled by a one-dimensional graph whose location in space is determined by the centerlines of the three-dimensional vessels. 

Despite eliminating a significant portion of complexity through this approach, efficiently solving the resulting model remains a common challenge.
Homogenization and multiscale methods are commonly used techniques for implementing upscaling to manage the extensive computational complexity associated with large network models.



Multiscale problems arise in many areas of science and engineering and typically involve multiple spacial length scales.   Traditional numerical methods, such as finite element or finite volume methods, can become prohibitively expensive when the number of degrees of freedom required to capture all relevant scales becomes large.
To address this issue, multiscale and homogenization methods have been developed that seek to efficiently capture the essential features of the problem at each length scale.  

In homogenization methods, the fine-scale problem is replaced by an equivalent coarse-scale problem that captures the essential features of the original problem at a lower computational cost <cit.>.  The effective properties for the coarse grid approximation are computed by analyzing the behavior of unit cells. 
The most common homogenization method is based on periodicity, which assumes that the fine-scale problem can be represented as a periodic array of repeating unit cells. 

Compared to the homogenization techniques, a multiscale methods form a broader class of numerical techniques, for example,  the Multiscale Finite Element Method  <cit.>,  the Heterogeneous Multiscale Method <cit.>, the Local Orthogonal Decomposition method <cit.>,  the Variational Multiscale Method <cit.>, the Generalized Multiscale Finite Element Method <cit.>,  the Multiscale Finite Volume Method  <cit.> and many others. 
Most multiscale methods are based on constructing the multiscale basis functions in the local domains to capture fine-scale behavior. 



There are various applications of multiscale and upscaling methods for network problems that have been studied extensively.

Upscaling traffic flows in complex networks is considered in <cit.>. The problem is considered in two-dimensional regions whose size (macroscale) is much greater than the characteristic size of the network arcs (microscale).  

The numerical upscaling method is presented in <cit.>, where a finite element approximation with a  localized orthogonal decomposition method represents the macroscale model. Moreover, the application to a two-dimensional network model of paper-based materials in the form of fiber networks is considered <cit.>.   
In <cit.>, the nodal displacement in a fiber network model is analyzed using a multiscale method based on the LOD. 

In <cit.>, the heterogeneous multiscale method (HMM) is proposed to couple a network model on the microscale with a continuum scale over the same physical domain.

The coarsening procedures for graph Laplacian problems written in a mixed saddle point form were presented in <cit.>. 

The numerical methods for computing the effective heat conductivity of fibrous insulation materials are presented in  <cit.>. The fast algorithm is constructed based on the upscaling procedure. It contains the solution of the auxiliary boundary value problems of the steady-state heat equation in a representative elementary volume occupied by fibers and air. The presented approach ignores air and is further simplified by taking advantage of the slender shape of the fibers and assuming that they form a network. 

A multiscale method for networks representing flows in a porous medium is presented in <cit.>. 
In <cit.>,  the mortar coupling is presented to couple pore-scale network models to additional pore-scale or continuum-scale models using mortars. Mortars are finite-element spaces in two dimensions connecting distinct subdomains by ensuring pressure and flux continuity at their shared interfaces.

In <cit.>, a cerebral blood flow is modeled as fluid flow in a complex network.   The authors construct an upscaling algorithm that significantly reduces the computational cost. Furthermore,  the upscaled model no longer requires extensive information regarding the topology of the capillary bed.
The reduction of the computational complexity of the simulation of blood flow in microvascular networks and the surrounding tissue is considered in  <cit.>. The authors employ homogenization to the microvascular network's fine-scale structures in the study, leading to a new hybrid approach. This approach models the fine-scale structures as a heterogeneous porous medium, while the larger vessels' flow is modeled using one-dimensional flow equations.



This paper introduces the novel approach for upscaling the complex network model based on the Generalized Multiscale Finite Element Method <cit.>. The GMsFEM has a significant advantage in incorporating small-scale features from heterogeneities into coarse-grid basis functions. The multiscale basis functions are constructed by solving local eigenvalue problems. The online solutions can be calculated for any suitable boundary condition or forcing by these greatly reduced-dimension multiscale basis functions.
In this work, we present the construction of the reduced-order model for complex, highly heterogeneous networks. The network model represents the fine-scale model. In contrast, the coarse-scale approximation is represented by a much coarser finite element mesh than the fine-scale network. We design multiscale basis functions to account for the networks' microscale features by solving the local spectral problems in the primary local network cluster.   The constructed multiscale approximation can handle highly varying connectivity and random network structure with a huge system size reduction. Convergence analysis of the proposed method is presented for semi-discrete and discrete network models. We establish the stability of the multiscale discrete network for implicit approximation.   
Numerical implementation of the network model is performed based on the DMNetwork framework. DMNetwork is a part of PETSc library  <cit.> for high-resolution multiphysics simulations on the large-scale complex network  <cit.>. DMNetwork provides data and topology management, parallelization for multiphysics systems over a network, and hierarchical solvers. 
We present the main components of the multiscale method for constructing the coarse-scale system.   
To test the presented upscaling approach, we consider regular cubic lattice networks with and without random elimination processes and random networks <cit.>. We show that the multiscale method can provide an accurate solution with a huge system size reduction with fine-scale solution reconstruction.  


The paper is organized as follows. In Section 2, we consider the fine-scale model of the complex network in semi-discrete and discrete forms with stability analysis. In Section 3, we present the construction of the coarse-scale model using the multiscale method and provide a priory estimate.   The numerical results are presented in Section 4 for structured and random heterogeneous three-dimensional network models. The conclusion and future works discussion are given in Section 5.  




§ PROBLEM FORMULATION


We consider network represented as a undirected graph G = (P, E), where 
P = {v_1, v_2, ..., v_N_v} is a set of vertices (nodes v_i ∈ℛ^d, d=2,3) and 
E is a set of two-elements subsets of P (connections or edges e_ij = {v_i, v_j } that connect vertices v_i (head) and v_j(tail) and i ≠ j).  
Here N_v is the total number of vertices/nodes, and N_e is the total number of edges/connections <cit.>.  We suppose that the graph is connected and weighted. 

Furthermore, we assume that the network is embedded into the rectangular cuboid

    Ω  = [0, L_1] × ... × [0, L_d].

In Figures <ref>, we illustrate the two and three-dimensional networks embedded into the square and cube.  For network generation, we use OpenPNM library <cit.>.   
We note that the method can be constructed for a general case without embedding it into a hyper-rectangle.   However, it will affect the coarse-grid construction, and a more general way should be considered based on the graph partitioning. 






We assign a heterogeneous property c_i for each vertex v_i  that can be associated with the volume in the pore-network model <cit.>.  Then, we associate a weight w_ij  to each edge proportional to the area of the edge and inversely proportional to the distance between nodes <cit.> (Figure <ref>).  
Additionally, we label nodes associated with the top and bottom boundaries (Γ⊂∂Ω)  to set Dirichlet boundary conditions. 



Let D(G) be a degree matrix, W = {w_ij}_i,j=1^N_v be a symmetric weight matrix (w_ij≥ 0,  w_ii = 0 and w_ij > 0 if e_ij is an edge) and L(G) = D(G) - W be a graph Laplacian of G, where  D(G) = diag(d_1,  …, d_N_v) with d_i = ∑_j=1^N_v w_ij <cit.>.  Note that we will write L instead of L(G) for simplicity.  

For all u ∈ℛ^N_v, we have 

    (Lu)_i = ∑_j ∼  i w_ij (u_i - u_j),

where we write j ∼  i if  nodes x_i and x_j are connected. 

We suppose that w_ij are bounded weights (0 < w≤ w_ij≤w < ∞) then   

    u^T L u
       = u^T D u - u^T W u
    = ∑_i=1^N_v d_i u_i^2 - ∑_i, j=1^N_v w_ij u_i u_j
       = 1/2( 
     ∑_i=1^N_v d_i u_i^2  
     - 2 ∑_i, j=1^N_v  w_ij u_i u_j 
     + ∑_i=1^N_v d_i u_i^2
    )
     = 1/2∑_i, j=1^N_v w_ij (u_i - u_j)^2 ≥ 0,

and operator L is positive semidefinite. 


Furthermore, we can write the following representation of the graph Laplacian L 

    L = B M B^T,

where matrix M={m_ij}_i,j=1^N_e is the diagonal matrix filled with the edge weights and B = {b_ij} is the vertex-edge incident matrix with size N_v × N_e and element b_ij = 1 if v_i is a head of edge e_j,  b_ij = -1 if v_i is a tail of edge e_j and zero otherwise.  
Them, we have  (Lu, u) = u^T L u = (B^T u)^T M (B^T u) = 1/2∑_i, j=1^N_v w_ij (u_i - u_j)^2.






 §.§ Semi-discrete network model


On the weighted undirected graph G, we consider the following time-dependent  problem 

    c_i ∂ u_i/∂ t  + ∑_j w_ij (u_i - u_j) = f_i,  
      ∀ i  ∈𝒩,   
       0 < t ≤ T,

where 𝒩={1,…,N_v} is the set of indexes,  u_i is defined on on the node x_i, f_i is the given nodal source term and c_i is the bounded coefficient, 0 < c≤ c_i ≤c < ∞.

We consider equations (<ref>) with initial conditions

    u_i = u_i,0, 
      ∀ i  ∈𝒩, 
       t = 0,

and boundary conditions 

    u_i = g_i, 
       i ∈𝒩_D, 
        0 < t ≤ T,

where 𝒩_D is the set of node corresponded to the boundary <cit.>.


Let u = (u_1, …, u_N_v)^T ∈ℛ^N_v be the vector defined on the set of nodes,  f = (f_1, …, f_N_v)^T ∈ℛ^N_v be the given source term vector, and C be the diagonal matrix, C=diag(c_1,  …, c_N_v). Then we can write problem (<ref>) in the following matrix form

    C ∂ u/∂ t  + L u = f,    0 < t ≤ T.

with initial conditions u = u_0 for t=0. 
Note that, in this formulation, the matrices and right-hand side vector are modified to incorporate boundary conditions. For instance, to have a symmetric matrix, we do not include boundary nodes and impose boundary conditions by setting flux on the nodes connected to the boundary. 


Let  V ⊂ℛ^N_v be the subspace of real-valued functions defined on the set of nodes satisfying Dirichlet boundary conditions  <cit.>.  Then we can write problem formulation in the following equivalent form: find u ∈  V such that 

    ( C ∂ u/∂ t, v )  + (L u, v) = (f, v), 
      ∀ v ∈ V, 
       0 < t ≤ T,

where 
(u, v) = u^T v and (Lu, v) = u^T L v is scalar products. 


Next, we introduce the following two norms ||u|| = √((u,u)) and ||u||_L = √( (L u, u))  and  show a stability estimate for a semi-discrete fine-scale scheme (<ref>) <cit.>. 



The  solution of the problem (<ref>) satisfies the following a priory estimate

    ||u(t)||_L^2
     ≼ 
     ||u_0||_L^2 
    +∫_0^t ||f||^2_C^-1 ds.



Let v = du/dt in  (<ref>) then we have

    (C du/dt, du/dt) 
    + ( L u,   du/dt) = (f, du/dt).

Using Young's inequality 
and 
1/2d/dt(L u,  u) =  ( Lu,   du/dt), we obtain the following estimate

    du/dt^2_C  + 1/2d/dt( Lu,  u ) 
    =
     (f, du/dt)
    ≤du/dt^2_C + 1/4 f ^2_C^-1,

or

    d/dt u _L^2
    ≤1/2 f ^2_C^-1.

Finally,  after integration by time, we obtain inequality (<ref>) . 





 §.§ Time-approximation and discrete network model


Let u^n = u(t_n) and u^n-1 = u(t_n-1), where t_n = n τ,  n=1,2, ... and τ > 0 be the uniform time step size.  
For approximation by time, we use backward Euler's  to obtain an implicit fully discrete scheme: find u^n ∈ V such that 

    ( C  u^n - u^n-1/τ, v) + (L u^n, v) = (f^n, v), 
      ∀ v ∈ V, 
       n = 1,2,...

with initial conditions u^0 = u_0. 

Similarly to the semi-discrete problem (<ref>), we have the following stability estimate <cit.>. 



The solution of the  problem (<ref>) is unconditionally stable and satisfies the following estimate

    ||u^n||_L^2 
    ≼ 
    ||u_0||_L^2 
    + τ∑_k=1^n ||f^k||^2_(C + τ/2 L )^-1.



The equation (<ref>) can be written as follows

    ( C  u^n - u^n-1/τ, v) + (L u^n, v) 
    = 
    ( (C + τ/2 L )  u^n - u^n-1/τ, v  ) 
    + 1/2 (L (u^n  + u^n-1), v) = (f^n, v).

Let v =  u^n - u^n-1/τ then we have

    ( (C + τ/2 L )  u^n - u^n-1/τ ,  u^n - u^n-1/τ) 
    + 1/2( L (u^n + u^n-1),    u^n - u^n-1/τ) = (f^n,  u^n - u^n-1/τ).

Using  the following inequality for the right-hand side

    (f^n,  u^n - u^n-1/τ) ≤ u^n - u^n-1/τ_(C + τ/2 L ) ^2 
    + 1/4 f^n _(C + τ/2 L )^-1^2,

we obtain the following estimate for L = L^T

    1/2 τ( L (u^n + u^n-1),   (u^n - u^n-1)  )
    =
    1/2 τ (L u^n,u^n) + 
    1/2 τ (L u^n-1, u^n-1) 
    ≤1/4 f^n _(C + τ/2 L )^-1^2,

or

    ||u^n||^2_L ≤ ||u^n-1||^2_L  + τ/2 f^n _(C + τ/2 L )^-1^2.




We note that the presented discrete network model with implicit approximation by time leads to solving the large system of equations on each time step. To reduce the size of the system, we use a homogenization approach and construct a coarse-scale approximation by introducing local spectral multiscale basis functions. 
 



§ MULTISCALE MODEL ORDER REDUCTION


Multiscale methods form a broad class of numerical techniques where most multiscale methods are based on constructing the multiscale basis functions in the local domains to capture fine-scale behavior.  In this work, we construct a local spectral multiscale basis function for the network model described above and follow the procedure defined in the Generalized Multiscale Finite Element Method (GMsFEM) <cit.>.  




Let 𝒯_H be a coarse mesh in Ω

    𝒯_H = ⋃_i=1^N_c K_i,

where K_i is the coarse cell, and N_c is the number of coarse cells. 

We denote {y_i}_i=1^N the nodes of the coarse mesh 𝒯_H, where N is the number of coarse grid nodes <cit.>.  We define ω_i as a neighborhood of the node y_i

    ω_i = ⋃{K_j ∈𝒯_H:  y_i ∈K̅_̅j̅},

and ω^K_j as a the neighborhood of the element K_j

    ω^K_j = ⋃{ω_i : K_j ∈ω_i  }.

In this work we suppose that Ω  = [0, L_1] × ... × [0, L_d], where d = 2,3 is the dimension.  Furthermore,  we consider a uniform mesh with square and cubic cells for simplicity

    K_i = [y_i, y_i+H e_1] × ... × [y_i, y_i+H e_d],

where H is the coarse grid size.

Next, we associate a sub-network G^ω_i to each ω_i. 
In Figure <ref>, we depict a coarse 5 × 5 mesh for 2D network and 5 × 5 × 5 coarse grid for 3D network with local domain ω_i and corresponded sub-network G^ω_i.  
Note that in implementing the sub-network extraction, we ensure that any node from the global graph G is contained on exactly one coarse element. 

The objective of this paper is to develop and analyze a multiscale finite element method for discrete network models defined above.  The main idea of the multiscale method is in construction of the coarse-scale approximation of the network by construction of the accurate multiscale basis functions {ψ^ms_r,i}_r=1^M_i in each local domain G^ω_i (sub-network), where M_i is the number of local basis functions.  
In the multiscale method, we define offline and online stages:


  * Offline stage:


  * Coarse grid (𝒯_H) and local domains construction (sub-network, G^ω_i).

  * Multiscale space construction by the solution of the local eigenvalue in each G^ω_i

    V_ms = span{ψ^ms_r,i: 1 ≤ i ≤ N,  1 ≤ r ≤ M_i  }.

 

  * Online stage:   


  * Solution of the coarse-scale problem (Galerkin approximation): find u ∈ V_ms such that 

    ( C ∂ u_ms/∂ t, v_ms)  + (L u_ms, v_ms) = (f, v_ms), 
      ∀ v_ms∈ V_ms.





We use a same-time approximation as for fine-scale problem with u_ms^n = u_ms(t_n), where t_n = n τ,  n=1,2, ... and τ > 0.  
By backward Euler's scheme, we obtain a fully discrete scheme: find u_ms^n ∈ V_ms such that 

    ( C  u_ms^n - u_ms^n-1/τ, v_ms) + (L u_ms^n, v_ms) = (f^n, v_ms), 
      ∀ v_ms∈ V_ms, 
       n = 1,2,...


Next, we discuss details of the multiscale space construction and coarse-scale model generation with some implementation aspects. 




 §.§ Multiscale basis functions and algorithm


In this work, we use the OpenPNM <cit.> to construct a network and assign properties. The constructed network is saved in the text files with the corresponding properties. To implement the numerical solution of the fine-scale discrete network model, we use a PETSc library with DMNetwork framework <cit.>. DMNetwork provides data and topology management, parallelization for multiphysics systems over a large network, and hierarchical solvers. 


For the local network G^ω_i, we have

    L(ω_i) = D(ω_i) - W(ω_i),     
    D(ω_i) = diag(d_1, …,d_N_v^ω_i),    
    d_i = ∑_j=1^N_v^ω_i w_ij,

where W(ω_i) is the local weight matrix and  N_v^ω_i is the number of nodes of local network G^ω_i associated with the local domain ω_i. 

For the  sub-network G^ω_i,  the local matrix L^ω_i is positive semidefinite for w_ij > 0.  Therefore (1) the eigenvalues 0 = λ^ω_i_1 ≤λ^ω_i_2 ≤…≤λ^ω_i_N_v^ω_i  of L^ω_i are real and nonnegative,  and there is an orthonormal basis of eigenvectors of  L^ω_i;  and (2) the smallest eigenvalue λ^ω_i_1 = 0 and corresponded eigenvector is one  <cit.>. 

Therefore the dimension of the nullspace of L^ω_i (dimension of the multiplicity of the eigenvalue 0) is equal to the number of connected components C_1, …,C_C (clusters) of the local network G^ω_i.  Therefore the vectors f_l (f_l = (f_l,1,…,f_l, N_v^ω_i) such that  f_l,j = 1 if x_j ∈ C_l and f_j = 0 otherwise, l=1,…,C) are form a basis of the nullspace. 

In the numerical implementation of the eigenvalue problem solution, we use the SLEPc library <cit.>. We preprocess the local network G^ω_i by extracting the primary cluster on which we solve a local eigenvalue problem. We suppose that the size of the primary cluster is much larger than the others. In order to incorporate disconnected clusters into the local multiscale space, we use indicator vectors f_l. Note that we can solve an eigenvalue problem in each sufficiently big cluster to incorporate local fine-scale behavior and construct an accurate approximation for the general case. For network preprocessing, we use existing tools from the OpenPNM library <cit.>. 

Then, we suppose that G^ω_i has no isolated vertices, and therefore the local degree matrix D^ω_i contains positive entities and is invertible.   


Then similarly to the  GMsFEM for standard finite element-based definition, we use a generalized eigenvalue problem that gives a good approximation space <cit.>. 

To construct a multiscale basis functions in each sub-network G^ω_i, we solve a local generalized eigenvalue problem 

    L^ω_iϕ_r^ω_i
    = λ_r^ω_i D^ω_iϕ_r^ω_i

where L^ω_i = L(ω_i) and  D^ω_i = D(ω_i) are the matrices defined on the sub-network G^ω_i. 





Next, we choose first M_i eigenvectors corresponded to the smallest eigenvalues,  ϕ_r^ω_i, r = 1,…,M_i in each G^ω_i with i = 1,…,N.  
In order to form a continuous space, we multiply eigenvalues to the linear partition of unity functions and form a multiscale space

    V_ms = span{ψ^ms_r,i: 1 ≤ i ≤ N,  1 ≤ r ≤ M_i  }.

where
ψ^ms_r,i = χ_i ϕ_r^ω_i and χ_i is the projection of the linear partition of unity functions in ω_i to the nodes of corresponded sub-network G^ω_i.


We implement the presented multiscale algorithm based on the algebraic way by forming the projection matrix 

    R =[  
    ψ^ms_1,1, …, ψ^ms_M_1,1, … ,  
    ψ^ms_1,N, …,  ψ^ms_M_N,N ]^T.

To define a coarse grid approximation, we use a projection matrix and form a coarse-scale matrix and right-hand side vector

    C_H = R C R^T,     L_H = R L R^T,    F_H= R F.

Finally,  we solve a coarse-scale system each time t_n

    C_H u_H^n - u_H^n-1/τ + L_H u_H^n = F_H,

and reconstruction of the fine-scale solution u^n_ms = R^T u^n_H. 
The illustration of the main steps of the algorithm is depicted in Figure <ref>. 
We note that all offline calculations are independent for each local domain ω_i and can be done entirely parallelly. All coupling is done on the coarse level and generally does not require parallelization. However, for large systems, the size of the coarse-scale approximation can still lead to the large system, and further parallelization of the linear system solution at time t_n can be done based on the PETSc implementation of the parallel linear or nonlinear solvers (KSP and SNES classes). Furthermore, the more advanced models can be upscaled similarly with advanced time-stepping techniques available in the TS framework of PETSc.   We will consider such advanced coupling in future works. 




 §.§ Convergence analysis


In this part of the work, we present a priory error estimate of the multiscale method for network models.  We note that the analysis presented below is closely related to the analysis in <cit.>.  
We start with the definition of global norms used in the presented analysis. 



Let d(u, u) = (Du,u) =  u^T Du and A = D^-1 L  then 

    d(u,u) = ||u||^2_D, 
        d(Au, u) = (Au)^T D u = u^T A^T D u  
    
    = u^T (D^-1 L)^T D u 
    = u^T L u = a(u,u) = ||u||^2_L,
        d(Au, Au) = (Au)^T D (Au) = u^T (D^-1 L)^T D (D^-1 L) u 
    =  u^T (L^T  D^-1 L) u  = ||Au||^2_D,
 
for D= D^T and L = L^T <cit.>.



For each ω_i (i=1,…,N) we define the local projection

    P_M_i^ω_i v 
    = ∑_r=1^M_i ( v^T D^ω_iϕ_r^ω_i) ϕ_r^ω_i,

for v ∈ V  and given M_i in the local eigenvalue problem (<ref>).



For P_M_i^ω_i the following inequalities holds <cit.>

    ||v - P_M_i^ω_i v ||^2_D^ω_i≤1/λ_M_i+1^ω_i ||v - P_M_i^ω_i v ||^2_L^ω_i≤1/λ_M_i+1^ω_i ||v||^2_L^ω_i, 
    
       ||v - P_M_i^ω_i v ||^2_D^ω_i≤1/ (λ_M_i+1^ω_i)^2  ||A (v - P_M_i^ω_i v) ||^2_D^ω_i≤1/ (λ_M_i+1^ω_i)^2  ||Av||^2_D^ω_i,
    
       ||v - P_M_i^ω_i v ||^2_L^ω_i≤1/λ_M_i+1^ω_i ||A (v - P_M_i^ω_i v) ||^2_D^ω_i≤1/λ_M_i+1^ω_i ||Av||^2_D^ω_i.


Next, we define the coarse projection Π: V → V_ms by

    Π v = ∑_l=1^N  χ_l (P_M_l^ω_l v).

and v - Π v  = ∑_l=1^N χ_l (v - P_M_l^ω_l  v). 




Assume that  u ∈ V is the fine scale solution of (<ref>) then the  following estimate  holds

    ||u - Π u ||_D^2 ≤1/λ_M+1^2  ||A u||^2_D

where 
λ_M+1 = min_K  λ_K, M+1 and  
λ_K, M+1 = min_y_l ∈ Kλ_M_l+1^ω_l.


Using that χ_l ≤ 1 we have

    ||u - Π u ||_D^K^2 
    ≤∑_y_l ∈ K  ||χ_l (v - P_M_l^ω_l  u)||^2_D^K≤∑_y_l ∈ K ||u - P_M_l^ω_l  u||_D^ω_l^2.

By combing with estimate (<ref>), we obtain the result.




Assume that  u ∈ V is the fine scale solution of (<ref>) then the  following estimate  holds

    ||u - Π u ||_L^2 ≤( 1/H^2 λ_M+1^2 + 1/λ_M+1) ||A u||^2_D

where 
λ_M+1 = min_K  λ_K, M+1 and  
λ_K, M+1 = min_y_l ∈ Kλ_M_l+1^ω_l.



For z = {z_j}_j=1^N_v with z_j = χ_l,j v_j,  we have

    ||z||^2_L^ω_l     = z^T L^ω_l z
    = 1/2∑_i, j=1^N^ω_l_v 
    w_ij (z_i - z_j)^2 
    = 1/2∑_i, j=1^N^ω_l_v 
    w_ij (χ_l,i v_i - χ_l,j v_j)^2 
        = 1/2∑_i, j=1^N^ω_l_v 
    w_ij (χ_l,i (v_i - v_j) - v_j( χ_l,j - χ_l,i))^2
    ≤1/2∑_i, j=1^N^ω_l_v w_ijχ_l,i^2 (v_i - v_j)^2
    + 
    1/2∑_i, j=1^N^ω_l_v w_ij v_j^2 (χ_l,i - χ_l,j)^2.

Then for z = χ_l (v - P_M_l^ω_l  v)  using properties of partition of unity function <cit.>,   
we have

    ||u - Π u ||_L^K^2 
    ≤∑_y_l ∈ K  ||χ_l (v - P_M_l^ω_l  u)||^2_L^ω_l≤∑_y_l ∈ K1/H^2 ||u - P_M_l^ω_l  u||_D^ω_l^2  
    +  
    ∑_y_l ∈ K ||u - P_M_l^ω_l  u||_L^ω_l^2.

By combing with estimate (<ref>), we obtain  

    ||u - Π u ||_L^K^2 
    ≤
    
    
    ( 1/H^2 λ_K, M+1^2 + 1/λ_K, M+1) ∑_y_l ∈ K 
    ||Au||^2_D^ω_l

with λ_K, M+1 = min_y_l ∈ Kλ_M_l+1^ω_l. 

Then, we have

    ||u - Π u ||_L^2   
     = ∑_K ∈𝒯_H  ||u - Π u||_L^K^2 
       ≤∑_K ∈𝒯_H( 1/H^2 λ_K, M+1^2 + 1/λ_K, M+1) ∑_y_l ∈ K  ||A (u - P_M_l^ω_l  u) ||^2_D^ω_l
       ≤( 1/H^2 λ_M+1^2 + 1/λ_M+1) ∑_y_l ∈ K||Au||_D^ω_l^2,

with λ_M+1 = min_K  λ_K, M+1. 





  §.§.§ Multiscale semi-discrete network


In this subsection, we consider semi-discrete network (<ref>) and give a priory estimates for multiscale approximation. 

Stability estimate derivation is similar to <ref>. We let v_ms = du_ms/dt = (u_ms)_t in  (<ref>) then we have

    (C (u_ms)_t,  (u_ms)_t ) 
    + ( L u_ms,   (u_ms)_t ) = (f, (u_ms)_t ).

Then by Young's inequality, we obtain the following estimate 

    (u_ms)_t ^2_C  + 1/2d/dt( Lu_ms,  u_ms) 
    ≤(u_ms)_t ^2_C + 1/4 f ^2_C^-1.

Then after integration by time, we obtain the following lemma. 



The  solution of the problem (<ref>) satisfies the following a priory estimate

    ||u_ms(t)||_L^2 ≤ ||u_ms(0)||_L^2 
    + 1/2∫_0^t ||f||^2_C^-1 ds.




Next, we can obtain the error estimate of the GMsFEM for the semi-discrete network model. 



Assume that  u ∈ V is the fine scale solution of (<ref>) and 
u_ms∈ V_ms is the multiscale solution of (<ref>)  then the  following error estimate  holds

    ||(u - u_ms)(t)||_C^2 + ∫_0^t || u - u_ms||_L^2 ds
       ≼
    ||(w - u_ms)(0)||_C^2 
    + 1/λ_M+1^2 ||A u||^2_D 
       + ∫_0^t  (  
    ( 1/H^2 λ_M+1^2 + 1/λ_M+1)  || A u||^2_D 
    +  1/λ_M+1^2 || (A u)_t ||^2_D 
    ) ds



Subtracting (<ref>)  from  (<ref>) , we have 

    ( C (u -u_ms)_t, v )  + (L (u - u_ms), v)  = 0
      ∀ v ∈ V_ms.


We write u - u_ms = (u  - w) + (w - u_ms) with w = Π u and since w ∈ V_ms we take v = w - u_ms and obtain

    ( C (w -u_ms)_t, w - u_ms)  + (L (w - u_ms), w - u_ms)  
    =
    ( C (w - u)_t, w - u_ms)  + (L (w - u), w - u_ms).


For the left part of (<ref>), we have

    ( C(w -u_ms)_t, w - u_ms) = 1/2d/dt||w - u_ms||_C^2,   
    (L (w - u_ms), w - u_ms)   = ||w - u_ms||^2_L.


For the right part of (<ref>) by Young's inequality, we obtain

    ( C(w - u)_t,  w - u_ms) ≤1/4 δ_1|| (w - u)_t ||^2_C + δ_1 ||w - u_ms ||^2_C,

and

    (L (w - u), w - u_ms) ≤1/4 δ_2 ||w - u||^2_L + δ_2 ||w - u_ms||^2_L.


Combing estimates for right and left parts of (<ref>) with δ_1 = δ_2 = 1/2, we can obtain 

    d/dt||w - u_ms||_C^2 + ||w - u_ms||^2_L 
    
    ≤
    || (w - u)_t ||^2_C + ||w - u_ms ||^2_C + ||w - u||^2_L.

By the Gronwall Lemma 
we obtain


    || (w - u_ms)(t)||_C^2 
    + ∫_0^t ||w - u_ms||^2_L  ds
    
    ≤
    || (w - u_ms)(0)||_C^2
    + ∫_0^t ( || (w - u)_t ||^2_C + ||w - u||^2_L  ) ds.

Then using triangle inequality we have







    || (u - u_ms)(t)||_C^2 
    + ∫_0^t ||u - u_ms||^2_L  ds
       ≤
    || (w - u_ms)(0)||_C^2 + || (w-u)(t) ||_C^2 
       + ∫_0^t ( ||(w-u)_t ||^2_C + ||w-u||^2_L ) ds

Finally with ||u-v||^2_C ≼ ||u-v||^2_D and using Lemmas <ref> and <ref>, we obtain the result.




Note that, we have Au = D^-1 Lu = D^-1 (f - C u_t) then

    ||Au||_D^2 = ||f- C u_t||^2_D^-1.
 
Therefore the error estimate in Theorem <ref> can be written as follows

    ||(u - u_ms)(t)||_C^2 + ∫_0^t || u - u_ms||_L^2 ds 
       ≼
    ||(w - u_ms)(0)||_C^2 
    + 1/λ_M+1^2 ||f- C u_t||^2_D^-1
       + ∫_0^t  (  
    ( 1/H^2 λ_M+1^2 + 1/λ_M+1) ||f- C u_t||^2_D^-1
    +  1/λ_M+1^2||(f- C u_t)_t||^2_D^-1) ds

If we map the local eigenvalue problem to size one domain <cit.>, then the eigenvalues scale with H^-2 and we obtain 

    ||(u - u_ms)(t)||_C^2 + ∫_0^t || u - u_ms||_L^2 ds 
    ≼  H^2 λ_M+1^-1

with accurate choose of u_ms(0).




  §.§.§ Multiscale discrete network


Finally, we present a priory estimate for a fully discrete problem (<ref>). 

Stability estimate derivation is similar to Lemma <ref>. We let v_ms = ( u_ms^n - u_ms^n-1)/ τ then 

    ( C  u_ms^n - u_ms^n-1/τ,   u_ms^n - u_ms^n-1/τ)     + (L u_ms^n,  u_ms^n - u_ms^n-1/τ) 
       = 
    ( (C + τ/2 L )  u_ms^n - u_ms^n-1/τ ,  u_ms^n - u_ms^n-1/τ) 
    + 1/2( L (u_ms^n + u_ms^n-1),    u_ms^n - u_ms^n-1/τ) 
       = 
    (f^n,  u_ms^n - u_ms^n-1/τ)
    ≤ u_ms^n - u_ms^n-1/τ_(C + τ/2 L ) ^2 
    + 1/4 f^n _(C + τ/2 L )^-1^2,

Then we have

    1/2 τ( L (u_ms^n + u_ms^n-1),   (u_ms^n - u_ms^n-1)  )
    =
    1/2 τ (L u_ms^n,u_ms^n) + 
    1/2 τ (L u_ms^n-1, u_ms^n-1) 
    ≤1/4 f^n _(C + τ/2 L )^-1^2,

or

    ||u_ms^n||^2_L ≤ ||u_ms^n-1||^2_L  + τ/2 f^n _(C + τ/2 L )^-1^2.




The solution of the  problem (<ref>) is unconditionally stable and satisfies the following estimate

    ||u_ms^n||_L^2 
    ≼ 
    ||u_ms^0||_L^2 
    + τ∑_k=1^n ||f^k||^2_(C + τ/2 L )^-1.



Next, we present the error estimate of the multiscale method  for the discrete  network model. 



Assume that  u^n ∈ V is the fine scale solution of (<ref>) and 
u^n_ms∈ V_ms is the multiscale solution of (<ref>)  then the  following error estimate  holds

    ||u^n - u^n_ms||_C^2 + τ∑_k=1^n || u^k - u^k_ms||_L^2
       ≼
    ||u^0 - u_ms^0||_C^2 
    
    + 
       ∑_k=1^n ( 
    τ( 1/H^2 λ_M+1^2 + 1/λ_M+1)  || A u^k||^2_D
    +
     1/λ_M+1^2 || A u^k||^2_D
     ).



Subtracting (<ref>)  from  (<ref>) , we have 

    ( C  (u^n -u^n_ms) /τ, v ) - ( C  (u^n-1 -u^n-1_ms) /τ, v )  + (L (u^n - u^n_ms), v)  = 0
      ∀ v ∈ V_ms.


For u^n - u^n_ms = (w^n - u^n_ms) + (u^n - w^n) with w^n = Π u^n and v = w^n - u^n_ms, we have

    ( C  (w^n -u^n_ms) /τ, w^n - u^n_ms) 
    
    +  (L (w^n - u^n_ms), w^n - u^n_ms)  
        =  
    ( C (w^n - u^n) /τ, w^n - u^n_ms) 
    
    + ( C (u^n-1 - u_ms^n-1) /τ, w^n - u^n_ms)  
    + (L (w^n - u^n), w^n - u^n_ms).


For the left part of (<ref>), we have

    ( C (w^n -u^n_ms), w^n - u^n_ms) = ||w^n - u^n_ms||_C^2,   
    (L (w^n - u^n_ms), w^n - u^n_ms)   = ||w^n - u^n_ms||^2_L.


For the right part of (<ref>) by Young's inequality, we have

    ( C (w^n - u^n), w^n - u^n_ms) ≤1/4 δ_1||w^n - u^n ||^2_C + δ_1 ||w^n - u^n_ms ||^2_C,


    ( C (u^n-1 - u_ms^n-1), w^n - u^n_ms)  ≤1/4 δ_2||u^n-1 - u_ms^n-1 ||^2_C + δ_2 ||w^n - u^n_ms ||^2_C,

and

    (L (w^n - u^n), w^n - u^n_ms) ≤1/4 δ_3 ||w^n - u^n||^2_L + δ_3 ||w^n - u^n_ms||^2_L.


Then with δ_1 = δ_2 = 1/4 and δ_3 = 1/2, we obtain 

    ||w^n - u^n_ms||_C^2 + τ ||w^n - u^n_ms||^2_L 
    ≼ 
    ||u^n-1 - u_ms^n-1 ||^2_C 
    +
    ||w^n - u^n||^2_C +  τ ||w^n - u^n||^2_L.


Using triangle inequality we have





    ||u^n - u^n_ms||_C^2 
    + τ ||u^n - u^n_ms||^2_L 
    ≼
    ||u^n-1 - u_ms^n-1 ||^2_C 
    +
    ||w^n - u^n||^2_C + τ ||w^n - u^n||^2_L.

Therefore

    ||u^n - u^n_ms||_C^2 
    + τ∑_k=1^n ||u^n - u^n_ms||^2_L 
    ≼
    || u^0 - u_ms^0 ||^2_C 
    + ∑_k=1^n  (
    ||w^n - u^n||^2_C + τ ||w^n - u^n||^2_L
    ).

Finally using Lemmas <ref> and <ref>, we obtain the result.





§ NUMERICAL RESULTS


In this section, we present numerical results for the proposed multiscale method. 

To test our approach, we consider three networks: 


  * Network-1: Structured network with  15625 nodes and 45000 edges.  
This network is related to the finite volume approximation in the heterogeneous domain on the 25 × 25 × 25 structured grid. 


  * Network-2: Structured network with dropping out of  12765 nodes and 24718 edges.  
This network is related to the previous 25 × 25 × 25 structured grid with an additional dropout procedure, where we randomly remove approximately  20% of nodes and connections with the additional removal of disconnected pores and broken connections to make a 'healthy' network.  


  * Network-3: Unstructured network with  15625 nodes and 119713 edges.  
This random network is formed by the Delaunay tessellation of arbitrary base points implemented in openpnm. 


In OpenPNM, we use a Cubic with and without throats, nodes dropping, and random structures (from left to right in Figure <ref>).  Random elimination is applied for some nodes and throats in the second network.  After that, the isolated elements and clusters are removed to prevent numerical problems in simulating flow within the network <cit.>.  

Next, to generate a heterogeneous property, we randomly assign node volumes, and throat area, find the length of each throat,  and calculate the connectivity coefficient based on the area and length of the throat.  Additionally, we set labels for nodes on the top and bottom boundaries. 



The first network is related to the convenient finite volume approximation of the parabolic equation in the heterogeneous domain. In this network, we assign a heterogeneous diffusion coefficient to each node (cell) and calculate connection weight as a harmonic average of the cell's coefficients that it connect. 
 For the Network-1, we have k_min=3.49 and k_max = 148.41 as minimum and maximum for the diffusion coefficient. We have  c_min=0.05 and c_max = 0.2 as minimum and maximum values for the time derivative coefficient. 

To set heterogeneous properties on the second and third networks,  we randomly assign node (pore) size diameters, assign edge (throat) diameters based on the pores it connects,  and calculate connection weight based on the throat diameter and distance between pores. This algorithm is convenient in pore-scale simulations using pore network models.  
For the Network-2, we have k_min=0.0055 and k_max = 0.116 as minimum and maximum for the connection's coefficients. For the time derivative coefficient, we have  c_min=0.0027 and c_max = 0.03 as a minimum and maximum for the node's coefficients. 
For the Network-3, we have k_min=0.00029 and k_max = 1.45 as minimum and maximum for the connection's coefficients. For the time derivative coefficient, we have  c_min=0.0027 and c_max = 0.03 as a minimum and maximum for the node's coefficients. 














 
 






































 




















 
We apply the presented multiscale method for each network and investigate accuracy for different number numbers of basis functions.  
The networks are embedded into domain (Ω=[0,1] × [0,1]× [0,1].
We set zero source terms and set Dirichlet boundary conditions on top and bottom boundaries with one and zero values, respectively.   The time used in simulation was T = 0.5 (Network-1) and T = 50 (Network-2 and 3) with 50 time steps. 
The coarse grid is 5 × 5 × 5.   In Figure <ref>, we represent eigenvectors corresponding to the first eight smallest eigenvalues for each network. 
Numerical implementation of the network model is performed based on the PETSc DMNetwork framework <cit.>. The basis functions are constructed using a generalized eigenvalue solver from SLEPc  <cit.>. 


To compare the accuracy of the presented multiscale method, we use the relative L_2 error in percentage.  
The errors are calculated using the following formula on the fine grid: 

    e^n = ||u^n - u_ms^n||_L_2/||u^n||_L_2× 100 %,
 
where n is the time layer, u^n are the reference  (fine grid) solution  and  u_ms^n are the solution using the multiscale method. 
We use the corresponding fine-grid solution for each test problem as a reference solution. 

In Figures <ref>,<ref>, and <ref>, we depict the reference (fine grid) and multiscale solution for Network-1, Network-2, and Network-3, respectively. 
The simulation results are presented for three time layers t_5, t_20 and t_50  (from left to right).   The reference solution is presented on the first row, and a multiscale solution using 16 multiscale basis functions on the second row.  
From the fine-scale solution, we observe a significant influence of the heterogeneity properties and network structure on the solution.  
The fine-scale network model leads to the system of equations with  DOF_h = 15625 for Network-1, DOF_h = 12765 for Network-2 and DOF_h = 15625 for Network-3.   By the multiscale method, we reduce the size of the system to DOF_H = 1728 for eight multiscale basis functions with 0.4 %, 2.06 %, and 2.8 % of errors for Network-1, Network-2, and Network-3.  We have DOF_H = 3456 for 16 multiscale basis functions with 0.2 %, 1.4 % and 1.9 % of errors for Network-1, Network-2 and Network-3.  


In Table <ref>, we present errors at the final time for Network-1, Network-2, and Network-3.  
We have good accuracy for the multiscale method with a sufficient number of multiscale basis functions. For example, we have 3.07 %, 16.1 %, and 4.3 % of relative error for one basis function per local domain for Network-1, Network-2, and Network-3, respectively. When we use eight multiscale basis functions, errors reduce to 0.4 %, 2.06 %, and 2.8 %. We observe a significant influence of the network structure on the multiscale method solution, where for the first structure network, we obtain a better accuracy for a multiscale method using a smaller number of basis functions. Unstructured networks give a more significant error using one basis function, but accurate approximation can be performed with a larger number of basis functions. 
The dynamic of the errors is presented in Figure <ref> for Network-1, Network-2, and Network-3. We observe that the errors reduce over time and give an accurate solution for structured and unstructured networks.  




§ CONCLUSION


We considered a time-dependent model on structured and random networks.  The time approximation was performed using an implicit scheme.  The stability of the semi-discrete and discrete networks was presented. 
The multiscale method for the network model was developed and analyzed for network models.  The proposed approach is based on the generalized multiscale finite element method.  To find a multiscale basis function, we solve local eigenvalue problems in sub-networks.  
Convergence analysis of the proposed method was presented for semi-discrete and discrete network models with stability estimates.   
Numerical results were presented for structured and random heterogeneous networks to confirm the theory.  



unsrt




