数理解析研究所講究録
1129 巻 2000 年 108-114

108

球船内 Stokes 問題の並列計算
Parallel computation of a Stokes problem in a spherical shell
鈴木厚 (Atsushi Suzuki)

九州大学大学院数理学研究科
1. はじめに遅い流れの非圧縮流体を記述する Stokes 方程式について考える. Stokes 方程式は, 例えば, 地球マントル対流の数学モデルである無限 Prandtl 数流体の Rayleigh-Benard 方程式の流速と圧力に関する支配方程式をなす. 3 次元領域の問題に対し効率よく数値解を求めるアルゴリズムを開発することが必要である.
偏微分方程式に対し, 有限要素法による離散化を行い CG 法に代表される Krylov 部分空間反復法などを用いて) 連立方程式を解く場合, 剛性行列の記憶に多くのメモリ $-$ が必要である.
取り扱う領域が球殻であることに着目し, 領域の対称性を利用する. ある部分領域の剛性行列のみを記憶することにより, メモリー使用量の削減を図る. 反復法の前処理に部分領域毎の不完全 Cholesky 分解を用いることで, 並列化が容易に実現できる.

2. 支配方程式
3 次元球殻領域で, 滑り境界条件を課す Stokes 方程式について考える. 3 次元球殻領域$\{(x_{1}, x_{2}, x_{s})\in \mathbb{R}^{3} ; R_{1}<\sqrt{x_{1}^{2}+x_{23}2+X2}<R_{2}\}$ で定義された流速, 圧力 に関する次の Stokes 方程式を考える.

