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


\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{\ }]


\def\udcs{519.612} %Здесь автор определяет УДК своей работы
\def\mscs{65F08} %Здесь автор определяет классификаторы AMS своей работы
\setcounter{page}{144}

\newtheorem{theorem}{Теорема}
\newtheorem{remark}{Замечание}

\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}
            \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
}

\providecommand{\doi}[1]{doi:\discretionary{}{}{}#1}

\begin{document}
\renewcommand{\refname}{References}
\thispagestyle{empty}

\title[\sc On the iterative solution of the Stokes problem]{\large On the iterative solution of the Stokes problem}
\author{
    А.М.~Гурин$^{1}$,
    В.П.~Ильин$^{2}$,
    Д.И.~Козлов$^{2, 3}$,
    Е.А.~Кузьмин$^{2, 3}$
}

\thanks{\rm Работа поддержана РНФ (грант 24-21-00402)}

\semrtop \vspace{1cm}
\maketitle 
\begin{center}
    $^{1)}$ Институт гидродинамики имени М. А. Лаврентьева СО РАН, \\
    $^{2)}$ Институт вычислительной математики и математической геофизики СО РАН, \\
    $^{3)}$ Новосибирский государственный университет\\
    \email{guralisk@mail.ru}
\end{center}
{\small
    \begin{quote}
        \noindent{\bf Abstract:} Block preconditioned iterative conjugate gradient methods for solving the three-dimensional Stokes problem are investigated. A formulation with a computational domain in the form of a parallelepiped is considered. “Standard” approximations on a cubic staggered grid, seven-point and two-point for the Laplace and derivative operators are used. The resulting saddle-type SLAE is regularized to ensure the uniqueness of the solution. The block preconditioner is constructed by the method of incomplete factorization with diagonal compensation and using band approximations for matrices inverse to the Schur complement and the grid Laplace operator, as well as using an algebraic multigrid algorithm. Examples of numerical experiments on a representative series of methodological problems using parallel algorithms on different numbers of processors are given. The issues of generalizing the proposed approaches to broader classes of problems are considered.\medskip
        
        \noindent{\bf Keywords:} Stokes problem, large sparse SLAEs, iterative preconditioned methods, Krylov subspace, Schur complement.
    \end{quote}
}

\bigskip

% % \begin{comment}
% % \begin{abstract}
% % Исследуются блочные предобусловленные итерационные методы сопряжённых градиентов для решения трёхмерной задачи Стокса. Рассматривается постановка с расчётной областью в форме параллелепипеда. Применяются ``стандартные'' аппроксимации на кубической разнесённой сетке, семиточечные и двухточечные для операторов Лапласа и дифференцирования. Получаемая СЛАУ седлового типа регуляризуется с обеспечением единственности решения. Построение блочного предобуславливателя производится методом неполной факторизации с диагональной компенсацией и с использованием ленточных приближений для матриц, обратных к дополнению Шура и сеточному оператору Лапласа, а также с применением алгебраического многосеточного алгоритма. Приводятся примеры численных экспериментов на представительной серии методических задач с использованием параллельных алгоритмов на различном числе процессоров. Рассматриваются вопросы обобщения предложенных подходов на более широкие классы задач.
% % \end{abstract}
% % \end{comment}

% % {\bf Ключевые слова:} задача Стокса, большие разрежённые СЛАУ, итерационные предобусловленные методы, подпространство Крылова, дополнение Шура.

\section{Введение}

Задача Стокса, получаемая из системы уравнений Навье-Стокса при малых скоростях флюида, является актуальной во многих приложениях гидромеханики, решению которых посвящено огромное количество исследований, см. монографии и обзорные работы \cite{1} - \cite{7}. Мы рассмотрим её в классической бездивергентной постановке для трёхмерного случая:

\begin{equation} \label{non-divergent statement}
    \frac{\partial \rho u}{\partial t}-\mu\Delta \overline{u} + \nabla p = f(\overrightarrow{r}), \quad div\ \overline{u} = 0, \quad \overline{r} \in \Omega, \quad \overline{u} = (u_1, u_2, u_3),
\end{equation}
где $\rho, \mu$ -- предполагаемые для простоты постоянными плотность и динамическая вязкость, а для вектора скорости $\overline{u}$ на границе $\Gamma$ и для давления в области $\Omega$ выполняются условия

\begin{equation} \label{boundary conditions}
    \overline{u}|_\Gamma = \overline{g}(\overline{r}), \quad \int\limits_{\Omega} p d \overline{r} = 0,
\end{equation}
обеспечивающие существование единственного решения в пространствах $\overline{u} \in \left( H^1_0 (\Omega) \right)^3, \ p \in L_0^2(\Omega)$ \cite{1}. Отметим, в частности, что на границе жидкости с твердой неподвижной стенкой задается условие прилипания $\overline{g}(\overline{r})=0$.

Для простоты расчётную область предполагаем параллелепипедом. Исходную задачу \eqref{non-divergent statement} - \eqref{boundary conditions} аппроксимируем на кубической смещённой сетке с шагом $h$, в которой компоненты вектора скорости определены в центрах перпендикулярных им сеточных граней, а давление -- в центрах сеточных ячеек ($x_i=ih, y_j=jh, z_k=kh$):
\begin{equation*}
    \begin{gathered}
        u_x = u_1(x_i,\ y_{j+1/2},\ z_{k+1/2}), \quad u_y = u_2(x_{i+1/2},\ y_j,\ z_{k+1/2}),
        \\
        u_z = u_3(x_{i+1/2},\ y_{j+1/2},\ z_k), \quad p = p(x_{i+1/2},\ y_{j+1/2},\ z_{k+1/2}),
        \\
        i = 0, 1, \ldots, N_1; \quad j = 0, 1, \ldots, N_2; \quad k = 0, 1, \ldots, N_3.
    \end{gathered}
\end{equation*}

После применения ``стандартных'' неявных конечно-разностных (или конечно-объёмных, или простейших конечно-элементных \cite{8}) аппроксимаций задачи \eqref{non-divergent statement} - \eqref{boundary conditions} мы получаем симметричную систему линейных алгебраических уравнений (СЛАУ) седлового типа, которая в матричном представлении имеет вид
\begin{equation} \label{saddle type system}
    Au = 
    \begin{bmatrix}
		\Delta & \nabla \\ 
		\nabla^{\top} & 0
	\end{bmatrix}
    \begin{bmatrix}
		\overline{u} \\ 
		p
	\end{bmatrix}
    =
    \begin{bmatrix}
		\overline{f} \\ 
		0
	\end{bmatrix}
    = f, \quad \overline{u} = (u_1, u_2, u_3), \quad \overline{f} = (f_1, f_2, f_3).
\end{equation}

Здесь и далее символы $\Delta$ и $\nabla$ обозначают матрицы, получаемые в результате сеточных аппроксимаций соответствующих дифференциальных операторов (для Лапласа -- с учётом граничных условий Дирихле \eqref{boundary conditions}, фактически этот оператор включает учёт коэффициента $\mu$ и аппроксимацию производной по времени, которая усиливает его диагональное преобладание). Отметим, что в общем случае сеточные функции, или подвекторы, $u_1, u_2, u_3, p$ имеют различные размерности, т. е.
\begin{equation*}
    \begin{gathered}
        u_1 \in \mathbb{R}^{N_x}, \quad
        u_2 \in \mathbb{R}^{N_y}, \quad
        u_3 \in \mathbb{R}^{N_z}, \quad
        p \in \mathbb{R}^{N_p},
        \\
        N_x = (N_1 - 1) N_2 N_3, \quad
        N_y = N_1 (N_2 - 1) N_3, \quad
        N_z = N_1 N_2 (N_3 - 1),
        \\
        N_p = N_1 N_2 N_3,
    \end{gathered}
\end{equation*}
общий порядок матрицы $A$ при этом равен $N = N_1 + N_2 + N_3 + N_p \equiv N_u + N_p$. Заметим также, что вектор 
\begin{equation} \label{eigenvector}
    u_0 = \left(0 \atop e_p\right), \quad e_p = \{1\} \in \mathbb{R}^{N_p}
\end{equation}
есть собственный вектор матрицы $A$, соответствующий её нулевому собственному числу, причём оно не является кратным, см. \cite{2} - \cite{6}. Это означает, что для согласования СЛАУ \eqref{saddle type system} необходимо условие ортогональности $(f,\ e) = 0, \ e = \{1\} \in \mathbb{R}^N$, a если вектор $u = \left(\overline{u} \atop p\right)$ есть решение СЛАУ \eqref{saddle type system}, и его подвектор $p$ ортогонален соответствующему подпространству ядра матрицы $A$, т. е. $(p,\ e_p) = 0$, то такое решение единственно.

