\documentclass[11pt,reqno]{amsart}
\usepackage[cp1251]{inputenc}
\usepackage[english,russian]{babel}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsfonts}
\usepackage{graphicx}
\usepackage{xcolor}
\usepackage[colorlinks]{hyperref}


\newcommand{\grad}{\mathop{\rm grad}\nolimits}
\renewcommand{\div}{\mathop{\rm div}\nolimits}
\renewcommand{\exp}{\mathop{\rm exp}\nolimits}
\newcommand{\erf}{\mathop{\rm erf}\nolimits}
\usepackage{bm}
\graphicspath{ {./pics/} }


\usepackage{titlesec}
\titlespacing{\section}{0pt}{2.5ex}{1.5ex}
\titlespacing{\subsection}{0pt}{1.5ex}{1ex}
\titlespacing{\subsubsection}{0pt}{1.5ex}{1ex}
\titleformat{\section}{\large\bfseries\centering}{\thesection}{1em}{}
\titleformat{\subsection}[runin]{\bfseries}{\thesubsection.}{0.5em}{}[.\mbox{\ }]
\titleformat{\subsubsection}[runin]{\bfseries}{\thesubsubsection.}{0.4em}{}[.\mbox{\ }]





%%%%%---------- ORCiD
\newcommand\orcidicon[1]{\href{https://orcid.org/#1}{\includegraphics[scale=0.02]{orcid.pdf}}}


\def\udcs{519.6} %Here the author places classificators of the paper according to Russian classification system
\def\mscs{65N30} %Here the author places  classificators according to the AMS classification list.
\setcounter{page}{144}

\newtheorem{lemma}{Lemma}
\newtheorem{theorem}{Theorem}
\newtheorem{definition}{Definition}
\newtheorem{corollary}{Corollary}
\newtheorem{proposition}{Proposition}



\hypersetup{linkcolor=blue,filecolor=black,urlcolor=blue, citecolor=blue}

\def\logo{\textcolor{blue}{{\bf\Huge S\raisebox{0.2ex}{\hspace{0.55ex}\raisebox{0.05ex}e\hspace{-1.65ex}$\bigcirc$}MR}}}

\def\semrtop
     {
        \vbox{
            \noindent\logo \hspace{-3mm}  \parbox{12cm}{     \begin{center}
            {\LARGE \textsf{СИБИРСКИЕ ЭЛЕКТРОННЫЕ}} \\[2mm]
            {\LARGE \textsf{МАТЕМАТИЧЕСКИЕ ИЗВЕСТИЯ}} \\[2mm]
            {\Large \textsf{Siberian Electronic Mathematical Reports}} \\[1mm]
            {\Large\tt{\textcolor{blue}{http://semr.math.nsc.ru}}}\\[1mm]
            \raisebox{1ex}{ISSN 1813-3304 }
            \vspace*{-3mm}
       %     {\small Институт математики им. С.Л.Соболева СО РАН}\\[-1mm]
       %     {\small Sobolev Institute of Mathematics SB RAS}
            \end{center}
}
     \vspace{-1mm}
     \noindent
      \begin{tabular}{c}
     \hphantom{aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa!} \\
     \hline\hline
     \end{tabular}

     \vspace{3mm}
     {\flushleft\it Том 20, № 2, стр. 144--145 (2023) \hspace{\fill}{\rm\small УДК \udcs \hspace{-4mm}}}
     \newline  %Volume, pages in Russian, filled by Editorial board
     	{\rm\small \textcolor{blue}{https://doi.org/10.33048/semi.2023.20.???}}\hphantom{aaaaaaaaaaaaaaaaaaaaaaaa\;}{\rm\small MSC\ \ \mscs }
  }%\vbox
}




\begin{document}
\renewcommand{\refname}{References}
\renewcommand{\proofname}{Proof.}
\renewcommand{\figurename}{Fig.}




\thispagestyle{empty}

\title[\sc DG-GMsFEM for Poroelasticity Equations in Het.]{\large Discontinuous Galerkin Generalized Multiscale Finite Element Method for Poroelasticity Equations in Heterogeneous Media}

\author[V.N. ALEKSEEV] 
{
{\bf \large \textsc{V.N. Alekseev}}\protect\orcidicon{0000-0002-6952-3139}\\ \\ 
}% 

\address{Alekseev Valentin Nikolaevich
\newline\hphantom{iii} Laboratory of Computational Technologies and Artificial Intelligence, 
\newline\hphantom{iii} North-Eastern Federal University, Belinsky st., 58,
\newline\hphantom{iii} 677000, Yakutsk, Russia}

\email{\textcolor{blue}{alekseev.valen@mail.ru}}


\thanks{\sc Alekseev V.N., 
Discontinuous Galerkin Generalized Multiscale Finite Element Method for Poroelasticity Equations in Heterogeneous Media}
\thanks{\copyright \ Alekseev V.N.}
\thanks{\rm The work is supported by the Russian government project Science and Universities (project FSRG-2024-0003) aimed at supporting junior laboratories and the Russian Science Foundation grant 23-71-30013}
\thanks{\it Received September, 30, 2024, Published September, 30, 2024}%


\semrtop \vspace{1cm}
\maketitle {\small
\vspace{-12pt}
\centerline{{\it Communicated by} {\sc P.P. Petrov}
%({\it filled by Editorial board})
}

\bigskip
\bigskip

\begin{quote}
\noindent{\bf Abstract:} This paper presents the Discontinuous Galerkin Generalized Multiscale Finite Element Method (DG-GMsFEM) for solving poroelasticity equations. The method is designed to capture fine-scale variations in pressure and displacement by using multiscale basis functions, enabling efficient computation on coarse grids. Two test cases are considered: one with homogeneous elasticity and permeability coefficients, and the other with heterogeneous coefficients. The accuracy of the method is validated by calculating the relative $L_2$ errors for both pressure and displacement, which decrease as the number of multiscale basis functions increases. The results demonstrate that the DG-GMsFEM method provides accurate approximations for both pressure and displacement, significantly reducing computational cost compared to the fine-scale solutions. Numerical experiments show that the DG-GMsFEM approach is robust and accurate for complex poroelasticity problems, even in highly heterogeneous media.

\noindent{\bf Keywords:} poroelasticity, DG-GMsFEM, multiscale methods, finite element method, pressure, displacement, permeability, heterogeneous media, coarse grid. 
 \end{quote}
}
\bigskip

\section{Introduction}

Poroelasticity, the study of fluid-saturated porous media, plays a critical role in many engineering and geophysical applications, such as subsurface flow, reservoir simulation, and biomechanics. The poroelasticity equations describe the interaction between mechanical deformations and fluid flow within porous materials, which often requires resolving fine-scale heterogeneities in material properties, such as permeability and elasticity. These fine-scale variations introduce a significant computational challenge, as capturing them across large domains using traditional numerical methods on fine grids leads to a substantial increase in computational cost \cite{1,2,3,4}.

To address this challenge, multiscale methods have been developed to efficiently capture fine-scale effects while performing computations on a coarse grid. One such method is the Generalized Multiscale Finite Element Method (GMsFEM), which constructs localized multiscale basis functions to account for subgrid variations in material properties \cite{5,6,7,8}. Unlike traditional methods, GMsFEM enables computational savings by reducing the need for fine-scale resolutions across the entire domain. The Discontinuous Galerkin (DG) version of GMsFEM (DG-GMsFEM) further enhances this approach by combining the flexibility of discontinuous Galerkin methods with the efficiency of multiscale modeling, allowing for the accurate resolution of complex phenomena, such as jumps in material properties and fluxes across element boundaries \cite{9,10,11,12}.

In recent years, several studies have extended the GMsFEM to poroelasticity problems. Fu, Chung and Mai \cite{13} developed a constraint energy minimizing GMsFEM for nonlinear poroelasticity and elasticity, demonstrating its effectiveness in capturing complex material behaviors. Tyrylgin et al. \cite{14} introduced a GMsFEM framework for poroelasticity in multicontinuum media, addressing the challenge of modeling interactions between multiple porous media types. Brown and Vasilyeva \cite{15} proposed a GMsFEM approach specifically for linear poroelasticity, showing its applicability in heterogeneous media. Their results indicate that multiscale methods can significantly reduce computational costs while preserving solution accuracy.

Moreover, the method has been extended to fractured and heterogeneous media. Ammosov et al. \cite{16} applied GMsFEM to thermoporoelasticity problems, handling the additional complexity of thermal effects in fractured media. Similarly, Tyrylgin et al. \cite{17,18,19} have developed GMsFEM approaches for coupled poroelasticity problems, focusing on efficient modeling of fractured and heterogeneous media using embedded fracture models.

In this paper, we extend the DG-GMsFEM framework to solve the poroelasticity equations. The DG formulation is chosen for its flexibility in handling discontinuities in material properties and ensuring stability across coarse grids, which is particularly beneficial for heterogeneous media. Our goal is to develop a method that efficiently solves for both displacement and pressure in poroelastic materials, while reducing computational cost compared to fine-scale methods. The key idea behind DG-GMsFEM is to construct multiscale basis functions for both pressure and displacement by solving local problems on fine grids, and then using these basis functions to approximate the global solution on a coarse grid.

We consider two test cases to evaluate the performance of our method: one with homogeneous coefficients for elasticity and permeability, and another with heterogeneous coefficients. The accuracy of the method is validated by comparing the coarse-grid multiscale solutions with reference fine-grid solutions. 

The remainder of this paper is organized as follows: Section 2 formulates the problem of poroelasticity and outlines the governing equations. Section 3 describes the fine grid and mixed finite element formulation. Section 4 presents the multiscale basis construction using DG-GMsFEM for pressure and elasticity. Section 5 provides numerical results, and finally, Section 6 concludes the paper with a discussion of the method's advantages.


\section{Problem formulation}
Let the domain $\Omega \subset \mathbb{R}^d$ represent a poroelastic medium where $d=2$ for two-dimensional problems. The fluid flow and mechanical behavior in this poroelastic medium are governed by a system of coupled equations. The fluid's movement is described by the mass conservation equation along with Darcy's law:

\begin{equation}
\label{eq:mass_balance}
\begin{split}
    \alpha \frac{\partial (\nabla \cdot \mathbf{u})}{\partial t} + \frac{1}{M} \frac{\partial p}{\partial t} + \nabla \cdot \mathbf{q} = f, \quad x \in \Omega, \\
    \mathbf{q} = -k \nabla p, \quad x \in \Omega,
\end{split}
\end{equation}
where $\mathbf{u}$ is the displacement vector, $\mathbf{q}$ represents the fluid velocity, $p$ is the fluid pressure, and $k$ is the effective permeability (related to the intrinsic permeability $\kappa$ via $k = \frac{\kappa}{\mu}$, where $\mu$ is the dynamic viscosity of the fluid). The term $f$ accounts for external sources or sinks of fluid, $\alpha$ is the Biot coefficient, and $M$ is the Biot modulus, which relates fluid pressure to volumetric strain.

The mechanical response of the poroelastic material is described by a stress balance equation, where the forces due to fluid pressure and solid deformation are combined. The equation for the displacement field $\mathbf{u}$ is given by:
\begin{equation}
\label{eq:mechanics}
\begin{split}
    - \nabla \cdot \boldsymbol{\sigma}_T(\mathbf{u}, p) = 0, \quad x \in \Omega, \\
    \boldsymbol{\sigma}_T(\mathbf{u}, p) = \boldsymbol{\sigma}(\mathbf{u}) - \alpha p \mathbf{I}, \quad x \in \Omega,
\end{split}
\end{equation}
where $\boldsymbol{\sigma}_T$ is the total stress tensor, $\boldsymbol{\sigma}(\mathbf{u})$ is the mechanical stress tensor, and $\mathbf{I}$ is the identity matrix. The term $\alpha p \mathbf{I}$ represents the contribution of pore pressure to the total stress.

The coupled system of poroelasticity equations then becomes:

\begin{equation}
\label{eq:coupled_system}
\begin{split}
    \alpha \frac{\partial (\nabla \cdot \mathbf{u})}{\partial t} + \frac{1}{M} \frac{\partial p}{\partial t} - \nabla \cdot (k \nabla p) = f, \quad x \in \Omega, \\
    - \nabla \cdot \boldsymbol{\sigma}(\mathbf{u}) + \alpha \nabla p = 0, \quad x \in \Omega.
\end{split}
\end{equation}

For a linear elastic material, the stress-strain relationship is given by:

\[
\boldsymbol{\sigma}(\mathbf{u}) = 2 \mu \boldsymbol{\varepsilon}(\mathbf{u}) + \lambda (\nabla \cdot \mathbf{u}) \mathbf{I}, \quad \boldsymbol{\varepsilon}(\mathbf{u}) = \frac{1}{2} (\nabla \mathbf{u} + \nabla \mathbf{u}^T),
\]
where $\boldsymbol{\varepsilon}(\mathbf{u})$ is the strain tensor, and $\lambda$ and $\mu$ are the Lame's parameters.

The system of equations \eqref{eq:coupled_system} is solved with the following initial conditions:
\begin{equation}
\label{eq:initial_conditions}
p(x, 0) = p_0(x), \quad \mathbf{u}(x, 0) = \mathbf{u}_0(x), \quad x \in \Omega,
\end{equation}
where $p_0(x)$ and $\mathbf{u}_0(x)$ are the initial pressure and displacement fields, respectively.

The boundary conditions are defined as follows:
\begin{equation}
\label{eq:boundary_conditions}
\begin{split}
    & -k \nabla p \cdot \mathbf{n} = 0, \quad x \in \partial \Omega, \quad p = g, \quad x \in \Gamma_T, \\
    & u_x = 0, \ \sigma_y = 0, \quad x \in \Gamma_L \cup \Gamma_R, \quad u_y = 0, \ \sigma_x = 0, \quad x \in \Gamma_T \cup \Gamma_B,
\end{split}
\end{equation}
where $\Gamma_L$, $\Gamma_R$, $\Gamma_T$, and $\Gamma_B$ represent the left, right, top, and bottom boundaries of the domain $\Omega$, respectively, and $\mathbf{n}$ is the outward-pointing normal vector on the boundary $\partial \Omega$.

\section{Fine grid approximation}

In this section, we present the fine grid approximation for the poroelasticity system using the mixed finite element method. The computational domain $\Omega$ is discretized into a fine grid partition $T_h$, where $h$ is the characteristic element size of the fine grid (see Figure \ref{gfine}).

Let $T_h$ be the fine grid partition, defined as the union of all cells $K_i$ such that:
\[
T_h = \bigcup_{i=1}^{N_{\text{cell}}} K_i,
\]
where $N_{\text{cell}}$ is the total number of cells in the fine grid. We use
$\mathcal{E}^h$ to denote the set of facets in $\mathcal{T}^h$ with  
$\mathcal{E}^h = \mathcal{E}^h_{o} \cup \mathcal{E}^h_{b}$, 
where 
$\mathcal{E}^h_{o}$ and $\mathcal{E}^h_{b}$ are the set of interior and boundary facets (see Figure \ref{gfine}).

\begin{figure}[h!]
\centering
\includegraphics[width=0.5\textwidth]{finemesh}
\caption{Illustration of fine grid $\mathcal{T}^h$ (right).}
\label{gfine}
\end{figure}

In the mixed finite element formulation, the displacement field $\mathbf{u}_h$ and the pressure field $p_h$ are solved simultaneously. The mixed finite element space is defined as a product of two spaces: one for the displacement and one for the pressure.

We define the discrete spaces for displacement and pressure as:
\[
V_h = \{ \mathbf{v}_h \in [L^2(\Omega)]^d : \mathbf{v}_h|_K \in [P_1(K)]^d, \forall K \in T_h \},
\]
\[
Q_h = \{ q_h \in L^2(\Omega) : q_h|_K \in P_1(K), \forall K \in T_h \}.
\]

The mixed finite element space is then the product space $V_h \times Q_h$, where we approximate the displacement as $\mathbf{u}_h \in V_h$ and the pressure as $p_h \in Q_h$.

The discretized unknowns are \(\mathbf{u}_h\) for the displacement field and \(p_h\) for the pressure field, defined in the spaces \(V_h\) and \(Q_h\), respectively.

The weak form of the governing system \eqref{eq:coupled_system} is obtained by multiplying the first equation by a test function \(v_h \in Q_h\) and the second equation by a test function \(\mathbf{v}_h \in V_h\). After integrating by parts over each element \(K \in T_h\), we obtain the following weak form:

\begin{equation}
\begin{split}
    \int_K \left( \alpha \frac{\partial (\nabla \cdot \mathbf{u}_h)}{\partial t} + \frac{1}{M} \frac{\partial p_h}{\partial t} - \nabla \cdot (k_h \nabla p_h) \right) v_h \, dx &= \int_K f v_h \, dx, \\
    \int_K \boldsymbol{\sigma}_h (\mathbf{u}_h) : \nabla \mathbf{v}_h \, dx &= \int_K \alpha \nabla p_h \cdot \mathbf{v}_h \, dx.
\end{split}
\end{equation}

To handle discontinuities between neighboring elements, numerical flux terms are introduced. The flux for the pressure \(p_h\) across an interior face \(E\) between two elements \(K^+\) and \(K^-\) is defined as:
\[
\{ \nabla p_h \} = \frac{1}{2} \left( \nabla p_h^+ + \nabla p_h^- \right), \quad [p_h] = p_h^+ - p_h^-.
\]

The full weak form of the system in the mixed space, including the flux terms for the pressure and displacement, becomes:

\begin{equation}
\begin{split}
    &\sum_{K \in T_h} \int_K \left( \alpha \frac{\partial (\nabla \cdot \mathbf{u}_h)}{\partial t} + \frac{1}{M} \frac{\partial p_h}{\partial t} - \nabla \cdot (k_h \nabla p_h) \right) v_h \, dx \\
    &+ \sum_{E \in \mathcal{E}_h} \int_E \left( \{ k_h \nabla p_h \cdot \mathbf{n} \} \cdot [v_h] + \{ k_h \nabla v_h \cdot \mathbf{n} \} \cdot [p_h] \right) \, ds \\
    &= \sum_{K \in T_h} \int_K f v_h \, dx + \sum_{E \in \mathcal{E}_h} \int_E \frac{\gamma_f}{h} k_h [p_h] [v_h] \, ds, \\
    &\sum_{K \in T_h} \int_K \boldsymbol{\sigma}_h (\mathbf{u}_h) : \nabla \mathbf{v}_h \, dx \\
    &+ \sum_{E \in \mathcal{E}_h} \int_E \left( \{ \boldsymbol{\sigma}_h \} : [\mathbf{v}_h] + \{ \mathbf{v}_h \} : [\boldsymbol{\sigma}_h] \right) \, ds  \\ 
    &= \sum_{K \in T_h} \int_K \alpha \nabla p_h \cdot \mathbf{v}_h \, dx + \sum_{E \in \mathcal{E}_h} \int_E \frac{\gamma_f}{h} [\boldsymbol{\sigma}_h] [\mathbf{v}_h] \, ds.
\end{split}
\end{equation}

For the time discretization, we use an implicit backward Euler scheme. Let \(t_n\) denote the \(n\)-th time step, and \(\Delta t\) be the time step size. The time derivatives of displacement and pressure are approximated as:

\[
\frac{\partial p_h}{\partial t} \approx \frac{p_h^{n+1} - p_h^n}{\Delta t}, \quad \frac{\partial \mathbf{u}_h}{\partial t} \approx \frac{\mathbf{u}_h^{n+1} - \mathbf{u}_h^n}{\Delta t}.
\]

Substituting these into the weak form of the governing equations, we get:
\begin{equation}
\begin{split}
    &\sum_{K \in T_h} \int_K \left( \alpha \frac{\nabla \cdot \mathbf{u}_h^{n+1} - \nabla \cdot \mathbf{u}_h^n}{\Delta t} + \frac{1}{M} \frac{p_h^{n+1} - p_h^n}{\Delta t} - \nabla \cdot (k_h \nabla p_h^{n+1}) \right) v_h \, dx \\
    &+ \sum_{E \in \mathcal{E}_h} \int_E \left( \{ k_h \nabla p_h^{n+1} \cdot \mathbf{n} \} \cdot [v_h] + \{ k_h \nabla v_h \cdot \mathbf{n} \} \cdot [p_h^{n+1}] \right) \, ds \\
    &= \sum_{K \in T_h} \int_K f^{n+1} v_h \, dx + \sum_{E \in \mathcal{E}_h} \int_E \frac{\gamma_f}{h} k_h [p_h^{n+1}] [v_h] \, ds, \\
    &\sum_{K \in T_h} \int_K \boldsymbol{\sigma}_h (\mathbf{u}_h^{n+1}) : \nabla \mathbf{v}_h \, dx \\
    &+ \sum_{E \in \mathcal{E}_h} \int_E \left( \{ \boldsymbol{\sigma}_h^{n+1} \} : [\mathbf{v}_h] + \{ \mathbf{v}_h \} : [\boldsymbol{\sigma}_h^{n+1}] \right) \, ds \\
    &= \sum_{K \in T_h} \int_K \alpha \nabla p_h^{n+1} \cdot \mathbf{v}_h \, dx + \sum_{E \in \mathcal{E}_h} \int_E \frac{\gamma_f}{h} [\boldsymbol{\sigma}_h^{n+1}] [\mathbf{v}_h] \, ds.
\end{split}
\end{equation}

The resulting system of equations can be written in matrix form as:
\[
A_h \mathbf{U}_h^{n+1} = \mathbf{F}_h^{n+1},
\]
where \(\mathbf{U}_h^{n+1}\) is the vector of unknowns, consisting of displacement \(\mathbf{u}_h^{n+1}\) and pressure \(p_h^{n+1}\) at the current time step, and \(A_h\) is the system matrix derived from the weak form.

The system matrix incorporates both the mass balance and elasticity equations in the mixed form, with flux terms ensuring weak continuity across element boundaries.

\section{Coarse Grid Approximation using DG-GMsFEM}

In this section, we describe the coarse grid approximation using the Discontinuous Galerkin Generalized Multiscale Finite Element Method (DG-GMsFEM) for the system of poroelasticity. 
The computational domain is discretized into a coarse grid, and multiscale basis functions are constructed to solve the system of equations on this grid. 
We present the algorithm separately for pressure and elasticity, followed by details of the construction of multiscale basis functions.


Let $\Omega$ be the computational domain, and let $T_H$ denote the partition of $\Omega$ into a coarse grid, where $H$ is the characteristic size of the coarse grid elements. The domain $\Omega$ is divided into a set of non-overlapping coarse elements:
\[
T_H = \bigcup_{i=1}^{N_{\text{cell}}} K_i,
\]
where $N_{\text{cell}}$ is the number of coarse grid elements. In Figure \ref{coarse}, we present illustration of the coarse grid.

\begin{figure}[t!]
\centering
\includegraphics[width=0.5\textwidth]{pics/coarsegrid.png}
\caption{Coarse grid $T_H$ with coarse cell $K_i$}
\label{coarse}
\end{figure}

Defining $V_H$ and $Q_H$ as the multiscale space for displacement and pressure:
\[
V_H = \text{span} \{ \phi_i \}_{i=1}^{N_u}, \quad Q_H = \text{span} \{ \psi_i \}_{i=1}^{N_p}
\]
where  
$N_u =\text{dim}(V_H)$ is the number of basis functions for displacement and $N_p =\text{dim}(Q_H)$ is the number of basis functions for pressure.
  
For the coarse grid approximation, we use a DG approach and have the following variational formulation: 
\begin{itemize}
\item Elasticity problem: find $u_H \in V_H$ such that
\begin{equation}
\label{ms}
a(u_H, v) = l(v), \quad \forall v \in V_H.
\end{equation} 
\item Pressure problem: find $p_H \in W_H$ such that
\begin{equation}
\label{ms}
a(p_H, \mathbf{v}) = l(\mathbf{v}), \quad \forall \mathbf{v} \in Q_H.
\end{equation} 
\end{itemize}

Note that the coarse-scale system can be formed by projecting the fine-scale system onto the coarse grid.  The projection matrix can be assembled using the multiscale basis functions. 

We aim to solve the poroelasticity system on this coarse grid by constructing multiscale basis functions within each coarse element $K_i$. The system is solved for pressure and displacement separately, and the details for each are provided below.

To construct the multiscale space, we start with the construction of the snapshot space that contains a set of basis functions formed by the solution of local problems with all possible boundary conditions up to the fine grid resolution in each coarse cell $K_i$ (local domain)  for $i=1, \cdots, N$, where $N$ is the number of coarse blocks in $\Omega$. After that, we solve a spectral problem to select dominant modes of the snapshot space.  


\subsection{Multiscale Basis Functions for Pressure}

The construction of multiscale basis functions for pressure is a crucial step in the DG-GMsFEM framework. These basis functions are designed to capture fine-scale features of the pressure within each coarse element and represent the pressure solution at the coarse scale.

The generation of the multiscale basis functions for pressure consists of the following key steps:
\begin{enumerate}
    \item Generation of a local snapshot space.
    \item Reduction of the dimension of the snapshot space by solving a local spectral problem.
\end{enumerate}

For each coarse element \(K_i \in T_H\), we generate a local snapshot space, which is used to capture the fine-scale variations of the pressure field within \(K_i\). The local snapshot space is constructed by solving local fine-scale problems with different boundary conditions.

The local snapshot functions are solutions of the following local problem:

\[
\begin{aligned}
    -\nabla \cdot (k \nabla p^i_l) &= 0, \quad x \in \, K_i, \\
    p^i_l &= \delta^i_j, \quad x \in \Gamma^i,
\end{aligned}
\]
where $l= 1, \cdots,  J^g_i$, $J^g_i$ is the number 
of fine grid facets on $\Gamma^i$, and \(\delta^i_j\) represents the Kronecker delta function that has value $1$ if $i = l$ and value $0$ otherwise.

The local snapshot space is denoted by \(V_{\text{snap}}^{i}\):
\[
Q_{\text{snap}}^{i} = \{p^i_l:1 \leq l \leq J_i\},
\]
and the snapshot space projection matrix is defined as:
\[
R_{\text{snap}}^{i} = \left[ p_1^{i}, \ldots, p_{J^g_i}^{i} \right]^{T}.
\]


To reduce the dimension of the local snapshot space, we solve a local spectral problem within each coarse element. This step allows us to extract the most significant modes that capture the fine-scale behavior while reducing the computational cost.

The local spectral problem is formulated as follows:
\begin{equation}
\label{eq:spectral-out}
\tilde{A}^{K_i}_g \tilde{\psi}^{i} =  
\lambda_j \tilde{S}^{K_i} \tilde{\psi}^{i}, 
\end{equation}
where:
\begin{itemize}
    \item $\tilde{A}^{K_i}$ is the stiffness matrix defined on the snapshot space, given by $\tilde{A}^{K_i} = R^{i, \text{snap}} A_h^{K_i} (R_{\text{snap}}^{i})^{T}$.
    \item $\tilde{S}^{K_i}$ is the mass matrix defined on the snapshot space, given by $\tilde{S}^{K_i} = R_{\text{snap}}^{i} S_h^{K_i} (R_{\text{snap}}^{i})^{T}$.
    \item \(\lambda_j\) and $ \tilde{\psi}^{i}$ are the eigenvalues and eigenvectors, respectively, representing the solution to the spectral problem.
\end{itemize}

Here $A_h^{K_i}$ and $S_h^{K_i}$ are matrix representations of the bilinear forms $a^{K_i}(p,\mathbf{v})$ and $s^{K_i}(p,\mathbf{v})$.


In the DG-GMsFEM framework, the bilinear forms include terms that account for jumps and averages across the boundaries of fine grid elements inside each coarse element \(K_i\). The bilinear forms are written as:

\[
a^{K_i}(p, \mathbf{v}) = \sum_{K \in T_h(K_i)} \int_K k_h \nabla p \cdot \nabla \mathbf{v} \, dx 
\]
\[
- \sum_{E \in \mathcal{E}^h(K_i)} \int_E \left( \{ k_h \nabla p \cdot \mathbf{n} \} \cdot [\mathbf{v}] + \{ k_h \nabla \mathbf{v} \cdot \mathbf{n} \} \cdot [p] - \frac{\gamma_f}{h} \{k_h\} [p] \cdot [\mathbf{v}] \right) ds,
\]

\[
s^{K_i}(p, \mathbf{v}) = \sum_{E \in \mathcal{E}^h(K_i)} \int_E k_h p \mathbf{v} \, ds.
\]

We remark that the integral in $s^{K_i}(p,\mathbf{v})$ is defined on the boundary of the coarse block, and $T_h$ is the fine grid for the local domain $K_i$.


After solving the local spectral problem, we select the dominant eigenfunctions corresponding to the smallest eigenvalues \(\lambda_j\). These eigenfunctions form the multiscale basis functions \(\psi_j^p\) for pressure in each coarse element \(K_i\):

\[
Q_H = \text{span} \{ 
\psi_{k}^i: \, 1 \leq i \leq N_{cell}^H, \,  1 \leq k \leq M^i_p \},
\]
where $\psi_{k}^i= (R_{\text{snap}}^{i})^{T} \tilde{\psi}^i_{k}$.

We also add one interior basis function, constructed by solving the following local problem:
\begin{equation}
\label{eq:snap-int}
-\nabla \cdot (k \nabla p^i_l) = f, \quad x \in  K_i,
\end{equation}
with $f = 1$ and homogeneous boundary conditions:
\[
p^i = 0, \quad x \in \Gamma^i.
\]

Thus, the final multiscale space for pressure is:
\[
Q_H 
= \text{span} \{ 
\psi_{k}^i, \psi^{i}: \, 
1 \leq i \leq N_{cell}^H, \,   
1 \leq k \leq M^i_p \}.
\]

To construct the coarse grid system, we generate a projection matrix using the multiscale basis functions:
\[
R_p = \left[
\psi_{1}^{1},  \ldots, \psi_{M^i_p}^{N_{cell}^H}, 
\psi_1 , \ldots, \psi_{N_{cell}^H}
\right]^{T}.
\]

Using the projection matrix, we obtain the following coarse grid system in matrix form for pressure:
\[
A_H p_H  =  F_H,
\]
where:
\[
A^p_H = R A^p_h R^{T}, \quad 
F^p_H = R F^p_h.
\]
After solving the coarse-scale system, the fine-scale solution for pressure is recovered as \(p_{ms}= R^{T} p_H\).

\subsection{Multiscale Basis Functions for Displacement}

The construction of multiscale basis functions for displacement is similar to the process for pressure, but now we consider the displacement field \( \mathbf{u} \) as a vector. These basis functions are designed to capture fine-scale features of the displacement within each coarse element and represent the displacement solution on the coarse grid.

The construction of the multiscale basis functions for displacement consists of the following key steps:
\begin{enumerate}
    \item Generation of a local snapshot space for the displacement vector \( \mathbf{u} \).
    \item Reduction of the dimension of the snapshot space by solving a local spectral problem.
\end{enumerate}

For each coarse element \(K_i \in T_H\), we generate a local snapshot space that captures the fine-scale variations of the displacement field \( \mathbf{u} \) within \(K_i\). The local snapshot space is constructed by solving local fine-scale elasticity problems with different boundary conditions. The snapshot functions are solutions of the following local problems:

\[
\begin{aligned}
    -\nabla \cdot \boldsymbol{\sigma}(\mathbf{u}^i_l) &= 0, \quad x \in K_i, \\
    \mathbf{u}^i_l &= g_i^l, \quad x \in \Gamma^i,
\end{aligned}
\]
where \( l = 1, \dots, d \cdot J^g_i \), $d$ is the dimension, \( J^g_i \) is the number of fine grid facets on \(\Gamma^i\), $g_i^l = (\delta_i^l, 0)$ and $(0, \delta_i^l)$ for $d = 2$.

The local snapshot space is denoted as \( V_{\text{snap}}^{i} \):
\[
V_{\text{snap}}^{i} = \{\mathbf{u}^i_l : 1 \leq l \leq 2 \cdot J^g_i\},
\]
and the snapshot space projection matrix is defined as:
\[
R_{\text{snap}}^{i} = \left[ \mathbf{u}_1^{i}, \ldots, \mathbf{u}_{2 \cdot J^g_i}^{i} \right]^{T}.
\]

To reduce the dimension of the local snapshot space, we solve a local spectral problem within each coarse element for the displacement vector \( \mathbf{u} \). The local spectral problem is formulated as follows:

\[
\tilde{A}^{K_i}_g \tilde{\psi}^{i} = \lambda_j \tilde{S}^{K_i} \tilde{\psi}^{i},
\]
where:
\begin{itemize}
    \item \(\tilde{A}^{K_i}\) is the stiffness matrix defined on the snapshot space and is given by:
    \[
    \tilde{A}^{K_i} = R_{\text{snap}}^{i} A_h^{K_i} (R_{\text{snap}}^{i})^{T}.
    \]
    \item \(\tilde{S}^{K_i}\) is the mass matrix defined on the snapshot space and is given by:
    \[
    \tilde{S}^{K_i} = R_{\text{snap}}^{i} S_h^{K_i} (R_{\text{snap}}^{i})^{T}.
    \]
    \item \(\lambda_j\) and \(\tilde{\psi}^{i}\) are the eigenvalues and eigenvectors representing the solution to the spectral problem.
\end{itemize}

The bilinear forms for the displacement vector \( \mathbf{u} \) are written as:

\[
a^{K_i}(\mathbf{u}, \mathbf{v}) = \sum_{K \in T_h(K_i)} \int_K \boldsymbol{\sigma}(\mathbf{u}) \cdot \nabla \mathbf{v} \, dx
\]
\[
- \sum_{E \in \mathcal{E}^h(K_i)} \int_E \left( \{ \boldsymbol{\sigma}(\mathbf{u}) \cdot \mathbf{n} \} \cdot [\mathbf{v}] + \{ \boldsymbol{\sigma}(\mathbf{v}) \cdot \mathbf{n} \} \cdot [\mathbf{u}] - \frac{\gamma_f}{h} \{ \boldsymbol{\sigma} \} [\mathbf{u}] \cdot [\mathbf{v}] \right) ds,
\]

\[
s^{K_i}(\mathbf{u}, \mathbf{v}) = \sum_{E \in \mathcal{E}^h(K_i)} \int_E (\lambda + 2 \mu) \mathbf{u} \cdot \mathbf{v} \, ds.
\]

After solving the local spectral problem, we select the dominant eigenfunctions corresponding to the smallest eigenvalues \(\lambda_j\). These eigenfunctions form the multiscale basis functions \(\phi_j^{\mathbf{u}}\) for displacement in each coarse element \(K_i\):

\[
V_H = \text{span} \{ \phi_{k}^i : 1 \leq i \leq N_{cell}^H, 1 \leq k \leq M^i_u \},
\]
where $\phi_{k}^i = (R_{\text{snap}}^{i})^{T} \tilde{\psi}^i_{k}$.

We also add an interior basis function, constructed by solving the following local problem:
\begin{equation}
\label{eq:snap-int-elasticity}
-\nabla \cdot \boldsymbol{\sigma}(\mathbf{u}^i_l) = \mathbf{f}, \quad x \in K_i,
\end{equation}
with \( \mathbf{f} = (1, 1) \) and homogeneous boundary conditions:
\[
\mathbf{u}^i = 0, \quad x \in \Gamma^i.
\]

Thus, the final multiscale space for elasticity is:
\[
V_H = \text{span} \{ \psi_{k}^i, \psi^{i} : 1 \leq i \leq N_{cell}^H, 1 \leq k \leq M^i_u \}.
\]


To construct the coarse grid system for elasticity, we generate a projection matrix using the multiscale basis functions:
\[
R_u = \left[ \psi_{1}^{1}, \ldots, \psi_{M^i_u}^{N_{cell}^H}, \psi_1, \ldots, \psi_{N_{cell}^H} \right]^{T}.
\]

Using the projection matrix, we obtain the following coarse grid system in matrix form for elasticity:
\[
A_H \mathbf{u}_H = F_H,
\]
where:
\[
A^u_H = R_u A^u_h R_u^{T}, \quad F^u_H = R_u F^u_h.
\]
After solving the coarse-scale system, the fine-scale solution for displacement can be recovered as \( \mathbf{u}_{ms} = R_u^{T} \mathbf{u}_H \).

\section{Numerical Results}

In this section, we present numerical results for two test cases to demonstrate the effectiveness of the DG-GMsFEM method for solving the poroelasticity equations. In both cases, we compare the reference fine-scale solution with the multiscale solution obtained using the DG-GMsFEM approach. We consider two scenarios:
\begin{itemize}
    \item Test 1: Homogeneous coefficients for elasticity and permeability.
    \item Test 2: Heterogeneous coefficients for elasticity and permeability.
\end{itemize}

The computational domain is defined as \( \Omega = [0, L_x] \times [0, L_y] \), where \( L_x = L_y = 1 \). Both the coarse and fine grids were generated using the gmsh software \cite{20}. The coarse grid contains 121 vertices and 100 cells, while the fine grid consists of 12,699 vertices and 25,396 cells.

For Test 1, the material properties are homogeneous with elasticity coefficient \( E = 1 \) and Poisson's ratio \( \nu = 0.3 \). The Biot coefficient is \( \alpha = 1 \). The boundary conditions for the displacement are defined as follows:
\[ 
u_x = 0, \, \sigma_y = 0 \quad x \in \Gamma_1 \cup \Gamma_3, 
\]
\[ 
u_y = 0, \, \sigma_x = 0 \quad x \in \Gamma_2 \cup \Gamma_4. 
\]
For the pressure, the boundary conditions are set as:
\[ 
p = 1.0 \quad \text{on} \, \Gamma_2, \quad \nabla p \cdot \mathbf{n} = 0 \quad \text{on} \, \Gamma_1 \cup \Gamma_3 \cup \Gamma_4. 
\]

The numerical simulations were carried out using the FEniCS library implemented in C++ \cite{21}. The multiscale coarse grid solution was obtained using \( m = 20 \) multiscale basis functions, and the fine-scale reference solution was used for comparison.

\begin{figure}[t!]
\centering
\includegraphics[width=0.5\textwidth]{pics/grid.png}
\caption{Coarse grid (red color) and fine grid (green color).}
\label{fig1}
\end{figure}

We first solve for displacement and pressure in the homogeneous case. Figure \ref{fig1} shows the fine and coarse grids used for the simulations. The fine grid is used for the reference solution, while the multiscale coarse grid is used for the DG-GMsFEM solution.

Figures \ref{fig2}, \ref{fig3}, and \ref{fig4} show the distribution of displacement along the $X$ and $Y$ directions, and the pressure distribution at three different time steps, \( n = 1, 10, 20 \). The time steps \( n \) represent different moments during the simulation. The top row in each figure corresponds to the reference fine-scale solution, and the bottom row corresponds to the multiscale solution with \( m = 20 \) multiscale basis functions.

\begin{figure}[t!]
\centering
\includegraphics[width=0.8\linewidth]{pics/ux}
\caption{Distribution of displacement along $X$ at $n = 1, 10, 20$ for test 1 (from left to right). Top: reference solution, $DOF_u^h=52050$. Bottom: multiscale solution with 20 multiscale basis functions, $DOF_u^H=4100$.}
\label{fig2}
\end{figure}

\begin{figure}[t!]
\centering
\includegraphics[width=0.8\linewidth]{pics/uy}
\caption{Distribution of displacement along $Y$ at $n = 1, 10, 20$ for test 1 (from left to right). Top: reference solution, $DOF_u^h=52050$. Bottom: multiscale solution with 20 multiscale basis functions, $DOF_u^H=4100$.}
\label{fig3}
\end{figure}

\begin{figure}[t!]
\centering
\includegraphics[width=0.8\linewidth]{pics/p}
\caption{Distribution of pressure at $n = 1, 10, 20$ for test 1 (from left to right). Top: reference solution, $DOF_p^h=26025$. Bottom: multiscale solution with 20 multiscale basis functions, $DOF_p^H=2100$.}
\label{fig4}
\end{figure}

We also evaluate the accuracy of the multiscale solution by calculating the relative \( L_2 \) error for both displacement and pressure. The relative \( L_2 \) error is defined as:
\[ 
e_{L_2} = \frac{|| u_{\text{ref}} - u_{\text{ms}} ||_{L_2}}{|| u_{\text{ref}} ||_{L_2}},
\]
where \( u_{\text{ref}} \) is the fine-scale reference solution and \( u_{\text{ms}} \) is the multiscale solution.

Figure \ref{error1} shows the relative \( L_2 \) errors for displacement (left) and pressure (right) as a function of the number of multiscale basis functions \( m \). The error decreases as the number of basis functions increases, demonstrating the effectiveness of the multiscale method. For \( n=20, m=40 \), the relative error in displacement reduces to around $2.38\%$. Similarly, for pressure, the relative error reduces to $11.45\%$ for \( n=20, m=40 \).

\begin{figure}[t!]
\centering
\includegraphics[width=0.45\linewidth]{pics/test1-error-Elasticity-L2}
\includegraphics[width=0.45\linewidth]{pics/test1-error-Pressure-L2}
\caption{Thermoelasticity problem for test 1. Relative $L_2$ error for displacement (left) and pressure (right) with a different number of multiscale basis functions \( m \).}
\label{error1}
\end{figure}

In the second test, the elasticity coefficient \(E\) and permeability \(k\) are heterogeneous. The distribution of these coefficients is shown in Figure \ref{kdmuE}. As in the first test, we compare the reference fine-scale solution with the multiscale solution using \( m = 20 \) multiscale basis functions.

\begin{figure}[h!]
\centering
\includegraphics[width=0.4\linewidth]{pics/E.png}
\includegraphics[width=0.4\linewidth]{pics/kdmu.png}
\caption{Elasticity coefficient $E$ (left), heterogeneous permeability $k$.}
\label{kdmuE}
\end{figure}

Figures \ref{fig5}, \ref{fig6}, and \ref{fig7} show the distribution of displacement along the $X$ and $Y$ directions, and the pressure distribution at three different time steps, \( n = 1, 10, 20 \), for Test 2. The top row in each figure corresponds to the reference fine-scale solution, and the bottom row corresponds to the multiscale solution with \( m = 20 \) multiscale basis functions.

\begin{figure}[t!]
\centering
\includegraphics[width=0.8\linewidth]{pics/ux2}
\caption{Distribution of displacement along $X$ at $n = 1, 10, 20$ for test 2 (from left to right). Top: reference solution, $DOF_u^h=52050$. Bottom: multiscale solution with 20 multiscale basis functions, $DOF_u^H=4100$.}
\label{fig5}
\end{figure}

\begin{figure}[t!]
\centering
\includegraphics[width=0.8\linewidth]{pics/uy2}
\caption{Distribution of displacement along $Y$ at $n = 1, 10, 20$ for test 2 (from left to right). Top: reference solution, $DOF_u^h=52050$. Bottom: multiscale solution with 20 multiscale basis functions, $DOF_u^H=4100$.}
\label{fig6}
\end{figure}

\begin{figure}[t!]
\centering
\includegraphics[width=0.8\linewidth]{pics/p2}
\caption{Distribution of pressure at $n = 1, 10, 20$ for test 2 (from left to right). Top: reference solution, $DOF_p^h=26025$. Bottom: multiscale solution with 20 multiscale basis functions, $DOF_p^H=2100$.}
\label{fig7}
\end{figure}

The displacement and pressure distributions for Test 2 are similar to those shown for Test 1. The relative \( L_2 \) error for both displacement and pressure is shown in Figure \ref{error2}. As with the homogeneous case, the error decreases as the number of multiscale basis functions \( m \) increases, demonstrating the accuracy of the method even for heterogeneous coefficients. For \( n = 20, m=40 \), the relative error reduces to around $4.34\%$. Similarly, for pressure, the relative error reduces to $5.35\%$.

\begin{figure}[t!]
\centering
\includegraphics[width=0.45\linewidth]{pics/test2-error-Elasticity-L2}
\includegraphics[width=0.45\linewidth]{pics/test2-error-Pressure-L2}
\caption{Thermoelasticity problem for test 2. Relative $L_2$ error for displacement (left) and pressure (right) with a different number of multiscale basis functions \( m \).}
\label{error2}
\end{figure}

The numerical results demonstrate the robustness and accuracy of the DG-GMsFEM method for solving the poroelasticity equations. In both homogeneous and heterogeneous cases, the method achieves significant computational efficiency by reducing the number of degrees of freedom, while maintaining acceptable accuracy levels as the number of multiscale basis functions \( m \) increases. This shows the potential of DG-GMsFEM for efficient and reliable modeling of complex poroelasticity problems, even in heterogeneous media.


\section{Conclusions}
In this paper, we applied the Discontinuous Galerkin Generalized Multiscale Finite Element Method (DG-GMsFEM) to solve the poroelasticity equations. By constructing multiscale basis functions for pressure and displacement, the method effectively captures fine-scale features while reducing the computational cost associated with solving on a fine grid. We considered two test cases: one with homogeneous coefficients and one with heterogeneous coefficients for elasticity and permeability. In both cases, the multiscale method demonstrated high accuracy, with relative $L_2$ errors for pressure and displacement decreasing as the number of multiscale basis functions increased. The computational efficiency of the DG-GMsFEM method allows for significant cost savings while maintaining accuracy, making it a promising approach for solving complex poroelasticity problems in heterogeneous media. 
\bigskip
\begin{thebibliography}{1}

\bibitem{1} M. A. Biot, 
\href{https://doi.org/10.1063/1.1712886}{\it General theory of three-dimensional consolidation}, Journal of Applied Physics, \textbf{12}:2 (1941), 155--164.

\bibitem{2} J. M. Carcione, C. Morency, J. E. Santos, 
\href{https://doi.org/10.1190/1.3474602}{\it Computational poroelasticity-A review}, Geophysics, \textbf{75}:5 (2010), 75A229--75A243.

\bibitem{3} A. H. D. Cheng, 
\href{https://doi.org/10.1016/B978-0-08-040615-2.50011-3}{\it Poroelasticity}, Switzerland: Springer International Publishing, \textbf{27} (2016).

\bibitem{4} E. Detournay, A. H. D. Cheng, 
{\it Fundamentals of poroelasticity}, Analysis and Design Methods, Pergamon (1993), 113--171.

\bibitem{5} U. Kalachikova, A. Tyrylgin, M. Vasilyeva, D. Ammosov, 
\href{https://doi.org/10.1016/j.cam.2022.114977}{\it Generalized Multiscale Finite Element Method for scattering problem in heterogeneous media}, Journal of Computational and Applied Mathematics, \textbf{424} (2023), 114977.

\bibitem{6} Y. Efendiev, J. Galvis, T.Y. Hou,
\href{https://10.1016/j.jcp.2013.04.045}{\it Generalized multiscale finite element methods}, Journal of computational physics, {\bf 251} (2013), 116--135.

\bibitem{7} E. T. Chung, Y. Efendiev, G. Li, 
\href{https://doi.org/10.1016/j.jcp.2014.05.007}{\it An adaptive GMsFEM for high-contrast flow problems}, Journal of Computational Physics, \textbf{273} (2014), 54--76.

\bibitem{8} V. Alekseev, U. Gavrileva, D. Spiridonov, A. Tyrylgin, M. Vasilyeva, 
\href{https://doi.org/10.1063/1.5064930}{\it Numerical simulation of the transport and flow problems in perforated domains using generalized multiscale finite element method}, AIP Conference Proceedings, AIP Publishing, \textbf{2025}:1 (2018).

\bibitem{9} V. Alekseev, M. Vasilyeva, U. Kalachikova, E. T. Chung,  
\href{https://doi.org/10.3390/computation9070075}{\it DG-GMsFEM for problems in perforated domains with non-homogeneous boundary conditions}, Computation, \textbf{9}:7 (2021), 75.9.

\bibitem{10} K. Gao, S. Fu, R. Gibson Jr, E. T. Chung, Y. Efendiev, 
\href{https://doi.org/10.1016/j.jcp.2015.03.068}{\it Generalized multiscale finite-element method (GMsFEM) for elastic wave propagation in heterogeneous, anisotropic media}, Journal of Computational Physics, \textbf{295} (2015), 161--188.

\bibitem{11} E. T. Chung, Y. Efendiev, W. T. Leung, 
\href{https://doi.org/10.1137/140986189}{\it An adaptive generalized multiscale discontinuous Galerkin method for high-contrast flow problems}, Multiscale Modeling \& Simulation, \textbf{16}:3 (2018), 1227--1257.

\bibitem{12} E. T. Chung, U. Kalachikova, M. Vasilyeva, V. Alekseev,
\href{https://doi.org/10.1016/j.matcom.2021.11.001}{\it Generalized multiscale discontinuous Galerkin method for convectionпїЅdiffusion equation in perforated media}, Mathematics and Computers in Simulation, \textbf{193} (2022), 666--688.

\bibitem{13} S. Fu, E.T. Chung, T. Mai, 
\href{https://doi.org/10.1016/j.jcp.2020.109569}{\it Constraint energy minimizing generalized multiscale finite element method for nonlinear poroelasticity and elasticity}, Journal of Computational Physics, \textbf{417} (2020), 109569.

\bibitem{14} A. Tyrylgin, M. Vasilyeva, D. Spiridonov, E.T. Chung, 
\href{https://doi.org/10.1016/j.cam.2020.112783}{\it Generalized multiscale finite element method for the poroelasticity problem in multicontinuum media}, Journal of Computational and Applied Mathematics, \textbf{374} (2020), 112783.

\bibitem{15} D. L. Brown, M. A. Vasilyeva, 
\href{https://doi.org/10.1016/j.cam.2015.08.007}{\it A generalized multiscale finite element method for poroelasticity problems I: Linear problems}, Journal of Computational and Applied Mathematics, \textbf{294} (2016), 372--388.

\bibitem{16} D. Ammosov, M. Vasilyeva, E. T. Chung, 
\href{https://doi.org/10.1016/j.cam.2021.113995}{\it Generalized Multiscale Finite Element Method for thermoporoelasticity problems in heterogeneous and fractured media}, Journal of Computational and Applied Mathematics, \textbf{407} (2022), 113995.

\bibitem{17} A. Tyrylgin, M. Vasilyeva, A. Alikhanov, 
\href{https://doi.org/10.1007/978-3-030-97020-8_14}{\it Multiscale Model Reduction for the Poroelasticity Problems Using Embedded Fracture Model}, International Conference on Mathematics and its Applications in New Computer Systems, Springer International Publishing, Cham (2021), 153--162.

\bibitem{18} A. Tyrylgin, J. Huang, 
\href{https://doi.org/110.1134/S1995080223100426}{\it Online Coupled Generalized Multiscale Finite Element Method for the Poroelasticity Problem in Three-Dimensional Media}, Lobachevskii Journal of Mathematics, \textbf{44}:10 (2023), 4183--4194.

\bibitem{19} A. Tyrylgin, M. Vasilyeva, D. Ammosov,  E. T. Chung,Y. Efendiev, 
\href{https://doi.org/10.3390/fluids6080298}{\it Online Coupled Generalized Multiscale Finite Element Method for the Poroelasticity Problem in Fractured and Heterogeneous Media}, Fluids, \textbf{6}:8 (2021), 298.

\bibitem{20} C. Geuzaine, J.-F. Remacle, 
\href{https://doi.org/10.1002/nme.2579}{\it Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities}, International Journal for Numerical Methods in Engineering, \textbf{79}:11 (2009), 1309--1331.

\bibitem{21} A. Logg, K.-A. Mardal, G. Wells, eds., 
{\it Automated solution of differential equations by the finite element method: The FEniCS book}, Vol. 84, Springer Science \& Business Media, 2012.


\end{thebibliography}
\end{document}
