時間 1, 2 次精度の圧力安定化・特性曲線法結合有限要素スキームによる 2, 3次元非定常 Navier-Stokes 方程式の数値計算
野津裕史*1 , 田端正久*2

1 はじめに
我々は非定常 Navier-Stokes 方程式のための時間 1 次精度圧力安定化・特性曲線法結合有限要素スキーム [6]を開発した. これに文献 [9] の結果を適用すれば, 時間 2 次精度圧力安定化・特性曲線法結合有限要素スキームを得る. スキームは P1/P1 要素を用いており, かつ連立 1 次方程式の係数行列が対称であるため, 大規模数値計算に有用であると考えられる. その有用性を, MPI 並列を用いた 2, 3 次元数値計算をおこなって実証するために本プロジェクトに応募した. 得られた研究成果の概要は以下のとおりである.
スキームに現れる連立一次方程式の係数行列は不定値対称行列である. 対称行列に対する代表的な反復解法
として, CG, CR, MINRES 法があり (例えば [1, 2, 4, 11] などを参照), 我々は以前 CG 法を用いていたが, CR,MINRES 法の方が有効であることを数値的に確認した. 線形解法を変更することで, より効率的な数値計算が可能となった.
線形解法ライブラリ Lis [3] の中には, MPI 並列の実装を容易にする便利な演算関数 (以後, Lis 並列演算関数と呼ぶ) が含まれている. MINRES 法の MPI 並列プログラムを Lis 並列演算関数を用いて実装した (本プログ
ラムは Lis に組み込まれた). 1 次精度スキームの 3 次元 MPI 並列プログラムを Lis 並列演算関数を用いて実装し, 約 130 万自由度の問題において 128 コアまでの MPI 並列計算の有効性が確認できた.
時間 1 次精度スキームによる 2, 3 次元キャビティ流れの数値計算結果は論文 [5] に掲載された. 計算機

FUJITSU PRIMEQUEST 580 はこの論文の結果を得るために非常に有用に利用したが, 計算結果の統一性の観点から, IBM eServer p5 モデル 595 で得られたものに絞って報告した. 実際, 本稿の 5 節を除く計算結果はFUJITSU PRIMEQUEST 580 によるものである.

2 圧力安定化・特性曲線法結合有限要素スキームΩ ⊂ Rd (d = 2, 3) を有界領域, Γ ≡ ∂ Ω を Ω の境界とする. T を正定数とする. 非定常 Navier-Stokes 方程式で支配される未知関数 (u, p) : Ω × (0, T ) → Rd × R を求める問題;

2
∂u


+ (u · ∇)u − ∇D(u) + ∇p = f ,∂tRe∇ · u = 0,u = g,
u = u0 ,(x,t) ∈ Ω × (0, T ),(x,t) ∈ Ω × (0, T ),(x,t) ∈ Γ × (0, T ),(1)x ∈ Ω , t = 0,を考える. ここで, u は流速, p は圧力, Re は Reynolds 数, f は外力, g は境界での流速, u0 は初期流速, D(u) は変形速度テンソル

Di j (u) ≡12µ∂ ui ∂ u j+
∂ x j ∂ xi¶
(i, j = 1, · · · , d)

*1 九州大学大学院数理学研究院, 現在, 産業技術総合研究所水素材料先端科学研究センター, h.notsu@aist.go.jp*2 九州大学大学院数理学研究院, tabata@math.kyushu-u.ac.jp

1

である.

2.1 有限要素スキーム
時間刻み ∆t と関数 u, w : Ω → Rd に対して, 関数 X1 (w, ∆t), X2 (u, w, ∆t) : Ω → Rd を