При выбираемой нами последовательно-естественной нумерации разных типов точек определения подвекторых компонент матрицы $\Delta, \nabla, \nabla^{\top}$ примут блочную форму, и матричное представление СЛАУ \eqref{saddle type system} записывается как
\begin{equation} \label{block form}
    Aw = 
    \begin{bmatrix}
		\Delta_1 & 0 & 0 & \nabla_1 \\ 
		  0 & \Delta_2 & 0 & \nabla_2 \\
        0 & 0 & \Delta_3 & \nabla_3 \\
        \nabla^{\top}_1 & \nabla^{\top}_2 & \nabla^{\top}_3 & 0\\
	\end{bmatrix}
    \begin{bmatrix}
		u_1 \\
        u_2 \\
        u_3 \\
		p
	\end{bmatrix}
    =
    \begin{bmatrix}
		f_1 \\
        f_2 \\
        f_3 \\
		0
	\end{bmatrix}.
\end{equation}

В рассматриваемых случаях матрицы $\Delta_l, \ l = 1, 2, 3$, являются семидиагональными, симметричными и положительно определёнными (с.~п.~о), а $\nabla_l$ -- треугольными двухдиагональными. Современные практические требования к точности численных методов обуславливают высокие порядки алгебраических систем вида \eqref{block form}, до нескольких миллиардов $(N \approx 10^9)$, что делает востребованным построение высокопроизводительных итерационных алгоритмов, в том числе с их распараллеливанием на многопроцессорных вычислительных системах (МВС) с распределённой и иерархической памятью.

Среди большого числа современных публикаций на данную тему (см., например, \cite{9} - \cite{13}) общая тенденция исследований направлена на использование предобусловленных методов в подпространствах Крылова \cite{7}. Длительное время популярностью пользуются алгоритмы типа SIMPLE (semi-implicit method for pressure-linked equations) и варианты метода Узавы \cite{1}. Отметим, что для решения сеточных уравнений Стокса также эффективно применяются многосеточные подходы, см. \cite{Oosterlee, Notay}/

В данной работе для решения СЛАУ мы исследуем метод сопряжённых градиентов с применением дефляции \cite{14}, в целях обеспечения свойства ортогональности давления \eqref{boundary conditions} на каждой итерации. Рассматривается также подход с регуляризацией седловой матрицы, аналогичный разработанному в \cite{15} при решении вырожденной задачи Неймана. В качестве предобуславливателя используется матрица блочной неполной факторизации, в которой операторы Лапласа и дополнения Шура аппроксимируются или с помощью приближения редкой обратной матрицы (SPAI, \cite{16}), или ``чисто алгебраическим'' многосеточным алгоритмом неполной факторизации, предложенным в \cite{17}.

Настоящая работа построена следующим образом. В п. 2 излагаются предобусловленные методы сопряжённых градиентов, в том числе с ускорением при помощи регуляризации СЛАУ и применении дефляции к итерационному процессу. П. 3 посвящён описанию предобуславливающих матриц с помощью подходов SPAI и AMG, а в последнем пункте обсуждаются предварительные результаты численных экспериментов с распараллеливанием алгоритмов на различном числе процессов.

\section{Блочный предобусловленный метод сопряжённых градиентов для решения задачи Стокса}

Определим алгебраическую регуляризованную задачу Стокса как СЛАУ, являющуюся обобщением \eqref{saddle type system}:
\begin{equation} \label{generalized system}
    A_{\gamma} u = 
    \begin{bmatrix}
		A_{1,1} & \nabla \\ 
		\nabla^{\top} & \gamma A_{2,2}
	\end{bmatrix}
    \begin{bmatrix}
		\overline{u} \\ 
		p
	\end{bmatrix}
    =
    \begin{bmatrix}
		\overline{f} \\ 
		0
	\end{bmatrix},
\end{equation}
где $A_{1,1}=\Delta+\nabla K\nabla^{\top}$, $A_{2,2}=ME$,  $E = \{1\} \in \matbb{R}^{N_p, N_p}$ есть матрица с единичными компонентами, $K \in \mathcal{R}^{N_u, N_u}$ и $M\in \mathcal{R}^{N_p,N_p}$ -- некоторые с.п.о. - матрицы, а $\gamma$ -- вещественный параметр регуляризации. В силу свойств матрицы градиента $\quad \nabla e_p = 0$ можно показать следующие два результата.

\begin{theorem} \label{1}
    Матрица $A_{\gamma}$ при достаточно большом значении $\gamma \geqslant \gamma_0 > 0$ имеет положительные собственные числа, за исключением одного (нулевого), соответствующего собственному вектору $u_0$ из \eqref{eigenvector}. СЛАУ \eqref{generalized system} является согласованной при условии $(f,u_0)=0$.
\end{theorem}

\begin{theorem}
    Алгебраическая система \eqref{generalized system} в условиях теоремы \ref{1} имеет единственное решение, которое обладает условием ортогональности для давления $(p, e_p) = 0$, и это решение не зависит от $\gamma$.
\end{theorem}

\begin{remark}
    Вместо условия ортогональности давления $p$ ядру матрицы $A_{\gamma}$, однозначность решения СЛАУ \eqref{generalized system} можно обеспечить, фиксируя значение (например, нулевое) давления $p$ в одном из узлов сетки. Однако можно ожидать, по аналогии с вырожденной задачей Неймана для уравнения Пуассона, что при этом значительно вырастет число обусловленности матрицы $A_{\gamma}$. Более строгий анализ данного вопроса требует специальных исследований. 
\end{remark}

\begin{remark}
    Наличие регуляризующего параметра $\gamma$ в матрице $A_{\gamma}$ формально даёт возможность минимизировать её число обусловленности и количества итераций при решении СЛАУ \eqref{generalized system}. Теоретически этот вопрос пока не исследован, а приближённая оптимизация $\gamma$ проводится экспериментально.
\end{remark}

Метод сопряжённых градиентов с предобуславливающей матрицей $B$ для решения СЛАУ \eqref{block form}, запишем в следующем виде \cite{14}:
\begin{equation} \label{preconditioned CG}
    \begin{gathered}
        r^0 = f - A_{\gamma}u^0, \quad q^0 = B^{-1} r^0, \quad \alpha_n = \sigma_n / \rho_n,
        \\
        u^{n+1} = u^n + \alpha_n q^n, \quad r^{n+1} = r^n - \alpha_n A_{\gamma}q^n, \quad q^{n+1} = B^{-1} r^{n+1} + \beta_n q^n,
        \\
        \sigma_n = (r^n, q^n) = (B^{-1} r^n, r^n), \quad \rho_n = (A_{\gamma}q^n, q^n), \quad \beta_n = \sigma_{n+1} / \sigma_n,
    \end{gathered}
\end{equation}
где вектор начального приближения предполагается ортогональным ядру матрицы $A_{\gamma}$, т.~е. $(u^0, u_0) = 0$,  для чего достаточно обеспечить условие $\sum\limits_{l = 1}^{N_p} p_l^0 = 0$. При этом для достижения условия сходимости
\begin{equation*}
    ||r^n|| \leqslant \varepsilon ||f||, \quad \varepsilon \ll 1
\end{equation*}
достаточное число итераций оценивается величиной
\begin{equation*}
    n(\varepsilon) \leqslant \frac{1}{2} |ln \frac{\varepsilon}{2}| (cond (B^{-1} A_{\gamma})^{1/2} + 1).
\end{equation*}

Построение предобуславливающей матрицы $B$ определим на основе блочно-факторизованного представления матрицы $\overline{A}_{\gamma}$, приближённо равной матрице $A_{\gamma}$ из СЛАУ \eqref{generalized system}:
\begin{equation} \label{factorization A_gamma}
    B = \overline{A}_{\gamma} =
    \begin{bmatrix}
        \overline{A}_{1,1}  & \nabla \\
        \nabla^{\top} & \gamma A_{2,2}
    \end{bmatrix}=
    \begin{bmatrix}
        I & 0 \\
        \nabla^{\top}\overline{A}^{-1}_{1,1} & I
    \end{bmatrix}
    \begin{bmatrix}
        \overline{A}_{1,1} & 0 \\
        0 & \overline{S}_{\gamma}
    \end{bmatrix}
    \begin{bmatrix}
        I & \overline{A}^{-1}_{1,1}\nabla \\
        0 & I
    \end{bmatrix},
