第2回シミュレーションスクール - 第2日: 熱拡散方程式のプログラムをつくろう

Size: px

Start display at page:

Download "第2回シミュレーションスクール - 第2日: 熱拡散方程式のプログラムをつくろう"

えりかよしくに
7 years ago
Views:

1 第 2 回シミュレーションスクール第 2 日 : 熱拡散方程式のプログラムをつくろう陰山聡政田洋平神戸大学大学院システム情報学研究科計算科学専攻

2 具体的な例題を一つに絞る 2 次元熱伝導問題実習中心各自自分のコードを 1 から作るそのコードを少しづつ改良していく最後は並列化されたコードをスーパーコンピュータに

3 初歩からゆっくり丁寧に網羅的な紹介よりも基礎的な手法と考え方の習得を具体的な例を通じてシミュレーションプログラムの作り方を体験的に学ぶ

4 実習中心各自がプログラムを書く C/C++ 言語で書いてもOK file 名のつけ方は例題に習って ( 例 :001.f c )

5 今日と明日の講義ではMakefileは使わない使ったほうが便利実用的毎回コンパイルコマンドを自分で打つほうが教育的 Makefileの文法の問題にとらわれがち

6 コードプログラム時間発展 2 次元 2-D 2D ディレクトリフォルダ解析的に解析解

7 マシン scalar にログイン cd mkdir Tue mkdir Wed cd Tue mkdir data cp -r /tmp/ss2/tue/src-work.

8 実習用 ( スケルトンコードあり一部空白 ) 各自のデータを出力する場所実習用コードのファイル名 : 数字.f95 実習問題の解答例

9 シミュレーション研究の流れ & 1 # ! " # $ % & ' ( ) * +, -.! & " / % & ' ( 9 : -.! ; ( < = -

10 シミュレーション研究の流れ今日の目標 & 1 # ! " # $ % & ' ( ) * +, -.! & " / % & ' ( 9 : -.! ; ( < = -

11 シミュレーション研究の流れ今日の目標 & 1 # ! " # $ % & ' ( ) * +, -.! & " / % & ' ( 9 : -.! ; ( < = -

12 シミュレーション研究の流れ今日の目標 & 1 # ! " # $ % & ' ( ) * +, -.! & " / % & ' ( 9 : -.! ; ( < = -

13 シミュレーション研究の流れ今日の目標 & 1 # ! " # $ % & ' ( ) * +, -.! & " / % & ' ( 9 : -.! ; ( < = -

14 例題 ( 熱伝導問題 ) の設定差分法の基礎シミュレーションプログラミングの基礎

16 0

17 0

18 0

19 一様な加熱 0

20 " # Ly $ Lx #!

21 L x L y 平方メートルの長方形領域辺上の温度は常に0 ( 固定 ) 面内に熱源が分布面の熱拡散率 k 面内の温度分布は?

23 T (x, t) x

27 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

28 温度 T (x, t) に対する基本方程式 k: 熱拡散係数 T t = k 2 T x 2

29 T t = s s : 熱源 (heat source) ろうそくで温度計を熱している図

30 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 T (x, y, t) y 2 + s(x, y)

