まえがき

こちらの論文が興味深い題材だったので、使用されている数値手法について記事にしてみました。
実装はこちらです。

TL;DR

差分法で非圧縮のNavier-Stokes方程式を数値計算する場合に現れる、Poisson方程式の直接解法についてです。
全方向に等間隔格子、もしくは１方向だけ不等間隔格子の場合にはFFTベースの解き方がよく使われますが、格子が２方向以上に不等間隔な場合にも応用が効くよ（本質は同じだよ）というお話です。

背景

差分法で非圧縮流体を考える際に解かれる方程式についてかなり簡単に触れます。
運動量の式を積分しながら非圧縮条件を満たすため、Fractional-step法と呼ばれる方法が標準的に利用されます。
すなわち、まず運動量の式を時間方向に数値積分し：

u_i^*
=
u_i^n
+
\int
\left(
    -
    u_j \frac{\partial u_i}{\partial x_j}
    -
    \frac{\partial p}{\partial x_i}
    +
    \frac{1}{Re} \frac{\partial}{\partial x_j} \frac{\partial u_i}{\partial x_j}
    +
    a_i
\right)
dt,

その後速度場が非圧縮条件を満たすように修正する：

\frac{
    u_i^{n + 1}
    -
    u_i^*
}{
    \Delta t
}
=
-
\frac{\partial \phi}{\partial x_i}

というアプローチです。
ここでスカラー$\phi$は以下のPoisson方程式を満たすように決められます：

\frac{\partial}{\partial x_i} \frac{\partial \phi}{\partial x_i}
=
\frac{1}{\Delta t}
\frac{\partial u_i^*}{\partial x_i}.

これを不等間隔格子の上でどのように解くかが今回のテーマです。

導出は有名ですので省略します。
ここまでの添字はダミーインデックスで、項に複数現れる場合には総和をとります。
簡単のためニュートン流体で均一密度・粘性を仮定しましたが、そうでなくても基本は同じです。

不等間隔格子

この記事の主眼は不等間隔格子についてですが、なぜそんなものを考える必要があるかという点についても簡単に触れます。

現実世界では理想的な形の系は稀ですが、流体の基礎研究分野では矩形領域（２次元・３次元）を考えることが多いです。
流体の動きをよく記述できるとされる非圧縮Navier-Stokes方程式は空間微分を含む方程式ですので、各方向に適切な境界条件を設定する必要があります。
上で示したPoisson方程式も同じです。

（私見ですが）最も数値計算で扱いやすいのは周期境界条件で、等方性乱流やその亜種の計算に多く用いられる印象です。
周期的ですので、どこか特定の箇所にメッシュを細かく設定するということは少なく、等間隔の格子を用いるのが自然な選択です。

流れ場の状態に応じてメッシュを逐次調整するAdaptive Mesh Refinementのことはここでは忘れることにします。
また周期境界であればスペクトル法が有力な選択肢ですが、今回は差分法、特に二次精度中心差分に絞ります。

基礎研究では壁面を伴う流れも頻繁に扱われます。
この場合には滑りなし条件を課すような壁をドメイン端に設定するのが最も一般的だと思いますが、壁面付近では粘性の影響が強く現れる領域（境界層）が存在します。
特にReynolds数が大きい場合には壁方向の速度変化が急峻になるため、壁面付近での速度場を正しく捉えるには、壁に垂直な方向のメッシュを細かく設定する必要があります。
どの程度細かくするかは色々議論がありますが、モデルを導入しない直接数値計算ではwall unitで１未満に必ず１点以上を置く（かなり粗めの基準）と思います。

後述しますが、数値計算アルゴリズムの観点からはメッシュが等間隔であれば扱いやすいです。
全領域で等幅でwall unit以下の細かいメッシュを設定できれば良いのですが、現実的な計算コストを軽く超えてしまうことがあります。
従って壁面付近のみを細かくし、壁から離れるにつれて徐々に（乱流であればKolmogorovスケール程度まで）粗くしていくのが一般的で、これの実現のために不等間隔格子の利用が必要です。
記事先頭の画像左のような形ですね。

離散Poisson方程式

Poisson方程式の離散化とその数値的扱いを考えます。
前述した２次元矩形領域でのPoisson方程式を二次精度中心差分で離散化すると以下のようになります：