\end{equation}
где $\overline{A}_{1,1}$ -- некоторая аппроксимация матрицы $A_{1,1}$, $\overline{S}_{\gamma} = \gamma A_{2,2} - \nabla^{\top} \overline{A}^{-1}_{1,1} \nabla$ есть дополнение Шура для матрицы $\overline{A}_{\gamma}$. Предобуславливатель конструируем с помощью одной из двух аппроксимаций матрицы $\overline{A}^{-1}_{1,1}$. Первый способ заключается в приближённом ленточном представлении обратной матрицы SPAI (см. \cite{16} и приведённую там литературу), а второй -- в применении алгебраического многосеточного метода (AMG, см. обзор в \cite{17}). В обоих случаях предобуславливатель строится путём замены точной матрицы $A_{1,1}$ из \eqref{generalized system} на приближённую. На каждой итерации требуется решать вспомогательную СЛАУ вида $Bq = r$, что экономично реализуется с помощью прямого и обратного хода факторизации:
\begin{equation} \label{forward and backward factorization}
    Bq =
    \begin{bmatrix}
        \overline{A}_{1,1} & 0 \\
        \nabla^{\top} & \overline{S}_{\gamma}
    \end{bmatrix}
    \begin{bmatrix}
        v_1 \\
        v_2
    \end{bmatrix}
    =
    \begin{bmatrix}
        r_1 \\
        r_2
    \end{bmatrix}, \quad
    \begin{bmatrix}
        I & \overline{A}_{1,1}^{-1} \nabla \\
        0 & I
    \end{bmatrix}
    \begin{bmatrix}
        q_1 \\
        q_2
    \end{bmatrix}
    =
    \begin{bmatrix}
        v_1 \\
        v_2
    \end{bmatrix}.
\end{equation}
В компонентной форме эти соотношения записываются следующим образом:
\begin{equation} \label{component form}
    \begin{gathered}
        v_1 = \overline{A}^{-1}_{1,1} r, \quad \overline{S}_{\gamma} v_2 = (\gamma E - \nabla^{\top} \overline{A}^{-1}_{1,1} \nabla) v_2 = r_2,
        \\
        q_2 = v_2, \quad q_1 = v_1 - \overline{A}^{-1}_{1,1} \nabla q_2,
    \end{gathered}
\end{equation}
где $\overline{A}_{1,1}$ означает какой-то из рассматриваемых предобуславливателей к матрице $A_{1,1}$, на которых остановимся в следующем разделе. Отметим, что в соотношениях \eqref{component form} при вычислениях $v_1$ и $q_1$ необходимо дважды умножать вектор на ленточную матрицу $\overline{A^{-1}_{1,1}}$, однако нахождение $v_2$ требует решения СЛАУ с плотной матрицей $\overline{S}_{\gamma}$. Здесь, тем не менее, можно заметить два обнадёживающих фактора. Во-первых, умножение на неё делается очень экономично, поскольку $Ev = (v, e_p) e_p$. Во-вторых, по аналогии с работой \cite{15}, следует ожидать хорошей обусловленности (не зависящей от $h$) матрицы $\overline{S}_{\gamma}$, за счёт выбора параметра $\gamma$ и матриц $K$, $M$. В таком случае целесообразно было бы применить эффективный алгоритм чебышёвского ускорения, не требующий плохо распараллеливаемых векторных скалярных произведений, но при этом требуется определять границы спектра матрицы $\overline{S}_{\gamma}$.

Отметим, что для обеспечения условия ортогональности давления $p$ к ядру матрицы $A$ из \eqref{saddle type system} на всех итерациях (которые выполняются с неизбежными машинными округлениями) предобусловленный метод сопряжённых градиентов следует применять с дефляцией, заключающейся в добавлении к базисам подпространств Крылова ``дефляционного'' вектора $v$, ортогональному ядру матрицы $A_{\gamma}$, т.е. $(v, u_0) = 0$ что приводит к следующим формулам (при этом предварительно осуществляется ортогонализация начального приближения и направляющего вектора):

\begin{equation*}
    \begin{gathered}
        u^0 = u^{-1} + vc = u^{-1} + Hr^{-1},\quad H = v \hat{A}_{\gamma}^{-1} v^{\top},
        \\
        c = \hat{A}^{-1}_{\gamma} v^{\top} r^{-1},\quad \hat{A}_{\gamma} = v^{\top} Av,\quad r^{-1} = f - A_{\gamma}u^{-1},
    \end{gathered}
\end{equation*}
где $u^{-1}$ -- произвольный вектор, удовлетворяющий условию $(u^{-1}, v) = 0$.

\begin{equation} \label{CG with deflation}
    \begin{gathered}
       r^0 = f - A_{\gamma}u^0, \quad q^0 = QB^{-1}r^0, \quad Q = I - HA_{\gamma}, n = 0, 1, \ldots :
       \\
       u^{n+1} = u^n + \alpha_n q^n, \quad r^{n+1} = r^n - \alpha_n A_{\gamma}q^n, \quad \alpha_n = \sigma_n / \rho_n,
       \\
       q^{n+1} = QB^{-1}r^{n+1} + \beta_n q^n, \quad \beta_n = \sigma_{n+1} / \sigma_n,
       \\
       \sigma_n = (B^{-1} r^n, r^n), \quad \rho_n = (Aq^n, q^n).
    \end{gathered}
\end{equation}
В этих соотношениях матрица $\hat{A}$ введена для общности (если дефляционных векторов несколько, то $v = (v_1, \ldots, v_m) \in \mathbb{R}^{N, m}$ есть дефляционная матрица, а в нашем случае $m = 1$ и $\hat{A}$ есть фактически скаляр).

Для решения блочного СЛАУ второго порядка на основе точной LDU-факторизации \eqref{factorization A_gamma} можно также предложить двустороннее экономичное предобуславливание следующего типа:
\begin{equation} \label{two-sided}
    \begin{gathered}
    A_{\gamma} = LDU, \quad L =
    \begin{bmatrix}
        I & 0\\
        A_{2,1}A^{-1}_{1,1} & I
    \end{bmatrix}, \quad
    U =
    \begin{bmatrix}
        I & A^{-1}_{1,1}A_{1,2}\\
        0 & I
    \end{bmatrix}, \quad
    D =
    \begin{bmatrix}
        A_{1,1} & 0\\
        0 & S
    \end{bmatrix},
    \\
    S_{\gamma} = \gamma A_{2,2} - A_{2,1} A_{1,1}^{-1}A_{1,2},
    L^{-1} =
    \begin{bmatrix}
        I & 0\\
        -A_{2,1}A^{-1}_{1,1} & I
    \end{bmatrix},
    U^{-1} =
    \begin{bmatrix}
        I & -A^{-1}_{1,1}A_{1,2}\\
        0 & I
    \end{bmatrix}, \\
    
    \overline{A}\overline{u} = L^{-1}A U^{-1}\overline{u} = L^{-1}f = \overline{f}, Uu, u = U^{-1} \overline{u},
    \\
    D\overline{u} =
    \begin{bmatrix}
        A_{1,1} & 0\\
        0 & S_{\gamma}\\
    \end{bmatrix}
    \begin{bmatrix}
        \overline{u}_{1}\\
        \overline{u}_{2}\\
    \end{bmatrix}
    =
    \begin{bmatrix}
        \overline{f}_1\\
        \overline{f}_2\\
    \end{bmatrix}.
    \end{gathered}
\end{equation}

При итерационном решении СЛАУ \eqref{two-sided} с матрицей $D$ на каждой итерации требуется однократное решение независимых подсистем с матрицами $A_{1,1}$ и $S_{\gamma}$. Векторы $\overline{f}$, $u$ вычисляются по одному разу, перед и после итераций соответственно.

В данном алгоритме требуется многократное решение алгебраических систем с матрицей $A_{1,1}$ в частности, при реализации двухуровневого итерационного процесса для СЛАУ с матрицей $S_{\gamma}$. В этих случаях наиболее целесообразно применять экономичный многосеточный метод AMG. Что касается внешних итераций, то здесь, вообще говоря, необходимо использовать предобусловленные алгоритмы с длинными рекурсиями типа FCG или GMRES. Предобуславливающую матрицу целесообразно брать на основе неполной факторизации или SSOR метода в применении к приближённому дополнению Шура
\begin{equation}\label{eq_diag_compens}
    \overline{S}_{\gamma} = \gamma M-A_{2,1}\overline{A}_{1,1}^{-1}A_{1,2},
