U - 日本海洋科学振興財団

海洋モデリング夏の学校 2013/8/23
非静力学海洋モデルの現状と課題
~個人開発の限界~
松村義正
北海道大学 低温科学研究所
自己紹介
●
2003年 数値モデル開発を志しこの道へ
●
東大気候システム研究センターで学位 (H21)
●
H23年度より北大低温研
●
昨年後半は南極へ...
●
●
●
大気も海洋も、一般に大循環モデル(GCM)と呼ばれるものはプリミ
ティブ方程式系に基づく
–
薄い球殻上の運動、水平運動が卓越
–
静水圧近似、Boussinesq近似、コリオリのcos項を無視
–
鉛直速度は予報変数から診断的に求める
例えばIPCCで引用されるような全球規模・数百年スケールの温暖化実
験等は（大気と海洋を結合した)大循環モデルによる
計算コスト優先、表現されない現象はパラメタリぜーションで導入
非静力学モデル 開発の動機
●
●
●
海洋には大小さまざまなスケールのプロセスが存在
数千年規模の大循環は乱流混合や局所的な沈み込みといったス
モールスケールプロセスによって駆動されている
現状のOGCMではスモールスケールプロセスをパラメタライズ
●
●
●
OGCMの信頼性はパラタメタリゼーションのチューニングによっ
て決まるといって過言ではない
計算機の発展に伴い, スモールスケールプロセスも陽に解像するこ
とが可能になりつつある
1kmスケールを境に鉛直運動が重要, 非静力項が無視できない
●
既存のモデルを高解像度するだけでは不十分
Non-hydrostatic
Hydrostatic
海洋非静力学モデルの現状
●
2次元モデルやスペクトルモデル(地形無し)は比較的簡単
–
個人レベルで作成・維持可能、混合層の乱流混合 etc.
●
Mahadevan et al., 1996 (JPO)
●
3次元フルモデルはMITgcm (Marshall et al., 1997)がリファレンス
●
●
–
Z座標、全球結合モデル
–
これをベースにしたモデルは多い
ROMS (Regional Ocean Model)
–
Terrain-following (σ座標), 領域モデル
–
Kanaska et al., 2006で非静力化
最近では非構造格子モデルも
–
FVCOM
–
Fringer et al. (2006), Kramer et al. (2010) etc.
日本の統合モデルを作るのであれば：
少なくとも領域モデルも対象とするなら、国際競争力の
観点から非静力オプションは絶対必要
海洋非静力学モデル
Navier-Stokes方程式
(ブジネスク近似)
圧力 p を静力学成分 phs とその他 pnh に分割
∂u
1 ∂p

=− u⋅∇  u−
fvV u

∂t
0 ∂x
∂v
1 ∂p

=− u⋅∇  v−
−fuV v

∂t
0 ∂y
∂w
⋅∇  w− 1 ∂ p −gV w 
=− u
0 ∂ z
∂t
非圧縮の連続の式
 =0
∇⋅u
水温・塩分（溶存物質)の移流拡散方程式
状態方程式 (非線形、多項式近似)
= , S , z 
phs : hydrostatic pressure
0
phs  x , y , z ≡∫z  x , y , z g dz
modified pressure
pnh
P≡
0
P以外の項は評価可能、Gにまとめる
(音速 ~1500m/s >> v)
∂
⋅∇  DQ
=− u
∂t
∂S
⋅∇  SDSQ S
=− u
∂t
p=phspnh
 
du
= G−∇ P
dt
2

