情報知能工学演習V (前半第6週) 政田洋平 システム情報学研究科計算科学専攻 TA : 菅 新 菅沼智史 水曜 新行紗弓 馬淵隼 木曜
演習 V( 前半 ) の内容 第 1 週 : 高性能計算 (High Performance Computing = HPC) 向けプログラミングの基礎 第 2 週 : シミュレーションの基礎 第 3 週 : 波の移流方程式のシミュレーション 第 4,5 週 : OpenMP による並列化 第 6 週 : シミュレーションコードの OpenMP による並列化 到達目標 :HPC シミュレーション 並列化の基礎を理解する 授業予定 水曜 4/23, 4/30, 5/7, 5/14, 5/21, 5/28 木曜 4/24, 5/1, 5/8, 5/15, 5/22, 5/29, 6/5 - 参考図書 - C/C++ プログラマーのためのOpenMP 並列プログラミング OpenMPによる並列プログラミングと数値計算法流体力学の数値計算法 菅原清文著 牛島省著 藤井孝蔵著
いよいよ最終回 : 移流方程式のシミュレーションコードを並列化
OpenMP の復習 :for 構文 ( 先々週の資料を参照 ) 並列リージョンの中で #pragma omp for を挿入する. 並列リージョン #include <stdio.h> #include <omp.h> int main(void){ int i; double x[100]; #pragma omp parallel { #pragma omp for for(i=0; i<100; i++){ x[i]=i; } } } 直後の for ループを複数のスレッドで分割 ループ変数の分散は 4 スレッドの場合 i スレッド 0 0~24 スレッド 1 25~49 スレッド 2 50~74 スレッド 3 75~99
OpenMP の復習 :for 構文 ( 先週の資料を参照 ) どの for ループでも並列化をしても良いわけではない. ループ内で変数 配列に依存関係がある場合には不可. 例 ) 再帰参照を含むループ for(i=0; i<100; i++){ x[i+1]=x[i]+a; } ループ内で再帰参照が行われていないか 十分確認しなが ら並列化を進める必要がある ( 先週の課題と関連 ).
OpenMP を使って波の移流プログラムを並列化 練習問題 :#pragma omp for 文を挿入し 3 週目に作った波の移流方程式のプログラムを並列化せよ. 諸注意 omp.hをincludeするのを忘れない. #pragma omp parallel 文を忘れない. どのループを並列化すべきかよく考える. 空間ループか時間ループか? -fopenmpのコンパイルオプション. 環境変数 : export OMP_NUM_THREADS =? 空間分割数 N=10000で移流方程式の並列シミュレーションを実行しなさい. 並列化による加速が実感できるか? ( 結果を統一するためにtend=50.0で実行すること )
サンプルプログラム :advection_parallel.c #include <stdio.h> #include <omp.h> #define N 10000 int main(void){ int i,n,nend; double xleft, xright; double x[n+1],dx; double u[n+1],uu[n+1]; double c,t,dt,tend; xleft = 0.0; xright = 100.0; tend = 50.0; c = 0.5; dx = (xright - xleft)/(double)n; dt = 0.5*dx; nend = (int)tend/dt; /*Initialization*/ for (i=0; i < N+1; i++){ x[i] = ((double)i)*dx; } for (i=0; i < N+1; i++){ u[i] = 0.0; if (i <= N/2) u[i] = 1.0; } /*--------------*/ } /*--- Start Update ---*/ /*----- End Update ----------*/ for (i=0; i < N+1; i++){ printf("%10.5f %10.5f\n", x[i], u[i]); } return 0;
並列加速効率を調べる : 練習問題 : 並列化による計算の加速を実感するために for ループの分割スレッド数を変化させて 下のような図を作成せよ. 環境変数 : export OMP_NUM_THREADS =? (bash) 又は setenv OMP_NUM_THREADS? (tcsh) ( 注 ) 実行時間の計測には time./a.out を用いよ. (sec) 1 2 3 4 5 6 7 8 Number of threads 8スレッドで何倍の加速が実現できたか? N=10 4, 5 10 4, 10 5 で調べよ.
サンプル : スレッド数 v.s., 計算時間 (imac で計測 ) ( 複数コアでの実行時間 )/(1 コアでの実行時間 ) RUNNNING TIME (Normalized by the time using ONE THREAD) 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 分割数が小さな計算は スレッド数が多くなると逆に計算時間がかかる [ 通信処理に時間をとられる ]. 分割数 = 10 5 分割数 = 10 4 N = 10000 N = 50000 N = 100000 分割数 = 5x10 4 1 2 3 4 5 6 7 8 NUMBER OF THREADS 分割数が十分多い時 8 並列で計算時間は約 3 分の 1(3 倍加速 ).
最終課題 : 並列シミュレーションコードの作成
最終課題 : 最終課題は以下の 2 つの課題から 1 つを選択しプログラムを作成. ソースコードと図をメールに添付して提出. 6.1 波の移流方程式の並列プログラムを空間 2 次元に拡張. - 空間を x[i] と y[j] で分割. - 離散化方程式も y 方向の空間微分項の寄与が追加される. ( 詳細は次ページ ) 6.2 1 次元熱伝導方程式の並列プログラムの作成. - 2 回微分の離散化が必要 ( 手法は後述のヒントを参照 ). - 基本的には移流方程式プログラムのアップデート部分を変更. - 境界条件は自由境界から固定境界に変更. ( 詳細は次次ページ以降 )
波の移流方程式の並列プログラムを空間 2 次元に拡張 (1) 離散化 u n+1 i,j u t + c 1 u x + c 2 u y =0, ( 空間 2 次元の移流方程式 ) = u n i,j c 1 t 2 x (un i,j u n i 1,j) c 2 t 2 y (un i,j u n i,j 1) ( 2 週目の資料参照 ) y5 y4 Δx { { Δy y ( ) y3 y2 y1 y0 x0 x1 x2 x3 x4 x5 x
波の移流方程式の並列プログラムを空間 2 次元に拡張 (2) 手順 1. x[i] とy[j] および2 次元配列 u[i][j], uu[i][j] を確保. 2. 2 次元面内で配列 u[i][j] を初期化. 3. アップデート部分にy 方向の空間微分項の寄与を追加. 4. アップデート終了後 ループを回して結果を出力. 諸注意 x 方向とy 方向の分割数は同じにせよ : Nx = Ny = 500. 方程式を解く範囲は0 x, y 100とせよ. 初期条件は x 50 かつ y 50でu=1.0 それ以外でu=0.0とせよ. パラメータとしてc1=c2=0.5, dt=0.5dxを用いよ. 上下左右の境界条件は自由境界とせよ. i = jの時の (x[i] 2 +y[j] 2 ) 1/2 を横軸に i=jの時のu[i][j] を縦軸にとった図を提出せよ. 時刻 tend = 0, 50, 100での出力を1つの図にまとめて提出せよ.
熱伝導方程式の並列プログラムの作成 (1) 初夏の工学部棟内の平均温度は一体何度なのか? 25 現実の全ての条件を考慮に入れるのは あまりにも難しすぎる. 簡略化 ( モデル化 ) された問題をシミュレート. 簡略化手順 1. 工学部棟を1 次元の鉄の棒だとみなしましょう. 2. 鉄の棒の両端が一定温度 30 の熱浴に接していると考えましょう. 3. 冷房が一定の割合で鉄の棒を一様に冷やしていると考えましょう.
熱伝導方程式の並列プログラムの作成 (2) 初夏の工学部棟内の平均温度は一体何度なのか? 木もいらない... 25 人間もいらない... モデル化 クーラーがないとさすがに暑くて勉強できない. 太陽 25 1 次元の鉄の棒 25 太陽 クーラー = 内部は一定の割合で一様に冷却されている
熱伝導方程式の並列プログラムの作成 (3) 1 次元の鉄の棒の温度 T(t,x) は熱伝導方程式に従って発展する. T (t, x) t = κ 2 T (t, x) x 2 C, 熱が棒の内側に伝わる効果 クーラーによる冷却効果 端の温度 : 固定 諸注意 25 1 25 N=1000, κ = 1, C = 0.01, 初期温度分布を T(t=0,x) = 15 と仮定せよ. 工学部棟の幅(= 棒の長さ ) を100メートル [= xright - xleft] とせよ. 両端の温度 T(x=xleft) とT(x=xright) は25 で常に一定とせよ ( 固定境界 ). Δt = 0.1Δx 2 を用いよ. 時刻 t = 0, 50, 100, 500での温度分布 T(x) を求め 出力結果を1つの図にまとめて提出せよ ( 横軸 : 位置 x 縦軸: 温度 T). 端の温度 : 固定
熱伝導方程式の並列プログラムの作成 (4) 1 回空間微分の離散化 (3つの離散化手法を紹介) 2 x x f(x) = f(x i+1) f(x i 1 ) f(x i+1 ) f(x i ) f(x i ) f(x i 1 ) x xi 2 x ( ) x x (( 中心差分 ) )( ) (( 前進差分 ) ( ) ) ( 後退差分 ) 2 回空間微分の離散化 単純に考えると以下の 2 通り : fi+1/2 緑丸での関数 f の 2 回微分 Δx = ( 青丸 での関数 f の 1 回微分 fi-2 fi-1 fi fi+1 fi+2 赤丸 での関数 f の 1 回微分 )/2Δx 緑丸での関数 f の 2 回微分 fi-1/2 Δx Δx = ( 青星 での関数 f の1 回微分 赤星 での関数 f の1 回微分 )/Δx どちらの方法を使っても良い. ただし f の 1 回微分は中心差分を用いよ.
課題提出方法 ソースファイルを圧縮せずにメールに添付して提出 メール本体には 学籍番号, 氏名を必ず明記 締め切りは 1 週間後の午後 1 時 ( 次回の授業の直前 ) 分からなかった点や授業の感想があれば書いてください 提出先メールアドレスは ymasada_at_harbor.kobe-u.ac.jp [ _at_ @ ] Subject( 件名 ) は ENSHU5 :KADAI6: 学籍番号 ファイル名は ENSHU5 _KADAI6_ _ 学籍番号.c ENSHU5 _KADAI6_ _ 学籍番号.png にはクラス (a または b), には課題番号 (1 または 2), 学籍番号は半角小文字 加点課題は最終ページ
余裕のある人のための問題
2 次元熱伝導方程式の並列プログラム ( 余裕のある人用 ) [ 2. T(t,x) 2. ) T t = κ ( 2 T x 2 + 2 T y 2 C y ytop (y=ytop) = 30 (T=30 ) xleft=ybottom=0, xright=ytop=100[m] [ ] (x=xleft) = 30 (x=xright) = 30 Nx=Ny=500 κ=1, c=0.01 T(t=0, x)=20 ybottom xleft (y=ybottom) = 30 xright x