\end{equation}
где $\overline{A}_{1,1}^{-1}$ -- ленточная аппроксимация обратной матрицы $A_{1,1}^{-1}$ (SPAI), в качестве $M$ проще взять единичную матрицу, а параметр $\gamma$ подбирается экспериментально. При этом матрица $\overline{S}=D_s + L_s + U_s$ разбивается на диагональную, нижнюю и верхнюю треугольные части, а предобуславливатель определяется формулой $B_s = (G_s+L_s)G_s^{-1}(G_s+U_s)$, подробнее о $G_s$ -- см. \cite{7}.

Для решения СЛАУ \eqref{generalized system} с матрицей $A_{\gamma}$ мы рассмотрим ещё построение предобуславливателя на основе классического блочного метода неполной факторизации с диагональной компенсацией, базирующегося на построении предобуславливателя следующего вида:
\begin{equation}\label{13}
    B_{\gamma} =
    \begin{bmatrix}
        A_{1,1} & 0 \\
        A_{2,1} & \gamma A_{2,2}
    \end{bmatrix}
    \begin{bmatrix}
        A_{1,1}^{-1} & 0 \\
        0 & G_{\gamma}^{-1} \\
    \end{bmatrix}
    \begin{bmatrix}
        A_{1,1} & A_{1,2} \\
        0 & G_{\gamma} \\
    \end{bmatrix},
\end{equation}
где матрица $G_{\gamma}$ есть фактически некоторая аппроксимация дополнения Шура, представляемая формулой
\begin{equation}\label{14}
    \begin{gathered}
    G_{\gamma} = \gamma A_{2,2} - A_{2,1}\overline{A}^{-1}_{1,1}A_{1,2} - \theta R, \quad \theta \in (0,1) \\
    Re_2 = (A_{2,1}A_{1,1}^{-1}A_{1,2} - A_{2,1}\overline{A}_{1,1}^{-1}A_{1,2})e_2 = \\
    = A_{2,1}(A^{-1}_{1,1} - \overline{A}^{-1}_{1,1})A_{1,2}e_2, \quad e_2=\{1\}\in \mathcal{R}^{N_p}.
    \end{gathered}
\end{equation}

Здесь $G \in (0,1)$ -- компенсирующий параметр, а $R$ -- диагональная матрица, а применение соотношений \eqref{13}, \eqref{14} (без регуляризации, т.е. $\gamma = 0$) обеспечивает при $G = 1$ условие полной компенсации, или согласования строчных сумм
\begin{equation}\label{15}
   B_{\gamma}e = A_{\gamma}e, \quad e = \{1\}\in \mathcal{R}^N,
\end{equation}
характеризующее в определённом смысле близость исходной и предобуславливающей матриц.

\section{Построение предобуславливателей \newline для оператора Лапласа}

Самый традиционный из универсальных способов аппроксимации симметричной матрицы заключается в приближённом разложении её на треугольные множители.
Если матрица $A_{1,1}$ представима в виде $A_{1,1} = D + L + U$, где $D, L, U$ -- в общем случае её блочные диагональная, нижняя и верхняя треугольные части, то предобуславливатель можно записать в виде
\begin{equation} \label{preconditioner}
    B_1 = (G + L)\ G^{-1} (G + U) = G + L + U + L\ G^{-1} U.
\end{equation}

Поскольку в формуле \eqref{preconditioner} представление $G = D - L\ G^{-1} U$, обеспечивающее точную факторизацию $(B_1 = A_{1,1})$, слишком ресурсоёмко, оно заменяется одной из двух возможных диагональных ``аппроксимаций'', реализующих неполную (приближённую) факторизацию. Первая основана на принципе компенсации, или согласования строчных сумм:
\begin{equation*}\label{theta, R}
    G = D - \overline{L\ G^{-1}U} - \theta R,\quad Re = (L\ G^{-1} U - \overline{L\ G^{-1} U}) e, \quad \theta \in [0,1],
\end{equation*}
где $\theta$ -- компенсирующий параметр, $e$ -- пробный вектор, обычно выбираемый из единичных компонент, а черта над матрицей означает взятие её диагональной части, т. е. $\overline{M} = diag\{m_{i,i}\}$. Здесь $R$ есть диагональная матрица, обеспечивающая выполнение условия ``полной компенсации'' $B_1 e = A_{1,1} e$ при $\theta = 1$, или согласование строчных сумм \cite{7}. Данный алгоритм будем называть CIF (Compensated Incomplete Factorization).

Второй подход заключается в использовании последовательной верхней релаксации с параметром релаксации $\omega$, который обозначается как SSOR (Symmetric Successive Over Relaxation):
\begin{equation} \label{SSOR}
    G = \omega^{-1}D, \quad \omega \in (0,2).
\end{equation}
По поводу выбора итерационных параметров $\theta$ или $\omega$ см. \cite{18}.

Для построения аппроксимации непосредственно обратной матрицы $\overline{A_{1,1}^{-1}}$ существует оригинальный подход SPAI, в котором приближение определяется в виде разрежённой матрицы с тем же портретом, что и у исходной матрицы $\Delta$. Конкретные алгоритмы при этом строятся на основе минимизации нормы Фробениуса для матричной невязки
\begin{equation}
    ||I - A_{1,1} \overline{A}_{1,1}^{-1}||_F = \sum_{j} (e_j - A_{1,1}m_j)^2,
\end{equation}
где $m_j$ и $e_j$ -- столбцы матриц $\overline{A}_{1,1}^{-1}$ и $I$ с номером $j$. Здесь ненулевые неизвестные $m_j$ находятся из независимого решения систем $Am_j = e_j$ методом наименьших квадратов \cite{16}, а получаемый при этом алгоритм естественным образом распараллеливается. Отметим, что поскольку в алгоритме приближённая обратная может получиться несимметричной даже для симметричной матрицы $A_{1,1}$, рекомендуется применять следующую процедуру симметризации:
\begin{equation}
    (\overline{A}^{-1}_{1,1})_{s} = \frac{1}{2}(\overline{A}^{-1}_{1,1} + (\overline{A}^{-1}_{1,1})^{\top}).
\end{equation}

В качестве дальнейшей модификации метода SPAI можно предложить использованный в \eqref{SSOR} алгоритм диагональной компенсации, в котором приближенная обратная матрица определяется как
\begin{equation} \label{17}
    (\overline{A}^{-1})_c = (\overline{A}^{-1})_s + \theta R, \quad Re_1 = A^{-1}e_1 - (\overline{A}^{-1})_s e_1, \quad e_1 \in R_1^{N_1}
\end{equation}
где $R$ -- диагональная матрица, обеспечивающая согласование строчных сумм $(\overline{A}^{-1})_c e_1 = A^{-1}e_1$.

Наиболее перспективным представляется подход, основанный на ``чисто алгебраической'' версии многосеточного метода AMG, предложенной в \cite{17} и основанного на неполной факторизации с диагональной компенсацией предобуславливателя блочно-треугольного вида, формируемого из иерархического определения матричного графа решаемой на вложенных сетках СЛАУ.

Данный подход базируется на построении рекурсивной сеточной и соответствующей алгебраической структур. Как обычно, многосеточные алгоритмы трактуются как рекурсивное применение двухсеточных, определение которых проведём на простом, но достаточно общем, примере следующим образом.

Пусть исходная СЛАУ \eqref{saddle type system} аппроксимируется на кубической сетке $\Omega^h$ с шагом $h$, содержащей последовательность вложенных сеток $\Omega_{l}^h, \ l = 1, \ldots, m$, каждая с числом узлов $N_l$, таких, что $\Omega^h = \Omega_1^h \supset \Omega_2^h \supset \ldots \supset \Omega_m^h$, причём граница $\Gamma$ проходит по линиям (или плоскостям) самой редкой сетки $\Omega_m^h$, а каждая более густая сетка 
$\Omega_{l-1}^h, \ l = 2, \ldots, m$, получается  из $\Omega_l$  делением шага предыдущей сетки пополам. Рассматриваемый далее алгоритм для решения систем с матрицей $A_{1,1}$ фактически основывается на применении классического алгебраического метода неполной факторизации при специальной перенумерации сеточных узлов и соответствующих векторных компонент.