∇⋅G=∇
P
∂P
=0
境界条件は:
∂n
連速の式より
タイムステップ毎に3D Poisson方程式を
解かなければならない
Poisson ソルバとは
2
q=∇ P =LP ,
●
●
●
●
●
●
直接法:
Gauss-Seidel:
SOR:
共役勾配(CG):
FFT:
Multigrid:
−1
P= L q
O(N2)
O(N2)
O(N3/2)
O(N3/2)
O(N log N)
O(N log N)
ソルバーの性能が超重要
離散化したPoisson方程式は単なる線形問題、逆行列を解くだけ
● 総格子数は > 109, 直接法は不可能
● 毎タイムステップ解かないといけない、計算コスト最優先
● 複雑な陸・海底地形形状 → FFTによるスペクトル法は無理
● 古典的反復法 (e.g., Gauss-Sidel, SOR, Conujugate Gradient) は
反復回数が格子数に依存 ＝＞ コストが指数関数的に悪化
● 自由表面 → 係数行列が時々刻々変化
● 並列化効率がとても重要 （ソルバは全域計算, 袖通信だけではダメ)
100
~
5000
m
100~10000 km
反復解法
MをLの近似逆行列として:
を残差r=q-Lpが十分小さくなるまで繰り返す
ランダムに分布している誤差を拡散・平滑化することに相当
任意の座標系・領域形状・境界条件を扱える。
●
通常の反復法:
→ 誤差の拡散, 高波数の誤差ほど効率良く平滑化される.
→ 長波長の誤差には時間がかかる
→ 格子数が増えると相対的に最大波長が増える
→ 収束までの反復回数が増加
→ 総計算コスト “一反復当たりの計算量 x 反復回数”
指数関数的に増大, 大規模実験は困難
Multi-Grid (多重格子) 法とは
●
長波長の誤差は荒い格子でも表現できる
→ 荒い格子での反復はコスト小
→ 誤差の各波長成分毎に最適な格子で拡散させる
→ 反復回数が総格子数によらない、スケーラブル
a prioriに誤差を波数毎に分解できるわけではない
→ 再帰的に格子間を移動しながら反復
●
●
マルチグッリドにしたら全て解決するわけではない:
–
地形が複雑(e.g. 急峻な海嶺、細い海峡)だと収束率悪化
–
並列化効率が必ずしもよくない(多数のレベルで何度も通信)
–
打ち切りεを甘くすると誤差が蓄積
色々工夫しました:
色々工夫しました
–
係数行列が対称なのでMGをCGの前処理として使う.
–
Using sparse approximate inverse (SAI) as multigrid
smoothers instead of Gauss-Sidel / Red-Black.
–
Add correction terms for unphysical divergence fields
induced by truncation error.
100億格子、数万並列(eg. 京)でも10反復くらいで収束!!
ソルバーコストは全体の3割程度
→ もはや「非静力」自体はコスト増にならない
B-Grid or C-Grid
C-Grid
B-Grid
P
P
v
P
P
v
u
P
u
P
P
P
P
P
P
P
Poisson eq.

∇ ⋅ ( ∂t u ) = ∇ 2 P = LP
2D : 9行対角
3D :11行対角
P
∗
2D : 5行対角
3D : 7行対角
計算効率に無視できない違い
P
鉛直差分は z-level, shaved-cell
z-level
σ-coordinate
s-coordiante
(σ+z hybrid)
shaved-cell
一般的には, 例えば底層水形成を扱うならσ座標のほうがいい,
(スムースな地形表現, 海底境界層がよく表現できる.)
標系が直交しない、圧力勾配に誤差
(高次の差分を使えば緩和できるが、ソルバーがつらくなる)
解像度が十分であれば z-level でも問題ない (Ezer2004, Legg2006, Wang2008)
ただし地形に埋まったmask領域が増えると無駄が多くなる
→ 3次元並列分割により、mask領域を効率的に計算対象から除去
Implicit Free Surface
●
●
海面昇降(外部重力波)の伝搬速度は~200m/s, 流速より2桁速い
海面昇降を陰解法により2次元Helmholtz方程式から予報
rigid­lidモデルの海面圧力診断Poisson方程式にηの時間微分が加わる
●
メリット:
モード分割が不要, タイムステップ間隔は傾圧成分を基準に
保存則はperfectに満たされる
●
●
デメリット:
–
外部モードがなまる
–
2次元とはいえ大規模な線形問題を解く必要, 並列化に不向き
採用した理由:
–
そもそもFree surface化は外部重力波を見るためではない,
–
どうせ非静水圧力に関して3次元のPoissonを解くので, Helmholtz
方程式を解くのは負担にならない, コードも流用
–
棚氷存在下では外部モードが ❝stiff­equationになる”、
粘性項の表現: サブグリッドスケールモデル
Navier-Stokes方程式
∂ ui
∂ ui
1 ∂ p ∂  ij
u j
=−
−
F i
0 ∂ x i ∂ x j
∂t
∂xj
Fi : 外力項 (コリオリ, 重力etc)
応力テンソル をどう評価するか?
→ そもそも何の時間発展を解いているのか ?
Direct Numerical Simulation (DNS)
●
流れの全成分を直接解く:
 ij ≡−2  m S ij (m: 分子粘性)