$\Omega=$$u=(u_{1}, u_{2}, u_{\mathrm{a}})$$-\triangle u+\nabla p=f$$(P)\{$$\nabla\cdot u=0$,$x\in\Omega$,$p$.

境界条件に滑り境界条件を課す :$u\cdot n=0$,$\{$$t^{(k)}\cdot D(u)n=0$$(k\backslash =1,2)$$x\in\partial\Omega$,.

は境界での外向き単位法線ベクトル, $t^{(1)(2)},$ は境界での独立な単位接ベクトル, $D(u)$ は変形速度テンソルである: $D_{ij}(u):= \frac{1}{2}(\partial_{j}u_{i}+\partial_{i}u_{j})(1\leq i,j\leq 3)$ . 粘性係数は, 定数 1 の場合を考え$n$$\theta$る.

3.

$\mathrm{P}1/\mathrm{P}1$

安定化有限要素法と離散化行列

Stokes 方程式を 3 次元領域で解くことを考慮し, 有限要素近似において最低次最小自由度の四面体要素, $P1/P1$ 要素を用いる. しかし, $P1/P1$ 要素は Stokes 方程式の混合型有限要素近似において必要となる下限上限条件を満たさない. この条件の克服のため, 最小 2 乗型 Galerkin 安定化有限要素法 [1] を用い離散化を行う.
はの多角形近似とする.の四面体要素による分割とする:.
五を$K$は省略してを最大要素直径とする. 煩雑さを避けるため, 領域に関するを各要素の直径,$\overline{\Omega}_{h}$$\overline{\Omega}=\bigcup_{K\in \mathcal{T}_{h}}\overline{K}$$h$$\Omega_{h}$$\Omega$$h_{K}$$h$109表記する.の有限要素空間を用いる:

を–次多項式の関数からなる空間とする. 流速, 圧力に対し次$S_{h}(\subset H^{1}(\Omega)\cap C0(\overline{\Omega}))$$W_{h}:=\{v_{h}\in S_{h}^{3} ; (v_{h}\cdot n_{\Omega})(P)=0(\forall P)\}$,$V_{h}:=\{v_{h}\in W_{h} ; (v_{h}, v^{(k)})=0(k=1,2,3)\}$$M_{h}$,:=&,$Q_{h}:=\{q_{h}\in S_{h} ; (q_{h}, 1)=0\}$.

ここで, $v^{(k)}=e^{(k)}\cross X(k=1,2,3)$ は剛体回転の自由度を表すベクトルである.
軸方向の単位ベクトルをとした. $P$ はは上の節点,の外向き単位法線ベクトルを表す.
$x_{k}$$e^{(k)}$$\partial\Omega$$\partial\Omega$$(\cdot, \cdot)$$n_{\Omega}$

を$L^{2}(\Omega)$または$L^{2}(\Omega)^{3}$での内積とする.

は球殻領域で滑り境界条件を課す Stokes 方$q\in S_{h}$ に対して, 双–次形式を次のよ$V_{h}\mathrm{x}Q_{h}$程式の–意可解性のために訴要な空間である.
うに定義する.

$u,$ $v\in S_{h}^{3},$$a(u, v):=2 \int_{\Omega}\sum_{1i,j\leq 3}Dij(u)D_{ij}(v)dx\leq$ ’$b(v, q):=- \int_{\Omega}\nabla\cdot vqdx$

安定化有限要素法によるスキームは$a(u_{h}, v_{h})+b(vh_{J}.

$(P_{h})$を満たすph)=(f, v_{h})$.

$(u_{h}, p_{h})\in V_{\iota},\cross Q_{h}$を求めることである.

,$v_{h}\in V_{h}$$(P_{h})\{$$b_{h}(u_{h}, q_{h})- \delta\sum h2K(\nabla ph, \nabla q_{h})_{K}=-\delta\sum h_{K}^{2}(f, \nabla q_{h})_{K_{)}}$$K\in \mathcal{T}_{h}$ここで$(\cdot, \cdot)_{K}$$q_{h}\in Q_{h}$.

$K\in \mathcal{T}_{h}$は要素 $K$ での内積を表す. 正定数$\delta$は安定化パラメータである.の有限要素基底を用いることは煩雑である. このため,の有限要素基底を用いるが, この場合係数行列は正則でない. そこで, 剛体回転の自由度と圧力の定数の任意性を取り除くための正射影を付加する.
$N_{W}=\dim W_{h},$ $N_{M}=\dim M_{h}$ , とし, 流速, 圧力に対する有限要素基底をそれぞれ,行列成分計算に$V_{h}\cross Q_{h}$$W_{h}\cross$$M_{h}$$W_{h}=\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}[\varphi_{1\varphi]},$$\cdots NW$ ’$M_{h}=\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}[\psi 1, \cdots\psi_{N_{M}}]$とし, 有限要素解の節点での値を$M_{h},$$p_{h}= \sum_{j}^{N_{M}}=1P\psi_{j}j$$V^{(k)}\in \mathbb{R}^{N_{W}}$$\{U_{j}\}_{j=1}^{N}W,$$\{\text{乃}\}_{j1}^{N_{M}}=$とする:$u_{h}\in W_{h},$$p_{h}\in$.

を v(りの節点での値からつくられるベクトル:$(1,1, \cdots, 1)^{T}\in \mathbb{R}^{N}\mathrm{n}f$$u_{h}--- \sum_{j=1}^{N},WU_{j\varphi}j,$$v^{(k)}= \sum_{j=1}N_{W}V^{()}\varphi_{j}jk$とする.
$N_{V}:=\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}[V(1), V^{(}2), V(3)]$$N_{Q}:=\mathrm{s}_{\mathrm{P}}\mathrm{a}\mathrm{n}[C]$,,$N:=\mathrm{s}_{\mathrm{P}}\mathrm{a}\mathrm{n}[(V(1), 0), (V^{(2)}, 0), (V^{(3)}, 0), (0, C)]$

, また $C=$

110に対し, 流速, 圧力に対する正射影と$P_{N}\perp$を次のように定義する:$P_{V}$:$\mathbb{R}N_{lV}arrow N_{V}^{\perp}$,$P_{V}U=U- \sum_{1k=}\mathrm{a}\frac{(U,V^{(k)})}{||V^{(k})||^{2}}V^{(k}$$P_{Q}$:$\mathbb{R}^{N_{M}}arrow N_{Q}^{\perp}$,$P_{Q}P=P- \frac{(P,C)}{||C||^{2}}C,$$\mathbb{R}^{N_{M}+N_{M}}arrow N^{\perp}$$P_{N^{\perp}}:$$(P_{h})$$P_{C}$$P_{V},$),$U\in \mathbb{R}^{N_{W}}$$P\in \mathbb{R}^{N_{M}}$,,$P_{N^{\perp}}:=$ .