Каждое из соответствующих множеств сеточных узлов и определённых в них векторных компонент разобьём на два подмножества, в результате чего получаем следующие рекурсивные представления исходной сетки и векторов:
\begin{equation}
    \begin{gathered}
        \Omega^h_{1} = \overline{\Omega}_1^h \cup \Omega^h_2 = \overline{\Omega}_1^h \cup \overline{\Omega}^h_2 \ldots \cup \overline{\Omega}_{m-1}^h \cup  \Omega^h_m,
        \\
        u = u^{(1)} = \left( (\overline{u}^{(1)})^\top, ({u}^{(2)})^\top \right)^\top =
        \\
        \left((\overline{u}^{(1)})^\top, (\overline{u}^{(2)})^\top, \ldots, \overline{u}^{(m-1)})^\top, u^{(m)})^\top \right)^\top.
    \end{gathered}
\end{equation}
Здесь каждая блочная векторная компонента $u^{(l)}$ уровня $l = 1, \ldots, m-1$ разбивается на две части $(\overline{u}^{(l)}$ и $u^{(l+1)})$, вторая из которых представляет следующий уровень и тоже делится надвое так, что вторая часть $(l+1)$-го уровня совпадает с целым $(l+2)$-ым уровнем, и т.д. Нумерация сеточных узлов и соответствующих векторных компонент осуществляется последовательно: сначала для $\overline{\Omega}^h$ и $\overline{u}^{(1)}$, затем для $\overline{\Omega}_2^h, \ \overline{u}^{(2)}$, и т. п. В этом случае для двухсеточного варианта матрица $A_{1,1}$ приводится к блочному виду второго порядка. В трёхсеточном методе правый нижний блок матрицы аналогично, в свою очередь, сводится к блочному виду второго порядка, а в многосеточном подходе эта процедура повторяется многократно.

При более детальной классификации узлов, которые мы будем использовать, матрица $A_{1,1}$ приводится к блочному четвёртому порядку. Для этого при $l = 1$ узлы $\Omega^h_l$ разобьём на четыре подмножества: $\Omega^h_l = \Omega_l^1 \cup \Omega_l^2 \cup \Omega_l^3 \cup \Omega_l^4, \quad \Omega_l^4 = \Omega^h_{l+1}$, которые относятся последовательно к центрам объёмов и граней, к серединам ребер и узлам вложенной редкой сетки $\Omega^h_{l+1}$. Для сокращения дальнейшего изложения рекурсивного подхода мы здесь, вплоть до формулы \eqref{preconditioner in factored form}, при описании исходной матрицы $A_{1,1} = A_{1,1}^{(1)}$ используем индекс $l$, подразумевая его значение $l = 1$. Обозначая через $\overline{u}_1^{(l)}, \overline{u}_2^{(l)}, \overline{u}_3^{(l)}, \overline{u}_4^{(l)}$ и $\overline{f}_1^{(l)}, \overline{f}_2^{(l)}, \overline{f}_3^{(l)}, \overline{f}_4^{(l)}$ ассоциированные подвекторы (с размерностями $N_1^{(l)}, N_2^{(l)}, N_3^{(l)}, N_4^{(l)})$ в $l$-й СЛАУ, запишем её в соответствующем блочном виде:
\begin{equation} \label{Block matrix of order 4}
    A_{1,1}^{(l)} u^{(l)} =
    \begin{bmatrix}
        D_{1,1}^{(l)} & A_{1,2}^{(l)} & 0   & 0 \\
        A_{2,1}^{(l)} & D_{2,2}^{(l)} & A_{2,3}^{(l)} & 0 \\
        0 & A_{3,2}^{(l)} & D_{3,3}^{(l)} & A_{3,4}^{(l)} \\
        0& 0 & A_{4,3}^{(l)} & D_{4,4}^{(l)}
    \end{bmatrix}
    \begin{bmatrix}
        \overline{u}_1^{(l)} \\
        \overline{u}_2^{(l)} \\
        \overline{u}_3^{(l)} \\
        \overline{u}_4^{(l)}
    \end{bmatrix} =
    \begin{bmatrix}
        \overline{f}_1^{(l)} \\
        \overline{f}_2^{(l)} \\
        \overline{f}_3^{(l)} \\
        \overline{f}_4^{(l)}
    \end{bmatrix}. 
\end{equation} 
Здесь можно обозначить $\overline{u}_4^{(l)} = u^{(l+1)}, \ \overline{f}_4^{(l)} = \overline{f}^{(l+1l)}$, так как соответствующие подвекторы представляют структуры следующего сеточного уровня. Заметим, что блоки $D^{l}_{1,1}, ..., D_{4,4}^{(l)}$ в \eqref{Block matrix of order 4} (напомним, что пока $l = 1$) являются диагональными матрицами. 

Отметим, что если рассмотреть уравнения узлового типа на вложенных трёхмерных сетках, составленных из тетраэдров, то множество узлов густой сетки $\Omega_1^h$ представляется состоящим только из двух подмножеств: $\overline{\Omega}_1^h$ -- подмножество середин рёбер редкой сетки и само подмножество узлов $\Omega_2^h$. Поэтому блочная структура СЛАУ в данном случае имеет второй порядок. В целом рассмотренные матричные представления двухсеточных, или двухуровневых, СЛАУ можно определить как $k$-блочные структуры, где $k = 2, 3, 4$ обозначает блочный порядок алгебраической системы \eqref{Block matrix of order 4} соответственно.

Для матрицы $A^{(l)}$ предобуславливатель определяем в факторизованной форме:
\begin{equation} \label{preconditioner in factored form}
    \begin{gathered}
        B^{(l)} =
        \begin{bmatrix}
            G_1^{(l)} & 0 & 0 & 0 \\
            A_{2, 1}^{(l)} & G_2^{(l)} & 0 & 0\\
            0 & A_{3, 2}^{(l)} & G_{3}^{(l)}& 0 \\
            0 & 0 & A_{4, 3}^{(l)} & G_4^{(l)}
        \end{bmatrix}
        \left( G^{(l)} \right)^{-1}
        \begin{bmatrix}
            G_1^{(l)} & A_{1, 2}^{(l)} & 0 & 0 \\
            0 & G_2^{(l)} & A_{2, 3}^{(l)} & 0 \\
            0 & 0 & G_3^{(l)} & A_{3, 4}^{(l)} \\
            0 & 0 & 0 & G_4^{(l)}
        \end{bmatrix},
        \\
        G_1^{(l)} = D^{(l)}_{1,1}, \quad G_2^{(l)} = D_{2, 2}^{(l)} -
        (A_{2, 1}^{(l)} (G_1^{(l)})^{-1} A_{1, 2}^{(l)})_1 - \theta S_2^{(l)},
        \\
        S_2^{(l)} e_2 = \left[A_{2, 1}^{(l)} (G_1^{(l)})^{-1} A_{1,2}^{(l)} - (A_{2,1}^{(l)} (G_1^{(l)})^{-1} A_{1,2}^{(l)})_1\right] e_2,
        \\
        G_3^{(l)} = D_{3, 3}^{(l)} - \left(A_{3, 2}^{(l)} (G_2^{(l)})^{-1} A_{2, 3}^{(l)}\right)_1 - \theta S_3^{(l)},
        \\
        S_3^{(l)} e_3 = \left[A_{3, 2}^{(l)} (G_3^{(l)})^{-1} A_{2, 3}^{(l)} - \left(A_{3, 2}^{(l)} (G_2^{(l)})^{-1} A_{3, 2}^{(l)}\right)_1 \right] e_3,
        \\
        G_4^{(l)} = D_{4, 4}^{(l)} - A_{4, 3}^{(l)} (G_3^{(l)})^{-1} A_{3, 4}^{(l)}.
    \end{gathered}
\end{equation}
Здесь  $A_{k, k}^{(l)}, \ k = 1, \ldots, 4$ и $G_1^{(l)}, G_2^{(l)}, G_3^{(l)}$ -- диагональные матрицы, а  $G_4^{(l)}$~-- семидиагональная с тем же матричным портретом, что и $A^{(l)}_{1,1}$, и именно её мы будем считать матрицей СЛАУ для сетки следующего уровня. Таким образом, при $l = 1$ и $m > 1$ в \eqref{preconditioner in factored form} следует положить $G_4^{(1)} = A^{(2)}_{1,1}$, определить для данной матрицы блочно-трёхдиагональное представление \eqref{Block matrix of order 4} и предобуславливатель $B^{(2)}$. Далее для $l = 2,..,m$ рекурсивное формирование $G^{(l)}_4$ продолжается.