\left( {l}_1 \right)_{i}
{f}_{i-1, j}
-
\left( {\left( {l}_1 \right)_{i}} + {\left( {u}_1 \right)_{i}} \right)
{f}_{i, j}
+
\left( {u}_1 \right)_{i}
{f}_{i+1, j}
+
\left( {l}_2 \right)_{j}
{f}_{i, j-1}
-
\left( {\left( {l}_2 \right)_{j}} + {\left( {u}_2 \right)_{j}} \right)
{f}_{i, j}
+
\left( {u}_2 \right)_{j}
{f}_{i, j+1}
=
{g}_{i, j}.

ここで添字は配列のインデックスを表していて、係数は以下の通り定義しています：

{{\left( {l}_1 \right)_{i}} = \frac{{1}}{{\left( {h}_{{\xi}^1} \right)_{i}}} \frac{{1}}{{\left( {h}_{{\xi}^1} \right)_{i-\frac{1}{2}}}}},

{{\left( {u}_1 \right)_{i}} = \frac{{1}}{{\left( {h}_{{\xi}^1} \right)_{i}}} \frac{{1}}{{\left( {h}_{{\xi}^1} \right)_{i+\frac{1}{2}}}}},

{{\left( {l}_2 \right)_{j}} = \frac{{1}}{{\left( {h}_{{\xi}^2} \right)_{j}}} \frac{{1}}{{\left( {h}_{{\xi}^2} \right)_{j-\frac{1}{2}}}}},

{{\left( {u}_2 \right)_{j}} = \frac{{1}}{{\left( {h}_{{\xi}^2} \right)_{j}}} \frac{{1}}{{\left( {h}_{{\xi}^2} \right)_{j+\frac{1}{2}}}}}.

離散化は例えばこちらに従い、一般直交座標での方程式に対して行なっています。
見た目はいかついですが、scale factor $h_{{\xi}^i}$はその場所での格子幅と同じですので、実際は格子幅を係数として自身と周囲４点を含む線形方程式です。

以降は全方向が滑りなし壁面の場合、つまりFractional-step法の文脈で現れるPoisson方程式の境界条件としてNeumann条件（勾配が０）を課す場合を考えます。
離散Poisson方程式においては端点で$f$の差分が$0$となるような処置を行います。
すなわち上の線型方程式の係数に補正が入り、例えば$x$方向マイナス側の端（$i = 1$）では$f_{0, j} = f_{1, j}$から

-
{\left( {u}_1 \right)_{1}}
{f}_{1, j}
+
\left( {u}_1 \right)_{1}
{f}_{2, j}
+
\left( {l}_2 \right)_{j}
{f}_{1, j-1}
-
\left( {\left( {l}_2 \right)_{j}} + {\left( {u}_2 \right)_{j}} \right)
{f}_{1, j}
+
\left( {u}_2 \right)_{j}
{f}_{1, j+1}
=
{g}_{1, j}

となります。
他の境界での扱いも同様です。

線型方程式の解法

線形方程式ですので、理屈の上では行列を愚直にひっくり返すことで解くことができますが、行列のサイズ（一辺が$N_x \times N_y$の正方行列）を考えると困難であることがわかります。
昨今の計算では$N$は1000超えが普通で、10000を超えることも珍しくないですので、行列だけでテラバイト・ペタバイトのサイズ感となり、そうそうRAMに乗りません。
データ量の問題が仮にどうにかなるとして、ほとんどが０の疎行列ですので、無駄な計算が非常に多くなってしまいます。
ということで疎行列に対する現実的な手法が色々と存在する訳ですが、ざっくり直接解法と反復解法に分類できます。

反復解法は誤差を評価し、何らかの方法で誤差を小さくするように解を修正していく手法全般を指し、共役勾配法やマルチグリッド法が実用に耐えうる手法の代表例です。
これらは直接解法比でより汎用的で、計算量の面で大きなメリットもありますが、閾値の設定・前処理・初期値への依存性・並列化などの考慮が必要です。

一方で直接解法は線形方程式と同値の問題に対して直接解を求めるものです。
反復解法比での大きな特徴としては、どんな問題に対しても（行列および右辺値によらず）一定の計算量であること、並列化が比較的容易であること、収束しない心配がないことが挙げられます。
（個人的には最後は実計算を行う上でありがたい性質だと思っています。）

今回以下でフォーカスするのは、直接解法の一つである固有値分解を利用した方法です。

固有値分解（等間隔格子の場合）

まずは等間隔格子の場合を簡単におさらいします。

等間隔格子であれば$h_{{\xi}^1}$は$\Delta x$（定数）となり、離散コサイン変換を利用した高速解法が適用できます（周期境界の場合は離散実フーリエ変換）。
有名な方法で各所に記事があるので詳細は割愛しますが、各$j$に対して$x$方向に変換を行うことで、上の離散Poisson方程式は以下のように変形できます：