,に対応する連立–次方程式は$A==$$(M)$となる. 行列 $A,$ $B$ および $D$ はそれぞれ双–次形式,.

および安定化項に対応する行法は破綻する可能性があるが, 反復列である. この行列は対称であるが正定値ではないため,解法として用いる. また $A$ の核は . $N$ である. CG 法の反復の各ステップにおいて, 行列とベクトルの乗算, ベクトルの加減算を行う際, 数値誤差の混入により, $N$ の成分が現れることを防ぐ$a(\cdot, \cdot)$$b(\cdot, \cdot)$$\mathrm{C}\mathrm{G}$ため, 正射影を用いる.

$P_{N}\perp$

4. 領域の対称性を利用したアルゴリズム
全体領域を, 8 個の部分領域に分ける.

また$\Omega^{(4)},$$\cdots$,.

$\Omega^{(7)}$:$\overline{\Omega}=\bigcup_{0\leq d\leq 7}\overline{\Omega}^{(d)}$$\Omega^{(0)}:=\{(x_{1}, x_{2}, x_{3})\in\Omega ; x_{1}>0, x_{2}>0, x_{3}>0\}$,$\Omega^{(1)}:=\{(x_{1}, x_{2,3}X)\in\Omega ; x_{1}<0., x_{2}>0, x_{3}>0\}$,$\Omega^{(2)}:=\{(X_{1}, x_{23}, X)\in\Omega ; x_{1}<0, x_{2}<0, x_{3}>0\}$,$\Omega^{(3)}:=\{(x_{1}, x_{23}, X)\in\Omega ; x_{1}>0, x_{2}<0, x_{3}>0\}$.

は $x_{3}<0$ である部分領域とする. この分割に対し, 直交変換 $R^{(d)}$ が存在し,$\Omega^{(d)}=R^{(d)}\Omega^{\langle 0})(0\leq d\leq 7)$が成り立つ.

次に行列の分割を定義する. 説明の簡単化のため, スカラー変数の圧力に対する行列 $D$ の分割を扱うが 流速に対する剛性行列に関しても同様に考えることができる.

添字の集合を次のように定義する:

A $:=\{1,2, \cdots, N_{M}\}$ ,$\Lambda_{I}$$:=\{1,2, \cdots, N_{M.I}\}$,$\Lambda_{B}:=\{N_{M},I+1, \cdots, N_{M}\}$

ここで $N_{M,I}<N_{M}$ であり, 添字に対応する節点に対して次が成り立つ:$i\in\Lambda_{I}$:$P_{i}$,

j\in AB: 乃,$\exists e,$$P_{i}\in\Omega^{(e)}$,$P_{j}\in\cup(\partial\Omega(e)_{\cap}\partial\Omega^{(f)})0\leq e<f\leq 7^{\cdot}$.

111は部分領域問の境界上の節点に対応する添字の集合である. 部分領域のように定義する:$A_{B}$$\Lambda^{(0)}:=\{1,2, \cdots, N_{M}^{(0)}\}$$\Lambda^{(d)}=\Lambda_{I}^{(d)}\cup\Lambda_{B}^{(d)},$,$\Lambda_{I}^{(0)}:=\{1,2, \cdots, N_{M,I}^{(0)}\}$$\Lambda_{I}^{(d)}\cap\Lambda_{B}^{(d)}=\emptyset(1\leq d\leq 7)$,$\overline{\Omega}^{(0)}$の添字の集合を次$A_{B}^{(0)}:=\{N_{M,I^{+1}}^{()}, \cdot*0.

はそれぞれ部分領域$\overline{\mathrm{f}}l^{(d)}$, N_{M}^{(0)}\}$.の添字の集合であり,次が成り立つものとする:$\Lambda_{I}=\bigcup_{d0\leq\leq 7}\Lambda^{(d}I)$$\exists \mathcal{T}_{I}^{(d)},$

