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 の配列関数機能を活用
最新版は spml-0.2.2
• http://www.gfddennou.org/arch/spmodel/spml.tar.gz
 対応環境・コンパイラ



x86 Linux: Fujitsu frt, Intel ifc (6.0, 7.0, 8.0)
Fujitsu VPP
NEC SX （移植作業中）
サンプルプログラム
 GFD の標準的な方程式系のモデルを整備

1 次元モデル
• 移流拡散方程式, 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 次元線形移流拡散方程式


 2
U
D 2
t
x
x
 コード例 (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
ソースコード例：時間積分（２）
 球面浅水方程式（渦度発散型）

1
      2 J  ,   0
t
R
D
1
      2 J  ,    2 g h  hs   E 
t
R
h
1
   h   2 J h,   0
t
R
 渦度方程式は 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 を用いる場合
多少複雑だが十分読めるコードが書ける
1 ght  h'
2
2
n1


 1  g ht 22 h'n1 2tF
n
1


n 1
F     h'    2 J (h' , )  hD
R


1

2 
 ht      2 J  ,    E 
R


n
ソースコード例：時間積分（２）
1 ght  h'
2
2
n1


 1  g ht 22 h'n1 2tF
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
ソースコード例：時間積分（２）
n
1


F     h'    2 J (h' , )  hD n 1
R


1

2 
 ht      2 J  ,    E 
R


n
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 日モデル時間計算
• Δt=1800 sec, 720 ステップ（T106 はΔt=900 sec)


表示は hh:mm:ss
まずまずの実行速度
Xeon
2.4GHz
T21
T42
T63
T106
00:00:03
00:00:33
00:01:09
00:13:52
00:00:03
00:00:09
----------
00:01:53
VPP800
SPMODEL のご利益
 数式のような数値コードが簡単に書ける
 書法が決まっている
• コード作成時のバグも出にくい

変数配列の添え字を管理しなくてよい
• do ループは時間積分ループだけ

改変も簡単
 ポイント
 すべての変数が同じ格子点上にある
 交互格子の差分モデルでは格子点情報も必要
• 中野さんの Grid Modeling System を参照
gt4f90io のご利益
 手間のかかる入出力手続きから解放


4 つのサブルーチンだけ覚えればよい
netCDF ライブラリを直に使うよりずっと簡単
SPMODEL の課題
 差分法と併用する場合
 ex). 球面プリミティブ方程式モデル
 Grid Modeling System を使う？
• 変数配列に構造体は使いたくないなあ…

配列添え字を陽に扱う？
• 安直だけど可読性は良くないだろう

差分しない
• 全部スペクトル, 使い物になるだろうか？
 コードを作成しながら検討を予定
今後の予定
 ドキュメントの整備

サンプルプログラムマニュアル
• 1 次元サンプルプログラムのマニュアルのみ完成

チュートリアル
• プログラム書法など
 3 次元モデルの開発と試験


球面プリミティブ方程式モデル
球面 MHD 方程式モデル
参考 URL
 SPMODEL

http://www.gfd-dennou.org/arch/spmodel/
 ISPACK

http://www.gfd-dennou.org/arch/ispack/
 gt4f90io

http://www.gfd-dennou.org/arch/gtool4/
おまけ
 debian パッケージあります
 /etc/apt/sources.list に以下の 2 行を追加
deb ftp://www.gfd-dennou.org/arch/spmodel/debian woody/
deb-src ftp://www.gfd-dennou.org/arch/spmodel/debian woody/
 apt でインストール
# apt-get install spml


ISPACK, netCDF のパッケージも用意
詳細は
http://www/gfddennou.org/arch/spmodel
メモ

Download Report