SPMODEL - ISPACK と gt4f90io による数値モデル開発 - 電脳 davis – gtool4 に関する ワークショップ SPMODEL - ISPACK と gt4f90io による数値モデル開発 - 小高正嗣 SPMODEL 開発グループ http://www.gfd-dennou.org/arch/spmodel/

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

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

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

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

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

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

ソースコード例：時間積分（１） 1 次元線形移流拡散方程式 コード例 (Euler 法で時間積分) 数式のように数値コードを書ける s_ZetaA = s_ZetaB + 　 & dt*( - U*s_Dx_s(s_Zeta) + D*s_Dx_s(s_Dx_s(s_Zeta)) http://www.gfd-dennou.org/arch/spmodel/1d-cyclic-e/advection-diffusion/sample/f90/advdiff1.f90

ソースコード例：時間積分（２） 球面浅水方程式（渦度発散型） 渦度方程式は 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) http://www.gfd-dennou.org/arch/spmodel/2d-sphere-w/shallow/shallow-zd/f90/shallow_zd.f90

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

ソースコード例：時間積分（２） 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

ソースコード例：時間積分（２） 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**2 & + w_Lapla_w(w_E)/R**2)

ソースコード例：データ入出力 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 ! 入出力終了

計算時の 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

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

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

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

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

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

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

メモ