階層的地球流体スペクトルモデル集
SPMODEL
佐々木洋平 (京大・数学), 竹広真一 (京大・数理研)
with地球流体電脳倶楽部
SPMODEL プロジェクト
dcmodel プロジェクト
davis プロジェクト

2019 年 3 月 29 日

佐々木 & 竹広 with 地球流体電脳倶楽部

dcmodel/davis チュートリアル

2019/03/29

1 / 24

本日のお題

数式を書くようにプログラミング!!
しかもスペクトル法で!!
さらに簡単に計算データの出力 & 描画を!!

佐々木 & 竹広 with 地球流体電脳倶楽部

dcmodel/davis チュートリアル

2019/03/29

2 / 24

例題: 一次元移流方程式∂ζ∂ζ= −c∂t∂x
一次元周期境界条件: ζ(x + L, t) = ζ(x, t)▶ 初期条件: ζ(x, t) = ζ0 (x)▶ ちなみに解析解は ζ(x, t) = ζ0 (x − ct)▶ サンプルプログラムその 1: advect.f90▶

佐々木 & 竹広 with 地球流体電脳倶楽部

dcmodel/davis チュートリアル

2019/03/29

3 / 24

まずは使ってみよう!
▶

コンパイル% spmfrt -o advect.out advect.f90⇒advect.out が生成される▶

実行% ./advect.out⇒advect.nc が出力される▶

結果表示% gplist advect.nc% gpview --anim t advect.nc@zeta% gpview --range 0:1.5e3 --anim t --Gaw advect.nc@zeta

佐々木 & 竹広 with 地球流体電脳倶楽部

dcmodel/davis チュートリアル

2019/03/29

4 / 24

用いているテクニックとライブラリ▶

Fortran(>=90) の配列演算記法: 配列を返す関数
DO I=1,IMA(I) = B(I) + C(I)
ENDDO

⇒

A = B + C

DO I=1,IMDATA(I) = EXP(-X(I)**2)
ENDDO

⇒

DATA = EXP(-X**2)

spmodel library(spml): 正/逆変換, 空間微分などを提供▶ 結果出力: gtool5 library▶ 結果表示: Dennou Ruby 製品群 (GPhys, gpview · · ·)▶

佐々木 & 竹広 with 地球流体電脳倶楽部

dcmodel/davis チュートリアル

2019/03/29

5 / 24

スペクトル法による数値計算
1. 境界条件を満たす関数系で展開
ζ(xj , t) =K∑iζ̃k (t)e

2πkxj
L

2πkxj
1∑
ζ̃k (t) =ζ(xj , t)e−i L
JJ−1,

j=0k=−K

2. 方程式に代入すると常微分方程式になる
dζk2πk= −ciζ̃kdtL3. 適当な初期値から離散化して時間積分ζ̃k (t + ∆t) = ζ̃k (t) − ci

2πkζ̃k ∆tL

4. 実空間の変数に戻す
ζ(xj , t) =K∑ζ̃k (t)ei

2πkxj
L,k=−K
佐々木 & 竹広 with 地球流体電脳倶楽部

dcmodel/davis チュートリアル

2019/03/29

6 / 24

数値計算の手順
1. 使用するモジュールの宣言
use ae_module

2. 変数宣言real(kind(1.0D0)) :: g_Zeta(0:im-1) ! 格子点データ
real(kind(1.0D0)) :: e_Zeta(-km:km) ! スペクトルデータ

3. スペクトル変換の初期化
call ae_Initial(im,km,xmin,xmax)

4. 初期値を与える
g_Zeta = U1*sech((g_X-X1)/sqrt(12/u1))**2& + U2*sech((g_X-X2)/sqrt(12/u2))**2&

5. スペクトルデータに変換
e_Zeta = e_g(g_Zeta)

6. スペクトルデータを時間積分
e_Zeta = e_Zeta + dt*( - c *e_Dx_e(e_Zeta))

7. 実空間データに戻して出力
g_Zeta = g_e(e_Zeta)

佐々木 & 竹広 with 地球流体電脳倶楽部

dcmodel/davis チュートリアル

2019/03/29

7 / 24

spml/ae_module
スペクトル計算のためのサブルーチン・関数を提供
初期化▶

subroutine ae_Initial(im, km, xmin, xmax)▶

スペクトル正変換/逆変換
e_g(g_Data)g_e(e_Data)▶

! 格子点 → スペクトルデータ! スペクトルデータ → 格子点データ

微分計算
e_Dx_e(e_Data) ! X 微分▶

積分・平均演算Int_g(g_Data) ! 全領域積分
Avr_g(g_Data) ! 全領域平均

佐々木 & 竹広 with 地球流体電脳倶楽部

dcmodel/davis チュートリアル

2019/03/29

8 / 24

spml/ae_module▶

