2018 年 8 月 27 日 (Ver.2)

第 103 回お試しアカウント付き
並列プログラミング講習会「OpenFOAM 初級」

今野 雅(東京大学情報基盤センター客員研究員，株式会社 OCAEL)目次
1

Ahmed body とは

5

2

解析モデル

8

3

基本的な解析条件

9

4

CFD 解析条件

10

5

解析領域

11

6

解析ケースの作成

12

7

解析手順

13

8

前処理
前処理ジョブの投入 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

14
16

9

mpi_core_setting.sh による MPI プロセスのピニング

1810

特徴辺抽出 (surfaceFeatureExtract)

19

8.1

11
ベース格子生成 (blockMesh)
11.1 ベース格子での解析領域の指定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
11.2 ベース格子での分割数・パッチの指定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

20
20
21

12

並列格子生成用領域分割 (decomposePar)2213
並列格子生成 (snappyHexMesh)
13.1 三角分割表面形状の定義 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13.2 直方体形状の定義 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13.3 特徴辺の設定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13.4 階段状格子の設定パラメータ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13.5 形状表面の細分割 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13.6 領域の細分割 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13.7 解析領域の内部点指定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13.8 境界適合 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13.9 境界層レイヤーの挿入 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13.10 前処理のログの確認 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
14

ファイルの転送34

15格子の可視化15.1 最終生成格子の可視化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
15.2 特徴辺の可視化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
16

流体解析用の領域分割

流体解析

35
35
37
38

17
初期値作成
17.1 初期値作成ジョブの投入 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
18

23
23
24
25
25
26
27
28
28
28
31

2

39
40
41

19

関数出力のモニター

42

20
20.1
20.2
20.3
20.4
20.5
20.6
20.7
20.8

流体解析の設定
流体物性値の設定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
乱流モデルの設定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
離散化スキームの設定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
解法の設定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
実行制御の設定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
計算履歴の出力 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
境界条件の設定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
解析結果の再構築 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

44
44
44
45
46
49
49
51
52

21

空力係数の実測値との比較

53

22

後処理

54

23

速度分布のプロットと実験値との比較

57

24

ParaView による計算結果の確認

59

25
流体解析の性能プロファイル実行
25.1 流体解析の性能プロファイル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
25.2 プロファイラ結果の表示 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
26

並列計算ベンチマークテスト

61
61
62
64

3

格子生成・空力解析演習
格子生成・空力解析演習手順 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
格子生成・空力解析演習設定変更例 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
格子生成・空力解析演習実行例 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

66
66
67
67

28
補足
28.1 Intel MPI ライブラリでの実行エラー . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
28.2 Oakforest-PACS での OpenFOAM のコンパイル . . . . . . . . . . . . . . . . . . . . . . . . . .

68
68
68

27
27.1
27.2
27.3

4

1 Ahmed body とは

図 1 Ahmed body 後流の特性 (左図：傾斜角度 θ = 25◦ ，右図：傾斜角度 θ = 35◦ )[4]• Ahmed body は，1984 年に当時のドイツ航空宇宙研究所 (DFVLR)，現在のドイツ航空宇宙センター (DLR) の Syed R. Ahmed が SAE(Society of Automotive Engineers) に投稿した論文 [2] において風洞実験を行なっている自動車モデルの略称である．
• このモデルは図 1 に示す通り，非常に単純化されたモデルであるため，風洞実験や CFD 解析での再現が容易である．
• 一方で，車体後部上面の傾斜角度が θ = 25◦ では，車体後部上面の両端に後引き渦が発生しているなど，複雑な流れ場となっている．

5

• 図 2 に示す通り，抗力係数 Cd 値が車体後部上面 (スラント) の傾斜角度 θ が 25◦ ∼ 30◦の範囲で急増し，θ = 30◦ を境に急減する特徴を持つことから，スラント角によって引き起こされる車体後方領域の流れを調査するのによく用いられた．

図 2 傾斜角度による空気抗力係数の特性 [2]6図 3 LSTM 低風速風洞でのAhmed body の実験)[3]

図 4 オイルと煤による Ahmed body 後部表面の流脈可視化 (左図：傾斜角度 θ = 25◦ ，右図：傾斜角度 θ = 35◦ )[3]• ヨーロッパの Models for Vehicle Aerodynamics(MOVA) プロジェクトとして，ドイツ，
エアランゲン大学の流体力学研究室 (通称 LSTM) の Lienhart ら [3] により，図 3 に示す
LSTM の低風速風洞において，スラント角が 25o，35o の Ahmed body について，2 成分
レーザードップラー風速計 (LDA) を用いた詳細な流速の計測や，図 4 に示すスラント表面の流脈可視等が行われた．
• これらの実験結果および Ahmed body の形状データは，ERCOFTAC(European ResearchCommunity On Flow, Turbulence And Combustion) の Classic ケース 8.2 として Web[6] で掲載されており，簡単に詳細なデータ入手できることから，自動車空力解析用のベンチ
マークモデルとして良く用いられる．

7

2 解析モデル• 車体の形状は車長は D = 1.044 m ，車幅は W = 0.389 m ，車高は h = 0.288 m である．
• 本演習では，車体後部上面の傾斜角度 (以下，スラント角) は ϕ = 25◦ とする．

図 5 Ahmed Body モデル [5]83 基本的な解析条件• RANS 乱流モデルを用いた非圧縮性流体の定常解析を行う．
• 基本的な解析条件は，ERCOFTAC Classic ケース 8.2 の WEB ページ [6] に記載されている
表 1 の条件を用いる．
表1

風洞高さ
空気の動粘性係数
流入速度
車高
レイノルズ数 (車高べース)

解析条件

H = 1.4 [m]ν = 15 × 10−6 [m2/s ]U = 40 [m/s] (風洞実験でのバルク風速)
h = 0.288 [m]Re = U h/ν = 768, 000 [-]9

4 CFD 解析条件
表 2 CFD 解析条件
上面・側面境界

すべり壁 (slip wall)

壁面境界

すべり無し壁 (no-slip wall)，滑面に対する壁関数

乱流モデル

標準 k − ω SST モデル

流入境界乱れ強さ

I = 0.25[%]l = 0.07H = 0.098[m]
全ての場について TVD(limited linear)
simpleFoam(定常乱流解析ソルバ)
SIMPLEC 法

流入境界乱流長さスケール
移流項の離散化スキーム
解析ソルバ
圧力・速度連成手法• 風洞実験では天井面と側面は開放されていたが，CFD 解析ではすべり壁とする．
• 床面と車体表面はすべり無し壁とし，壁関数を設定する．
• 乱流モデルは k - ω SST モデルを使用する．
• 風洞実験での測定を行なったドイツ・バイエルン州・エアランゲン大学の LSTM 低風速風
洞の流入風は，平均乱流強度が 0.25% 以下と記載されているので，乱流強度 I は 0.25%．
• 乱流長さスケール l は，風洞高さを H=1.4[m] として，l = 0.07 × H = 0.098[m]．
105 解析領域

水平面

鉛直面
図 6 解析領域の格子分割と Ahmedbody モデル• 本解析では，車体が対称面を持つため，解析領域も空間対称に設定することが可能である．
• また，RANS 乱流モデルを用いた定常解析を行うので，車体の対称面に対して，流れ場の定常解も対称であると仮定し，対称面を設けた半分の領域の解析を行う．
• 解析領域は 15[m] ×0.935[m] ×1.4[m] とする．
• 計算格子はベース格子の分割幅 ∆x が概ね 0.1[m] となるように分割数を指定する．
• そのベースメッシュを基準とし，車体周辺の流れを捉えるための再分割領域の設定と，形状を表現するための細分割の設定を行う．

11