必要な格子解像度は(L/)3  (Re3/4)3
海洋物理の対象では現実問題として不可能
Reynolds Averaged Navier-Stokes (RANS)
●
解いているのはensemble mean
 ij ≡u i  u j  Reynolds Stress: 非線形項によって生じる擾乱成分の相関
2次のmomentの導出に3次のmomentが..., closure problem
高次のmomentを低次のmomentにて近似 (Mellor and Ymada, 1974 etc.)
Large Eddy Simulation (LES)
Filteringされた物理量の時間発展を解く
G x− x xdx
∫
 x=
∣∫ G xdx∣
G(x): Filter Function
{
1 if ∣x i∣ i
G x ≡
0 otherwise
 ij ≡ui u j −ui u j −2 m Sij  : subgrid-scale stress tensor
1
 ij ≡ ij −  ij  kk : deviatoric part of stress tensor
3
p 1
P ≡   kk : modified pressure
0 3
1
1
1

=
u
u
−
kk
k
k
2
2
2 u k u k = {Total K.E.}− {Resolved K.E.}
有限体積法では
= {Subgrid Scale K.E.}
∂ ui
∂ ui ∂ P
∂  ij
⇒ N-S eq.:
u j

=−
F i
∂t
∂ x j ∂ xi
∂xj
(非圧縮モデルではPはPoisson方程式から診断)
Eddy Viscosity Model

1 ∂ ui ∂ u j
 ij =−2  Sij , Sij ≡

2 ∂ x j ∂ xi

: eddy viscosity coefficient, UNKNOWN
Smagorinsky's Model (Smagorinsky, 1964)
1
 ≃ C 2  2 ∣S∣ , ∣S∣≡
Sij Sij
2

慣性小領域まで解像し、かつ等方乱流ならば局所平衡仮説から解析解が求まりC~0.2
現実的にはtuning parameter:
等方性が前提, 海洋モデルでは通常水平のみに適用
Deardorff's Model (Deardorff, 1980)
e
1/3

 ≃ C K  e ,  ≡ min C l
, [ x  y  z ]
N


e: subgrid-scale kinetic energy
3/ 2
∂ ui
∂e
∂e
∂
e
e
=−u j
−  ij
 ∂ 
−  ,  ≡ C

∂t
∂ xj
∂ xj ∂ xj
∂ xj

成層に起因する非等方性を陽に扱う
Ck, Cl, Ce: 経験定数, 適切な値は不明
transientな問題では初期のeに依存

Dynamic LES (Germano et al., 1991)
Grid filterより大きなスケールのtest filter*を考える
∗
∗
G
 x− x xdx
∫
1
if
∣x
∣
∗
i
i
  x=

,
G

x
≡
0 otherwise
∣∫ G∗  xdx∣
{
Test Filtered N-S eq. :
i
i ∂ P ∗
∂u
∂u
∂  ∗ij
j
u

=−
 F ∗i
∂t
∂ x j ∂ xi
∂ xj
where subtest-scale deviatoric stress tensor:
1
∗





 ij ≡ ui u j − ui u j −  ij  u
k uk − uk u k 
3
∣ S
ij
≃−C 2  ∗ 2 ∣S
2つのスケールと*でSmagorinsky
constant C は等しいと仮定する.
Germano Identity :
1





Lij ≡ ∗ij − ij = u
u
−
u
u
−
 ij  u
i
j
i
j
k u k − uk u k 
3
∣ S
 ≡ 2 C2 M
= 2 C2 
 2 ∣S∣ S −  ∗2 ∣S
[
ij
ij
(ただし
]
ij
i
i = u
u
に注意する)
C以外は全て既知, ただし自由度5, 過剰決定系をなす
⇒ 最小二乗法によりCを決定 (Lily, 1992)
2
2
Q ≡∣Lij − 2 C M ij∣ ,
∂Q
2
=−2
L
−
2
C
Mij  Mij = 0

ij
2
∂C
Lij Mij
⇒C =
2 M ij Mij
2
Eddy Viscosity Coefficient
Lij Mij
2
2
2
 = C  ∣S∣=
 ∣S∣
2 Mij Mij
毎ステップCの最適値が解像されている予報変数から動的に診断される
Tracer Transport (Moin et al., 1991)
Grid Filtered equation:
2
∂fi
∂T
∂T
C
∂T
2
 ui
=−
, f i ≡ ui T − ui T ≃−
 ∣S∣
∂t
∂ xi
∂ xi
Pr
∂ xi
Test Filtered equation:
∗
2



∂
f
∂T
∂T
C
∂
T
i
∗
∗
2
∣


i
u
=−
, f i ≡u
 ∣S
i T − ui T ≃−
∂t
∂ xi
∂ xi
Pr
∂ xi
Germano Identity:
2 
2

C
∂
T
∂
T
C
∗
2
∗
2

 ≃
∣

Li ≡ f i − fi = u
T
−
u
T
 ∣S∣
−  ∣S
≡
Mi
i
i
Pr
∂ xi
∂ xi
Pr
C 2 Li M i
=

Pr Mi Mi
[
Eddy Diffusivity Coefficient:
Li Mi 2
C2 2
=
 ∣S∣=
 ∣S∣
Pr
M i Mi
]
開発しているモデルの概要
●
名前はkinaco
●
3次元非静力海洋モデル
●
非圧縮系, 圧力はPoisson方程式から診断 (following Marshall et al., MITgcm)
●
水平一般曲線直交座標, 構造化格子, C-grid
●
鉛直z-level, partial cell
●
●
一部shaved-cell的表現
●
解像度が十分であれば z-levelはデメリットにならない
Implicit free surface
●
●
Multigrid-SAI-CG Poisson solver
●
●
モード分解不要, 海面昇降は2次元Helmholtz方程式から予報
計算コストが格子点数に対してスケーラブル, かつ並列化効率高
Dynamic LES subgrid-scale model
●
経験的パラメータを極力排除
南極底層水形成のシミュレーション
Maqueda et al., 2004
冷却と海氷生成が高密度陸棚水を形成、大陸斜面を下る
● 粗い解像度のOGCMではパラメタ化せざるを得ない
●
地形(~1km)の効果が最も重要
本当に非静力が重要になるのはさらに
小さいスケール
Olbers and Visbeck, 2005
●
沈降する高密度水と周りの軽い水が混合
●
鉛直シアによるKelvin-Helmholtz不安定
●
Large Eddy Simulation による理想化実験
Super­high resolution idealized experiment:
Direct Numerical Simulation (DNS*)
for the reference
•
Lower layer: DSW (θ =­2℃, S =34.6 psu)
Upper Layer: MCDW (θ =1℃, S =34.7 psu)
∆U = 50cm/s
Hyperbolic Tangent­shaped vertical profile for θ, S, U Domain: 1280m x 320m x 400m (periodic boundary)
Without Rotation •
Resolution: isotropic grid with ∆x=∆y=∆z=1.25m
•
•
•
•
Subgrid viscosity: ν=10­4 m2/s (isotropic, constant)
• Subgrid diffusivity:κ=10­5 m2/s (isotropic, constant)
•
(note: it is not exactly DNS, Re ~ 106 due to computational limitations)
•
Integrated for 20­hours using non­hydrostatic model Result of DNS*
billows of ~ 200m develops and
then merged and collapsed.
• Highly coherent in the cross­flow direction in the early stage.
• Large EKE continues for 8 hours.
• Vertical shear has layered structure of ~20m after the event.
Potential Density
•
Eddy Kinetic Energy
640m
360m
Strain Rate (shear)
θ
U
σθ
N2
バルクで見積もった鉛直渦拡散は (200 m)2 / 6h ~ O(1) m2/s
質量輸送はキャベリングにより増える
Large Eddy Simulation (LES) subgrid model for coarser resolution experiment inertial
sub range
~ k ­5/3
Results of Smagorinsky Model
C=0.2
C=0.5
C=1.0
Dynamic LES result
DNS* (= 1.25m)
D-LES (=10m x 5m)
D­LES scheme well reproduce the DNS* solution.
Note again that D­LES does not use any tuning parameter.
Sensitivity to the horizontal resolution ⊿x=10m
⊿x=20m
⊿x=40m
⊿x=80m
⊿x=160m
⊿x=320m
Requires ∆x < 40 m, at least.
Consistent with the inertial­subrange in DNS results. The Thermobaric Effect:
0db
1000db
MCDW
DSW
2000db
MCDW
DSW
深く沈むほど密度差大
→ 安定化
MCDW
DSW
High-resolution idealized experiments on constant slope
grid length is 20m x 20 m x 10 m
turbulent kinetic energy
最近追加した事
●
棚氷
– 現在の気候研究界隈で最もホットな話題
– 今後の南極沿岸モデリングには必須
– 非静力・棚氷カップルモデルはそんなにない
– 底面境界層を解像し、熱交換・融解パラメタリ
ゼーション
既存のOGCMではかなり簡略化(定数x温度差)
– 低温研氷床観測グループとコラボ
最近追加した事
●
ラグランジアン粒子追跡 (オンライン)
– 各粒子は固有IDと複数のプロパティ(質量・浮上沈降速度・
年齢etc.)を保持
– 予報流速+粒子個別の速度を用いて4次Runge­Kuttaで移流
– 双方向リンクリストデータ構造を採用：
高速・完全並列化対応
任意の場所・時刻で粒子の生成消滅
各モデル格子中に存在する粒子毎にグループ化
→ 膨大な数の粒子を扱う場合も高い検索性,
粒子間相互作用も可能
Frazil Ice (氷晶)を陽に解像する
海氷生成シミュレーション
●
既存の海氷モデルでは表層で熱が奪われると過冷却が発
生しないよう薄氷を形成
●
海氷は断熱材、薄氷が成長するに従い熱交換を遮断
●
実際はまず微小なFrazil Iceが形成、集積して海氷に。
●
風が強いなど海面が乱流状態の場合グリースアイスに
→ Frazil particle (cluster)をラグランジュ粒子として
–
質量を保持、現場の温度が結氷点以下であれば
凍結、以上であれば融解
–
粘性は粒子集積度の関数:
Slurry Formula: ν* = (1+2.1φ+14.1φ2)ν
–
粒子が上層10cmに来たら固着し薄氷へ
–
結氷点は高圧・高塩ほど下がる
Frazil Ice (氷晶)を陽に解像する
海氷生成シミュレーション
●
格子幅1m x 1m x 1m
●
格子数 64 x 64 x 64, 周期境界
●
●
θ
初期値 -1.6°C一様, 60m深に塩分躍層
(混合層 30psu)
気温-20°, Haney型 (~800W/m2)
w
Frazil Distribution
Melt rate
frazil粒子が(負の)潜熱を下層に輸送し下層がpotential supercoolingに
オンライン粒子追跡の応用
●
南極底層水トレーサー, 棚氷融解水トレーサー
–
●
堆積物からの微量元素の巻き上げ
–
●
●
海域毎の寄与を同定, 年齢も判別可能
海氷にも取り込まれ、融解時に放出
ラグランジアン海氷モデル
–
既存のレオロジーモデルは格子幅が氷板のスケール
より小さいと破綻
–
氷板一つ一つをラグランジュ粒子として扱う
–
氷板間の相互作用がレオロジー
生態系モデル
–
生態系モデル屋は生物種を増やしたい
–
トレーサー数増大、しかし密度に大きな偏り
–
異なる種間の相互作用が重要
今後の計画：
●
とにかく使える計算機リソースに依存
●
京 10.51PFlops
●
東大FX10 1.13PFlops
●
ノードあたり~200GFlops,
●
100万ノード時間くらいは確保
今後の大規模計算はO(104)並列が前提
●
(1000node x 16core x 1000hr)
---
60億格子 x 200万ステップ
PLAN 1: 領域を拡大
1kmメッシュ x 20m で南極全周を覆う
(20000 x 1000 x 200 格子)
->どこでどれだけ南極底層水ができるか明らかに
PLAN 2: 超高解像度に
Mooring Observations: Feb. 2008 — Jan. 2009
Line III
Hydrographic Obs. : 21 Jan. — 27 Jan. 2009
by H/V Umitaka; R/V Hakuho
Line I
Line II
M1
M2
M4
M3
M1-M4: Mooring
stations
High ice production area
( >5m/yr)
Cape Darnley
Polynya
PLAN 2: 超高解像度に
LESが機能し、エントレインメントが表現できる
100ｍメッシュで500km x 500kmを覆う
Oxygen
現実設定のモデリングで南極底層水の
エントレインメントを陽に解像
まとめ
●
●
海洋非静力学モデルを開発
–
多重格子Poissonソルバ、Dynamic-LES、online粒子追跡
–
O(104）並列対応, 実効性能7%以上
–
コストは静力学モデルの高々2倍
今後は個人レベルでPeta Flops計算が可能な時代
–
●
今まで不可能だと思われていた計算が可能
南極底層水形成モデリングに関して、ブレイクスルーが見える
–
他にも沿岸域の防災・水産関連など応用は幅広い
やりたいこと・できることは沢山あるが、モデル開発・
チューニングばっかりやっていると実研究に手がつかない
→ 論文書けない、死活問題
非静力海洋モデル・ソルバー単体・粒子追跡単体でも
誰か是非共同研究しましょう、計算は私がやりますので
I/O performance:
●
●
●
Read/Write whole data on a single PE is impossible
due to memory limitation.
Dividing data files into each reagion of corresponding
PE is not good solution for ~10000 parallel run
Test 3 different parallel IO schemes
1) Each PE reads/writes the corresponding area of single file
non-sequencial access, too slow.
2) gather/scatter for each vertical layer
→ reads/writes on PE#0
sequential access
3) MPI-IO
I/O Test on HA8000 8192PE
Serial­IO [s]
MPIIO [s]
1024 x 1024 x 60
16.29
20.56
1024 x 1024 x 120
27.29
17.63
1024 x 1024 x 240
51.85
24.65
2048 x 2048 x 60
49.57
16.57
2048 x 2048 x 120
101.7
69.26
2048 x 2048 x 240
199.1
178.9 ???
4096 x 4096 x 60
241.0
45.93
4096 x 4096 x 120
431.2
29.06
4096 x 4096 x 240
N/A
93.89
8192 x 8192 x 60
N/A
235.8
MPI-IO is much faster when ~109 grids