X1 (w, ∆t)(x) ≡ x − w(x)∆t,no ∆t
X2 (u, w, ∆t)(x) ≡ x − u(x) + w(x − w(x)∆t),2で定義する. 以下では, 記号 ◦ は関数の合成を表し, Ω 上の関数 ψ に対して ψ ◦ X1 (w, ∆t)(x) ≡ ψ (X1 (w, ∆t)(x))とする. t n と NT は, スキームによって異なる定義を与える. すなわち, 1, 2 次精度スキームに対して, それぞれ,1 次精度: t n ≡ n∆t,(∆t0 + (n − 1)∆tn2 次精度: t ≡0NT ≡ [T /∆t],(n ≥ 1)(n = 0)NT ≡ [(T − ∆t0 )/∆t] + 1,,とする. ここに, ∆t0 は 2 次精度スキームの第 1 ステップに用いる (十分小さい) 時間刻みである. Th ≡ {K} を領域 Ω の三角形 (四面体) 分割とする. 近似領域を Ωh とし, Γh ≡ ∂ Ωh とする. g ∈ C0 (Γ )d に対して, P1/P1 有限要素空間を©
ª
Xh ≡ vh ∈ C0 (Ω h )d ; vh |K ∈ P1 (K)d , ∀K ∈ Th ,©
ª
Mh ≡ qh ∈ C0 (Ω h ); qh |K ∈ P1 (K), ∀K ∈ Th ,©
ªVh (g) ≡ vh ∈ Xh ; vh (P) = g(P), ∀P : Γh 上の節点 ,©
ªQh ≡ qh ∈ Mh ; (qh , 1) = 0 ,により定義し, Vh ≡ Vh (0) とする. u, w, ζ ∈ H 1 (Ωh )d に対して Vh 上の一次形式 Mh1 (u, w; ∆t), Mh2 (u, ζ , w; ∆t)を⟨Mh1 (u, w; ∆t), vh ⟩ ≡⟨Mh2 (u, ζ , w; ∆t), vh ⟩ ≡³ u − w ◦ X (w, ∆t)
1´, vh ,∆t³ u − w ◦ X (ζ , w, ∆t)
2∆t´, vh ,とする. H 1 (Ωh )d × H 1 (Ωh )d , H 1 (Ωh )d × L2 (Ωh ), H 1 (Ωh ) × H 1 (Ωh ) 上の双一次形式 ah , bh , Ch をそれぞれ,¢2¡
D(u), D(v) ,Re¡¢bh (v, q) ≡ − ∇ · v, q ,¡¢Ch (p, q) ≡ − δ ∑ h2K ∇p, ∇q K ,ah (u, v) ≡K∈Thとする. ここで, δ は正定数, hK は要素 K の直径, (·, ·)K は L2 (K)d 内積である. w, v ∈ H 1 (Ωh )d , q ∈ H 1 (Ωh ) に対して, 記号 ãh , b̃h をそれぞれ,ãh (w, v) ≡´´ 2 ³
2³D(w) ◦ X1 (w, ∆t), D(v) + ∆t D(w)J(w)T , J(v) ,Re³´³ Re
´b̃h (v, q, w) ≡ − ∇ · v, q ◦ X1 (w, ∆t) − ∆t qJ(w)T , J(v) ,とする. ここに, 記号 J は Jacobi 行列,Ji j (v) ≡

∂ vi
∂xj

(i, j = 1, · · · , d),2である.

u0h は u0 の近似関数とする. (1) のための時間 1 次精度圧力安定化特性曲線有限要素スキーム [5, 6] は,(nnn⟨Mh1 (unh , un−1h ; ∆t), vh ⟩ + ah (uh , vh ) + bh (vh , ph ) = ( f , vh ), ∀vh ∈ Vh ,bh (unh , qh ) + Ch (pnh , qh ) = 0,

∀qh ∈ Qh ,(2)NTを満たす {(unh , pnh )}n=1⊂ Vh (g) × Qh を求めよ, となる. また, (2) に対応する, 時間 2 次精度圧力安定化特性曲線有限要素スキームの一般ステップ (n ≥ 2) は,

o 1no1nnn−1nn n−1n n n−1
a(u,v)+ã(u,v)+b(v,p)+b̃(v,p,u)⟨M(u,u,u;∆t),v⟩+hhhhhhhhh2hhhhhhhhh22´1³ nn−1n−1=f+f◦
X(u,∆t),v, ∀vh ∈ Vh ,
1hh

2
nnbh (uh , qh ) + Ch (ph , qh ) = 0, ∀qh ∈ Qh ,(3)である. スキーム (2), (3) は P1/P1 要素を用いており, かつ現れる連立 1 次方程式の係数行列が対称となるため, 大規模数値計算に有用である,
注意 1. (i) スキーム (3) は文献 [9] のスキームに圧力安定化項を付加したものであり, 同文献で提案されているものと同様の内部反復法により非線形スキームの解を求めることができる. その際に現れる行列は対称である.

(ii) 問題 (1) には圧力の初期値がないため, スキーム (3) を n = 1 として適用できない. 従って, 第１ステップの解 (u1h , p1h ) は他のスキームで求める必要がある. 我々は, 十分小さな時間刻みを用いた時間１次精度スキーム (2) により求める.

3 線形解法の比較
スキーム (2), (3) に現れる連立一次方程式の係数行列は不定値対称行列である. CG 法は正定値行列に対して有用とされるが,破綻が起きなければ不定値対称行列においても収束すること [4], およびほとんどすべての初期ベクトルに対してその破綻は起きないこと [10] が知られている. 実際, 我々が行った計算では破綻が起きたことはなく, 以前は CG 法を用いていた. しかしながら, CG 法で並列計算を用いて計算を行った結果,コア数によって反復回数および計算時間が変化することが確認された. これは, 丸め誤差のみの影響である. そこで不定値対称行列の例として Helmholtz 方程式および定常 Stokes 方程式をとりあげて反復法選択の再考を行い, CG 法よりも CR, MINRES 法が適していることを数値的に確認し, その成果を報告した [7, 8]. 本稿ではスキーム (2) に対する, CG, CR, MINRES 法の丸め誤差の影響の大きさについて述べる.

Ω = (0, 1)3 , Re = 1 とする. 領域 Ω の一辺の分割数を 32 として四面体分割を行い, スキーム (2) において, δ = 0.5, ∆t = 1/32 とする. 連立一次方程式 (総自由度) の数は n = 167, 244 である. 現れる係数行列を
TnnA ∈ Rn×nsym , x∗ ≡ (1, · · · , 1) ∈ R , b ≡ Ax∗ ∈ R として,Ax = b(4)を満たす x ∈ Rn を求める問題を考える.
問題 (4) を 3 つの反復法, CG, CR, MINRES 法で, 並列計算 (OpenMP) を用いて解いた. すべての解法において点 Jacobi 前処理 [1] を用いた. コア数は 1, 2, · · · , 8 とした. 反復法の初期値は 0 ∈ Rn とした. k 回反復後の残差ベクトルを rk とする. 図 1 はコア数と反復回数のグラフおよびコア数と計算時間の両対数グラフである. コ3

ア数によって, 各反復法の反復回数が変化しており, その影響は計算時間に現れている. 3 つの解法を比較すると, CG 法の反復回数の変化は大きく, CR, MINRES 法は小さい. 計算時間は CR 法が最も短い結果となった.

CR, MINRES 法は, CG 法にはない 残差が単調減少する 性質をもっているため, その違いが, 図 2 の収束履
歴で確認できる. なお, 3 つの解法によって求めた解はすべて (ほぼ) 同じである. これらの結果から, CG 法よりも CR, MINRES 法が有用であることがわかる.
10,0001,200

CGCRMINRESCGCRMINRES

1,000
800

9,000

Elapsed time (sec.)# of iterations

600

8,000

400

200
7,000

100

6,000
1234567182345678# of cores# of cores

102

102

100

100

10-210-210-410-4||rk|| / ||r0||||rk|| / ||r0||

図 1 各解法の，コア数と反復回数のグラフ (左図) とコア数と経過時間の両対数グラフ (右図).

-610CG-810-10CG

10-6-8MINRES CR

10-101010-1210-1210MINRES

CR

10-1410-140

2000

4000

6000

8000

10000

02000

k (# of iteration)

4000

6000

8000

10000

k (# of iteration)

図 2 収束履歴, 6 コア (左図), 8 コア (右図).

4 Lis 並列演算関数を用いた MPI 並列計算
約 130 万自由度の問題に対する Lis 並列演算関数を用いた MPI 並列の効果を述べる. ある 3 次元問題において, スキーム (2) を用い, 線形解法に MINRES 法を採用して MPI 並列計算を 1 ステップおこなった. 図 3 左右の図はそれぞれ, コア数に対する, 高速化率と経過時間の両対数グラフである. ともに 4 コア利用時 (2,175秒) を基準とした結果である. 128 コアまでの MPI 並列の効果が観察された.

4128

2000

Ideal

Ideal

641000Elapsed time (sec.)

Speed up ratio

321684500

250

1252162.548163264# of cores

128256512

48163264# of cores

128256512

図 3 コア数に対する, 高速化率 (左図) と経過時間 (右図) の両対数グラフ.

5 数値計算結果
本節では, スキーム (2), (3) の数値解の, 厳密解への数値的収束精度, およびスキーム (2) によるキャビティ流れ問題の数値計算結果について述べる. 領域は Ω = (0, 1)d (= Ωh ) とする. 一辺の分割数を NΩ とし, 代表長さを h ≡ 1/NΩ とする.

5.1 収束精度の確認(u, p) を問題 (1) の解, (uh , ph ) をスキーム (2) による有限要素解とする. 誤差としてErr ≡∥Πh u − uh ∥l 2 (H 1 (Ω )d ) + ∥Πh p − ph ∥l 2 (L2 (Ω )d )∥uh ∥l 2 (H 1 (Ω )d ) + ∥ph ∥l 2 (L2 (Ω )d )T
⊂ X に対してを用いる. ノルム空間 X に対して, ∥ · ∥l 2 (X) は離散 L2 (X) ノルムを表し, {ψ n }n=1N∥ψ ∥l 2 (X)vuNu T≡ t (t n − t n−1 )∥ψ n ∥2∑Xn=1である.
以下のテスト問題を設定する.
問題 1 (2 次元, 斉次 Dirichlet 境界条件). 問題 (1) において, d = 2, T = 1, Re = 1 とする. 厳密解は

 2
sin (π x1 ) sin(2π x2 )µ ¶u(x,t) = {1 + sin(π t)}  − sin2 (π x2 ) sin(2π x1 )
pcos(π x1 ) cos(π x2 )である.

5

問題 2 (3 次元, 斉次 Dirichlet 境界条件). 問題 (1) において, d = 3, T = 1, Re = 1 とする. 厳密解は¡¢
sin2 (π x1 ) sin(π x2 ) sin(π x3 ) sin π (x2 − x3 )
¡¢

µ ¶u
 sin(π x1 ) sin2 (π x2 ) sin(π x3 ) sin π (x3 − x1 ) (x,t) = {1 + sin(π t)} ¡¢p sin(π x1 ) sin(π x2 ) sin2 (π x3 ) sin π (x1 − x2 ) cos(π x1 ) cos(π x2 ) cos(π x3 )である.
問題 3 (3 次元, 非斉次 Dirichlet 境界条件). 問題 (1) において, d = 3, T = 1, Re = 1 とする. 厳密解は

sin(x1 + 2x2 + x3 + t) − sin(x1 + x2 + 2x3 + t)



µ ¶u
 − sin(2x1 + x2 + x3 + t) + sin(x1 + x2 + 2x3 + t)
(x,t) = p sin(2x1 + x2 + x3 + t) − sin(x1 + 2x2 + x3 + t) sin(x1 + x2 + x3 + t) − 8 sin3 (1/2) sin(t + 3/2)である.

1 次精度スキーム (2) を用いて, 問題 1, 3 を解いた. スキーム (2) において, ∆t = h とした. 問題 1 では
NΩ = 8, 16, 32, 64, 128, 問題 3 では NΩ = 4, 8, 16, 32, 64 とした. 図 4 は Err と ∆t の両対数グラフであり,左図が問題 1, 右図が問題 3 の結果である. 時間に関して 1 次精度であることが確認された.
同様に, 2 次精度スキーム (3) を用いて, 問題 1, 2 を解いた. スキーム (3) において, ∆t = ch1/2 とした. ここに c = 0.2 (d = 2), 0.1 (d = 3) とした. 問題 1 では NΩ = 8, 16, 32, 64, 128, 問題 2 では NΩ = 4, 8, 16, 32 とした. 図 5 は Err と ∆t の両対数グラフであり, 左図が問題 1, 右図が問題 2 の結果である. 2 次元において時間2 次精度であることが確認された. 3 次元では, 精度 2 より大きな勾配が得られたが, h が小さくなれば, 2 に近づいていくと考えている.

0.1

Err

Err

0.1

0.01
110.01110.01∆t0.0010.010.1∆t0.1

図 4 1 次精度スキーム (2) による Err と ∆t の両対数グラフ (左: 問題 1, 右: 問題 3).

5.2 合法キャビティ流れ問題次の 2, 3 次元合法キャビティ流れ問題を設定する.

61

Err

Err

0.1

0.1

220.01110.010.1∆t0.010.010.1∆t

図 5 2 次精度スキーム (3) による Err と ∆t の両対数グラフ (左: 問題 1, 右: 問題 2).

問題 4 (2 次元, Re = 5, 000). (1) において Re = 5, 000, 境界流速は,(g1 (x,t) = g1 (x) =4x1 (1 − x1 ) (x2 = 1)
(その他)

0,

g2 = 0,(5)とし (図 6 左図および中央図), 初期流速は定常ストークス方程式の解とする.
問題 5 (3 次元, Re = 1, 000). (1) において Re = 1, 000, 境界流速は,(g1 (x,t) = g1 (x) =16x1 (1 − x1 )x2 (1 − x2 ) (x3 = 1)
(その他)

0,

g2 = g3 = 0,(6)とし (図 7 左図および中央図), 初期流速は定常ストークス方程式の解とする.
スキーム (2) を用いて問題 4, 5 を解く. 問題 4 に図 6 右図, 問題 5 に図 7 右図の非一様なメッシュを用いる.
総自由度はそれぞれ 34, 410, 298, 508 である. 問題 4 では ∆t = 1/32, δ = 0.2, 問題 5 では ∆t = 1/24, δ = 50として, 有限要素解が数値的に定常状態となるまで計算を行った. 図 8 は問題 4 の定常解の流線図および圧力等高線, 図 9 は問題 5 の定常解の流速ベクトルの各平面への射影と各平面における圧力等高線であり , 流れの特徴を捉えた解が得られている. スキーム (2) がこの問題に有効であることを示す結果といえる.

6 結び
本プロジェクト遂行中に, 適切な線形解法の選択および Lis 並列演算関数を用いた MPI 並列計算をおこなった. CG 法よりも CR, MINRES 法がスキームに現れる不定値対称行列に有用であることがわかった. MPI 並列による高速化は 128 コアまでの有効性が確認できた. それより多いコア数での高速化は今後の課題である. ス
キーム (2), (3) の数値解の, 厳密解への収束精度がそれぞれ 1, 2 次であることを確認した. 2, 3 次元合法キャビ
ティ流れ問題にスキーム (2) を適用し, 流れの特徴を捉える有用なスキームであることを観察した. P1/P1 要素を用い, かつ現れる行列は対称な本スキームが, 大規模数値計算に適した有力な数値解法として期待できる結果を得た.

7

g1 ( x1 ,1)x2u = ( g1 , 0)

1u=0

ΩOu=0T1u=01x1

Ox1

10.5

図 6 左図: 2 次元合法キャビティ流れ問題, 中央図: g1 (·, 1) のグラフ, 右図: メッシュ (NΩ = 256).

u = ( g1 , 0, 0)Tx3

g1 ( x1 , x2 ,1)

1(0.5, 0.5,1)

x2

Ω

x2

1

O11x1

O1x1

図 7 左図: 3 次元合法キャビティ流れ問題, 中央図: g1 (·, ·, 1) のグラフ, 右図: メッシュ (NΩ = 48).

図 8 流速ベクトルの各平面への射影 (上段) と各平面における圧力等高線 (下段), Re = 5, 000.

8

(a) x2=0.5x3x1(b) x1=0.5x3x2(c) x3=0.5x2x1

図 9 流速ベクトルの各平面への射影 (左側) と各平面における圧力等高線 (右側), Re = 1, 000.

9謝辞
この研究中に, 九州大学情報基盤研究開発センターの藤野清次教授には線形解法の選択について, 西田晃博士には MPI 並列プログラム作成に用いた Lis について, 有益な助言を頂いたことに感謝する. 本研究は, 日本学術振興会科学研究費補助金基盤研究 (S), No.16104001, から支援を受けた.

参考文献[1] R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo, C. Romine and H.
van der Vorst, Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, SIAM,Philadelphia, 1994.
[2] A. Greenbaum, Iterative Methods for Solving Linear Systems, SIAM, Philadelphia, 1997.
[3] Lis, http://www.ssisc.org/lis/.
[4] 森正武, 杉原正顕, 室田一雄, 線形計算, 岩波, 東京, 1994.
[5] H. Notsu, Numerical computations of cavity flow problems by a pressure stabilized characteristic-curve
finite element scheme, Transactions of Japan Society for Computational Engineering and Science, Vol.2008(2008), No.20080032, available at http://www.jstage.jst.go.jp/article/jsces/2008/0/20080032/ pdf.
[6] 野津裕史, 田端正久, Navier-Stokes 方程式のための圧力安定化・特性曲線法結合有限要素スキーム, 日本
応用数理学会論文誌, Vol.18, No.3 (2008), pp.427–445.

[7] H. Notsu and M. Tabata, On the influence of rounding errors of the CG method in indefinite problems,Proceedings of International Kyoto-Forum on Krylov Subspace method (2008), pp.59–64.
[8] 野津裕史, 田端正久, ヘルムホルツ・ストークス方程式に対する反復法の丸め誤差の影響について, 日本数学会秋季総合分科会 応用数学分科会講演アブストラクト (2008), pp.163–166.

[9] H. Notsu and M. Tabata, A single-step characteristic-curve finite element scheme of second order in time forthe incompressible Navier-Stokes equations, Journal of Scientific Computing, Vol.38, No.1 (2009), pp.1–14.
[10] 鈴木厚, 田端正久, 不定値対称行列に対する共役勾配法の収束について, 京都大学数理解析研究所講究録,
Vol.1265 (2002), pp.39–44.
[11] H. A. van der Vorst, Iterative Krylov Methods for Large Linear Systems, Cambridge University Press,
Cambridge, 2003.

10