行列 $D$ の$\overline{\Omega}^{(d)}$$\mathcal{T}_{B}.(d)$,$\Lambda_{I}(e)\cap\Lambda_{I}(f)\emptyset=(.0\leq.

:$\Lambda^{(0)}arrow\Lambda$,

での自由度に対する行列

ここに,$D^{(d)}$)$=\Lambda_{I}0(d)$$\mathcal{T}_{I}^{()}d\Lambda_{I}^{(}$$D^{(d)}$,,$\Lambda_{B}=\bigcup_{\tau 0\leq d\leq}\Lambda_{B}^{(d})$$\mathcal{T}_{B}^{()}d\Lambda_{B}^{(})=A_{B}0(d)$,.

を次のように表す:$D^{(d)}=(_{D_{B}^{(}}D_{II,d}^{(d)})I$全体の行列 $D$ は部分領域での行列e<f\leq 7)$$\Gamma J_{B}^{(\dot{a}))}D_{IB}^{()}dB^{\cdot}$の重ね合わせによりに表現できる.

$D=$

.

$T$

$D_{BB}= \sum_{0\leq d\leq 7}\tau_{B}^{(}d)D(d)\mathcal{T}^{()}BBBd$である.

剛性行列 $A$ と発散に対応する行列 $B$ に対しては, 成分が 1, $-1,0$ からなる行列$M_{A}^{(d)},$$M_{B}^{(d)}$をとなることに注意する必要がある.
. で行列を作成するのこれらの考察より, 領域全体の行列を作成する必要は無く, 部分領域みで良いことが分かる.

用いて$A^{(d)}=M_{A}(d)A(0),$$B^{(d)}=M_{B}(d)B(0),$$(0\leq d\leq 7)$$\Omega^{(0)}$アルゴリズム (a) 前処理行列に不完全 Cholesky 分解による行列を採用する.

$C=P_{N^{\perp}}$ .
アルゴリズム (b) 前処理行列に部分領域毎の不完全 Cholesky 分解による行列を採用する.

$C=P_{N^{\perp}}$ .

112

ここに,$\hat{D}^{-1}=$であり,$\Omega^{(d)}$内の節点に対しては Cholesky 分解による行列, 部分領域問境界上の節点に対しては対角スケーリングを採用する.

$\hat{A}^{-1}$も同様に定義される.

5. 数値実験を次のように与えた.の球殻領域において,$0\leq\varphi<2\pi)$ を極座標とするとき$R_{1}=11/9,$ $R_{2}=22/9$$R_{2}$,$0\leq\theta\leq\pi,\mu$$f$$(r, \theta, \varphi),$$f(r, \theta, \varphi)=e^{(r})(\frac{1}{r}\frac{R_{2}R_{1}}{R_{2}-R_{1}}-\frac{R_{1}}{R_{2}-R_{1}}+\epsilon\sin\pi(\frac{R_{2}-r}{R_{2}-R_{1}})Y^{2}3(\theta, \varphi))$