6 解析ケースの作成• OpenFOAM では，基本的に付属しているチュートリアルモデルの中から解きたい解析モ
デルをコピーして編集するのが一般的な解析ケースの作成方法である．
• 本解析は非圧縮性流体解析であるから，カテゴリーは incompressible であり，また定常乱流解析のため，ソルバーは simpleFoam となる．
• 対応するチュートリアルのディレクトリから，本解析に最も近いケースを探すと，バイクの空力解析 motorBike および車の 2 次元空力解析 simpleCar となるので，これらを参考に設定を作成する方法がある．
• 今回は，演習の時間が限られているので，既に作成済みの解析ケースをコピーする．

Oakforest-PACS にログインして，以下のコマンドを実行し，講習会用ファイルを展開する．
なお，講習会用ディレクトリの参照を容易にするためにシンボリックリンクを貼っている．
計算用ディレクトリへの移動と講習会用ファイルの展開
1 cd /work/gt00/$USER
2 tar xf /work/gt00/share/lecture20180828.tar.xz
3 ln -s /work/gt00/$USER/lecture ~/作成済の解析ケースのコピー
1 cd ~/lecture/Ahmed25
2 cp -a template case-1
3 cd case-1

12

7 解析手順1. 前処理 (スクリプト: pre.sh)（a）特徴辺抽出 (surfaceFeatureExtract)（b）ベース格子生成 (blockMesh)（c）並列格子生成用領域分割 (decomposePar -decomposeParDict 辞書名)（d）並列格子生成 (snappyHexMesh -decomposeParDict 辞書名 -overwrite -parallel)（e）格子の再構築 (注) (reconstructParMesh -constant)（f）格子の品質チェック (checkMesh -constant)
2. 初期値作成 (スクリプト: solveInit.sh)（a）格子順変更による行列バンド幅縮小 (renumberMesh -overwrite -parallel)（b）ポテンシャル流れを解いて初期値作成 (potentialFoam -parallel)
3. 流体解析用の領域分割 (スクリプト: decomposrPar.sh)（a）流体解析用の領域分割 (注) (decomposePar -force)
4. 流体解析 (スクリプト: solve.sh．性能プロファイル実行時は solveVTune.sh)（a）ソルバ実行 (simpleFoam -parallel)
5. 後処理 (スクリプト: post.sh)（a）解析結果のサンプリングなど (postProcess -func sample -parallel など)
6.結果プロット (スクリプト: plot.sh)

注) 並列格子生成による格子数が非常に多い場合等，格子の再構築と本解析用の領域再分割も行わない手順もあるが，本手順では格子生成と流体解析の領域分割が独立となり自由度が高く，計算効率も上げやすい．


13



8 前処理Code 1 pre.sh(前半)
1 #!/bin/bash
2 #PJM -S
3 #PJM -g gt00
4 #PJM -L elapse=0:15:00
5 #PJM -L rscgrp=tutorial-flat
6 #PJM -L node=1
7 #PJM --mpi proc=64
8 source share.sh # 共通の設定

ジョブ投入コマンド pjsub に与えるオプションを#PJM 以降に記述できる．また，これらの指定は，ジョブ投入コマンド pjsub のオプションにより上書きされる．

• -L rscgrp=tutorial-flat : リソースグループ=tutorial-flat(講習会の時間のみ有効)
講習会終了後はリソースグループを lecture-flat に書き変えるか，pjsub のオプションに-L rscgrp=lecture-flat を指定する．
• -L node=1 : ノード数=1• –mpi proc=64 : 総プロセス数=64• -L elapse=0:15:00 : 経過時間制限=15 分 (講習用リソースグループの最大値)• -g gt00 : 課題番号=gt00(講習用)• -S : ジョブ統計情報とノードごとの詳細情報をファイルに出力
14

Code 2 share.sh
1 #!/bin/sh
2 module purge # 標準で有効なmoduleに依存しないよう全moduleをunload
3 module load intel/2018.1.163 # intelコンパイラ環境をload
4 export HOME=/work/$(id -gn)/$USER # $HOMEを参照するOFのスクリプト用に設定
5 # OpenFOAM-v1712環境設定
6 # KNL版はログインノードやプリポストノードで動作せず，通常版に比べての高速化は僅か
7 source /work/gt00/share/OpenFOAM/OpenFOAM-v1712/etc/bashrc \8 WM_COMPILER=Icc2018_1_163 WM_MPLIB=INTELMPI2018_1_163 # 通常版
9 #source /work/gt00/share/OpenFOAM/OpenFOAM-v1712/etc/bashrc \10 #WM_COMPILER=Icc2018_1_163KNL WM_MPLIB=INTELMPI2018_1_163 # KNL版
11 source /usr/local/bin/mpi_core_setting.sh # MPIプロセスのピニング
12 export I_MPI_DEBUG=5 # MPIプロセスのピニング等のデバッグ情報取得
13 env # 念のため実行時の環境変数を記録しておく
14 numa="numactl␣-p␣1" # numactl -p 1 : MCDRAMをできるだけ使用
15 # numactl -p 1付きでのコマンド実行が失敗する場合はCacheモードなので，以降numactlを実行しない
16 $numa numactl -s 2> /dev/null || numa=""
17 mpirun="mpiexec.hydra␣-n␣$PJM_MPI_PROC" # MPI実行用コマンド($PJM_MPI_PROC=プロセス数)
18 log=$PJM_JOBNAME.l$PJM_JOBID # ログファイル名のベース=ジョブ名.lジョブID(任意)リソースグループが*-cache の時は Cache モードとなり，高バンド幅メ
モリの MCDRAM はキャッシュメモリとして利用されるので，特に設定の必要はないが，リソースグループが*-flat の時は Flat モードとなり，

MCDRAM を優先的に使うには，numactl などで設定する. 図引用元 [7]15Code 3 pre.sh(後半)
9 $numa surfaceFeatureExtract >& $log.surfaceFeatureExtract # 特徴辺抽出
10 $numa blockMesh >& $log.blockMesh # ベース格子生成
11 $numa decomposePar -decomposeParDict system/decomposeParDict.mesh \>& $log.decomposePar.mesh # 格子生成用領域分割．-decomposeParDictで領域分割の辞書指定
13 $mpirun $numa snappyHexMesh -decomposeParDict system/decomposeParDict.mesh \14-overwrite -parallel >& $log.snappyHexMesh # 並列格子生成
15 $numa reconstructParMesh -constant >& $log.reconstructParMesh # 格子の再構築
16 $numa checkMesh -constant >& $log.checkMesh # 格子の品質チェック
12

8.1 前処理ジョブの投入
ジョブスクリプト pre.sh を pjsub で投入する．
ジョブの投入1 pjsub pre.sh # 講習会終了後は pjsub -L rscgrp=lecture-flat pre.sh

以下のコマンドにより，投入したジョブの状況を確認する．
ジョブの状況確認
1 pjstat

本演習において，今後 pjsub でジョブを投入したら，pjstat で適宜ジョブの状態を確認する．

16

ジョブが実行状態になったら，tail コマンドを用いて，ログのトレースを行う．
ログファイルのトレース
1 tail -f *.sh.l*

ユーティリティ毎にログファイルが異なるので，End などのメッセージが出るなどしてログファイルの更新が止まったら，Ctrl+C(コントロールキーと C キー) を押して tail コマンドを一旦終了させて，カーソル↑ + リターン，または，!!で再度上記コマンドを実行する．




Tips

• 特徴辺抽出 (surfaceFeatureExtract)，ベース格子生成 (blockMesh)，領域分割 (decomposePar)，格子再構築 (reconstructParMesh) などの並列処理できないユーティリティ
は，後述するような複数ノードを使用するジョブ内ではなく，pre.sh のように 1 ノードで実行するジョブで実行させるほうがトークンを節約できる．
• 空いている隙間リソースに割りあてるバックフィルスケジューリング機能により実行開始が早まる可能性が高くなるので，経過時間を見積って，できるだけ適切な経過時間制
限値を-L elapse で設定する．




17

9 mpi_core_setting.sh による MPI プロセスのピニング

ピニングを適切に行わないと，プロセスが偏って割り当てられ，性能が劣化する可能性が高い．最初は必ずピニング結果を確認する

