SPMODEL - ISPACK と gt4f90io による数値モデル開発 -

Presentation on theme: "SPMODEL - ISPACK と gt4f90io による数値モデル開発 -"— Presentation transcript:

1 SPMODEL - ISPACK と gt4f90io による数値モデル開発 -
電脳 davis – gtool4 に関する ワークショップ SPMODEL - ISPACK と gt4f90io による数値モデル開発 - 小高正嗣 SPMODEL 開発グループ

2 我々の目標 理解のための数値モデルが欲しい そのような数値モデルの要件 数式のように数値モデルとその結果を理解
簡単モデルと複雑モデルの結果をつなげて理解 そのような数値モデルの要件 可読性が高い 拡張性が高い 階層化された数値モデル群の一つ

3 SPMODEL 階層的地球流体力学スペクトルモデル集 基本方針 GFD スペクトルモデルを階層的に整備
ライブラリ(spml)とサンプルプログラムから構成 基本方針 コードの可読性を重視 書法を統一 実行速度性能にはこだわらない 共通のライブラリを利用 プログラム構造を階層化 計算エンジンに ISPACK, データ I/O に gt4f90io

4 SPMODEL ライブラリ（spml） ISPACK の F90 インターフェース 対応環境・コンパイラ F90 の配列関数機能を活用
対応環境・コンパイラ x86 Linux: Fujitsu frt, Intel ifc (6.0, 7.0, 8.0) Fujitsu VPP NEC SX （移植作業中）

5 サンプルプログラム GFD の標準的な方程式系のモデルを整備 1 次元モデル 2 次元モデル 3 次元モデル
移流拡散方程式, KdV 方程式 2 次元モデル 水路領域ブシネスク方程式, 赤道β面浅水方程式 球面順圧渦度方程式, 球面浅水方程式 3 次元モデル 球殻ブシネスク方程式 球殻 MHD 方程式 プリミティブ方程式（開発予定）

6 読みやすいコードにするために 書法をきめる 手間のかかるデータ入出力は gt4f90io に 変数配列 配列関数
4 つのサブルーチンを call するだけ データの付加情報も渡すことができる netCDF　の大域属性を利用

7 SPMODEL 書法 変数 配列関数 変数にはデータの種類を示す接頭詞をつける (出力型)_(機能)_(入力型) のように書く
実空間データ g_Zeta スペクトルデータ s_Zeta 配列関数 (出力型)_(機能)_(入力型) のように書く スペクトル変換 s_Zeta=s_g(g_Zeta) 発散 s_Div_g(g_Zeta) 引数の型を間違えることが少ない

8 ソースコード例：時間積分（１） 1 次元線形移流拡散方程式 コード例 (Euler 法で時間積分) 数式のように数値コードを書ける
s_ZetaA = s_ZetaB 　 & dt*( - U*s_Dx_s(s_Zeta) + D*s_Dx_s(s_Dx_s(s_Zeta))

9 ソースコード例：時間積分（２） 球面浅水方程式（渦度発散型） 渦度方程式は leapfrog スキームで
w_EtaA=w_EtaB + 2*dt* & (- w_Div_xy_xy(xy_w(w_Eta)*xy_GradLon_w(w_Chi)/R, & xy_w(w_Eta)*xy_GradLat_w(w_Chi)/R)/R & + w_Jacobian_w_w(w_Eta, w_Psi)/R**2)

10 ソースコード例：時間積分（２） 質量保存の式 Leapforg + 台形 semi-implicit を用いる場合
多少複雑だが十分読めるコードが書ける

11 ソースコード例：時間積分（２） call get_dHsfc
w_HsfcA=((1 + dt**2*Grav*H0*n(:.1))*w_HsfcB　　　　　　　　&　 + 2*dt*w_dHsfc)/(1 - dt**2*Grav*H0*n(:.1)) … End

12 ソースコード例：時間積分（２） soubroutin get_dHsfc …
w_dHsfc=- (w_Div_xy_xy(xy_Hsfc*xy_GradLon_w(w_Chi)/R, & xy_Hsfc*xy_GradLat_w(w_Chi)/R)/R & - w_Jacobian_w_w(w_Hsfc,w_Psi)/R**2) & - H0*w_DivB & -H0*dt*(w_Div_xy_xy(xy_Eta*xy_GradLon_w(w_Psi)/R, & xy_Eta*xy_GradLat_w(w_Psi)/R)/R & + w_Jacobian_w_w(w_Eta,w_Chi)/R** & + w_Lapla_w(w_E)/R**2)

13 ソースコード例：データ入出力 use gt4history ! モジュール引用 …
call HistoryCreate( & ! 出力ファイル作成 　 file=“output_file”, title=“Run Number”, source=“shallow” & institution=“GFD Dennou Club SPMODEL project”, & dims=(/‘lon’,‘lat’,‘t ’/),dimsizes=(/im,jm,0/), & longnames=(/‘longitude’,‘latitude ’,‘time ’/), & units=(/'deg.','deg.','sec.'/), origin=real(tinit), & interval=real(intrst*delt) ) call HistoryAddVariable( varname=“eta”, & ! 出力変数定義 dims=(/’lon’,’lat’,’t ’/), logname=‘Vorticity’, units=‘1/s’,xtype=‘double’) do it =1, n ! ループ開始 call HistoryPut(‘eta’,xy_Eta) ! 変数出力 end do call HistoryClose ! 入出力終了

14 計算時の CPU 時間 球面浅水モデルの場合 15 日モデル時間計算 表示は hh:mm:ss まずまずの実行速度 T21 T42 T63
Δt=1800 sec, 720 ステップ（T106 はΔt=900 sec) 表示は hh:mm:ss まずまずの実行速度 T21 T42 T63 T106 Xeon 2.4GHz 00:00:03 00:00:33 00:01:09 00:13:52 VPP800 00:00:09 00:01:53

15 SPMODEL のご利益 数式のような数値コードが簡単に書ける ポイント 書法が決まっている 変数配列の添え字を管理しなくてよい 改変も簡単
コード作成時のバグも出にくい 変数配列の添え字を管理しなくてよい do ループは時間積分ループだけ 改変も簡単 ポイント すべての変数が同じ格子点上にある 交互格子の差分モデルでは格子点情報も必要 中野さんの Grid Modeling System を参照

16 gt4f90io のご利益 手間のかかる入出力手続きから解放 4 つのサブルーチンだけ覚えればよい
netCDF ライブラリを直に使うよりずっと簡単

17 SPMODEL の課題 差分法と併用する場合 コードを作成しながら検討を予定 ex). 球面プリミティブ方程式モデル
Grid Modeling System を使う？ 変数配列に構造体は使いたくないなあ… 配列添え字を陽に扱う？ 安直だけど可読性は良くないだろう 差分しない 全部スペクトル, 使い物になるだろうか？ コードを作成しながら検討を予定

18 今後の予定 ドキュメントの整備 3 次元モデルの開発と試験 サンプルプログラムマニュアル チュートリアル 球面プリミティブ方程式モデル
1 次元サンプルプログラムのマニュアルのみ完成 チュートリアル プログラム書法など 3 次元モデルの開発と試験 球面プリミティブ方程式モデル 球面 MHD 方程式モデル

19 参考 URL SPMODEL ISPACK gt4f90io http://www.gfd-dennou.org/arch/spmodel/
gt4f90io

20 おまけ debian パッケージあります apt でインストール /etc/apt/sources.list に以下の 2 行を追加
ISPACK, netCDF　のパッケージも用意 詳細は deb ftp:// woody/ deb-src ftp:// woody/ # apt-get install spml

21 メモ