第 2 回シミュレーションスクール 第 2 日 : 熱拡散方程式のプログラムをつくろう陰山聡 政田洋平 神戸大学大学院システム情報学研究科計算科学専攻 2010.12.07
具体的な例題を一つに絞る 2 次元熱伝導問題 実習中心 各自自分のコードを 1 から作る そのコードを少しづつ改良していく 最後は並列化されたコードをスーパーコンピュータに
初歩から ゆっくり 丁寧に 網羅的な紹介よりも 基礎的な手法と考え方の習得を 具体的な例を通じてシミュレーションプログラムの作り方を体験的に学ぶ
実習中心 各自がプログラムを書く C/C++ 言語で書いてもOK file 名のつけ方は例題に習って ( 例 :001.f95 001.c )
今日と明日の講義ではMakefileは使わない使ったほうが便利 実用的 毎回 コンパイルコマンドを自分で打つほうが教育的 Makefileの文法の問題にとらわれがち
コード プログラム時間発展 2 次元 2-D 2D ディレクトリ フォルダ解析的に 解析解
マシン scalar にログイン cd mkdir Tue mkdir Wed cd Tue mkdir data cp -r /tmp/ss2/tue/src-work.
実習用 ( スケルトンコードあり 一部空白 ) 各自のデータを出力する場所 実習用コードのファイル名 : 数字.f95 実習問題の解答例
シミュレーション研究の流れ & 1 # 0 2 3 -! " # $ % & ' ( ) * +, -.! & " / 0 4 5 6 7 8 - % & ' ( 9 : -.! ; ( < = -
シミュレーション研究の流れ今日の目標 & 1 # 0 2 3 -! " # $ % & ' ( ) * +, -.! & " / 0 4 5 6 7 8 - % & ' ( 9 : -.! ; ( < = -
シミュレーション研究の流れ今日の目標 & 1 # 0 2 3 -! " # $ % & ' ( ) * +, -.! & " / 0 4 5 6 7 8 - % & ' ( 9 : -.! ; ( < = -
シミュレーション研究の流れ今日の目標 & 1 # 0 2 3 -! " # $ % & ' ( ) * +, -.! & " / 0 4 5 6 7 8 - % & ' ( 9 : -.! ; ( < = -
シミュレーション研究の流れ今日の目標 & 1 # 0 2 3 -! " # $ % & ' ( ) * +, -.! & " / 0 4 5 6 7 8 - % & ' ( 9 : -.! ; ( < = -
例題 ( 熱伝導問題 ) の設定 差分法の基礎 シミュレーションプログラミングの基礎
0
0
0
一様な加熱 0
" # Ly $ Lx #!
L x L y 平方メートルの長方形領域辺上の温度は常に0 ( 固定 ) 面内に熱源が分布面の熱拡散率 k 面内の温度分布は?
T (x, t) x
T x = 2 T 2 = 0 x 2 T 2 x > 0 2 T 2 x < 0 T t = 0 T t > 0 T t < 0
温度 T (x, t) に対する基本方程式 k: 熱拡散係数 T t = k 2 T x 2
T t = s s : 熱源 (heat source) ろうそくで温度計を熱している図
1D 2D T (x, y, t) t T (x, t) t = k = k 2 T (x, t) x 2 + s(x) { 2 } T (x, y, t) x 2 + 2 T (x, y, t) y 2 + s(x, y)
T (x, y, t) t 2 : ラプラシアン ( 2 = k x 2 + 2 y 2 ) T (x, y, t) + s(x, y) T (x, y, t) = k 2 T (x, y, t) + s(x, y) t T t = k 2 T + s
(0, 0) (x, y) (L x, L y ) の長方形領域で T (0, y) = T (L x, y) = T (x, 0) = T (x, L y ) = 0 という境界条件のもと T (x, y, t) t という熱拡散方程式を解き = k 2 T (x, y, t) + s(x, y) 定常状態 ( T t T (x, y) を求めよ = 0) の温度分布
離散化とはなにか? なぜ離散化が必要か? 解くべき方程式が分かっているのなら 計算機に放りこめばそれを解いてくれるのではないか? そういう時もある
解ける? du dx = x2 + u 2
解ける
数式処理ソフトウェア
解析的に ( 紙とペンで ) 解けないか? 数式処理ソフトで解けないか? 数値的に解く あるいは計算機シミュレーションをする
計算機の 計算力 を活かして熱拡散方程式を数値的に解く T (x, y, t) t = k 2 T (x, y, t) + s(x, y)
スーパーコンピュータ 何が スーパー か? 計算機のスピード競争 Top500 リスト
単位はFLOPS Floating Point Operations Per Second 1 秒あたりの倍精度浮動小数点数演算の演算回数倍精度浮動小数点数 : 15 桁ほどの小数どんな演算? 四則演算 (* + / )
計算機の速度を活かす計算機の四則演算の速度を活かす対象問題を四則演算の問題に焼き直す 離散化
T (x, y, t) t = k 2 T (x, y, t) + s(x, y) をどう四則演算の問題として表現するか?
T t を四則演算で表現する ( 一階の微分 ) 2 T x 2 ( 二階の微分 )
微分は接線の傾き T (t) t
T (t +Δt) T (t) T (t) t t +Δt t 傾き T (t+ t) T (t) t
引き算 1 回 割り算 1 回で微分を近似 T t T (t + t) T (t) t 1/ t をあらかじめ計算しておけば 割り算の代わりに掛け算
関数 の x = 2 における微分 f(x) = 1 + x 1 + x2 2 + x3 3 + x4 4 df dx を差分で近似的に計算するコードを書く
src-work/001.f95 を編集 自分で一から書いても可 正解例は /tmp/ss2/tue/src/001.f95
ポイント df(x) dx (f(x + x) f(x)) / x
コンパイル pgf95 -o 001.exe 001.f95 実行./001.exe
001.f95
x = 2 での微分の値を求めているがこれを他の x に変更して試してみる x を小さくするとごさが小さくなることを確認せよ
t t
Ny 3 2 Δy 1 0 0 1 2 3 Δx Nx
T (x, t) T (x8) T (x7) T (x6) x6 x7 x8 x
T (x, t) T (x8) T (x7) T (x6) Δx x6 x7 x8 x
[ ] d 2 dx 2 T (x) x 7 = = [ ( )] d dt dx dx x ( 7 dt ) dx 7.5 ( ) dt dx x ( T (x8 ) T (x 7 ) x ) 6.5 ( T (x7 ) T (x 6 ) 7.5 x = x = T (x 8) 2 T (x 7 ) + T (x 6 ) ( x) 2 ) 6.5
d 2 T (x) dx 2 T (x + x) 2 T (x) + T (x x) x 2
001.f95 にならい 同じ関数 のx = 2における2 階微分グラムを作れ f(x) = 1 + x 1 + x2 2 + x3 3 + x4 4 df dx (x = 2) の値を差分で計算するプロ ファイル名は 002.f95 とすること コンパイルは pgf95 -o 002.exe 002.f95 実行は./002.exe x を小さくすると誤差が小さくなることを確認せよ
d 2 T (x) dx 2 T (x + x) 2 T (x) + T (x x) x 2 d 2 T (x i ) dx 2 T (x i+1) 2 T (x i ) + T (x i 1 ) x 2 2 T (x i, y j ) x 2 T (x i+1, y j ) 2 T (x i, y j ) + T (x i 1, y j ) x 2 2 T (x i, y j ) y 2 T (x i, y j+1 ) 2 T (x i, y j ) + T (x i, y j 1 ) y 2
(xi, yj+1) (xi-1, yj) (xi, yj) (xi+1, yj) (xi, yj-1)
T (x, y, t) t ( 2 = k x 2 + 2 y 2 ) T (x, y, t) + s(x, y) T (x i, y j, t) T i,j (t) と略記 ( T i,j (t + t) T i,j (t) Ti+1,j (t) 2 T i,j (t) + T i 1,j (t) = k t x 2 + T ) i,j+1(t) 2 T i,j (t) + T i,j 1 (t) y 2 +s i,j
T (x i+1,y j ) 2 T (x i,y j )+T (x i 1,y j ) x 2, T (x i,y j+1 ) 2 T (x i,y j )+T (x i,y j 1 ) y 2 Ti,j+1 Ti-1,j Ti,j Ti+1,j Ti,j-1
T i,j (t + t) = T i,j (t) + k t x 2 {T i+1,j(t) 2 T i,j (t) + T i 1,j (t)} + k t y 2 {T i,j+1(t) 2 T i,j (t) + T i,j 1 (t)} +s i,j t この式は t だけ未来の温度 = 現在の温度分布の四則演算 という形をしている
格子系の上に定義された温度分布の時間変化 t = t0 t = t0 + Δt t = t0 + 2Δt
格子系 具体的な実体 ( モノ ) として考える カプセル化 格子系 以外のモノは 格子系 にメッセージを送る 他のモノは 格子系 の中身を直接操作することも ( 見ることも ) できないようにするのが基本 バグ防止 コード開発 メンテナンスの容易さ 大規模なシミュレーションコード開発には必須
$HOME/Tue/src-work/003.f95 エディタで 003.f95 を確認 実行 pgf95 -o 003.exe 003.f95 &&./003.exe
003.f95
003.f95 を自由に編集し コンパイル 実行せよ nx=-2 など 負の値を入れてみよ
assert 前提条件 ( 絶対に正しいはずの条件 ) を念の為に確認する 万が一条件違反 プログラム停止 走り続けておかしな結果に悩むよりも 潔く止まってくれた方がマシ
ut モジュール ( ユーティリティ ) どんどん充実させる ut assert
module ut use constants implicit none contains subroutine ut assert(condition, message) logical, intent(in) :: condition character(len=*), intent(in) :: message if (.not.condition ) then stop message! print out the message and die. end if end subroutine ut assert end module ut
そのモジュールの public なルーチン ルーチン名にアンダースコア二つ ( ) をつける モジュール名 名前 アンダースコア二つ ( ) がないルーチンは いま見ているモジュール内部の private なルーチンとわかる ファイルを見つけやすい
解答例 : $HOME/Tue/src-work/004.f95 自由に編集 コンパイル 実行 Ny 3 2 Δy 1 0 0 1 2 3 Δx Nx
004.f95
T i,j (t + t) = T i,j (t) + k t x 2 {T i+1,j(t) 2 T i,j (t) + T i 1,j (t)} + k t y 2 {T i,j+1(t) 2 T i,j (t) + T i,j 1 (t)} +s i,j t Ti,j+1 Si,j+1 Ti-1,j Ti,j Ti+1,j Si-1,j Si,j Si+1,j Ti,j-1 Si,j-1
温度場 tmp 熱源場 src 解答例 : $HOME/src-work/005.f95
005.f95
ルーチン名 iset_tmp T (x, y) = sin (xπ/l x ) sin (yπ/l y ) T(x, t) y x
t = t0 t = t0 + Δt t = t0 + 2Δt
T i,j (t + t) = T i,j (t) + k t x 2 {T i+1,j(t) 2 T i,j (t) + T i 1,j (t)} + k t y 2 {T i,j+1(t) 2 T i,j (t) + T i,j 1 (t)} +s i,j t
T i,j (t + t) = T i,j (t) + k t x 2 {T i+1,j(t) 2 T i,j (t) + T i 1,j (t)} + k t y 2 {T i,j+1(t) 2 T i,j (t) + T i,j 1 (t)} +s i,j t { k = t x 2 (T i+1,j + T i 1,j ) + k } y 2 (T i,j+1 + T i,j 1 ) + s i,j { ( k + 1 t 2 x 2 + k )} y 2 T i,j
とする α x = α y = k x 2 k y 2 β = 2(α x + α y )
T i,j (t + t) = t {α x (T i+1,j (t) + T i 1,j (t)) +α y (T i,j+1 (t) + T i,j 1 (t)) +s i,j } + (1 β t) T i,j (t) = ( 重み1)x( 周囲の温度 ) +( 重み2)x( 自分の温度 ) +( 加熱の効果 )
重みつき平均値 f = (1 w)f(a) + wf(b) 0 w 1ならば f は f(a) と f(b) の間それ以外だと ギザギザの関数の場合 振動しながら発散
拡散 ( ラプラシアン ) の効果 = 周囲の値の平均値に近づこうとする効果 Ti,j+1 Ti-1,j Ti,j Ti+1,j Ti,j-1
t 時間刻み幅 tは大きいほうが望ましい必要な計算ステップ数 ( 計算時間 ) が小さくなるのでだがむやみに大きくとれるはずはない限界値 t critical がある CFL 条件 クーラン条件 今の場合 1 β t > 0
解いている基本方程式に依存 用いている計算手法 ( 今の場合は差分法 ) に依存
差分法で熱拡散方程式を解くコードを作る解答例 : src-work/006.f95 今 熱源分布はなし ; src(:,:) = 0 最終的には温度は全ての点でゼロ T(:,:) = 0 になるはず
pgf95 -Mbounds 006.f95 配列の境界違反をチェック tmp(i,j) で i>nx のときなど
006.f95
熱源分布を 0 以外にする 時間刻み幅 t を変える t critical よりも大きくしたらどうなるか? コンパイル & 実行 : pgf95 -o 006.exe 006.f95 &&./006.exe