先 頭 の 物 理 コ ア は，CPU ス ケ ー ジ ュ ー リ
ングクロック割り込
み抑止設定 (動的 Tickless) がされていないので，なるべく先頭タイ
ルは使用しない

1 ノード内プロセス数が 1〜33
L2 キャシュを 1 プロセスが専有

1 ノード内プロセス数が 34〜66
L1 キャシュを 1 プロセスが専有

18




一般的に，ノード当り
66 プロセス以下の L1
専有より，ノード当り
33 プロセス以下の L2専有のほうが，メモリ帯域幅も多く使用でき，ピーク性能が高くなるが，L1 専有のほうが費用 (トークン) 対性
能が高い事も多いので，事前のベンチマー
クテストが重要図引用元 [7]10 特徴辺抽出 (surfaceFeatureExtract)• 格子生成において特徴辺を再現するため，形状データから特徴辺を抽出する．
• 特徴辺の抽出する形状を system/surfaceFeatureExtractDict で定義する．
• 辺に隣接する 2 つの三角形の法線ベクトルの成す角度が includedAngle 以下であれば特徴辺として抽出される．
• ここでは，車体の先端の特徴辺を取り込み，かつ平面内の分割線を特徴辺として認識しな
いよう includedAngle を 179.9[o] としている．
Code 4 system/surfaceFeatureExtractDict
17 AhmedBody.stl // 三角分割表面形状ファイル名
18 {
19
20

// 特徴辺の抽出法
extractionMethod

extractFromSurface; // 表面形状ファイルから抽出

21
22
23
24
25
26

extractFromSurfaceCoeffs{// 隣接する2つの面の法線の角度が以下より小さい辺を特徴辺とする
includedAngle
179.9;}

27
28
29// 可視化用に特徴辺をobj形式で出力
writeObjyes;

1911 ベース格子生成 (blockMesh)
11.1 ベース格子での解析領域の指定• ベース格子の作成には blockMesh を用いる．
• 計算格子はベースメッシュの幅が約 ∆x = 0.1[m] となるように分割数を指定し，各境界面にもパッチ名を付与する．
• 車体後部が x = 0 ，車体の対称面が y = 0 ，床面が z = 0 となるように解析領域 ( 15[m]×0.935[m] ×1.4[m] ) を指定する．
Code 5 system/blockMeshDict
20 xmin -5;
21 xmax 10;
22 ymin -0.935;
23 ymax
24 zmin25 zmax

0;
0;
1.4;

26
27 vertices
28 (
29
30

( $xmin $ymin $zmin )( $xmax $ymin $zmin )20

11.2 ベース格子での分割数・パッチの指定
ベース格子の分割幅 ∆x が，約 0.1[m] となるように分割数を指定する．
Code 6 system/blockMeshDict
39 blocks
40 (41hex ( 0 1 2 3 4 5 6 7)(150 10 14) simpleGrading (1 1 1)

42 );

解析領域の各面に対して，以下のようにパッチを指定する (設定は省略)．
表 3 解析領域の境界面
境界条件

パッチ名

パッチ型

対称面

centerside
inlet
outletlowerWallupperWallsymmetry
wall
patch
patch
wall
wall

側壁
流入面
流出面
床面
天井

21

12 並列格子生成用領域分割 (decomposePar)• 生成する格子が複雑になるに従って，snappyHexMesh の計算時間が長くなる．
• 特に Xeon Phi ではコアの性能が低く時間がかかるので，できれば並列計算する．
• 今回格子生成用の領域分割の設定は，本解析用の領域再分割とは別に行うため，通常のsystem/decomposeParDict ではなく，system/decomposeParDict.mesh で行う．
• numberOfSubdomains で領域分割数 (並列数)，method で分割手法を指定する．
• ここでは，1 ノードあたり 68 物理コアのうち 64 コアを用いて並列計算を行う (2 の階乗のほうが simple や hierarchical での分割数指定が容易)．
Code 7 system/decomposeParDict.mesh
18 numberOfSubdomains 64; // 領域分割数(並列数)1920 // method: 分割手法
21 methodhierarchical;

22
23 hierarchicalCoeffs
24 {
25
26
27

n
delta
order

(4 4 4); // x,y,z方向の分割数
0.001;
xyz; // 分割順序

28 }

22

13 並列格子生成 (snappyHexMesh)blockMesh で解析領域を作成したので，次は snappyHexMesh で車両周りのメッシュ作成，
細分割領域の作成，境界層レイヤーの挿入を行う．ここでは，snappyHexMesh を用いて，細かなメッシュ操作について説明する．

13.1 三角分割表面形状の定義
snappyHexMesh の設定ファイル snappyHexMeshDict の geometry で形状の指定と細分割領域を定義する．基本的な記述方法としては，形状の種類，形状の識別名，形状の領域を指定する．形状の識別名は細分割や境界層の挿入のメッシュ操作の際に用いられる．
Code 8 system/snappyHexMeshDict
28 geometry
29 {
30
31
32
33
34
35

// CADデータのファイル名(constant/triSurfaceディレクトリ内)
AhmedBody.stl{
type triSurfaceMesh; // 形状の種類
name AhmedBody;// 形状の識別名}

23

13.2 直方体形状の定義
本解析モデルでは車両周りの流れを捉えるために 4 段階の細分割領域を定義した．車両に近づくにつれて細かくなるように領域 refinebox1-refibebox4 を定義した．
Code 9 system/snappyHexMeshDict
37
38
39
40
41
42

refinebox1 // 形状の識別名{
type searchableBox;// 形状の種類
min ( -2.610 -0.935 0.0); // 直方体の座標最小値
max ( 2.610 0.935 0.845 ); // 直方体の座標最大値}

図 7 領域の細分割

2413.3 特徴辺の設定
後述する surfaceFeatureExtract で抽出される特徴辺に対する分割レベルを設定する．
Code 10 system/snappyHexMeshDict
81

features({file "AhmedBody.eMesh";level 4;
}

82
83
84
85
86

// 特徴線ファイル// 分割レベル

13.4 階段状格子の設定パラメータ
階段状格子を作成するために必要なメッシュ数の制限や格子間のレイヤ数に関するパラメー
タを設定する．指定した細分割レベルまで細分割が進むように，maxLocalCells や
maxGlobalCells は大きくしておく．
Code 11 system/snappyHexMeshDict
67 castellatedMeshControls
68 {
69
70
71

maxLocalCells 100000000; // プロセッサあたりの最大メッシュ数
maxGlobalCells 200000000; // 全体の最大メッシュ数
minRefinementCells 10;// 細分割される最小メッシュ数

25

72
73
74

nCellsBetweenLevels 3;// 細分割レベルが異なる格子間のレイヤ数
maxLoadUnbalance 1;// 並列時ロードバランシングを行なわない割合// 本解析用の領域分割を後で行うので，ここではロードバランシングを行なわない

13.5 形状表面の細分割
refinementSurfaces で geometry で定義した形状表面に細かいメッシュを作成する．分割
レベルを level で指定し，ベースメッシュの分割レベルを 0 として、1 上がるごとに 格子の辺を 2 分割する細分割を行う．分割レベルは最小値と最大値を指定するが，なるべく同一にするほうが望ましい．形状再現性や格子分解能を考慮して，必要に応じて分割レベルを調整する．また，作成されるパッチ面の型定義やグループ化には patchInfo を用いる．
Code 12 system/snappyHexMeshDict
98
99
100
101
102

refinementSurfaces{
AhmedBody // geometryブロックで定義した形状の識別名{
level (4 4); // 形状表面の細分割レベル (最小値 最大値)

103
104
105
106
107

patchInfo{type wall; // パッチの境界条件
inGroups (AhmedBodyGroup); // グループ名

26

108

}

109
110

}
}

13.6 領域の細分割
refinementRegions で gemetry で定義した形状の領域内 (外も可能) のメッシュを細分化を行う．細分割する領域を mode で設定するが，指定した形状の内部 (inside)、外部(outside)、距離 (distance) を選択できる．また，levels に対して ((距離 分割レベル)) を指定する．ここで，距離は mode が distance 以外の時は，特に意味を持たない．
Code 13 system/snappyHexMeshDict
129
130
131
132
133
134
135
136

refinementRegions{
refinebox1{
mode inside; // 領域内が分割されるlevels ((1E15 1)); // 2番目の数字が細分割レベル// 最初の数字はmodeがdistanceの場合のみ，1番目が距離を表わす．
}

137
138

refinebox2

27

13.7 解析領域の内部点指定
locationInMesh で解析領域の内部点を指定する．この点により定義した形状の内外判定が行われる．この点は，格子生成の途中の段階でも格子の界面や節点に一致してはいけない．
Code 14 system/snappyHexMeshDict
166

locationInMesh ( -1e-10 -1e-10 1e-10 );13.8 境界適合
snapControls で境界適合過程の制御パラメータを指定する．今回はデフォルトのままである．

13.9 境界層レイヤーの挿入
addLayersControls の layers で壁面に対して，境界層レイヤーを挿入することができる．
境界層レイヤーを挿入するために必要な設定パラメータは「レイヤ厚さ」，「レイヤ数」，「拡大比」である．これらのパタメータを組み合わせて，解析に適したレイヤーを挿入する．
Code 15 system/snappyHexMeshDict
213 addLayersControls
214 {
215
216

// 以下で指定するレイヤのサイズが相対値(true)か絶対値(false)か？
// 相対値の場合，レイヤの外側の格子幅に対する比を指定する．

28

217

relativeSizes false;

218
219
220
221
222
223
224
225

// パッチ毎のレイヤー情報
layers{lowerWall // パッチ名{
nSurfaceLayers 3; // レイヤ数}

226

"AhmedBody_.*" // パッチ名{
nSurfaceLayers 5; // レイヤ数
}

227
228
229
230
231

}

232
233

expansionRatio 1.2; //レイヤー拡大比

234
235
236 //

firstLayerThickness 0.0025; //最初のレイヤー厚さ [m]finalLayerThickness 1; // 壁面最遠界層レイヤー厚さ [m]29図 8 境界層レイヤー (車体表面 5 層，地表面 3 層)30

13.10 前処理のログの確認
ジョブが終了したら，必ずジョブのエラー出力ファイル*.sh.e[ジョブ ID] にエラーが出力されていない無い事を確認する．
1 more *.sh.e* # または， ls -l *.sh.e* でファイルサイズが0である事を確認する

次に，more コマンドを用いて，snappyHexMesh ログの確認を行う．
1 more *.snappyHexMesh

最初に Intel MPI のデバッグ情報が出力されているので，プロセスのピニングが意図したものか確認する．
pre.sh.l*.snappyHexMesh
1 [0] MPI startup(): Multi-threaded optimized library
2 [0] MPI startup(): shm data transfer mode

# [ランク] 使用されるファブリック

3 # 略# [ランク] 使用されるファブリック
5 # Omni-Pathのデフォルトでは，ノード内はshm(共有メモリ)，
ノード間はtmi(Tag Matching Interface)
6 [0] MPI startup(): Rank
Pid
Node name Pin cpu # プロセスのピニング情報
7 [0] MPI startup(): 048661
c3847.ofp 2# Rank0 → Pin2
8 # 略
9 [0] MPI startup(): 63
48733
c3847.ofp 65# Rank63 → Pin65
10 [0] MPI startup(): I_MPI_DEBUG=5
11 [0] MPI startup(): I_MPI_FABRICS_LIST=tmi # 使用するファブリックリストの指定(注)12 [0] MPI startup(): I_MPI_FALLBACK=0# 指定したファブリックを必ず使用する(注)4 [63] MPI startup(): shm data transfer mode31

13 [0] MPI startup(): I_MPI_INFO_NUMA_NODE_MAP=hfi1_0:0
14 [0] MPI startup(): I_MPI_INFO_NUMA_NODE_NUM=2
15 [0] MPI startup(): I_MPI_JOB_FAST_STARTUP=1 # 高速プロセス起動アルゴリズムon(注)16 # 略
17 [0] MPI startup(): I_MPI_PIN_MAPPING=64:0 2,1 3,2 4,3 5,4 6,5 7,6 8,7 9,8 10,9 11,10

12,11 13,12 14,13 15,14 16,15 17,16 18,17 19,18 20,19 21,20 22,21 23,22 24,23 25,24
26,25 27,26 28,27 29,28 30,29 31,30 32,31 33,32 34,33 35,34 36,35 37,36 38,37 39,38
40,39 41,40 42,41 43,42 44,43 45,44 46,45 47,46 48,47 49,48 50,49 51,50 52,51 53,52
54,53 55,54 56,55 57,56 58,57 59,58 60,59 61,60 62,61 63,62 64,63 65 # ピニング設定(注)18 # (注) これらはmpi_core_setting.shで設定される環境変数

また，OpenFOAM のアプリケーションのログにおける nProcs や Hosts の出力が意図したプロセッサ数やノード数となっていない場合には，ジョブファイルの指定を確認する．
pre.sh.l*.snappyHexMesh: v1712 # ビルドバージョン
2 Arch: "LSB;label=32;scalar=64" # "バイト順;ラベルbit数;実数変数bit数"
3 Exec: snappyHexMesh -decomposeParDict system/decomposeParDict.mesh -overwrite parallel # 実行コマンド
4 Date: Jan 1 1970 # 開始日時
5 Time: 00:00:00# 開始時刻6 Host: "cxxxx.ofp" # ホスト名
7 PID
: xxxxx# プロセスID
8 I/O: uncollated # 領域分割ファイルを1ファイルにまとめない
9 Case: /work/0/gt00/txxxxx/lecture/Ahmed25/case-1 # ケースディレクトリ
1 Build32

10 nProcs : 64 # 計算で使用されているプロセス数
11 Hosts:

12 ((cxxxx.ofp 64) # (計算ノードのホスト名 MPIプロセス数)

13
14 )

パッチに付加されたレイヤの統計情報は，snappyHexMesh のログの最後に出力される．
Code 16 pre.sh.l*.snappyHexMesh
2281 patch

faces

layers----2284 lowerWall
38953
2285 AhmedBody_body 16595
2286 AhmedBody_head 2288
2287 AhmedBody_stilt 240
2288 AhmedBody_slant 1062-----2.99
4.85
4.65
1.3
4.58

2282
2283 -----

overall thickness[m][%]----0.00831
91.4
0.015583.1
0.016186.4
0.002714.5
0.0176
94.5

また，checkMesh のログにエラーが無い事を確認する．
checkMesh のログ確認
1 more *.checkMesh

33

14 ファイルの転送
今回は，データをクライアント側に転送して可視化するので，ユーザマシンの端末で以下のコマンドを実行して，講習用ファイル一式をユーザマシンに転送する．
(ユーザマシン) 講習会用ファイルの転送
1 cd
2 ls
3 # もし，
既にlectureが存在したら，mv lecture lecture.0 のようにして別名に変更する．
4 mkdir lecture
5 rsync txxxxx@ofp.jcahpc.jp:lecture/ ~/lecture/ -auv --exclude=processor*

ここで，txxxxx は利用者番号である．
なお，転送元と転送先どちらもに末尾に/(スラッシュ) を付ける．オプションの意味は，以下の通りである．

• -a (–archive) : ディレクトリを再帰的かつ，ファイル情報を保持したまま転送．
• -u (–update) : 新規・更新されたファイル・ディレクトリのみ転送．
• -v (–verbose) : 転送情報を表示．
• –exclude=processor* : processor ディレクトリを除外して転送．
転送元の講習会用ディレクトリの実体は/work/gt00/$USER/lecture であるが，講習会用ディレクトリの参照を容易にするために，ホームディレクトリにシンボリックリンクを貼っているので，上記の指定で参照可能である．

34

15 格子の可視化
15.1 最終生成格子の可視化
ユーザマシンで ParaView を立ち上げて，OpenFOAM データ可視化用のダミーファイル
pv.foam を読み込み，格子を可視化する．
表 4 ParaView によるベース格子の可視化

1. ‘File/Open‘メニューで，ケースのディレクトリの pv.foam を選択し，OK を押す．
2. Mesh Regions で center，lowerWall，AhmedBody_*を選択する (レイヤを可視化する場合にはcentar のみを選択する)．
3. ‘Apply‘ボタンを押す．
4. ‘Properties‘タブ内の ‘Representation‘で ‘Surface With Edges‘を選択する．
5. ‘Properties‘タブ内の ‘Coloring‘で ‘Solid Color‘を選択する．
6. 歯車状のボタンを押して，‘Advanced Properties‘を有効にする．

7. ‘Camera Parallel Projection‘をチェックして，格子の分割幅がわかりやすい平行投影モードにする．

35

図9

生成格子

36

15.2 特徴辺の可視化• surfaceFeatureExtract コマンドを実行すると，constant/triSurface に特徴辺ファイル.eMesh が作成される．
• 特徴辺ファイルは OpenFOAM の独自形式であるため ParaView 等で可視化できない．
• writeObj が yes の場合には，constant/extendedFeatureEdgeMesh に WaveFront OBJ 形式の特徴辺が AhmedBody_edgeMesh.obj として出力される．
• ParaView の ‘File/Open‘メニューから
constant/extendedFeatureEdgeMesh/AhmedBody_edgeMesh.obj を読み込む．

図 10

特徴辺

37

16 流体解析用の領域分割
Code 17 decomposePar.sh
1 #!/bin/bash
2 #PJM -S
3 #PJM -g gt00
4 #PJM -L elapse=0:15:00
5 #PJM -L rscgrp=tutorial-flat
6 #PJM -L node=1
7 source share.sh # 共通の設定
8 rm -rf 0;mkdir 0;cp -a 0.orig/* 0/ # 初期時刻ディレクトリを用意
9 $numa decomposePar -force >& $log.decomposePar # 流体解析用の領域分割
10 # -force：並列格子生成用のプロセッサディレクトリが残っているので消去する

Code 18 system/decomposeParDict
18 numberOfSubdomains 256; // 並列数
1920 // 分割手法
21 method

scotch;

今回は，1 ノードあたり 64 プロセスで 4 ノード使用し，計 256 プロセスの並列解析を行う．
流体解析用の領域分割
1 pjsub decomposePar.sh # 講習会終了後は pjsub -L rscgrp=lecture-flat decomposePar.sh38

17 初期値作成
Code 19 solveInit.sh
1 #!/bin/bash
2 #PJM -S
3 #PJM -g gt00
4 #PJM -L elapse=0:15:00
5 #PJM -L rscgrp=tutorial-flat
6 #PJM -L node=4
7 #PJM --mpi proc=256
8 source share.sh # 共通の設定
9 # 格子順変更による行列バンド幅縮小
10 $mpirun $numa renumberMesh -overwrite -parallel >& $log.renumberMesh
11 # ポテンシャル流れを解いて初期値とする
12 $mpirun $numa potentialFoam -parallel >& $log.potentialFoam• -L node=4 : ノード数=4• –mpi proc=256 : 総プロセス数=256=4×64(ノード毎のプロセス数ではない)• renumberMesh : 格子順変更→行列バンド幅縮小– 有限体積法では，圧力の Poisson 方程式および運動量 (速度) や乱流統計量の輸送方程式
は，各種モデル化や system/fvSchemes で設定した離散化スキームを用いて離散化され，最終的には線型一次方程式に帰着される．
– 線型一次方程式は，system/fvSolution で設定した線型ソルバを用い，許容残差の範39

囲で解かれる．
– 係数行列のバンド幅が小さいほうが一般に計算効率が良いため，バンド幅を縮小させる
ために，格子の順番を変更する．
• potentialFoam : ポテンシャル流れを解いて初期値とする– 流体の粘性を無視しているため短時間で解けるポテンシャル流れを解き，流体解析の初期値として用いる．
– ポテンシャル流れを初期値とすると，通常安定かつ早く収束できる．

17.1 初期値作成ジョブの投入
初期値作成
1 pjsub solveInit.sh # 講習会終了後は pjsub -L rscgrp=lecture-flat solveInit.sh40

18 流体解析
Code 20 solve.sh
1 #!/bin/bash
2 #PJM -S
3 #PJM -g gt00
4 #PJM -L elapse=0:15:00
5 #PJM -L rscgrp=tutorial-flat
6 #PJM -L node=4
7 #PJM --mpi proc=256
8 source share.sh # 共通の設定
9 $mpirun $numa simpleFoam -parallel >& $log.simpleFoam # 並列流体解析

流体解析実行
1 pjsub solve.sh # 講習会終了後は pjsub -L rscgrp=lecture-flat solve.sh41

19 関数出力のモニター
線型ソルバーの残差や空力係数をモニターする場合は foamMonitor コマンドを用いる．
Code 21 foamMonitor.sh
1 #!/bin/bash
2 module purge # 標準で有効なmoduleに依存しないよう全moduleをunload
3 module load intel/2018.1.163 # intelコンパイラ環境をload
4 source /work/gt00/share/OpenFOAM/OpenFOAM-v1712/etc/bashrc \5 WM_COMPILER=Icc2018_1_163 WM_MPLIB=INTELMPI2018_1_163 # OpenFOAM-v1712環境設定
6 foamMonitor -r 1 postProcessing/forceCoeffs/0/coefficient.dat & # -r :更新秒
7 foamMonitor -r 1 -l postProcessing/residuals/0/residuals.dat &solve.sh のジョブが開始されたら，以下を実行する．
関数出力のモニター
1 ./foamMonitor.sh

42# -l : 縦軸log

図 11 線型ソルバー残差のモニター

図 12

43

空力係数のモニター

20 流体解析の設定
20.1 流体物性値の設定
粘性係数などの流体物性値は constant/transportProperties で設定する．空気の動粘性
係数は ν = 1.5 × 10−5 m 2/s を設定する．
Code 22 constant/transportProperties
19 nu [0 2 -1 0 0 0 0] 1.5e-05; // 動粘性係数 [m^2/s]

20.2 乱流モデルの設定
乱流モデルは constant/turbulenceProperties で設定する．RANS の k - ω SST モデルの乱流モデルを設定する．
Code 23 system/turbulenceProperties
17 simulationType RAS;
18
19 RAS
20 {
21
22 //
23 //
24 //
25 //

RASModelRASModelRASModelRASModelRASModelkOmegaSST;// k-ω SSTモデル
kEpsilon;// 標準k-εモデル
RNGkEpsilon; // RNG k-εモデル
realizableKE; // Realizable k-εモデル
kOmega;// k-ωモデル

44

26
27

turbulenceon; // 乱流の計算：on=有効、off=無効

printCoeffs

on; // 乱流モデル係数の表示：on=有効、off=無効

28
29
30 }

RASModel で乱流モデルの種類を設定する．また，printCoeffs が on の場合は，計算初期に乱流モデル係数が表示される．バージョンによって，デフォルトの係数が変更される場合もあるので，計算条件保存のため on にしておくほうが良い．




Tips

例えば，OpenFOAM-5.0 および v1712 から，標準 k − ε モデル (kEpsion) のモデル係数 C3のデフォルト値が-0.33 から 0 に変更されている．



20.3 離散化スキームの設定
離散化スキームは system/fvSchemes で設定する．移流項のスキーム divSchemes(phi, 変数) は，速度について TVD 系スキーム limitedLinear のベクトル版 limitedLinearV を使用する．limitedLinear スキームのパラメータは 0 から 1 の範囲を取るが，最も安定となる 1 を指定した．乱流統計量については，最も安定である upwind(一次風上) を用いた．もちろん，他のスキームでの解析可能である．
Code 24 system/fvSchemes45

27 divSchemes
28 {29default

none;

30
31
32
33 //
34 //
35
36 //
37

// 移流項の離散化スキーム
div(phi,U)bounded Gauss limitedLinearV 1; // TVD系スキーム
div(phi,U)bounded Gauss upwind;div(phi,k)bounded Gauss limitedLinear 1;
div(phi,k)bounded Gauss upwind;// 風上(最安定・1次精度・数値拡散過大)
div(phi,omega)bounded Gauss limitedLinear 1;
div(phi,omega)bounded Gauss upwind;// 風上(最安定・1次精度・数値拡散過大)

20.4 解法の設定
Code 25 system/fvSolution
17 solvers
18 {
19
20
21
22
23
24

// 圧力のPoisson方程式の線型ソルバの設定
p{solver
PCG; // 線型ソルバ：PCG(前処理付き共役勾配法)
preconditioner DIC; // 前処理：DIC(diagonal incomplete-Cholesky)tolerance
1e-7; // 最終残差の許容誤差

46

25
26

relTol0.01; // 相対残差(最終残差/初期残差)の許容誤差

}

27
28
29
30
31
32

// potentialfoam用の線型ソルバの設定
Phi{$p; // p(圧力)の定義を参照する(pと同じ設定を用いる)
}

33
34
35
36
37
38
39
40
41
42

// U, k, omega, epsilonの輸送方程式の線型ソルバ"(U|k|omega|epsilon)"{solversmoothSolver; // 線型ソルバ：定常反復法
smoother
GaussSeidel; // スムーサ：ガウス-ザイデル法
tolerance
1e-8; // 最終残差の許容誤差
relTol0.1; // 相対残差(最終残差/初期残差)の許容誤差
nSweeps
1; // 反復内のスイープ数
}

43 }
44
45 SIMPLE
46 {
47
48

// 非直交補正反復回数// 格子の非直交性が悪く初期に発散する場合には非直交補正の反復を1以上にする

47

49

nNonOrthogonalCorrectors 1;

50
51
52

// SIMPLEC(SIMPLE consistent)型か？
consistent yes;

53 }54
55 // SIMPLE，
SIMPLEC法の緩和係数
56 // 発散する場合には小さくする
57 relaxationFactors
58 {
59
60
61
62
63

equations{U"(k|omega|epsilon)"
}

0.9;
0.9;

64 }
65
66 // potentialFoamの反復回数
67 potentialFlow
68 {69

nNonOrthogonalCorrectors 15;

70 }

48

20.5 実行制御の設定
Code 26 system/controlDict
simpleFoam;// アプリケーション名
18 startFrom
startTime; // 解析開始の設定法
19 startTime0;// 解析の開始時間(定常計算では，時間=反復番号)
20 stopAt
endTime;// 解析終了の設定法
21 endTime500;// 解析の終了時間
22 deltaT
1;// 解析時間刻み(定常計算では通常1)23 writeControl
timeStep; // 解析結果書き出しの決定法
24 writeInterval
500;// 書き出す間隔
25 runTimeModifiable true;// 各時間ステップで設定ファイルを再読み込みするか
26 writeFormat
binary;// データファイルのフォーマット(binary, ascii)
27 writePrecision
6;// データファイルの有効桁(上記がasciiの場合)
28 writeCompression off;// データファイルの圧縮(off, on)
29 timeFormat
general;// 時刻ディレクトリのフォーマット
30 timePrecision6;// 時刻ディレクトリのフォーマット有効桁
17 application

20.6 計算履歴の出力
抗力・揚力係数 ( Cd , Cl ) や，線型ソルバ残差の履歴を出力するには，OpenFOAM の関数
機能 functions を利用するため，system/controlDict の functions に，#include を用
いて関数を定義したファイルを設定する．関数定義のファイルは system に置く．

49

Code 27 system/controlDict
32 functions
33 {34#include "forceCoeffs"#include "forces"#include "residuals"

35
36
37 }

空力係数の関数を定義した system/forceCoeffs では，以下のように，車体のパッチ名，
自由流速度 (流入速度)，代表長さ，代表面積 (見付面積) などを定義する．空力の関数を定義した system/forces もほぼ同様である．
Code 28 system/forceCoeffs
9 forceCoeffs
10 {
11
12
13
14
15
16
17
18
19
20

type
forceCoeffs;// 関数型libs("libforces.so"); // 関数のライブラリ
writeControl timeStep;// 出力間隔手法(時刻ステップベース)
timeInterval 1;// 算出間隔(毎回)logyes;// ログ出力の有無
patches( "AhmedBody_.*"); // 車体のパッチ名
rhorhoInf;// 密度の定義
rhoInf
1.225;// 流体密度(非圧縮性解析の場合，ダミー)
liftDir(0 0 1);// 揚力方向dragDir(1 0 0);// 抗力方向50

21pitchAxis
CofR
magUInf
lRef
Aref

22
23
24
25

(0 1 0);// ピッチング軸方向(-0.5269 0 0.186); // モーメント中心
40;// 自由流速度(流入速度)0.47;// モーメント算出用代表長さ
0.0575;// 代表面積(見付面積．0.389*0.288/2+0.03*0.05)

26 }線型ソルバ残差の出力の関数定義は以下のように，残差を出力する場のリストなどを定義する．
Code 29 system/residuals
9 residuals
10 {

type
residuals;// 関数型libs("libutilityFunctionObjects.so"); // 関数のライブラリ
writeControl
timeStep;// 出力間隔手法(時刻ステップベース)
writeInterval1;// 算出間隔(毎回)
fields (p U k omega);// 残差を出力する場のリスト

11
12
13
14
15
16 }

20.7 境界条件の設定
本モデルの境界条件を表 5 に示す．境界条件は 0.orig ディレクトリに置かれた物理量ファイル (U, p, k, omega, nut...) で設定する．

51

境界名

境界型

U

inlet

patch

outletlowerWall,
Ahmed_.*
upperWall,sidecenterpatch

fixedValuepressure
InletOutlet
Velocitywall

noSlip表5

境界条件
εturbulent
MixingLength
Dissipation
Rate
Inlet

ω

νt

zero
Gradientkturbulent
Intensity
Kinetic
EnergyInletturbulent
MixingLength
Frequency
Inletcalculated

fixedValuezero
Gradient

inletOutlet
kqRWall
Function

inletOutlet
epsilonWall
Function

inletOutlet
omegaWall
Function

calculated
nutkWallFunctionpwall
symmetry

slip
symmetry

20.8 解析結果の再構築
Code 30 reconstructPar.sh
7 source share.sh # 共通の設定
8 $numa reconstructPar -latestTime >& $log.reconstructPar # 最終時刻の解析結果の再構築

解析結果の再構築
1 pjsub reconstructPar.sh # 講習会終了後は pjsub -L rscgrp=lecture-flat reconstructPar.sh52

21 空力係数の実測値との比較
Ercoftac ベンチマークでは車両の空気抗力係数 Cd と車両周りの対称面における主流方向速度分布について検証されているので，実験値との比較検証を行う．
空力係数の関数の算出結果は，流体解析実行時にpostProcessing/forceCoeffs/0/coefficient.dat 出力されているので，以下のようにして値を確認し，風洞実験値と比較する．結果を表 6 に示す．
空力係数の関数出力の表示
1 more postProcessing/forceCoeffs/0/coefficient.dat

Code 31 空力係数の関数出力
10 # Time
11 500Cm-2.446843e-01

Cd
3.247248e-01

Cl
3.799807e-01

-01

表6

抗力係数 Cd
揚力係数 Cl

空力係数の比較実験

CFD(template)

0.299
0.345

0.325
0.380

53

Cl(f)
-5.469401e-02

Cl(r)4.346747e

22 後処理
Code 32 post.sh
1 #!/bin/bash
2 #PJM -S
3 #PJM -g gt00
4 #PJM -L elapse=0:15:00
5 #PJM -L rscgrp=tutorial-flat
6 #PJM -L node=4
7 #PJM --mpi proc=256
8 source share.sh # 共通の設定
9 # サンプリング
10 $mpirun $numa postProcess -parallel -latestTime -func sample >& $log.sample
11 # 無次元化壁座標y+算出
12 $mpirun $numa simpleFoam -parallel -latestTime -postProcess -func yPlus >& $log.yPlus

解析結果をサンプリングしたり，無次元化壁座標 y + を算出するなどの後処理についても並列実行が可能である．
後処理ジョブの投入
1 pjsub post.sh # 講習会終了後は pjsub -L rscgrp=lecture-flat post.sh

サンプリングの設定は system/sample で行う．

54

Code 33 system/sample
18 type sets;
19 libs ("libsampling.so");// 関数型// 関数のライブラリ

20
21 // 補間方法: 格子中心と節点，界面における値を用いて補間
22 interpolationScheme cellPointFace;
23
24 // 集合サンプリング出力形式 raw:ASCII生データ形式
25 setFormat raw;
26// サンプルするフィールドのリスト
28 fields
29 (
30
U
31 );
27

32
33 // 集合サンプリングの定義
34 sets
35 (
36
37
38
39
40

x-243 // 鉛直プロファイルラインのサンプリング名{// サンプルタイプ: 線分と格子界面との交点の中点および界面
type
midPointAndFace;// サンプリング座標値の出力軸55また，無次元化壁座標 y + は，post.sh.l*.yPlus で確認できる．
Code 34 post.sh.l*.yPlus
601 yPlus yPlus write:
602
603
604
605
606
607
608
609

writing field yPlus
patch side y+ : min = 0.4074, max = 1883.17, average = 74.6756
patch lowerWall y+ : min = 24.8497, max = 392.265, average = 108.25
patch upperWall y+ : min = 3.88533, max = 174.989, average = 102.439
patch AhmedBody_body y+ : min = 27.5751, max = 462.126, average = 116.744
patch AhmedBody_head y+ : min = 13.1871, max = 143.468, average = 59.8327
patch AhmedBody_stilt y+ : min = 33.0879, max = 580.702, average = 126.926
patch AhmedBody_slant y+ : min = 68.1991, max = 370.604, average = 139.328

56

23 速度分布のプロットと実験値との比較• 実験値と解析結果のデータをプロットしてグラフ化するが，ここでは gnuplot を用いる．
• 今回は車体後部上面の各点についての速度分布を確認するために車体後部上面の形状もプロットした．他にも，線型ソルバーの残差，空力係数，空力のプロット用 gnuplot 入力
ファイルが plot ディレクリに収めてあるので，これらを gnuplot でプロットする (もちろん手動で gnuplot を実行しても良い)．
Code 35 plot.sh
7 for file in plot/*.gp
8 do
9

gnuplot $file

10 done

上記を実行し，画像ビューワ evince で表示するには以下のようにする．
プロット実行・表示
1 ./plot.sh
2 evince *.pdf

57

0.5

0.4

z [m]

0.3

0.2

Exp.
CFD

0.1

0
-0.3

-0.2

-0.1

x [m]

0

0.1

図 13 速度の鉛直プロファイルの実験値との比較

58

0.2

24 ParaView による計算結果の確認
解析結果をユーザマシンに転送し，ParaView で流線を描き，後引き渦の存在を確認する．
(ユーザマシン) 計算結果のファイル転送
1 rsync txxxxx@ofp.jcahpc.jp:lecture/ ~/lecture/ -auv --exclude=processor*
2 または，
カーソル↑で上記を呼び出してリターン(エンター)キーを押す

1. Refresh ボタンを押して，転送した解析結果を読み込む．
2. Mesh Regions で internalMesh, center，lowerWall および AhmedBody_*を選択する．
3. Cell Arrays で U のみを選択する．
4. Last Frame ボタンを押して，Apply ボタンを押す．
5. Pipeline Browser で pv.foam を選択し，Extract Block フィルターを選択後，internalMesh 以外の全パッチを選択し，Apply ボタンを押す．
6. Coloring を ◦U にする．
7. Pipeline Browser で pv.foam を選択し，Extract Block フィルターを選択後，internalMesh を選択し，Apply ボタンを押す．
8. さらに Stream Tracer フィルターを選択し，Seed type に High resolution line
source を選択し，Point1=(-5,-0.935,0.2)，Point2=(-5,0,0.2)，
Resolution=100 を入力する．
9. Apply ボタンを押して流線を表示する．
10. Coloring を ◦U にする．
59

図 14

流線と対称面速度分布 (Color Map や Color legend の位置等は適宜変更する)

60

25 流体解析の性能プロファイル実行
25.1 流体解析の性能プロファイル• 性能プロファイラにより，計算負荷が高い部分 (hotspots) など計算効率改善のためのデー
タを，ソースの変更や再コンパイル無しに取得可能である．
• ソースレベルの詳細なプロファイルには，デバッグオプションを付けた再コンパイルが必要である．
• VTune Amplifier を用いて，流体解析の性能プロファイル実行を行うスクリプトを示す．
Code 36 solveVtune.sh
9 module load vtune/2018.1.0.535340 # VTune Amplifier XEのmoduleをload
10 # VTune Amplifier XEのコマンドamplxe-clをgtoolオプションで指定
11 # -collect hotspots : hotspot(計算負荷が高い部分)の解析を行う
12 # -r ディレトリ名:MPIランク範囲 : MPIランク範囲に関する解析結果をディレトリ名に保存
13 $mpirun -gtool "amplxe-cl␣-collect␣hotspots␣-r␣$PJM_JOBNAME.v$PJM_JOBID:0" \14 $numa simpleFoam -parallel >& $log.simpleFoam # 並列流体解析

上記のジョブを投入することで，流体解析の性能サンプリングを行うことができる．
流体解析の性能サンプリング
1 vi system/controlDict # endTimeを101と短くする(計算負荷の概要把握には十分)
2 pjsub solveVtune.sh # 講習会終了後は pjsub -L rscgrp=lecture-flat solveVtune.sh61

25.2 プロファイラ結果の表示
Code 37 vtune.sh
7 module load vtune/2018.1.0.535340 # VTune Amplifier XEのmoduleをload
8 for dir in solveVtune.sh.v*/
9 do# amplxe-cl : VTune Amplifier XEのコマンドライン版コマンド# -r データディレトリ名# -show-as=percent : 関数の実行時間を秒ではなく，割合で示す# -format=csv : テキスト形式ではなく，CSV形式で出力する(区切記号はコンマではなくタブ)
amplxe-cl -R hotspots -r $dir -show-as=percent -format=csv > ${dir%/}.csv

10
11
12
13
14
15 done• ここでは，GUI ツールを用いず，CLI ツールによりプロファイル結果を CSV 形式に変換して表示する (GUI は [8] 参照)．．
• ジョブの実行終了後，上記のスクリプトを実行して，プロファイルデータをテキスト変換
して後，表示する．
プロファイルデータの変換と表示
1 ./vtune.sh
2 more *.csv

以下が CPU 時間の割合を示した CSV 形式の結果ファイルの一部である．

62

Code 38 simpleFoam-vtune.sh.v*.csv
1 Function

CPU Time
2 Foam::DICPreconditioner::precondition
3 PMPI_Allreduce 13.681334
4 Foam::lduMatrix::Amul
9.519370
5 PMPI_Waitall
4.315902
6 Foam::GaussSeidelSmoother::smooth)
7 Foam::PCG::solve
3.067522

15.049757 # 圧力PCG線型ソルバのDIC前処理# MPIのAllreduce(線型ソルバのベクトル内積)# 線型ソルバーでの疎行列ベクトル積SpMV# MPIの非同期通信の終了待ち
3.337809 # 速度と乱流量統計量線型ソルバ(GaussSeidel# 圧力PCG線型ソルバ• 計算時間の割合は，線型ソルバ関連の
Foam::DICPreconditioner::precondition(PCG 線型ソルバの DIC 前処理) が 1 位，Foam::lduMatrix::Amul(線型ソルバーでの疎行列ベクトル積 SpMV) が 3 位であるが，MPI 通信関連の PMPI_Allreduce が 2 位，PMPI_Waitall が 4 位となっている．
• MPI プロセス数が増加するに伴ない，線型ソルバにおける MPI の Allreduce が計算時間の多くを占めるようになるので，超大規模ノードの解析では，MPI 通信時間を削減・隠蔽する必要がある [9]．

63

26 並列計算ベンチマークテスト• 領域分割数や領域分割方法，メモリモード，圧力の線型ソルバ等を変更して流体解析を行い，解析時間を比較する．
– 領域分割数 (system/decomposeParDict の numberOfSubdomains)– 領域分割法: scotch, multiLevel, simple, hierarchical など(system/decomposeParDict の method)()– 圧力の線型ソルバ，前処理，スムーサ (system/fvSolution)– ノード数 (pjsub -L node=ノード数 –mpi proc=領域分割数)– メモリモード (pjsub -L rscgrp={tutorial,lecture}-{flat,cache}))• 予備検討で行う並列計算のベンチマークテストでは，時間 (反復) ステップ数を本計算より短くして検討する事が多い．
• 解析の初期化にかかる時間は，ステップ数が多い本計算では相対的に小さくなる．
• ここでは，初期ステップ完了から最終ステップ完了までの解析時間のステップあたりの平均値を比較する．
• P ノード並列時の初期と最終のステップの解析時間 (ExecutionTime) の差を t(P ) として，以下を求める．
– スピードアップ率：SP = t(1)/t(P )– 並列化効率 [%]：EP = SP /P × 100 = (t(1)/t(P ))/P × 100• 生成に時間がかる格子データ，領域分割データなどは，なるべく再生成せずに済むようサ
ブディレクトリを作成して，共有データはシンボリックリンクする．

64

並列計算ベンチマークテスト実行例
1 cd ~/lecture/Ahmed25/case-1# 領域分割ケースのディレクトリ作成
3 cd mpi_000064-method_scotch# cd !$ でも良い
4 ln -s ../{0.orig,constant,*.sh} ./# 親ディレクトリ内の必要なファイルにリンク
5 cp -a ../system ./# systemディレクトリ以下は変更するのでコピー
6 vi system/decomposeParDict# 領域分割の変更(ここではnumberOfSubdomainsを64)
7 pjsub decomposePar.sh# 領域分割
8 pjsub -L node=1 --mpi proc=64 solveInit.sh # decomposePar.shが終了後
9 mkdir solver_PCG-DIC# 圧力の線型ソルバのケース
10 cd solver_PCG-DIC# cd !$ でも良い
11 ln -s ../{processor*,constant,*.sh} ./# 親ディレクトリ内の必要なファイルにリンク
12 cp -a ../system ./# systemディレクトリ以下は変更するのでコピー
13 vi system/fvSolution# 圧力の線型ソルバを必要に応じて変更
14 # system/fvSolution.benchに使用可能な線型ソルバ，
前処理，
スムーサの定義を記述しているので参照する
15 pjsub -L node=1 --mpi proc=64 -L rscgrp=tutorial-flat solveVtune.sh # 64プロセス/ノード
16 pjsub -L node=1 --mpi proc=64 -L rscgrp=tutorial-cache solveVtune.sh # 上記のCacheモード
17 pjsub -L node=2 --mpi proc=64 -L rscgrp=tutorial-flat solveVtune.sh # 32プロセス/ノード
18 pjsub -L node=2 --mpi proc=64 -L rscgrp=tutorial-cache solveVtune.sh # 上記のCacheモード
19 # 以下同様に別の領域分割ケース/線型ソルバのディレクトリを作成し，解析ジョブ実行する
2 mkdir mpi_000064-method_scotch

解析時間の集計結果表示
1 cd ~/lecture/Ahmed25/case-1
2 ./averageExecutionTime.sh65

27 格子生成・空力解析演習
27.1 格子生成・空力解析演習手順• template ディレクトリをコピーし，設定の変更を行なう• 流体解析を実行し，空力係数，速度分布の実験値との比較を行う．
• ./Allrun.batch を実行すると，前処理から結果プロットまでの全ジョブを，依存関係があるサブジョブとするステップジョブとして投入する．

図 15 ステップジョブイメージ (図引用元:[7])6627.2 格子生成・空力解析演習設定変更例• 生成格子を変更 (設定: system/snappyHexMeshDict)– 床面や車体表面のレイヤ数・レイヤ最小厚み変更– 車体表面の分割レベルの変更• 流体解析の初期値の境界条件を変更– 乱流動粘性係数の nut の壁関数変更 (設定: 0.orig/nut)• 流体解析の計算条件を変更– 乱流モデルの変更 (設定:constant/turbulenceProperties)– 移流項離散化スキームの変更 (設定:system/fvSchemes)

27.3 格子生成・空力解析演習実行例
1 cd ~/lecture/Ahmed25
2 cp -a template case-2 # コピー先のディレクトリ名は任意
3 # 既往の生成格子を用いる場合には以下のようにしてコピーする
4 # cp -a case-1/constant/polyMesh case-2/constant/
5 cd case-2
6 ./Allrun.batch
7 pjstat -E
8 ./foamMonitor.sh# この後，適宜設定変更# 格子が生成済の場合pre.shジョブを実行しない(オプションは-hで確認)# -E: ステップジョブ内のサブジョブを展開して表示# solve.shがRUNNINGになったら実行する

9 evince *.pdf

67

28 補足28.1 Intel MPI ライブラリでの実行エラー
Intel MPI ライブラリを用いて OpenFOAM を実行すると，多ノードで実行エラーになる場合があるが，以下のようにして RDMA translation キャッシュ機能を OFF にすると，実行エ
ラーを回避できる可能性がある．
RDMA translation キャッシュ機能 OFF
1 export I_MPI_DAPL_TRANSLATION_CACHE=0
2 export I_MPI_DAPL_UD_TRANSLATION_CACHE=028.2 Oakforest-PACS での OpenFOAM のコンパイル
OpenFOAM 自動ビルドスクリプト installOpenFOAM の Web ページを参考にコンパイル
する．

URL: https://gitlab.com/OpenCAE/installOpenFOAM/blob/master/README.md

68

参考文献
[1] The OpenFOAM Foundation[2] S.R. Ahmed, G. Ramm, Some Salient Features of the Time-Averaged Ground Vehicle
Wake, SAE-Paper 840300, 1984
[3] H.Lienhart, C.Stoots, S.Becker: Flow and Turbulence Structures int the Wake of a
Simplified Car Model (Ahmed Model), DGLR Fach Symp. der AG STAB, Stuttgart
University, pp.15-17, November 2000
[4] Minguez, M., Pasquetti, R., Serre, E.: Spectral vanishing viscosity stabilized les of theahmed body turbulent wake, Communications in computational physics, 5(2-4),
pp.635-648, 2009[5] C. Hinterberger, M. Garcia-Villalba, W. Rodi: Large Eddy Simulation of flow aroundthe Ahmed body, "Lecture Notes in Applied and Computational Mechanics / The
Aerodynamics of Heavy Vehicles: Trucks, Buses, and Trains", R. McCallen, F.
Browand, J. Ross (Eds.), Springer Verlag, ISBN: 3-540-22088-7, 2004
[6] Flow Around a Simplified Car Body (Ahmed Body) , ERCOFTAC Classic Database
C82[7] Oakforest-PACS 利用支援ポータル「Oakforest-PACS システム 利用手引書」[8] Oakforest-PACS 利用支援ポータル「インテル VTune Amplifier XE 2018」ドキュメント[9] 今野 雅: OpenFOAM における Communication Avoiding CG 法の実装と性能評価，
オープン CAE シンポジウム 2017, 2017 年 12 月
69

Copyright
本ドキュメントは，クリエイティブ・コモンズ 表示 -非営利 4.0 国際 ライセンスの下に提供されています．

This work is licensed under a Creative Commons Attribution-NonCommercial 4.0
International License.

70