ここで,$e^{(r)}$は半径方向の単位ベクトル$e^{(r)}(x)= \frac{x}{|x|},$$Y_{3}^{2}$,$(R_{1}\leq r\leq$$\epsilon=0.1$.

は正規化された 3 次 2 階球面調和関数である. 図 1 に球殻領域の四面体による要素分割 (メッシュ B) を示す.

での計
表 1 に離散化条件を示す. メッシュ A での計算コストを表 2, より大規模なメッシュでの残差収束状況を図
算コストを表 3 に示す. 残差収束判定は $\in=10^{-}10$ とした. メッシュ
2 に示す. 横軸は経過時間 (秒) である. 表 2, 表 3 と異なり, 前処理行列の作成時間を含んでいる.
共有メモリー型の並列計算機 Origin 2000(MIPS R10000@250MHz,バイトキャッシュ) をにより並列処理を記述した.
法のアルゴリズム中, 内積計算を行う
用い, POSIXが, 処理を複数プロセッサー (スレッド) に分割して行う場合, 同期処理が必要となる. 同期処理のため, barrier-wait 関数 [51 を用いたが, スレッドの同期の際のオーバーヘッドが大きい. このため, 並列化が容易なアルゴリズムにも関わらず, 小規模な問題 (メッシュ A) では高い並列効率は得られなかった. プログラムのチューニングが必要である.
アルゴリズムは領域の対称性を利用し, メモリ一の節約を図っているため a と比較しメモ$\mathrm{B}$$\mathrm{B}$$4\mathrm{M}$$\mathrm{T}\mathrm{h}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{d}\mathrm{s}[5]$$\mathrm{C}\mathrm{G}$$\mathrm{b}$リー使用量が約 1/2 になっていることが分かる.
謝辞 : 科学技術庁高度計算科学技術共同推進制度 「地球深部ダイナミクスの数値シミュレーション」 の援助により, 防災科学技術研究所所有の並列計算機を使用した. ここに, 謝意を表する.
参考文献[1] L. P. Franca, S. L. Frey, and T. J. R. Hughes. Stabilized finite element methods: I.
Application to the advective-diffusive model. Computer Methods in Applied Mechanics
and Engineering, 95:253-276, 1992.

113[2] E. F. Kassechieter. Preconditioned conjugate gradients for solving singular systems.
J. Comput. Appl. Math., 24:265-275, 1988

[3] M. Tabata. Finite element approximation to infinite Prandtl number Boussinesq equations
with slip boundary condition. Computational Fluid Dynamics ’98, 22-27, John Wiley&Sons, Chichester, 1998.
[4] 鈴木厚. 有限要素離散化 Stokes 方程式に対する反復解法. 京都大学数理解析研究所講究録1084 (数値計算における前処理の研究): 103-110, 1999 年 2 月.

[5] D. R. Butenhof. Programming with POSIX Threads. Addison-Wesley, Massachusetts,

1997.
表 1. 離散化パラメータ

節点数

メッシュ

A25,538
225,838$\mathrm{B}$アルゴリズムa$\mathrm{b}$$\mathrm{b}$$\mathrm{b}$$\mathrm{b}$$\mathrm{P}\mathrm{E}$数

11
2
4
8

要素数

$h$

$N_{W}$$\delta$

141,168 0.30681 0.5 72,842
1,293,312 0.15466 0.5 657,542

表 2. メッシュ A での計算効率
計算時間 (秒) 速度倍率 反復回数395886
460.810310.800
215.757
249.143

残差$N_{M}$

25,538205,866$(\log_{1}\mathrm{o})$使用メモリ-10.2952$\overline{-250}$$(\mathrm{K}\text{バイト})1,4$

1

1.483
2.136
1.850

282
282
282281

表 3. メッシュでの計算効率計算時間 (秒) 速度倍率 反復回数

-10.1539-10.1090-10.1081-10.1105

82,94484,06486,33690.976$\mathrm{B}$アルゴリズム$\mathrm{P}\mathrm{E}$数残差$(\log_{1}\mathrm{o})$使用メモリ(a$\mathrm{b}$$\mathrm{b}$$\mathrm{b}$$\mathrm{b}$

1
1
2
48

3,602.274
4,462.921
2,676.887
2,044.857
1,932.546$\mathrm{K}$バイト)

896,704451-10.1490
1458 -10.6498
444,304
1.667
449
446,736-10.0143
2.183
457 -10.2025
457,424
2.309
457 -10.2207
478,512
SGI Origin2000, MIPS R10000@250MHz
MIPSpro Compilers: Version 7.2.1–$\mathrm{C}$114

図 1. 球殻領域の四面体要素分割 (メッシュ B)\={o}$\vee\underline{\mathrm{o}\mathrm{o}})$$\overline{\varpi\supset}$$. \frac{\mathfrak{D}}{\mathrm{q})(D,\underline}$

図 2. 残走」又宋腹歴