その中身は?
ISPACK(石岡, 2015) を Fortran90 の関数でくるんだもの
FFT 用の変換テーブル, 領域の大きさ → 微分計算に使用• 座標変数の設定: g_X, g_X_weight• 微分演算 → 波数のかけ算: e_Dx_e(k) = - 2*pi*(k/L)*e_Data(-k)••▶

御利益: 数式の如くプログラムが書ける
f=

∂ζ
∂x⇒

e_f = e_Dx_e(e_Zeta)

f=

∂ 2ζ
∂x2

⇒

e_f = e_Dx_e(e_Dx_e(e_Zeta))

佐々木 & 竹広 with 地球流体電脳倶楽部

dcmodel/davis チュートリアル

2019/03/29

9 / 24

spml プログラミング書法

先頭の文字で変数の種類を区別:▶ 実空間格子点データ: g_ で始める
▶ スペクトルデータ: e_ で始める
▶ spml の関数名の命名規則:(出力の型)_(機能)_(入力の型)▶ 縮約の如くプログラムを書けば型を間違えにくい▶

g_Data = g_e(e_Dx_e(e_g(g_Data1)))

佐々木 & 竹広 with 地球流体電脳倶楽部

dcmodel/davis チュートリアル

2019/03/29

10 / 24

練習問題 (1)▶

移流方程式に拡散項を付け加えてみよう∂ζ∂ζ ∂ 2ζ= −c + 2∂t∂x ∂x:::::▶

Euler scheme だとこんな感じ..?
e_Zeta = e_Zeta&& + dt * ( - c * e_Dx_e(e_Zeta)&&+ e_Dx_e(e_Dx_e(e_Zeta)) )佐々木 & 竹広 with 地球流体電脳倶楽部

dcmodel/davis チュートリアル

2019/03/29

11 / 24

非線型項の取り扱い (変換法)∂ζである場合には∂xスペクトル変数を一度実空間に戻してから積を取る▶

移流項が非線型 ζ▶

spml の関数を用いると一行で書ける
g_e(e_Zeta) * g_e(e_Dx_e(e_Zeta))

佐々木 & 竹広 with 地球流体電脳倶楽部

dcmodel/davis チュートリアル

2019/03/29

12 / 24

練習問題 (2)▶

advect_lf.f90 を元に KdV 方程式のプログラムを書いてみよう∂ζ∂ζ
∂ 3ζ
= −ζ
− 3
∂t
∂x ∂x:::▶

Leap frog scheme だとこんな感じ..?

e_Zeta = e_Zeta0&& + 2 * dt * ( - e_g(g_e(e_Zeta1) * g_e(e_Dx_e(e_Zeta1))) &&- e_Dx_e(e_Dx_e(e_Dx_e(e_Zeta1)))) &e_Zeta0 = e_Zeta1 ; e_Zeta1 = e_Zeta▶

切断波数の取り方に注意
変数の積から高波数成分が生成される⇒J∼2K+1 の格子点数では十分に表現できない• 2 次の非線型項なら格子点数は J ≥ 3K+1 に増やしておく•佐々木 & 竹広 with 地球流体電脳倶楽部

dcmodel/davis チュートリアル

2019/03/29

13 / 24SPMODEL プログラミング
1. 領域と境界条件に適合したモジュールを選択
2. 支配方程式を形式的にスペクトル変換[] []∂ ζ̃∂ 3ζ̃m∂ ζ̃m−=− ζ∂t∂x∂x3
3. 適当なスキームで時間に関して差分化{ [ττ +1+ ∆t ×ζ̃m= ζ̃m] []}∂ ζ̃∂ 3ζ̃m− ζ−∂x∂x3

4. 2,3 に従ってプログラムを書く
e_Zeta = e_Zeta& + dt * ( - e_g(g_e(e_Zeta) * g_e(e_Dx_e(e_Zeta)))- e_Dx_e(e_Dx_e(e_Dx_e(e_Zeta)))&佐々木 & 竹広 with 地球流体電脳倶楽部

dcmodel/davis チュートリアル&&)

2019/03/29

14 / 24

もっとおすすめの方法▶

サンプル・デモプログラムを元に改造
計算したい領域形状に即したモジュールのサンプル(たいてい) 拡散方程式用のプログラムが用意されている• 必要に応じて変数を追加, 演算を変更• 時間積分部分を改造••▶

サンプルの置き場所:•Web: http://www.gfd-dennou.org/library/spmodel/

佐々木 & 竹広 with 地球流体電脳倶楽部dcmodel/davis チュートリアル

2019/03/29

15 / 24

一般的注意▶

ライブラリ・サンプルプログラムは無保証••▶

必ず自分でテストしてみることできれば解析解と比較してみる

マニュアルを見よう••

気分でプログラムを書くのは危険マニュアルで確認しつつ書くこと
spml のマニュアルを見てみよう!
モジュールと計算領域・境界条件に注目!
−−▶