{\lambda_k} {{F}_{k, j}}
+
{\left( {l}_2 \right)_{j}} {{F}_{k, j-1}}
-
\left( {\left( {l}_2 \right)_{j}} + {\left( {u}_2 \right)_{j}} \right) {{F}_{k, j}}
+
{\left( {u}_2 \right)_{j}} {{F}_{k, j+1}}
=
{{G}_{k, j}}.

ここで$\lambda_k$は波数（固有値）で、

-
\left(
    \frac{2}{\Delta x}
    \sin \left( \frac{\pi k}{2 N_x} \right)
\right)^2

で与えられます。
なお固有値０が含まれるのは、Neumann境界条件の設定によるものです（定数分だけ解に自由がある）。

$x$方向がスペクトル空間に変換されたことで$x$方向の差分（$i - 1$や$i + 1$）が消えました。
$y$方向には引き続き３点のステンシルが残っていますが、これは三重対角行列ですので、Thomas algorithmなどで比較的簡単に解くことが可能です（$y$方向には等間隔格子である必要もありません）。
この場合の必要計算量は、$x$方向への離散コサイン変換を各$j$に対して行うことから、$O \left( N_x N_y \log N_x \right)$ と見積もれます。
フーリエ変換の効率の高さを最大限に活かせる手法です。

三重対角行列を解くコストは$O \left( N_y \right)$で、各固有値$\lambda_k$に対して解くため全体では$O \left( N_x N_y \right)$となります。
フーリエ変換比でBig-Oとしては小さいと言えるため、ここでは無視しています。

固有値分解（不等間隔格子の場合）

残念ながら２方向の格子が不等間隔の場合にはフーリエ変換ベースの方法は使用できませんが、同じアイデア（固有値分解を利用した$x$方向の対角化）が適用できます。

１次元

まずは簡単のため$x$方向のみを取り出して、１次元での離散Poisson方程式：

{{\left( {l}_1 \right)_{i}} {{f}_{i-1}} - \left( {\left( {l}_1 \right)_{i}} + {\left( {u}_1 \right)_{i}} \right) {{f}_{i}} + {\left( {u}_1 \right)_{i}} {{f}_{i+1}} = {{g}_{i}}}

に対する固有値分解を考えます。
これは線形方程式ですので、以下のように書けます：

\boldsymbol{L} \vec{f} = \vec{g}.

$\boldsymbol{L}$は非対称の実三重対角行列ですが、これを固有値分解することを考えます。
任意の実対称行列$\boldsymbol{S}$は直交行列$\boldsymbol{Q}$（およびその転置$\boldsymbol{Q}^T$）と対角行列$\boldsymbol{\Lambda} = \text{diag} \left( \vec{\lambda} \right)$を用いて

\boldsymbol{S} = \boldsymbol{Q} \boldsymbol{\Lambda} \boldsymbol{Q}^T

と書けますので、これを利用します。
$\boldsymbol{L}$は対称行列ではないので、何かしらの手段で対称行列に変換する必要がありますが、元の離散Poisson方程式の構造から、

\boldsymbol{D}_{ii} = \sqrt{{\left( {h}_{{\xi}^1} \right)_{i}}}

を作用させればよい、すなわち

\boldsymbol{S}
=
\boldsymbol{D}
\boldsymbol{L}
\boldsymbol{D}^{-1}

とわかります。
なお$\boldsymbol{D}$は対角行列ですので、逆行列$\boldsymbol{D}^{-1}$は単に各要素の逆数を取ったものとなります。
これらをまとめると、目的の方程式は以下のように記述できます：

\vec{f} 
=
\boldsymbol{D}^{-1} \boldsymbol{Q} \boldsymbol{\Lambda}^{-1} \boldsymbol{Q}^T \boldsymbol{D} \vec{g}.

なお式変形の途中で$\boldsymbol{Q}$が直交行列であること、すなわち

\boldsymbol{Q}^{-1}
=
\boldsymbol{Q}^T

を利用しました。

以下で実際の数値計算で行う手順として整理します。

ステップ１：対称化

$\boldsymbol{S}$を計算します。
対角成分は変わらず、その上下（対称なので同じ）の成分は以下のようになります：

\boldsymbol{S}_{i, i + 1}
\equiv
\boldsymbol{S}_{i + 1, i}
=
\frac{
    1
}{
    \sqrt{\left( {h}_{{\xi}^1} \right)_{i}}
}
\frac{
    1
}{
    \sqrt{\left( {h}_{{\xi}^1} \right)_{i + 1}}
}
\frac{
    1
}{
    \left( {h}_{{\xi}^1} \right)_{i + \frac{1}{2}}
}