31 T (x, y, t) t 2 : ラプラシアン ( 2 = k x 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

32 (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) の温度分布

33 離散化とはなにか? なぜ離散化が必要か? 解くべき方程式が分かっているのなら計算機に放りこめばそれを解いてくれるのではないか? そういう時もある

34 解ける? du dx = x2 + u 2

35 解ける

36 数式処理ソフトウェア

37 解析的に ( 紙とペンで ) 解けないか? 数式処理ソフトで解けないか? 数値的に解くあるいは計算機シミュレーションをする

38 計算機の計算力を活かして熱拡散方程式を数値的に解く T (x, y, t) t = k 2 T (x, y, t) + s(x, y)

39 スーパーコンピュータ何がスーパーか? 計算機のスピード競争 Top500 リスト

41 単位はFLOPS Floating Point Operations Per Second 1 秒あたりの倍精度浮動小数点数演算の演算回数倍精度浮動小数点数 : 15 桁ほどの小数どんな演算? 四則演算 (* + / )

42 計算機の速度を活かす計算機の四則演算の速度を活かす対象問題を四則演算の問題に焼き直す離散化

43 T (x, y, t) t = k 2 T (x, y, t) + s(x, y) をどう四則演算の問題として表現するか?

44 T t を四則演算で表現する ( 一階の微分 ) 2 T x 2 ( 二階の微分 )

45 微分は接線の傾き T (t) t

46 T (t +Δt) T (t) T (t) t t +Δt t 傾き T (t+ t) T (t) t

47 引き算 1 回割り算 1 回で微分を近似 T t T (t + t) T (t) t 1/ t をあらかじめ計算しておけば割り算の代わりに掛け算

48 関数の x = 2 における微分 f(x) = 1 + x 1 + x2 2 + x3 3 + x4 4 df dx を差分で近似的に計算するコードを書く

49 src-work/001.f95 を編集自分で一から書いても可正解例は /tmp/ss2/tue/src/001.f95

50 ポイント df(x) dx (f(x + x) f(x)) / x

51 コンパイル pgf95 -o 001.exe 001.f95 実行./001.exe

52 001.f95

53 x = 2 での微分の値を求めているがこれを他の x に変更して試してみる x を小さくするとごさが小さくなることを確認せよ

54 t t

55 Ny 3 2 Δy Δx Nx

56 T (x, t) T (x8) T (x7) T (x6) x6 x7 x8 x

57 T (x, t) T (x8) T (x7) T (x6) Δx x6 x7 x8 x

58 [ ] 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

59 d 2 T (x) dx 2 T (x + x) 2 T (x) + T (x x) x 2

60 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 を小さくすると誤差が小さくなることを確認せよ

61 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

62 (xi, yj+1) (xi-1, yj) (xi, yj) (xi+1, yj) (xi, yj-1)

63 T (x, y, t) t ( 2 = k x 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

64 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

65 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 だけ未来の温度 = 現在の温度分布の四則演算という形をしている

66 格子系の上に定義された温度分布の時間変化 t = t0 t = t0 + Δt t = t0 + 2Δt

67 格子系具体的な実体 ( モノ ) として考えるカプセル化格子系以外のモノは格子系にメッセージを送る他のモノは格子系の中身を直接操作することも ( 見ることも ) できないようにするのが基本バグ防止コード開発メンテナンスの容易さ大規模なシミュレーションコード開発には必須

68 $HOME/Tue/src-work/003.f95 エディタで 003.f95 を確認実行 pgf95 -o 003.exe 003.f95 &&./003.exe

69 003.f95

70 003.f95 を自由に編集しコンパイル実行せよ nx=-2 など負の値を入れてみよ

71 assert 前提条件 ( 絶対に正しいはずの条件 ) を念の為に確認する万が一条件違反プログラム停止走り続けておかしな結果に悩むよりも潔く止まってくれた方がマシ

72 ut モジュール ( ユーティリティ ) どんどん充実させる ut assert

73 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

74 そのモジュールの public なルーチンルーチン名にアンダースコア二つ ( ) をつけるモジュール名名前アンダースコア二つ ( ) がないルーチンはいま見ているモジュール内部の private なルーチンとわかるファイルを見つけやすい

75 解答例 : $HOME/Tue/src-work/004.f95 自由に編集コンパイル実行 Ny 3 2 Δy Δx Nx

76 004.f95

77 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

78 温度場 tmp 熱源場 src 解答例 : $HOME/src-work/005.f95

79 005.f95

80 ルーチン名 iset_tmp T (x, y) = sin (xπ/l x ) sin (yπ/l y ) T(x, t) y x

81 t = t0 t = t0 + Δt t = t0 + 2Δt

82 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

83 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

84 とする α x = α y = k x 2 k y 2 β = 2(α x + α y )

85 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( 自分の温度 ) +( 加熱の効果 )

86 重みつき平均値 f = (1 w)f(a) + wf(b) 0 w 1ならば f は f(a) と f(b) の間それ以外だとギザギザの関数の場合振動しながら発散

87 拡散 ( ラプラシアン ) の効果 = 周囲の値の平均値に近づこうとする効果 Ti,j+1 Ti-1,j Ti,j Ti+1,j Ti,j-1

88 t 時間刻み幅 tは大きいほうが望ましい必要な計算ステップ数 ( 計算時間 ) が小さくなるのでだがむやみに大きくとれるはずはない限界値 t critical がある CFL 条件クーラン条件今の場合 1 β t > 0

89 解いている基本方程式に依存用いている計算手法 ( 今の場合は差分法 ) に依存

90 差分法で熱拡散方程式を解くコードを作る解答例 : src-work/006.f95 今熱源分布はなし ; src(:,:) = 0 最終的には温度は全ての点でゼロ T(:,:) = 0 になるはず

91 pgf95 -Mbounds 006.f95 配列の境界違反をチェック tmp(i,j) で i>nx のときなど

92 006.f95

93 熱源分布を 0 以外にする時間刻み幅 t を変える t critical よりも大きくしたらどうなるか? コンパイル & 実行 : pgf95 -o 006.exe 006.f95 &&./006.exe

解析力学B - 第11回: 正準変換

解析力学B - 第11回: 正準変換解析力学 B 第 11 回 : 正準変換神戸大 : 陰山聡ホームページ ( 第 6 回から今回までの講義ノート ) http://tinyurl.com/kage2010 2011.01.27 正準変換バネ問題 ( あえて下手に座標をとった ) ハミルトニアンを考える q 正準方程式は H = p2 2m + k 2 (q l 0) 2 q = H p = p m ṗ = H q = k(q