Запись вида $(C)_1$ в \eqref{preconditioner in factored form} обозначает диагональную часть матрицы $C, \ \theta \in [0, 1]$ -- параметр компенсации, а пробные векторы (обычно их компоненты берут равными единице) $e_2^{(l)}, e_3^{(l)}$ имеют размерности $N_2^{(l)}, N_3^{(l)}$. Для конкретизации предобуславливающей матрицы $B^{(l)}$ нам ещё необходимо представить матрицу $G_4^{(l)}$ в факторизованной форме, и при этом будем рассматривать два варианта. Если $l + 1 = \overline{m}$ есть номер последней применяемой сетки, то используется точное треугольное разложение $G_4^{(l)} = A^{(l+1)} = L^{(l+1)} U^{(l+1)}$. В противном случае матрица $G_4^{(l)}$ в \eqref{preconditioner in factored form} приближённо заменяется факторизованным предобуславливателем $B^{(l+1)}$, определяемым в соответствии с этой же формулой. При использовании факторизованного предобуславливателя $B$ из \eqref{17} в методе сопряжённых градиентов \eqref{preconditioned CG} решение вспомогательной СЛАУ вида $Bq = r$ легко реализуется по следующим формулам (точнее говоря, ниже представлен псевдокод):
\begin{equation}
    \begin{gathered}
        G_1^{(l)}q_1^{(1)} = r_1^{(1)}, \quad G_2^{l}q_2^{(l)} = r_2^{(l)}-A_{2,1}^{(l)}q_1^{(1)}, \quad G_3^{(l)}q_3^{(l)} = r_3^{(l)}-A^{(l)}_{3,2}q_2^{(l)},\\
        G_4^{(l)}q_4^{(l)} = r_4^{(l)}-A_{4,3}^{(l)}q^{(l)}_3, \quad  q_3^{(l)} = q_3^{(l)} - (G_3^{(l)})^{-1}A^{(l)}_{3,4}q_4^{(l)},\\
        q_2^{(l)} = q_2^{(l)} - (G_2^{(l)})^{-1}A_{2,3}^{(l)}q^{(l)}_3, \quad q_1^{(l)} = q_1^{(l)}-(G_1^{(l)})^{-1} A^{(l)}_{1,2}q_2^{(l)}
    \end{gathered}
\end{equation}

Применение предобуславливателя в форме \eqref{preconditioner in factored form} формально соответствует методу неполной факторизации IFM (Incomplete Factorization Method, см.~\cite{18}) с согласованием строчных сумм, т. е. по условию $B^{(l)}e^{(l)} = A^{(l)}e^{(l)}$ при $\theta = 1$. Отметим, что реализация такого подхода фактически включает этапы редукции, грубосеточной коррекции и продолжения, обычных для ``классических'' многосеточных алгоритмов, но без применения сглаживания, см. обзор в \cite{21}.

\section{Примеры численных экспериментов}

Мы проиллюстрируем эффективность предложенных подходов результатами предварительных экспериментальных исследований на представительной серии методических задач. Исходные уравнения решаются в кубической расчетной области с однородной средой на регулярных равномерных сетках с числом ячеек $32^3, 64^3, 128^3$. 

Также расчёты выполнены на искусственном керне размером $120*120*124$ ячейки, который построен как случайная разномасштабная пористая среда. Для генерации модели создан набор структурированных кубических сеток с размером $14^3$, $20^3$, $40^3$, $60^3$. На каждой сетке определяются случайные величины в центрах соответствующих ячеек в диапазоне от 0 до 1. Затем была произведена трикубическая интерполяция, как описано в статье \cite{19}, с этих сеток на конечную сетку $120*120*124$, после чего значения, полученные с грубых сеток, были просуммированы и на конечной сетке получено распределение значений от 0 до 4. Ячейки, в которых значение превышало определённое пороговое значение, были обозначены как проницаемые. Таким образом, была получена пористая структура, напоминающая реальный керн. Пороговое значение было выбрано такое, чтобы получить пористость керна 50\%. 

Единственное принципиальное отличие от реальных кернов со связной непроницаемой средой -- это наличие непроницаемых островков, окружённых проницаемыми ячейками. Однако это отличие не мешает тестированию решателей. Для построения окончательной расчетной сетки непроницаемые ячейки были удалены. В приведенных ниже таблицах расчеты на этой сетке обозначены как ``kern''. Внешний вид расчетной сетки показан на рисунке \ref{mesh_image}.

\begin{figure}[h]
   \centering
   \includegraphics[width=0.8\textwidth]{mesh.png}
   \caption{Искусственный керн размером 120*120*124 и центральное сечение перпендикулярное оси z.}
   \label{mesh_image}
\end{figure}

Для всех сеток величина шага постоянна, $h = 1мкм$, так что фактически при различных дискретизациях размеры расчётных областей отличаются. Нижняя и верхняя грани куба, перпендикулярные оси $z$, являются точками входа и выхода потока флюида, с заданными одинаковыми вертикальными скоростями $u_z = 1мм/с$. Боковые поверхности куба рассматриваются как непроницаемые стенки с условиями прилипания $u_x = u_y = u_z = 0$. Вязкость среды задается постоянной и равной $\mu = 5\cdot10^{-3} \text{Па}\cdot\text{с}$, а плотность -- $\rho = 1000 кг/м^3$. Внутри расчётной области начальная скорость и давление полагаются нулевыми. Во всех таблицах результаты для сеток $32^3$, $64^4$, $128^3$ относятся к задачам с однородной средой, а слово kern соответствует модели случайной пористой среды на сетке размером 120*120*124.

Вычисления проводились в Yandex Cloud на узле Intel Xeon Gold 6338 с 48-ю ядрами, 96-ю потоками и с частотой 2 Ghz. В нижеследующих таблицах приводятся результаты предварительных экспериментальных исследований для модельных задач. Целью расчётов было сравнение эффективности итерационных процессов и анализ масштабируемости распараллеливания при использовании различных временных шагов $\Delta t$ и разного числа параллельных процессов $P = 2, 4, 8, 16, 32, 48$. Во всех таблицах результаты приводятся для расчёта одного (первого) шага по времени. Разработанный авторами программный код для решения задач на неструктурированных сетках реализован на платформе INMOST \cite{INMOST_citation}, с существенным использованием функциональных инструментов библиотеки PETSc \cite{PETSc_citation}.

\begin{table}[h]
    \centering
    \caption{Числа итераций и времени расчёта для $\Delta t = 0.5$.}
    \begin{tabular}{|c|c|c|c|c|c|c|}
         \hline
         N\backslash P & 2 & 4 & 8 & 16 & 32 & 48
         \\ \hline
         32^3 & 40/10.1 & 41/6.2 & 41/4.1 & 41/3.3 & 41/2.9 & 41/2.9
         \\ \hline
         64^3 & 81/215.8 & 82/123.3 & 82/79.7 & 82/60.6 & 82/45.7 & 82/41.2
         \\ \hline
         128^3 & - & - & - & 167/1136.6 & 166/790.2 & 167/684.3
         \\ \hline
         kern & 42/236.8 & 42/122.2 & 42/63.2 & 42/37.2 & 42/22.8 & 42/18.0
         \\ \hline
     \end{tabular}
    \label{table_base_LDU_0_5}
\end{table}

\begin{table}[h]
    \centering
    \caption{Результаты расчётов для $\Delta t = 0.5\cdot 10^{-3}$.}
    \begin{tabular}{|c|c|c|c|c|c|c|}
         \hline
         N\backslash P & 2 & 4 & 8 & 16 & 32 & 48
         \\ \hline
         32^3 & 40/9.4 & 39/5.5 & 40/3.7 & 40/2.8 & 40/2.5 & 40/2.4
         \\ \hline
         64^3 & 72/183.0 & 73/100.4 & 73/62.2 & 74/47.1 & 74/32.1 & 74/30.8
         \\ \hline
         128^3 & - & - & - & 131/816.9 & 131/571.0 & 131/489.4
         \\ \hline
         kern & 44/251.0 & 44/121.9 & 44/65.0 & 44/38.8 & 44/24.5 & 44/19.3
         \\ \hline
     \end{tabular}
    \label{table_base_LDU_5e_4}
\end{table}

\begin{table}[h]
    \centering
    \caption{Результаты расчётов для $\Delta t = 10^{-7}$.}
    \begin{tabular}{|c|c|c|c|c|c|c|}
         \hline
         N\backslash P & 2 & 4 & 8 & 16 & 32 & 48
         \\ \hline         
         32^3 & 10/1.9 & 10/1.1 & 10/0.6 & 10/0.5 & 10/0.4 & 10/0.4
         \\ \hline         
         64^3 & 10/20.9 & 10/10.7 & 10/6.0 & 10/3.9 & 10/2.5 & 10/2.2
         \\ \hline         
         128^3 & 10/205.4 & 10/105.6 & 10/60.0 & 10/37.8 & 10/25.7 & 10/21.1
         \\ \hline         
         kern & 11/62.0 & 11/31.1 & 11/16.5 & 11/9.8 & 11/5.8 & 11/4.6
         \\ \hline
     \end{tabular}
    \label{table_base_LDU_1e_7}