ステップ２：固有値分解

実対称（三重対角）行列$\boldsymbol{S}$を固有値分解します。
結果として固有値$\boldsymbol{\Lambda}$（対角行列）と直交行列$\boldsymbol{Q}$（固有ベクトルを列ベクトルとして集めたもの）が得られます。

固有値分解には例えば以下が利用できます：

LAPACK：stemr
Python NumPy：linalg.eigh

ステップ３：必要な行列を計算

以下のサイズ$N_x$の正方行列２つを計算し、保持しておきます：

\boldsymbol{Q}^T \boldsymbol{D},

\boldsymbol{D}^{-1} \boldsymbol{Q}.

$\boldsymbol{D}$は対角行列ですので、それぞれ$\boldsymbol{Q}^T$・$\boldsymbol{Q}$の列・行への掛け算・割り算を行うだけです。

ステップ４：右辺値を変換する

変換行列を作用させることで、右辺値$\vec{g}$をスペクトル空間$\vec{\gamma}$に持っていきます：

\vec{\gamma} 
=
\left( \boldsymbol{Q}^T \boldsymbol{D} \right) \vec{g}.

右辺の行列はステップ３で計算済みです。

ステップ５：Poisson方程式を解く

スペクトル空間でのPoisson方程式は、単なる固有値による除算です：

\vec{\phi}
=
\boldsymbol{\Lambda}^{-1} \vec{\gamma},

すなわち

\vec{\phi}_i
=
\frac{\vec{\gamma}_i}{\vec{\lambda}_i}.

ステップ６：解を物理空間に戻す

変換行列（ステップ４で利用したものの逆行列）を作用させることで実施します：

\vec{f} 
=
\left( \boldsymbol{D}^{-1} \boldsymbol{Q} \right) \vec{\phi}.

こちらも右辺の行列はステップ３で計算済みです。

まとめ

ステップ１から３までは左辺の行列のみで完結しており、右辺値を一切利用していないのがポイントです。
固有値分解は比較的重めの計算ですが、これは格子生成を行なった後に一度だけ実施すれば良い操作であり、流体方程式を積分する間には行う必要がありません。
流体ソルバー内でPoisson方程式を解く際に必要になるのはステップ４から６のみであり、毎ステップ行う操作は非常に単純です。

２次元

行列$\boldsymbol{Q}^T \boldsymbol{D}$を対角ブロックに置いたブロック行列を作用させることで、２次元へも直接拡張が可能です。
行列表記すると大層な計算に見えますが、数値計算コードとしては各$j$に対して別々に固有値分解を実施することと同値です。
よって３次元に対する拡張も一直線です。

サイズ$N_x$の正方行列とベクトルの積を各$j$に対して行う（実装上は$N_x$の正方行列と$N_x$行$N_y$列の行列積の）計算となるので、計算量は$O \left( N_x^2 N_y \right)$と見積もれます。
ここまで来ると、FFTベースの手法も核はあくまで固有値分解であり、行列積を（計算量ベースで）効率的に求めるためにFFTを利用しているだけ（固有ベクトルの成分が三角関数由来の周期性を持つため可能）、ということが明らかです。

結果と考察

GitHubのように実装して検証を行いました。
時間精度と拡散項の取り扱いでいくらか工夫していますが、いずれも一般的な内容です。

Reynolds数無限の極限で二乗量が保存する程度のことは確認していますが、基本的にPoisson方程式が解けているかにしか着目していないので、verificationは割と適当です。。

直接解法ですので、理論上は修正速度場の発散は０になります。
速度場の発散をモニターすると、確かにFFTベースと同等の誤差程度であることが確認でき、直接解法として成立していることが示唆されます。
GPU環境が手元にないためパフォーマンスに関しては明言できませんが、行列積はGPUの得意領域ですから、実用上（Nが1000程度）十分な効率が出るのではないかと予想できます。

（パフォーマンスに関して実測なしで議論するのはタブーですが）FFTベースの手法と比べてどうなのよ、というのは気になる点です。
計算量ベースで同じ自由度で比較しますと、大きめのNであればFFTに軍配が上がります。
ただし同じ物理現象をシミュレートする（境界層をきちんと解像する）という観点で比較した場合、等間隔格子の設定が必須なFFTでは、行列積の方法に比べて格段に自由度を大きく取る必要があります。
これを加味すると、原則として行列積ベースの手法の方が好ましいという結論になりそうです。