local: http://localhost/spmodel/reference/
Web: http://www.gfd-dennou.org/library/spmodel/

できればソースも見てみよう

佐々木 & 竹広 with 地球流体電脳倶楽部

dcmodel/davis チュートリアル

2019/03/29

16 / 24

出力操作: gtool5 library▶

モジュールの宣言
use gtool_history▶

出力ファイルの作成、次元の定義
call HistoryCreate(...

▶

変数の定義
call HistoryAddVariable(...

▶

出力call HistoryPut(...

▶

終了call HistoryClose

佐々木 & 竹広 with 地球流体電脳倶楽部

dcmodel/davis チュートリアル

2019/03/29

17 / 24

練習問題 (3): 出力変数の追加▶

advect.f90 を使っててみる▶ 変数定義の追加∂ζを x, t の 2 次元データとして出力し∂xcall HistoryAddVariable(&
varname='dzetadx', dims=(/'x','t'/), &longname='derivative of displacement',&units='1', xtype='double')▶

出力ルーチンの追加
call HistoryPut('dzetadx',g_e(e_Dx_e(e_Zeta)))

佐々木 & 竹広 with 地球流体電脳倶楽部

dcmodel/davis チュートリアル

2019/03/29

18 / 24

練習問題 (4): 出力変数の追加∫ L
ζ 2 dx を t の 1 次元データとして出力してみる▶
0•

変数定義の追加 (次元に注意)
call HistoryAddVariable(&varname='z2int', dims=(/'t'/),&longname='Integral of square of zeta', &units='1', xtype='double')•

出力ルーチンの追加
call HistoryPut('z2int',Int_g(g_e(e_Zeta)*g_e(e_Zeta)))•gpprint コマンドで数字を出力してみる% gpprint advect.nc@z2int

佐々木 & 竹広 with 地球流体電脳倶楽部

dcmodel/davis チュートリアル

2019/03/29

19 / 24

例題: 2 次元拡散方程式▶

2 次元周期境界条件の下で解いてみる( 2)∂ζ∂ζ∂ 2ζ=ν+∂t∂x2 ∂y2▶

サンプルプログラム: diffuse_2d.f90
コンパイルして実行してみよう
ソースファイルを眺めてみよう

佐々木 & 竹広 with 地球流体電脳倶楽部

dcmodel/davis チュートリアル

2019/03/29

20 / 24

gpview の便利なオプション▶

範囲指定% gpview --range 0:1 --anim t diffuse_2d.nc@zeta▶

データ切出し% gpview --anim t diffuse_2d.nc@zeta,x=0.5% gpview diffuse_2d.nc@zeta,x=0.5,y=0.5▶

平均操作% gpview --mean x --anim t diffuse_2d.nc@zeta▶

その他のオプション% gpview --help% man gpview

佐々木 & 竹広 with 地球流体電脳倶楽部

dcmodel/davis チュートリアル

2019/03/29

21 / 24

その他の例題

2 次元ベータ面: モドン▶ 2 次元チャネル領域: 熱対流問題▶ 2 次元球面領域: ロスビー波···▶ Web を見てみよう
▶

•

Web:
http://www.nagare.or.jp/mm/2006/spmodel/
http://www.gfd-dennou.org/library/spmodel/

佐々木 & 竹広 with 地球流体電脳倶楽部

dcmodel/davis チュートリアル

2019/03/29

22 / 24

練習問題 (5)▶

移流方程式に変えてみよう∂ζ∂ζ∂ζ= −cx− cy∂t
∂x
∂y

▶

Euler scheme だとこんな感じ?
ee_Zeta = ee_Zeta& + dt *( - cx * ee_Dx_ee(ee_Zeta)&- cy * ee_Dy_ee(ee_Zeta)▶

&&)

Leap frog scheme の方が安定性がよろしい
ee_Zeta = ee_Zeta0&& + dt * ( - cx * ee_Dx_ee(ee_Zeta1)&- cy * ee_Dy_ee(ee_Zeta1) )&ee_Zeta0 = ee_Zeta1; ee_Zeta1=ee_Zeta2

佐々木 & 竹広 with 地球流体電脳倶楽部

dcmodel/davis チュートリアル

2019/03/29

23 / 24

練習問題 (6)▶

線型β 面順圧方程式に変えてみよう∂ψ
∂ζ
= −β,
∂t
∂x

▶

∇2 ψ = ζEuler scheme だとこんな感じ?
ee_Zeta = ee_Zeta + dt * ( - beta*ee_Dx_ee(ee_Psi) )ee_Psi = ee_LaplaInv_ee(ee_Zeta)▶

(余力があれば) 非線型ベータ面順圧方程式に挑戦∂ζ∂ψ= −J(ψ, ζ) − β,∂t∂x

佐々木 & 竹広 with 地球流体電脳倶楽部

dcmodel/davis チュートリアル∇2 ψ = ζ

2019/03/29

24 / 24