\end{table}

\begin{table}[h]
    \centering
    \caption{Зависимость числа итераций от параметра регуляризации $\gamma$ (сетка $64^{3}$, см. табл.\ref{table_base_LDU_5e_4}).}
    \begin{tabular}{|c|c|c|c|c|c|c|c|c|c|}
        \hline
        \gamma & 0.05 & 0.01 & 0 & -0.1 & -0.2 & -1 & -2 & -4 & -8
        \\ \hline
        n & \infty & 83 & 74 & 44 & 35 & 20 & 17 & 16 & 10
        \\ \hline
    \end{tabular}
    \label{table_LDU_gamma1}
\end{table}


\begin{table}[h]
    \centering
    \caption{Зависимость числа итераций от параметра регуляризации $\gamma$ (сеточная модель kern, см. табл.\ref{table_base_LDU_5e_4}).}
    \begin{tabular}{|c|c|c|c|c|c|c|c|c|c|}
        \hline
        \gamma & 0.005 & 0.001 & 0 & -0.005 & -0.05 & -0.01 & -0.1
        \\ \hline
        n & \infty & 59 & 44 & 78 & 186 & 100 & 251
        \\ \hline
    \end{tabular}
    \label{table_LDU_gamma2}
\end{table}

\begin{table}[h]
    \centering
    \caption{Результаты расчёта для $\Delta t = 0.5$, при использовании диагональной компенсации \eqref{17} для аппроксимации дополнения Шура.}
    \begin{tabular}{|c|c|c|c|c|c|c|}
         \hline
         N\backslash P & 2 & 4 & 8 & 16 & 32 & 48
         \\ \hline
         32^3 & 23/5.4 & 23/3.9 & 23/2.7 & 24/2.3 & 23/2.0 & 23/1.9
         \\ \hline
         64^3 & 28/65.3 & 28/41.6 & 28/27.7 & 28/22.3 & 28/17.7 & 28/16.7
         \\ \hline
         128^3 & 31/635.4 & 31/387.1 & 31/284.8 & 31/229.4 & 31/163.1 & 31/136.6
         \\ \hline
         kern & 65/376.7 & 65/189.3 & 65/98.6 & 65/58.1 & 65/36.5 & 65/29.0
         \\ \hline
     \end{tabular}
    \label{table_LDU_diag_compens}
\end{table}

\begin{table}[h]
    \centering
    \caption{Результаты расчёта для $\Delta t = 0.5\cdot 10^{-3}$, алгоритм без предобуславливателя дополнения Шура.}
    \begin{tabular}{|c|c|c|c|c|c|c|}
         \hline
         N\backslash P & 2 & 4 & 8 & 16 & 32 & 48
         \\ \hline         
         32^3 & 14/2.1 & 14/1.4 & 14/1.0  & 14/0.81 & 14/0.77 & 14/0.74
         \\ \hline         
         64^3 & 15/21.4 & 15/13.3  & 15/9.2 & 15/7.6 & 15/5.7 & 15/5.3 
         \\ \hline         
         128^3 & - & 17/122.5 & 17/84.0 & 17/73.9 & 17/53.4 & 17/46.6 
         \\ \hline         
         kern & - & - & 490/354.2 & 490/208.6 & 490/126.1 & 489/104.2 
         \\ \hline
     \end{tabular}
    \label{wo_precond}
\end{table}

Исследования осуществлялись для итерационного процесса с предобуславливанием СЛАУ по формуле \eqref{two-sided}. При этом решение вспомогательной системы с дополнением Шура осуществлялось итерационными методами FGMRES или CG, с применением предобуславливателя $\overline{S}_{\gamma}$ (с разными значениями $\gamma$) или без него. Матрица $A^{-1}_{1,1}$ ``обращалась'' с помощью алгоритма BoomerAMG \cite{23} из библиотеки Hypre, а её аппроксимация $\overline{A}^{-1}_{1,1}$ определялась методом SPAI. После этого матрица $\overline{S}_{\gamma}$ ``обращалась'' (т.е. решалась соответствующая СЛАУ) также с помощью AMG.

В таблицах \ref{table_base_LDU_0_5}–\ref{table_base_LDU_1e_7} приведены результаты расчётов для $\gamma = 0$ при различных значениях временного шага: $\Delta t = 0.5$, $0.5 \cdot 10^{-3}$ и $10^{-7}$. Можно заметить, что число итераций и время решения растут с увеличением временного шага. Сетка ``kern'' ведёт себя нетипично: число итераций немного уменьшается при увеличении шага с 0.0005 до 0.5 секунд. Причина такого поведения неясна. 

Таблицы~\ref{table_LDU_gamma1} и \ref{table_LDU_gamma2} демонстрируют возможность уменьшения числа итераций за счёт оптимизации параметра регуляризации $\gamma$ (при $M = I$) из выражения \eqref{eq_diag_compens}. Можно заметить, что на расчетной сетке с однородным кубом $64^3$ число итераций уменьшается с отрицательной $\gamma$. Однако на сетке с искусственным керном никакого улучшения не наблюдается и параметр $\gamma$ выгоднее задавать равным нулю.

В таблице~\ref{table_LDU_diag_compens} даны результаты расчётов при $\Delta t = 0.5$, получаемые с использованием диагональной компенсации \eqref{17} при $\theta = 1$ для аппроксимации дополнения Шура. Получаемый при этом предобуславливатель $\overline{S}_{\gamma}$ строится без регуляризации, т.е. $\gamma = 0$, а его обращение, как и в предыдущих случаях, осуществляется с помощью алгоритма BoomerAMG.

Наконец, в таблице \ref{wo_precond} приведены результаты полученные при решении СЛАУ без предобуславливания, но с точной безматричной аппроксимацией дополнения Шура. На сетках с пустыми кубами $32^3$, $64^3$, $128^3$ наблюдается наивысшая эффективность расчета. При этом на сетке ``kern'' скорость счета падает во много раз.

Анализ представленных выше результатов позволяет сделать следующие предварительные выводы:
\begin{enumerate}
    \item[\textbullet] Решение сеточных начально-краевых задач для систем уравнений Стокса представляет собой достаточно трудоёмкую проблему, особенно при больших значениях временных шагов. В частности, значение $\Delta t=0.5$ в таблицах~\ref{table_base_LDU_0_5}, \ref{table_LDU_diag_compens} фактически соответствует стационарной задаче.

    \item[\textbullet] Эффекты, наблюдаемые в таблице \ref{wo_precond}, связаны с особенностью дополнения Шура для исследуемой матрицы. Оно представляет собой дивергенцию от обратного лапласиана, примененного к градиенту, что близко к единичной матрице I. Поэтому на однородных средах наблюдается высокая эффективность. На искусственном керне большую роль играют граничные условия прилипания, т.к. по всему керну распределены непроницаемые стенки и дополнение Шура далеко от матрицы I. 
    
    \item[\textbullet] Применение регуляризации показывает положительный эффект только на пустых кубах. Это связано с приближением аппроксимации дополнения Шура к единичной матрице. Можно сделать вывод, что применение такой регуляризации неэффективно для задачи о течении в порах.
    
    \item[\textbullet] Использование диагональной компенсации при ленточной аппроксимации обратной матрицы (алгоритм SPAI) в дополнении Шура позволяет существенно повысить производительность исследуемых методов на задачах с однородными средами, но неэффективно для пористой среды.
    
    \item[\textbullet] Эффективность распараллеливания рассматриваемых итерационных процессов можно считать вполне удовлетворительной, хотя здесь ещё есть большие резервы для ускорения за счёт учёта особенностей алгоритмов и архитектуры многопроцессорной вычислительной системы (МВС).

    \item[\textbullet] Для расчетной сетки ``kern'' с пористой средой из рассмотренных вариантов наибольшую производительность показывает предобуславливатель \eqref{two-sided} без дополнительных регуляризаций и диагональной компенсации.
\end{enumerate}

