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

\usepackage{tikz}
\usepackage{multirow}
\graphicspath{ {./figs/} }

\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} %Здесь автор определяет УДК своей работы
\def\mscs{65M60,} %Здесь автор определяет классификаторы AMS своей работы
\setcounter{page}{144}

\newtheorem{lemma}{Лемма}
\newtheorem{theorem}{Теорема}
\newtheorem{definition}{Определение}
\newtheorem{corollary}{Следствие}
\newtheorem{proposition}{Предложение}


\definecolor{myblue}{rgb}{0.03, 0.27, 0.79}

\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
     	{\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{\figurename}{Fig.}
\renewcommand{\tablename}{Table}

\thispagestyle{empty}

\title[\sc Online Meshfree Multiscale Method]{\large Online Meshfree Multiscale Method for flow with Forchheimer law in fractured media}
\author[D.Y. Nikiforov, S.P. Stepanov]{
{\bf \large \textsc{D.Y. Nikiforov\protect\orcidicon{0000-0002-5415-6284}, S.P. Stepanov\protect\orcidicon{0000-0001-9445-2726}}}
\\ \\
\textit{Представлено} \textsc{П.П. Петровым}
%\textit{(заполняется редактором)}
}%
\address{Djulustan Yakovlevich Nikiforov  % Контактные данные всех авторов и место работы указываются только на английском языке
\newline\hphantom{iii} North-Eastern Federal University,
\newline\hphantom{iii} 58 Belinsky str.,
\newline\hphantom{iii} 677000, Yakutsk, Russia}%
\email{\textcolor{blue}{dju92@mail.ru}}%

\address{Sergei Pavlovich Stepanov  % Контактные данные всех авторов и место работы указываются только на английском языке
\newline\hphantom{iii} North-Eastern Federal University,
\newline\hphantom{iii} 58 Belinsky str.,
\newline\hphantom{iii} 677000, Yakutsk, Russia}%
\email{\textcolor{blue}{cepe2a@inbox.ru}}%

\thanks{\sc Nikiforov D.Y., Stepanov S.P., 
Online Meshfree Multiscale Method}
\thanks{\copyright \ 2024 Nikiforov D.Y, Stepanov S.P}
\thanks{\rm Работа поддержана грантом РНФ (№ 24-71-00125)}
\thanks{\it Поступила 1 января 2023 г., опубликована 31 декабря 2023 г.}%


\semrtop \vspace{1cm}
\maketitle {\small
\begin{quote}
\noindent{\bf Abstract:} This paper presents a model of filtration process in a fractured porous medium. The filtration velocity in a porous medium is described by Darcy law, and in fractures is described by Forchheimer law. To describe filtration through fractures on a fine grid, Discrete Fracture Model (DFM) is used. To reduce the size of the discrete problem, Online Meshfree Generalized Multiscale Finite Element Method (Online MFGMsFEM) is developed. The results of numerical calculations showed higher accuracy relative to the offline multiscale method.\medskip

\noindent{\bf Keywords:} filtration process, fractured porous media, Forchheimer law, multiscale method.
 \end{quote}
}

\bigskip

\section{Introduction}
Traditional approaches to modeling fluid flow through porous media are based on Darcy's law, which describes laminar flow conditions. However, in geological formations, the presence of fractures can significantly alter flow dynamics, often resulting in high-velocity fluid movement that exceeds the assumptions of Darcy's law. In such cases, Forchheimer's law becomes relevant, accounting for the inertial effects that arise in non-laminar flow conditions \cite{Forchheimer1901}.

In this study, we present an approach that integrates Darcy's law for flow through porous media with Forchheimer's law for flow through fractures. By addressing the distinct flow laws in these two domains, we aim to enhance the understanding of filtration processes in heterogeneous media. Our approach involves the development of a mathematical model that combines Darcy and Forchheimer flow regimes, allowing for a more correct prediction of fluid behavior in complex fractured porous structures.

The Forchheimer model is a complex and nonlinear, which makes direct simulation computationally demanding, particularly for large-scale applications. Efficient solvers are therefore essential for accelerating forward simulations of these nonlinear systems. This study aims to address this need by developing solvers to enhance computational feasibility for applications.

In fractured media, high-resolution grids are critical for capturing multiscale variations, leading to increased computational costs due to the high number of degrees of freedom. Model reduction techniques are essential for managing these demands. Various local multiscale model reduction methods have been proposed, including upscaling techniques \cite{Wu2002,Efendiev} and multiscale approaches \cite{Chung2010,Chung2018}, which have shown promise in reducing computational cost.

To address nonlinear problems, the GMsFEM incorporates an efficient online basis strategy, as exemplified by the online GMsFEM \cite{Chung2015}. Instead of frequently updating computationally expensive offline basis functions, this approach adaptively introduces basis functions that minimize the residual in each local domain. By generating local online multiscale basis functions based on residual data, this method reduces dependence on an extensive offline basis while achieving faster error reduction. This approach improves accuracy with fewer basis functions and effectively accommodates temporal variations in coefficients, enhancing computational efficiency.

Based on GMsFEM \cite{ms,gmsfem1,Chung2023,step}, this study applies the meshfree approach MFGMsFEM \cite{mgmsfem, Nikiforov2023, Nikiforov2024, Nikiforov2024_2, Nikiforov2024_3}, facilitating efficient handling of complex geometries and heterogeneous domains to enhance simulation performance and computational efficiency. Using a finite element method approximation, we present numerical results for a two-dimensional fractured domain, specifically examining the impact of nonlinearity on Online MFGMsFEM accuracy. While multiple online strategies can address nonlinear problems, our focus here is on a simple, efficient approach.

The paper is organized as follows: Section 2 presents the problem formulation, while Sections 3 and 4 outline the proposed methodology. Section 5 discusses the numerical results, and Section 6 concludes with a summary.

\section{Problem formulation}
In this paper we consider the problem of modeling the flow of a single-phase fluid in fractured porous media. The mathematical model is described by a coupled system of parabolic differential equations for the pressure $p_m$ in the porous matrix $\Omega$ and the pressure $p_f$ in the fractures $\gamma$ (Fig. \ref{ris:domain})

\begin{equation}\label{eq:system1}\begin{split}
c_m \frac{\partial p_m}{\partial t} + \nabla \cdot \left( v_m \right) + T_{mf} (p_m-p_f) = g_m, \quad x \in \Omega, \\
c_f \frac{\partial p_f}{\partial t} + \nabla \cdot \left( v_f \right) + T_{mf} (p_f-p_m) = g_f, \quad x \in \gamma, 
\end{split}\end{equation} 
where $c_i$ is the media compressibility coefficient, $v_i$ is the velocity ($i = m, f$), $g_i$ is the source term and $T_{mf}$ is a mass transfer coefficient. 

\begin{figure}[!ht]
\begin{center}
\includegraphics[width=0.6\linewidth]{domain}
\end{center}
\caption{Computational domain $\Omega$ with fractures $\gamma$.}
\label{ris:domain}
\end{figure}

In the system \eqref{eq:system1}, the fluid velocity $v_m$ in porous medium $\Omega$ is described by Darcy's law
\begin{equation}\label{eq:eq2}
v_m = -\frac{k_m}{\mu} \nabla p_m,
\end{equation} 
where $k_m$ is a porous medium permeability and $\mu$ is a viscosity.

We assume that the fluid is filtrating at high velocity through the fractures, so we will use Forchheimer's law (for more details, see \cite{Auriault2007})
\begin{equation}\label{eq:eq3}
v_f = -\frac{K^F(p_f)}{\mu} \nabla p_f,
\end{equation} 
where $K^F(p_f)$ is a Forchheimer permeability, which is described by the following formula
\begin{equation}\label{eq:eq4}
K^F(p_f) = \frac{2k_f}{1 + \sqrt{1 + \frac{4 \beta \rho k_f^2}{\mu^2} |\nabla p_f|}}, \quad |\nabla p_f| = \sqrt{\nabla p_f \cdot \nabla p_f}.
\end{equation}
Here $k_f$ is the fracture permeability, $\rho$ is the fluid density and $\beta \geqslant 0$ an inertial coefficient.

By putting \eqref{eq:eq2} and \eqref{eq:eq3} into the system \eqref{eq:system1} taking into account \eqref{eq:eq4} we obtain the following system
\begin{equation}\label{eq:system2}\begin{split}
c_m \frac{\partial p_m}{\partial t} - \nabla \cdot \left( \frac{k_m}{\mu} \nabla p_m \right) + T_{mf} (p_m-p_f) = g_m, \quad x \in \Omega, \\
c_f \frac{\partial p_f}{\partial t} - \nabla \cdot \left( \frac{K^F(p_f)}{\mu} \nabla p_f \right) + T_{mf} (p_f-p_m) = g_f, \quad x \in \gamma. 
\end{split}\end{equation} 

The system of equations \eqref{eq:system2} is supplemented by initial and boundary conditions
\begin{equation}\label{eq:eq5}\begin{split}
p_i = p_0, \quad &t = 0,    
\\
p_i = p_D, \quad &x \in \xi.
\end{split}\end{equation} 

Note that the coefficient $\beta$ is an indicator of the nonlinearity of the filtration rate $v_f$. When $\beta = 0$ we obtain Darcy's law with $K^F = k_f$.

\section{Fine-grid discretization}
For spatial approximation on a fine grid, we use the continuous finite element method with the usual linear basis functions (FEM). For the fractures accounting we use the discrete fracture model (DFM) \cite{dfm1, dfm2}. Thus the first equation in the system \eqref{eq:system2} will be discretized in two-dimensional space for the porous matrix and the second equation in one-dimensional space for the fractures. Next, we will make the transition to a model with one pressure. Let the characteristic time of mass exchange between the pressures in the pores and fractures significantly exceed the characteristic times of the filtration processes. This means that $T_{mf} \rightarrow \infty, \ p_m \rightarrow p, \ p_f \rightarrow p$, where $p$ is the total pressure of the system. Adding up the equations of the system \eqref{eq:system2} and setting $p_m = p_f = p$ (see \cite{Laevsky}), we obtain the following variational problem: find $p \in H^1(\Omega)$ such that
\begin{equation}\label{eq:app1}
m(\frac{\partial p}{\partial t}, v) + a(p, v) = f(v), \quad \forall v \in H^1(\Omega),
\end{equation} 
where
\begin{equation}\label{eq:app2}\begin{split}
m = m_1 + m_2, \quad a = a_1 + a_2, \quad f = f_1 + f_2,\\
m_1(\frac{\partial p}{\partial t}, v) = \int\limits_{\Omega} c_m \frac{\partial p}{\partial t}vdx, \quad m_2(\frac{\partial p}{\partial t}, v) = \alpha \int\limits_{\gamma} c_f \frac{\partial p}{\partial t}vds, \\
a_1(p, v) = \int\limits_{\Omega} \left( \frac{k_m}{\mu} \nabla p, \nabla v \right)dx, \quad a_2(p, v) = \alpha\int\limits_{\gamma} \left( \frac{K^F(p)}{\mu} \nabla p, \nabla v \right)ds, \\
f_1(v) = \int\limits_{\Omega} g_m v dx, \quad f_2(v) = \alpha\int\limits_{\gamma} g_f v ds. 
\end{split}
\end{equation}
Here, when integrating over fractures $\gamma$, we take into account their aperture $\alpha$.

Let us define an uniform, for simplicity, time grid $\omega_t = \{ t^n = n \tau, \ n = 0, 1, \ldots, N_t-1, \ \tau N_t = t_{max} \}$, where $\tau$ - time step, $N_t$ - number of time steps and $t_{max}$ - solution time. Here we use implicit time discretization for the problem \eqref{eq:app1}. Then we obtain the following variational problem: find $p^{n+1} \in H^1(\Omega)$ such that 
\begin{equation}\label{eq:app11}
m(\frac{p^{n+1}-p^{n}}{\tau}, v) + a(p^{n+1}, v) = f(v), \quad \forall v \in H^1(\Omega), \quad n = 0, 1, \ldots, N_t-1,
\end{equation} 
where $p^{n} \approx p(t^n)$. Here in $a_2(p, v)$ for nonlinear coefficient $K^F(p)$ we use Picard iterations method
\begin{equation}
a_2(p^{n+1}, v) = \alpha\int\limits_{\gamma} \left( \frac{K^F(p^j)}{\mu} \nabla p^{n+1}, \nabla v \right)ds,
\nonumber
\end{equation}
where $j$ is the index of iterations. Thus, at each time point $n+1$ we also solve the problem \eqref{eq:app11} until the completion of the iterations according to the following criterion
\begin{equation}
\frac{||p^{j+1}-p^j||_{H^1}}{||p^{j+1}||_{H^1}} \cdot 100 \le \epsilon,
\nonumber
\end{equation}
where ${||u||}_{H^1} = \sqrt{\int_{\Omega} (\nabla u,\nabla u) dx}$. In our implementation, we took $\epsilon = 1\%$. At the beginning of iterations $j=0$, we took $p^j=p^n$, and at the end of the iterations we set $p^{n+1}=p^{j+1}$. 

We build an unstructured mesh $\mathcal{T}_h$ and construct discrete function space $V_h \subset V$. Then we use $p_h = (p_{h, 1}, p_{h, 2}, \ldots, p_{h, Nf})^T$ to denote the vector of the required unknowns, where $N_f$ -- the number of fine mesh vertices. Next we write the following matrix form for the fully discrete system 
\begin{equation}\label{eq:app3}
M \frac{p_h^{n+1}-p_h^{n}}{\tau} + A p_h^{n+1} = F,
\end{equation}
where
\begin{equation}\label{eq:app4}\begin{split}
M = \{m_{ij}\}, \quad m_{ij} = m_1(\phi_j, \phi_i) + m_2(\hat{\phi}_j, \hat{\phi}_i), \\ 
A = \{a_{ij}\}, \quad a_{ij} = a_1(\phi_j, \phi_i) + a_2(\hat{\phi}_j, \hat{\phi}_i), \\
F = \{f_{i}\}, \quad f_i = f_1(\phi_i) + f_2(\hat{\phi}_i),
\end{split}\end{equation}
where $\phi_i$ -- two-dimensional and $\hat{\phi}_i$ -- one-dimensional fine-scale piecewise linear basis functions for the porous medium and fractures.

\section{Meshfree GMsFEM}
The algorithm of MFGMsFEM can be divided into two stages: offline and online. In the offline stage, we  construct multiscale basis functions that take into account fractures. The basis functions are computed in local domains $S_i$. In offline stage we have the following steps
\begin{enumerate}
\item Generation of meshfree coarse-scale partitioning (point cloud).
\item Computation multiscale basis functions and construction a multiscale offline space.
\item Building a projection matrix into the multiscale offline space.
\end{enumerate}

In the online stage, the fine-grid system is converted to a coarse grid system using the resulting projection matrix. As a result, the problem \eqref{eq:app2} is solved as a coarse-scale problem. In online stage we have the following steps
\begin{enumerate}
\item For the current time step, project the fine-scale system to the coarse-scale using the multiscale projection matrix.
\item Solving the coarse-scale problem.
\item Reconstruction the solution on the fine grid also using the multiscale projection matrix.
\item Moving to the next time step.
\end{enumerate}


\subsection{Multiscale basis functions}
In MFGMsFEM we use a point cloud instead of a structured coarse grid. The point cloud is built on top of a fine mesh, so it is convenient to construct it relative to the most meshfree methods.

\begin{figure}[!ht]
\begin{center}
\includegraphics[width=0.6\linewidth]{point_cloud}
\end{center}
\caption{Point cloud illustration obtained with the centroidal Voronoi tessellations method.}
\label{ris:point_cloud}
\end{figure}

Let $\mathcal{S}_H$ be a partition of the computational domain $\Omega$ to the point cloud so that $\Omega \subset \bigcup_{i=1}^{N} S_i$ and suppose that each coarse element $S_i$ is partitioned into a connected union of elements of a fine grid (see Fig. \ref{ris:online}). Let ${\{ x_i \}}^{N}_{i=1}$ is the coarse-scale nodes, where $N$ denotes the number of coarse nodes. Here, the coarse elements $S_i$ are the supports of basis functions
\begin{equation}
S_i = \{ y \in \mathcal{R}^N : \ \parallel y - x_i \parallel \le r_i \},
\nonumber
\end{equation} 
where $r_i$ is radius of coarse element $S_i$.

To calculate the location of nodes $x_i$, we use the version of CVT (centroidal Voronoi tessellations) proposed in \cite{cvt1}. An important point here is the choice of the distribution density function. Here we calculate the distribution density function $\rho(x)$ of a random variable taking into account the fractures. Next, we use the algorithm proposed in \cite{cvt2} for finding the radii $r_i$ for $S_i$.

We denote the basis functions by $\psi_{i, k}$, which is supported in $S_i$, and the index $k$ represents the numbering of these basis functions. In turn, the solution will be sought as
\begin{equation}\label{eq:ms_sol}
p_{\text{c}}(x) = \sum_{i, k} c_{i, k} \psi_{i, k}(x).
\end{equation} 

Once the basis functions have been identified (see the next section), the global coupling of the CG is given through the variational form
\begin{equation}\label{eq:ms_problem}
m(\frac{p_{\text{c}}^{n+1}-p_{\text{c}}^{n}}{\tau}, v) + a(p_{\text{c}}^{n+1}, v) = f(v), \quad \forall v \in V_{\text{off}}, \quad n = 0, 1, \ldots, N_t-1,
\end{equation} 
where $V_{\text{off}}$ is multiscale function space.

\textbf{Local basis functions.} To construct the offline space $V^{S_i}_{\text{off}}$, we solve the following local eigenvalue problem in each $S_i$
\begin{equation}\label{eq:spectral_problem}
A \Psi^{\text{off}} = \lambda^{\text{off}} B \Psi^{\text{off}},
\nonumber
\end{equation}
where $A$ and $B$ denote similar fine-scale matrices as
\begin{equation}
\begin{split}
A = \{a_{ij}\}, \quad a_{ij} = \int_{S} \left(\frac{k_m}{\mu} \nabla \phi_j, \nabla \phi_i \right) dx + \alpha \int_{\gamma^S} \left( \frac{k_f}{\mu} \nabla \hat{\phi}_j, \nabla \hat{\phi}_i \right) ds,\\ 
B = \{b_{ij}\}, \quad b_{ij} = \int_{S} \left(\frac{k_m}{\mu} \phi_j, \phi_i \right) dx 
+ \alpha \int_{\gamma^S} \left(\frac{k_f}{\mu} \hat{\phi}_j, \hat{\phi}_i \right) ds,
\end{split} 
\nonumber
\end{equation}
where $\phi_i$ and $\hat{\phi}_i$ are fine-scale basis functions for the matrix $S_i$ and fractures $\gamma^{S_i}$ subdomains. Note that we use $k_f$ here instead of $K^F$, since offline basis functions only take into account spatial heterogeneity without the nonlinear part.

For definition of the multiscale basis functions, we select the first $M$ eigenvectors $\psi_k^{S_i}, k = 0, \ldots M$ corresponding to the first $M$ smallest eigenvalues.

\textbf{Global formulation.} In the meshfree multiscale method the shape functions $W_i(x)$ defined in $S_i$ form the initial coarse space
\begin{equation} 
V^{\text{init}}_0 = \text{span} \{ W_i(x) \}_{i=1}^{N}.
\nonumber
\end{equation}
Here the shape functions $W_i(x)$ are defined as
\begin{equation}
W_i(x) = \frac{\phi_i(x)}{\sum_{j=1}^N \phi_j(x)},
\nonumber
\end{equation}
where $\phi_i(x)$ are kernel functions, which here are cubic splines
\begin{equation}
\phi(r) = 2 \left\{
\begin{array}{lr}
2/3 + 4 \ (r-1) \ r^2,& r \le 0.5,\\
4/3 \ (1-r)^3,& 0.5 \le r \le 1,\\
0,& 1 \le r,
\end{array}
\right.
\nonumber
\end{equation}
where $r$ is the normalized distance.

Accordingly, by multiplying the form function $W_i(x)$ by the eigenvectors $\psi_k^{S_i}$, we obtain the basis functions in the offline space $V_{\text{off}}$ 
\begin{equation}
\psi_{i, k} = W_i \psi_k^{S_i}, \ 1 \le i \le N, \ 1 \le k \le M,
\nonumber
\end{equation}
where $k$ is an index of the eigenvector. Next, we construct the spectral multiscale space
\begin{equation}
V_{\text{off}} = \text{span} \{ \psi_{i, k}: 1 \le i \le N, \ 1 \le k \le M \}.
\nonumber
\end{equation}
Also, for further use, an operator matrix is constructed
\begin{equation}
R^T = [ \psi_1, \ldots, \psi_{N_c} ],
\nonumber
\end{equation}
where $N_c = N \cdot M$ and $\psi_i$ are used to denote the nodal values of each basis function defined on the fine grid.

The transition matrix $R$ is then used to solve the following system on the coarse scale
\begin{equation}\label{eq:coarse_problem} 
M_\text{c} \frac{p_\text{c}^{n+1}-p_\text{c}^{n}}{\tau} + A_\text{c} p_\text{c}^{n+1} = F_\text{c},
\end{equation}
where 
\begin{equation}
M_\text{c} = R M R^T, \quad A_\text{c} = R A R^T, \quad F_\text{c} = R F.
\nonumber
\end{equation}

After solving the system \eqref{eq:coarse_problem}, we can go from a coarse-scale solution to a fine-scale solution using also the transition operator $R$ and the solution $p_\text{c}$
\begin{equation}\label{eq:coarse_problem2} 
p^{n+1}_\text{ms} = R^T p^{n+1}_{c}, 
\end{equation}
where $p^{n+1}_\text{ms}$ is a fine-grid projection of the coarse-grid solution.

\begin{figure}[h!]
\begin{center}
\begin{minipage}[h]{0.325\linewidth}
\center{\includegraphics[width=\linewidth]{pou}}
\end{minipage}
\begin{minipage}[h]{0.325\linewidth}
\center{\includegraphics[width=\linewidth]{evec}}
\end{minipage}
\begin{minipage}[h]{0.325\linewidth}
\center{\includegraphics[width=\linewidth]{basis}}
\end{minipage}
\end{center}
\caption{Illustration of offline basis function. Left: shape function $W_i(x)$. Center: eigen vector $\psi_k^{S_i}$. Right: offline multiscale basis function $\psi_{i,k}$.}
\label{ris:offline}
\end{figure}

\subsection{Online enrichment of multiscale space}
In the online stage of GMsFEM, the projection matrix $R_n$ is updated at the $n$-th time step by incorporating online residual-based multiscale basis functions. 

For predefined heterogeneities and fracture distributions, a fixed projection matrix $R$ is used initially ($n = 0$), and the system is solved without updating $R_l$ at every time step. However, from the first time step ($n = 1, 2, \ldots$), the multiscale space is enriched by adding online residual-based multiscale basis functions $\vartheta^i_1$ calculated locally in $S_i$ using current residuals. This enrichment updates the projection matrix to:
\[
(R_l)^T = [ \psi_1, \ldots, \psi_{N_c}, \vartheta_1^l, \ldots, \vartheta_N^l].
\]
where $N_c = N \cdot M$, $N$ denotes the number of coarse nodes, $M$ the number of offline basis functions on each $S_i$ and $l$ denotes the index of online iterations. The process can be repeated by changing the online basis functions at every $n$-th time step for $l = n$ or at selected intervals (e.g., every 5th time step) to control the computational load. This enrichment of the multiscale space with online basis functions allows achieving adaptive accuracy and computational efficiency, especially for nonlinear problems with time variability.

Construction of the local residual based online multiscale basis functions is based on the solution of the following local problem in $S_i$:
\begin{equation}\label{eq:online}
\begin{split}
m(\frac{\Phi_i^l}{\tau}, v) + a(\Phi_i^l, v) = r(v), \\
r(v) = f(v) + m(\frac{p^n - p^l}{\tau}, v) - a(p^l, v),
\end{split}
\end{equation}
where $r(v)$ is the residual, $p^l$ is the solution obtained by solving the system \eqref{eq:coarse_problem}--\eqref{eq:coarse_problem2} with $R$.

The problem \eqref{eq:online} is supplemented by zero Dirichlet boundary condition $\Phi_i^l = 0$ on the $\partial S_i$. Finally, we obtain an online multiscale basis functions by multiplication on the partition of unity function (shape function) $ \vartheta^l_i = \chi_i \Phi_i^l$ (Fig. \ref{ris:online}). 

Thus, the general algorithm of the online stage at each time step is as follows (we will update $R_l$ at every 5th time layer)
\begin{itemize}
\item If $n\%5 = 0$: 
\begin{itemize}
\item Сalculation $p^l$ by solving \eqref{eq:coarse_problem}--\eqref{eq:coarse_problem2} with $R$.
\item Contruction $R_l$ by solving problem \eqref{eq:online} using $p^l$.
\item Using $R_l$ solve \eqref{eq:coarse_problem}--\eqref{eq:coarse_problem2}.
\end{itemize}
\item Else: Using $R_l$ solve \eqref{eq:coarse_problem}--\eqref{eq:coarse_problem2}.
\item Transition to the next time step.
\end{itemize}
Note that we use Picard iterations only when solving problem \eqref{eq:coarse_problem}--\eqref{eq:coarse_problem2} using $R_l$.

\begin{figure}[h!]
\begin{center}
\begin{minipage}[h]{0.325\linewidth}
\center{\includegraphics[width=\linewidth]{frac}}
\end{minipage}
\begin{minipage}[h]{0.325\linewidth}
\center{\includegraphics[width=\linewidth]{phi}}
\end{minipage}
\begin{minipage}[h]{0.325\linewidth}
\center{\includegraphics[width=\linewidth]{onbasis}}
\end{minipage}
\end{center}
\caption{Illustration of online basis function. Left: fine grid of local $S_i$ with fractures (red). Center: local function $\Phi_i^l$. Right: local online multiscale basis function $\vartheta^l_i$.}
\label{ris:online}
\end{figure}

\section{Numerical results}
In this section, numerical results are presented to test the proposed meshfree multiscale approach to solve a nonlinear 2D problem of single-phase filtration problem with discrete fractures. For Forchheimer law, we test an inertial coefficient with different values $\beta$ = 0, 1 and 1000. Note that, for $\beta = 0$ we obtain Darcy's law with linear permeability $K^F = k_f$.

For numerical experiments of the model problem, we use the following parameters $c_m = 1.0e-8$, \ $c_f = 1.0e-8$, \ $k_m = 1.0e-12$, \ $k_f = 1.0e-5$, \ $\alpha = 0.01$, \ $\rho = 850$, \ $g_m = 0$, \ $g_f = 0$, \ $p_0 = 40.0e6$, \ $p_D = 30.0e6$, \ $\tau = 1~\text{day}$, \ $T = 360~\text{days}$. The fine grid has 53875 vertices and 106947 elements. For multiscale method we use $N = 256$ coarse-scale vertices. 

In the numerical results, we call the solutions obtained by the meshfree multiscale method using only offline basis functions the offline solutions, and the solutions obtained by the meshfree multiscale method with offline and online basis functions the online solutions. Also, we call the solutions obtained by the finite element method on a fine grid the reference solutions.

Figures \ref{ris:sol1}, \ref{ris:sol2}, and \ref{ris:sol3} show the solutions at the final time for $\beta = 0$, $\beta = 1$, and $\beta = 1000$, respectively. Note that for a large value of $\beta = 1000$, Forchheimer's law applies and the flow velocity in fractures decreases significantly. Here, for the meshfree mutiscale method, we use 2 offline basis functions $M = 2$ on each local $S_i$. In the offline multiscale method, the degrees of freedom are 512 ($DOF_c = 512$). In the online multiscale method, we use $DOF_c = 768$, since we use 2 offline and 1 online basis fucntions on each $S_i$. As we can see, the online solutions are similar to the reference solutions.

\begin{figure}[h!]
\begin{center}
\begin{minipage}[h]{0.325\linewidth}
\center{\includegraphics[width=\linewidth]{fine0}}
\end{minipage}
\begin{minipage}[h]{0.325\linewidth}
\center{\includegraphics[width=\linewidth]{offline0}}
\end{minipage}
\begin{minipage}[h]{0.325\linewidth}
\center{\includegraphics[width=\linewidth]{online0}}
\end{minipage}
\end{center}
\caption{Pressure at the final time for $\beta = 0$. Left: reference solution. Center: offline solution with $M = 2$. Right: online solution with $M = 2$ and 1 online basis functions.}
\label{ris:sol1}
\end{figure}

\begin{figure}[h!]
\begin{center}
\begin{minipage}[h]{0.325\linewidth}
\center{\includegraphics[width=\linewidth]{fine1}}
\end{minipage}
\begin{minipage}[h]{0.325\linewidth}
\center{\includegraphics[width=\linewidth]{offline1}}
\end{minipage}
\begin{minipage}[h]{0.325\linewidth}
\center{\includegraphics[width=\linewidth]{online1}}
\end{minipage}
\end{center}
\caption{Pressure at the final time for $\beta = 1$. Left: reference solution. Center: offline solution with $M = 2$. Right: online solution with $M = 2$ and 1 online basis functions.}
\label{ris:sol2}
\end{figure}

\begin{figure}[h!]
\begin{center}
\begin{minipage}[h]{0.325\linewidth}
\center{\includegraphics[width=\linewidth]{fine1000}}
\end{minipage}
\begin{minipage}[h]{0.325\linewidth}
\center{\includegraphics[width=\linewidth]{offline1000}}
\end{minipage}
\begin{minipage}[h]{0.325\linewidth}
\center{\includegraphics[width=\linewidth]{online1000}}
\end{minipage}
\end{center}
\caption{Pressure at the final time for $\beta = 1000$. Left: reference solution. Center: offline solution with $M = 2$. Right: online solution with $M = 2$ and 1 online basis functions.}
\label{ris:sol3}
\end{figure}

To numerically evaluate the accuracy of the meshfree multiscale method with and without online basis functions, we compare solutions by multiscale methods with reference solutions using the following norms
\[ 
{||e||}_{L^2} = \frac{{||u_1-u_2||}_{L^2}}{{||u_1||}_{L^2}}\cdot 100, \quad
{||e||}_{H^1} = \frac{{||u_1-u_2||}_{H^1}}{{||u_1||}_{H^1}}\cdot 100,
\]
where ${||u||}_{L^2} = \sqrt{\int_{\Omega} u^2 dx}, \ {||u||}_{H^1} = \sqrt{\int_{\Omega} (\nabla u,\nabla u) dx}$, $u_1$ -- reference solution and $u_2$ -- multiscale solution. In Table \ref{tab:table}, we present the $L_2$ and $H_1$ errors between the reference solution and the multiscale solutions at the final time $T$. Here, we see a good agreement between the multiscale solutions and the reference solutions. We also observe that the errors decay as basis functions are added to each local $S_i$. 

For the online method, we update of the online basis functions every five time steps. Adding online basis functions can significantly improve accuracy. For example, for $\beta = 1$ if one adds online basis functions to four offline basis functions, one obtains the following result the relative $H^1$ error decreases from about 7.24\% to 1.32\%. Thus, one obtains better accuracy with fewer degrees of freedom $DOF_c$.

\begin{table}[h!]
\begin{center}
\begin{tabular}{*{8}{ c| @{\hspace{3.5mm}}}}
& \multirow{2}{*}{$M$} & \multicolumn{3}{c}{Offline} & \multicolumn{3}{c}{Online}\\
& & $DOF_c$ & ${||e||}_{L^2}$ & ${||e||}_{H^1}$ & $DOF_c$ & ${||e||}_{L^2}$ & ${||e||}_{H^1}$\\
\hline
\multirow{4}{*}{$\beta = 0$} & 1 & 256 & 7.17 & 96.53 & 512 & 2.27 & 47.64 \\
& 2 & 512  & 0.98  & 31.83 & 768  & 1.4e-02 & 3.06 \\
& 4 & 1024 & 5.6e-02 & 7.25 & 1280 & 4.4e-04 & 0.25 \\
& 8 & 2048 & 1.6e-03 & 1.03 & 2304 & 3.9e-05 & 4.6e-02 \\
& 16 & 4096 & 1.8e-04 & 0.23 & 4352 & 7.02e-06 & 1.4e-02 \\
\hline
\multirow{4}{*}{$\beta = 1$} & 1 & 256 & 2.26 & 60.34 & 512 & 0.49 & 32.48 \\
& 2 & 512  & 0.43 & 28.83 & 768 & 1.75e-02 & 4.45 \\
& 4 & 1024 & 4.6e-02 & 7.24 & 1280 & 3.13e-03 & 1.32 \\
& 8 & 2048 & 4.1e-03 & 1.67 & 2304 & 3.02e-04 & 0.39 \\
& 16 & 4096 & 1.3e-03 & 0.83 & 4352 & 1.2e-04 & 0.22 \\
\hline
\multirow{4}{*}{$\beta = 10^3$} & 1 & 256 & 0.93 & 54.1 & 512 & 0.63 & 47.45 \\
& 2 & 512  & 0.53  & 41.08 & 768  & 0.16 & 30.62 \\
& 4 & 1024 & 0.36  & 32.63  & 1280 & 3.24e-02 & 17.92 \\
& 8 & 2048 & 9.4e-02 & 24.32  & 2304 & 1.26e-02 & 10.3 \\
& 16 & 4096 & 1.8e-02 & 16.01 & 4352 & 2.5e-03 & 6.66
\end{tabular}
\end{center}
\caption{Error at the final time.}
\label{tab:table}
\end{table}

\section{Conclusion}
This paper introduces an effective model for simulating the filtration process in fractured porous media by integrating Darcy's law with Forchheimer's law. We developed an online multiscale approach based on the Generalized Multiscale Finite Element Method (GMsFEM) to enhance solution accuracy. In this approach, online multiscale basis functions are computed by solving local problems derived from the residual, allowing adaptive refinement. The use of the Online Meshfree GMsFEM significantly improves calculation accuracy, as demonstrated by numerical experiments showcasing the efficiency of our algorithm, particularly in the online enrichment of the multiscale space. This method has proven effective for solving flow problems in complex, fractured, and heterogeneous media, offering valuable insights for analyzing filtration processes in intricate geological settings with discrete fractures.

\bigskip
\begin{thebibliography}{1} % Ссылки на литературу приводятся только на английском языке (даже если у источника нет англоязычной версии).

\bibitem{Forchheimer1901} P.Z. Forchheimer, {\it Wasserbewegung durch Boden}, Zeitschrift des Vereins Deutscher Ingenieure, \textbf{45} (1901), 1781--1788. 

\bibitem{Wu2002} X.H. Wu, Y. Efendiev and T.Y. Hou, \href{https://doi.org/10.3934/dcdsb.2002.2.185}{\it Analysis of upscaling absolute permeability}, Discrete and Continuous Dynamical Systems --B, \textbf{2}:2 (2002), 158--204.

\bibitem{Efendiev} Y. Efendiev, V. Ginting and T.Y. Hou, {\it Multiscale Finite Element Methods for Nonlinear Problems and Their Applications}. Communications in Mathematical Sciences, \textbf{2}:4 (2004), 553--589.

\bibitem{Chung2010} E.T. Chung and Y. Efendiev, \href{https://doi.org/10.1137/09075620X}{\it Reduced-Contrast Approximations for High-Contrast Multiscale Flow Problems}, Multiscale Modeling \& Simulation, \textbf{8}:4 (2010), 1128--1153.

\bibitem{Chung2018} E.T. Chung, W. Leung, M. Vasilyeva and Y. Wang, \href{https://doi.org/10.1016/j.cam.2017.08.017}{\it Multiscale model reduction for transport and flow problems in perforated domains}, Journal of Computational and Applied Mathematics, \textbf{330} (2018), 519--535.

\bibitem{Chung2015} E.T. Chung, Y. Efendiev and W.T. Leung, \href{https://doi.org/10.1016/j.jcp.2015.07.068}{\it Residual-driven online generalized multiscale finite element methods}, Journal of Computational Physics, \textbf{302} (2015), 176--190.

\bibitem{ms} Y. Efendiev and T.Y. Hou, 
\href{https://doi.org/10.1007/978-0-387-09496-0}{\it Multiscale finite element methods: theory and applications}, Springer Science \& Business Media, \textbf{4} (2009).

\bibitem{gmsfem1} Y. Efendiev, J. Galvis and T.Y. Hou, \href{https://doi.org/10.1016/j.jcp.2013.04.045}{\it Generalized multiscale finite element methods}, Journal of Computational Physics, \textbf{251} (2013), 116--135.

\bibitem{Chung2023} E.T. Chung, Y. Efendiev and T.Y.Hou,  \href{https://doi.org/10.1007/978-3-031-20409-8}{\it Multiscale Model Reduction: Multiscale Finite Element Methods and Their Generalizations}, Springer Nature, \textbf{212} (2023).

\bibitem{step} S. Stepanov, D. Nikiforov and A. Grigorev, 
\href{https://doi.org/10.3390/math9202545}{\it Multiscale Multiphysics Modeling of the Infiltration Process in the Permafrost}, Mathematic, \textbf{9}:20 (2021), 2545.

\bibitem{mgmsfem} D.~Nikiforov, 
\href{https://doi.org/10.1016/j.jcp.2022.111798}{\it Meshfree Generalized Multiscale Finite Element Method}, Journal of Computational Physics, \textbf{474} (2023), 111798.

\bibitem{Nikiforov2023} D.Y. Nikiforov and S.P. Stepanov, \href{https://doi.org/10.1134/S1995080223030241}{\it Modeling of Artificial Ground Freezing Using a Meshfree Multiscale Method}, Lobachevskii Journal of Mathematics, \textbf{44}:3 (2023), 1206--1214.

\bibitem{Nikiforov2024} D. Nikiforov, S. Stepanov and N. Lazarev, \href{https://doi.org/10.1016/j.cam.2024.115988}{\it Meshfree multiscale method for the infiltration problem in permafrost}, Journal of Computational and Applied Mathematics, \textbf{449} (2024), 115988.

\bibitem{Nikiforov2024_2} D. Nikiforov and S. Stepanov, \href{https://doi.org/10.1016/j.cam.2024.116020}{\it Meshfree multiscale method with partially explicit time discretization for nonlinear Stefan problem}, Journal of Computational and Applied Mathematics, \textbf{451} (2024), 116020.

\bibitem{Nikiforov2024_3} D. Nikiforov, L.A. Poveda, D. Ammosov, Y. Sarmiento and J. Galvis, \href{https://doi.org/10.1016/j.cam.2024.116367}{\it Meshfree Generalized Multiscale Exponential Integration Method for Parabolic Problems}, Journal of Computational and Applied Mathematics, \textbf{459} (2025), 116367.

\bibitem{Auriault2007} J.L. Auriault, C. Geindreau and L. Org\'eas, 
\href{https://doi.org/10.1007/s11242-006-9096-x}{\it Upscaling Forchheimer law}, Transport in Porous Media, \textbf{70}:2 (2007), 213--229.

\bibitem{dfm1} M.V. Vasilyeva, V.I. Vasilyev, A.A. Krasnikov and D.Ya. Nikiforov, 
\href{}{\it Numerical simulation of single-phase fluid flow in fractured porous media}, Uchenye Zapiski Kazanskogo Universiteta. Seriya Fiziko-Matematicheskie Nauki, \textbf{159}:1 (2017), 100--115. 

\bibitem{dfm2} V.I. Vasil'ev, M.V. Vasil'eva, V.S. Gladkikh, V.P. Ilin, D.Ya. Nikiforov, D.V. Perevozkin and G.A. Prokop'ev, 
\href{https://doi.org/10.1134/S199047891804018X}{\it Numerical solution of a fluid filtration problem in a fractured medium by using the domain decomposition method}, Journal of Applied and Industrial Mathematics, \textbf{12}:4 (2018), 785--796. 

\bibitem{Laevsky} M.I. Ivanov, I.A. Kremer, Y.M. Laevsky,  
\href{https://doi.org/10.1016/j.cam.2023.115639}{\it Explicit-implicit schemes for non-isothermal filtration problem: Single-temperature model}, Journal of Computational and Applied Mathematics, \textbf{440} (2024), 115639.

\bibitem{cvt1} L. Ju, Q. Du and M. Gunzburger, 
\href{https://doi.org/10.1016/S0167-8191(02)00151-5}{\it Probabilistic methods for centroidal Voronoi tessellations and their parallel implementations}, Parallel Computing, \textbf{28}:10 (2002), 1477--1500.

\bibitem{cvt2} M. Griebel and M.A. Schweitzer, 
\href{https://doi.org/10.1137/S1064827599355840}{\it A particle-partition of unity method for the solution of elliptic, parabolic, and hyperbolic PDEs}, SIAM Journal on Scientific Computing, \textbf{22}:3 (2000), 853--890.



\end{thebibliography}
\end{document}