\section{Заключение}
Рассматриваемые в работе блочные методы неполной факторизации для системы сеточных уравнений Стокса, в комбинации с применением регуляризации, многосеточных подходов, алгоритмов аппроксимации обратных матриц с диагональной компенсацией, предобусловленных крыловских итерационных процессов и технологий распараллеливания вычислений на МВС с распределённой и общей памятью представляют собой широкое поле деятельности для перспективных разработок высокопроизводительного математического и программного обеспечения для решения актуальных проблем суперкомпьютерного моделирования. Следует отметить, что задача Стокса обычно представляет собой только часть больших междисциплинарных проблем. Однако рассматриваемые в работе вопросы имеют самодостаточных характер, и полученные предварительные результаты дают хорошую перспективу для дальнейших исследований в различных направлениях. Отметим, что в данной работе по техническим причинам использовался широко распространенный решатель BoomerAMG, однако предварительные расчеты показывают возможность ускорения вычислений при переходе на алгебро-геометрический многосеточный метод описанный в п. 3. Существенные перспективы видятся при использовании многосеточных методов непосредственно для задачи Стокса \cite{Oosterlee}, а также при распараллеливании алгоритмов с помощью декомпозиции областей \cite{7}.

\bigskip
\begin{thebibliography}{1}

\bibitem{1} M.A. Olshanskii, E.E. Tyrtyshnikov,
\href{https://doi.org/10.1137/1.9781611973464}{\it Iterative Methods for Linear Systems: Theory and Applications}, SIAM Philadelphia, 2014.

\bibitem{2} H.C. Elman, D.J. Silvester, A.J. Wathen,
\href{https://doi.org/10.1017/S0022112006220561}{\it Finite Elements and Fast Iterative Solvers: with Applications in Incompressible Fluid Dynamics}, Oxford University Press, New-York, 2005.

\bibitem{3} Yu.V. Bychenkov, E.V. Chizhonkov,
\href{https://www.researchgate.net/publication/316715243_Iterative_methods_for_saddle_point_problems_in_Russian}{\it Iterative Methods for Solving Saddle-Point Problems}, BINOM, Moscow, 2010.

\bibitem{4} F. Brezzi, M. Fortin,
\href{https://doi.org/10.1007/978-1-4612-3172-1}{\it Mixed and Hybrid Finite Element Method}, Springer Series in Computational Mathematics, 1991.

\bibitem{5} D. Boffi, F. Brezzi, M. Fortin,
\href{https://doi.org/10.1007/978-3-540-78319-0\_2}{\it Finite Elements for the Stokes Problem}, Lecture Notes in Mathematics. Mixed Finite Elements, Compatibility Conditions, and Applications, Springer, Berlin, Heidelberg, {\bf 1939} (2008), 45--100.

\bibitem{6} M. Benzi, G.H. Golub, J. Liesen,
\href{https://doi.org/10.1017/S0962492904000212}{\it Numerical Solution of Saddle Point Problems}, Acta Numerica, {\bf 14} (2005), 1--137.

\bibitem{7} V.P. Il’in,
\href{https://doi.org/10.1134/S0965542521110099}{\it Iterative Preconditioned Methods in Krylov Spaces: Trends of the 21st Century}, Computational Mathematics and Mathematical Physics, {\bf 61} (2021), 1750--1775.

\bibitem{8} V.P. Il'in,
{\it Mathematical modeling Part 1. Continuous and dicrete models}, SB RAS, Novosibirsk, 2017.

\bibitem{9} R. Estrin, C. Greif,
\href{https://doi.org/10.1137/16M1102410}{\it SPMR: A Family of Saddle-Point Minimum Residual Solvers}, SIAM Journal on Scientific Computing, {\bf 40}:3 (2018), 1884--1914.

\bibitem{10} F.J. Gaspar, Y. Notay, C.W. Oosterlee, C. Rodrigo,
\href{https://doi.org/10.1137/130920630}{\it A Simple and Efficient Segregated Smoother for the Discrete Stokes Equations}, SIAM Journal on Scientific Computing, {\bf 36}:3 (2014), A1187--A1206.

\bibitem{11} N. Gould, D. Orban, T. Rees  
\href{https://doi.org/10.1137/130916394}{\it Projected Krylov Methods for Saddle-Point Systems}, SIAM Journal on Matrix Analysis and Applications. {\bf 35} (2014), 1329--1343. 

\bibitem{12} C. Greif, D. Schoetzau
\href{https://doi.org/10.1137/130916394}{\it Preconditioners for Saddle Point Linear Systems with Highly Singular (1,1) Blocks}, Electronic Transactions on Numerical Analysis,  {\bf 22}, (2006), 114--121.

\bibitem{Oosterlee} C.W. Oosterlee, F.J.G. Lorenz,
\href{https://doi.org/10.1109/MCSE.2006.115}{\it Multigrid Methods for the Stokes System. Computing in Science and Engineering}, {\bf 8}:6 (2006), 34--43.

\bibitem{Notay} Y. Notay,
\href{https://doi.org/10.1137/16M1071419}{\it Algebraic Multigrid for Stokes Equations}, SIAM Journal on Scientific Computing, {\bf 39}:5 (2017), S88--S111. 

\bibitem{Benzi} M. Benzi, G.H. Golub, J. Liesen,
\href{https://doi.org/10.1017/S0962492904000212 }{\it Numerical solution of saddle point problems} Acta Numerica, {\bf 14} (2005), 1–-137.

\bibitem{14} Ya. Gurieva, V. Il'in,
\href{https://www.mathnet.ru/php/archive.phtml?wshow=paper&jrnid=znsl&paperid=7012&option_lang=rus}{\it On methods of conjugate directions for multiple solution of SLAE} Zapiski POMI, {\bf 496} (2020), 26--42.

\bibitem{15} M.I. Ivanov, I.A. Kremer, M.V. Urev,
\href{https://doi.org/10.1134/S1995423919040049}{\it Solving the Pure Neumann Problem by a Finite Element Method}, Numerical Analysis and Applications, {\bf 12} (2019), 359–-371.

\bibitem{16} M. Benzi, C.D. Meyer, M. T?ma,
\href{https://doi.org/10.1137/S1064827594271421}{\it A Sparse Approximate Inverse Preconditioner for the Conjugate Gradient Method}, SIAM Journal on Scientific Computing, {\bf 17}:5 (1996), 1135--1149.

\bibitem{17} V.P. Il’in,
\href{https://doi.org/10.1007/s10958-023-06446-6}{\it Multigrid Incomplete Factorization Methods in Krylov Subspaces}, Journal of Mathematical Sciences, {\bf 272} (2023), 523–-532.

\bibitem{18} V.P. Il’in,
\href{https://icmmg.nsc.ru/sites/default/files/pubs/scanbook2pro.pdf}{\it Finite Element Methods and Technologies}, SB RAS, Novosibirsk, 2017.

\bibitem{20} K. St\"uben,
\href{https://doi.org/10.24406/publica-fhg-290270}{\it Algebraic multigrid (AMG). An introduction with applications}, Sankt Augustin, 1999.

\bibitem{19} F. Lekien, J.E. Marsden,
\href{https://doi.org/10.1002/nme.1296}{\it Tricubic interpolation in three dimensions} International Journal for Numerical Methods in Engineering, {\bf 63}:3 (2005), 455-–471.

\bibitem{21} K. St\"uben,
\href{https://doi.org/10.1016/s0377-0427(00)00516-1}{\it A review of algebraic multigrid}, Journal of Computational and Applied Mathematics, {\bf 128}:1-2 (2001), 281--309.

\bibitem{23} V. Henson, U. Yang, 
\href{https://doi.org/10.1016/S0168-9274(01)00115-5 }{\it BoomerAMG: A parallel algebraic multigrid solver and preconditioner}, Applied Numerical Mathematics. {\bf 41}:1 (2002), 155--177.

\bibitem{PETSc_citation}S. Balay, W.D. Gropp, L.C. McInnes, B.F. Smith, 
\href{https://doi.org/10.1007/978-1-4612-1986-6_8 }{\it Efficient Management of Parallelism in Object-Oriented Numerical Software Libraries}, Modern Software Tools for Scientific Computing. Birkh\"auser, Boston, MA. (1997), 163--202.

\bibitem{INMOST_citation}I. Konshin, I. Kapyrin, K. Nikitin, K. Terekhov,
\href{http://dx.doi.org/10.1007/978-3-319-32152-3_26 }{\it Application of the parallel INMOST platform to subsurface flow and transport modelling}, In: Proc. 11th International Conference on Parallel Processing and Applied Mathematics (2015), 277--286.


\end{thebibliography}

\end{document}