京都大学防災研究所年報 第 66 号 B
DPRI Annuals, No. 66 B, 2023

動径基底関数に基づく大気大循環モデルを用いた傾圧不安定波実験
Baroclinic Instability Experiments Using an Atmospheric General Circulation Model
Based on Radial Basis Functions

小笠原宏司 (1)・榎本剛
Koji OGASAWARA(1) and Takeshi ENOMOTO(1) 京都大学大学院理学研究科(1) Graduate School of Science, Kyoto University, Japan
Synopsis
A primitive equation model using radial basis function generated finite differences(RBF-FD) was developed. The new model can use any quasi-uniform nodes because RBFFD is a meshless method. In addition, RBF-FD can achieve a high-order convergence. A
test case of baroclinic instability was used to evaluate the primitive equation model. The
test case examines whether a new primitive equation model numerically converge to thesolutions of the test case. In the test case, three norms of vorticity (2-norm, infinity-norm,
and infinity norm of 2-norm of gradient) are used to check the convergence of vorticity
field. The results show that the RBF-FD model approximately converge in 2-norm andinfinity norm of 2-norm of gradient, but the infinity norm is smaller than the numericalconverged solutions. The smaller infinity norm seems to be due to a Gibbs phenomenon
in the reference solution or to a different formulation of hyper viscosity.

キーワード: 動径基底関数，プリミティブ方程式モデル，メッシュレス法
Keywords: radial basis functions, primitive equation model, meshless method

1.

はじめに

置は節点間距離ができるだけ広く取られているため，
CFL条件で制約されるタイムステップが長くとれる

Radial Basis Functions (RBF) を用いた球面上の偏強みがある．RBF-FDを球面上の偏微分方程式に適用

微分方程式の解法（RBF離散化手法）はメッシュレス

した研究事例は移流モデル（Fornberg and Lehto, 2011;

法であり，現行で用いられているスペクトル手法とGunderman et al., 2020）と浅水波方程式モデル (Flyer

同等の精度である（Flyer and Wright, 2006；Flyer andet al. 2012; Tillenius et al. 2015) があり，球面調和関数

Wright, 2009）．しかしRBF離散化手法は時間積分のを適用したスペクトル変換法や不連続ガラーキン法

際に，節点数を𝑁とした時に計算量がO(𝑁 2 )かかり，よりも少ない節点数で場を再現できることが報告さ

大規模計算に適していない．そこでRBFを用いて生れている．また大規模計算に向けた並列計算の研究成された有限差分法（RBF-generated Finite Differenceも行われている（Bollig et al., 2012; Tillenius et al.,

(RBF-FD) ）が大規模な計算に用いられている．RBF-

2015）．

FDは高次の代数的な収束が可能であり，メッシュレ

2次元球面モデルに関する研究が行われている一ス法である．メッシュレス法は球面において準一様方で，3次元プリミティブ方程式モデルの研究は行わな節点配置を用いることができる．準一様な節点配れていない．そこで本研究では水平方向の離散化に― 151 ―

Fig. 2 Node selection in a unit plane. Green andred nodes are all nodes in the plane. The red nodes
include the central node and its nearest 50 nodes.

Fig. 1 Gaussian radial basis function(GARBF). The horizontal and vertical axes are
the distance between nodes and the function
value of GARBF, respectively.

𝑓(𝑥) ≅ 𝑠(𝑥) = ∑ 𝑤𝑖 𝜙(𝑟)(1)

RBF-FDを，鉛直方向には差分法を適用したプリミティブ方程式 モデルを 開発 し た．開発し たモデル は

重み𝑤𝑖 は選点条件から決定される．選点条件とは，

Polvani et al. (2004)で提案された傾圧不安定のテスト

内挿する関数の値がわかっている点において内挿がケースで性能が評価される．

厳密に成り立つという条件である．選点条件を課す

本論文は第2節でRBF，第3節でプリミティブ方程と式，第4節で時間積分と離散化についてまとめる．第𝑓(𝒙) = 𝑨𝒘

5節ではテストケースについて述べ，第6節で実験結

𝑨𝒊,𝒋 = 𝜙(‖𝒙𝑖 − 𝒙𝑗 ‖)果の結果を示し，第7節スケーラビティーの検証を行う．第8節で結果を考察し，第9節で全体のまとめをとなり，連立方程式をといて𝒘を決定する．決定したする．

重みを(1)に用いると任意の点での内挿値が得られる．
GAでは𝜀を小さく取ることでより高精度な内挿が2.

できるが，𝜀を小さくし過ぎると条件数が高くなり，

Radial Basis Functions

連立方程式を正確に解くことが困難になり，精度がRadial basis functions (RBF) は距離のみに依存する

落ちる．

関数であり ，数種類 提案さ れている ． 本研究で はGaussian RBF（GA）を用いた(Fig.1)．GAは以下の式

2.2

RBF-FD

RBF-FDは微分を求める点と𝑛個の近傍の点からなで計算される．

る𝑛 + 1個の点（ステンシル）（Fig.2）を用いて微分を計算する手法である．微分は重み𝑐𝑖 とデータの積

𝜙(𝑟) = exp(−𝜀 2 𝑟 2 )の和で計算する（Wright and Fornberg, 2006）．
𝑟はユークリッド距離がよく用いられる．𝜀は形状を決める任意のパラメーターである．𝜀を大きくとる

𝑛+1

𝜕𝑓≅ ∑ 𝑐𝑖 𝑓𝑖
𝜕𝑥

ほど短い距離でおよそ0の値になる。

𝑖=1

2.1

RBF内挿

RBF-FDはRBF内挿から導出されているため，RBF内

重みは３つのステップで求める．1)微分を求めた

挿の性質を多く引き継いでいる．そのためRBF内挿かい点（中心点）の近傍のn点（Fig.2）を求める．2)中ら順に説明する．RBF内挿は，格子点が不要なメッシ

心点と近傍の点における節点間距離から行列𝑨を計ュレス法に基づく内挿法である．内挿値は基底関数

算する．3)以下の連立方程式をといて重み𝒄を決定すとするRBFと重みとの積で定義される．

る．
𝑹𝒄 = 𝒃― 152 ―

𝑨
𝑹=[ 𝑇
𝒆
𝒃=

𝑼 = [𝑢𝑝 , 𝑣𝑝 ]
𝜅 = 𝑅d /𝑐p

𝒆
]
0

𝜕𝑨1
𝜕𝑥

𝜎 = 𝑝/𝑝𝑠(𝟐)

𝑨1 = [𝜙1,1 , … , 𝜙1,𝑛+1 , 1]

𝑢𝑝 は東西風，𝑣𝑝 は南北風，𝜓はジオポテンシャル，

𝑒 = [1,1, … ,1]𝑇

𝑝は気圧，𝑝s は地表面気圧，𝑓はコリオリパラメータ𝒆は𝑛 + 1個の１を持つ列ベクトルである．𝑨1 の微ー，𝑅d は乾燥大気の気体定数，𝑐p は定圧比熱．ラグラ

分は連鎖律を用いて行う．RBFは距離で微分し，距離ンジュ微分は以下の式で示される．球標系の水平方

𝑟を𝑥で微分する．行列やベクトルの末尾の1は精度向向の成分には添字hをつける．

上のために入れ，計算される𝒄の最後の要素𝑐𝑛+1 は重みとして用いない．このように計算した𝑑𝜕𝜕
= + 𝑼 ∙ ∇h + 𝜎̇
𝑑𝑡 𝜕𝑡
𝜕𝜎
𝑑𝜎
𝜎̇ =
𝑑𝑡

𝑛個の要素を持つ𝒄を行ベクトルとして積み上げて行列𝑫を作り，微分は以下のように計算する．
𝜕𝒇

≅ 𝑫𝒇
𝜕𝑥
𝒇 = [𝑓1 , 𝑓2 , … , 𝑓𝑁 ]4.

時間積分と離散化

本節では時間積分と離散化手法に関して述べる．
RBF-FDは線形安定性解析において，正の実固有値

この節以降の添字のない水平方向の変数は全て絶を持つ．この正の実固有値を数値積分の安定領域に対デカルト座標系の成分である．またベクトルと行押し込むために超粘性を偏微分方程式に加える必要

列の表記に以下のルールを適用する．

がある（Fornberg and Lehto, 2011）．超粘性は式 (2)ベクトルをBold体の小文字で表しの微分をΔk に置き換えて計算する．方程式には以下のようにして加える．
𝑑𝒇dt𝒗 = (𝑢, 𝑣, 𝑤)＝左辺 + (−1)k+1

節点ご との 変数 を並 べた ベ クトル は下 線付 き の

𝛾
𝑯(𝒇)
𝑁k

Bold体の小文字で表し
𝒗 = (𝑣1 , … , 𝑣𝑁 )𝑯は超粘性，𝛾は超粘性係数である．

3.

プリミティブ方程式

行列をBold体の大文字で示す．

プリミティブ方程式は予報方程式と診断方程式から構成される．予報方程式には運動方程式，連続の式，熱力学の式が含まれる．水平方向に球座標系，鉛

4.1

時間積分

時間積分には4次精度の古典的Runge-Kutta法を適用した．

直方向に地形に準拠した𝜎座標をそれぞれ適用した

4.2 水平・鉛直方向の節点配置

プリミティブ方程式は以下の式で表される．

水平方向には正20面体格子（Fig.3)を用いた．鉛𝑑𝑼= −𝜵h 𝜓 − 𝑅d 𝑇𝜵h ln(𝑝s ) − 𝑓𝐤 × 𝑼
𝑑𝑡

(3)

𝑑𝑇
𝜅𝑇
=
𝜔
𝑑𝑡 𝜎𝑝s(4)1

𝜕ln 𝑝s = − ∫ ∇h ∙ (𝑝s 𝑼) 𝑑𝜎
𝜕𝑡
0

(5)

𝜎

𝜓 = −𝑅d ∫ 𝑇 𝑑(ln 𝜎)(6)0𝜔 = 𝜎𝑼 ∙ ∇h 𝑝s +

𝜕
ln 𝑝s
𝜕𝑡

(7)

直方向にはスタッガード格子を適用する(Fig.4)．層の中心をフルレベル，境界面をハーフレベルと定義する．層数を𝐾とすると，天頂が𝜎1/2 であり，地表面が𝜎𝐾+1/2 である．フルレベルには予報変数とジオポテ
ンシャル，ハーフレベルには鉛直𝜎速度𝜎̇ が配置される．天頂と地表面における鉛直流の境界条件は𝜎̇ = 0
である．
フルレベルの𝜎はハーフレベルから以下の式で求められる．

― 153 ―𝜕 𝜕 𝜕
∇c = ( , , )𝜕𝑥 𝜕𝑦 𝜕𝑧
発散も同様に投影行列を用いて計算できる．
RBFを用いて離散化する際には投影行列を𝒃行列の計算に適用する(Flyer and Wright, 2009)．

𝜕Fig. 3 Icosahedral nodes. The number of grids
is 642. The Mollweide projection is applied.

𝑃

𝜕𝑥
𝜕

𝜙 = (𝑥𝒙𝑡 𝒙𝒌 − 𝑥𝑘 )

𝜕𝑦
𝜕

𝜕𝜙
𝜕𝑟

= 𝑷𝒃

[ 𝜕𝑧 ]
𝒙 = (𝑥, 𝑦, 𝑧)𝑡 , 𝒙𝑘 = (𝑥𝑘 , 𝑦𝑘 , 𝑧𝑘 )𝑡
𝑫𝑥 = (𝑷𝒃)𝑥 𝑹−𝟏(𝑷𝒃)𝑥 = [(𝑃𝑏)1,1 , (𝑃𝑏)1,2 , (𝑃𝑏)1,3 ]
𝑦成分は2行目，𝑧成分は３行目の行ベクトルである．
投影行列を適用した微分行列を用いて勾配と発散は
以下のように計算される．
∇h ≅ [𝑫x , 𝑫y , 𝑫z ]∇h ∙ 𝑈 ≅ 𝑫x 𝒖 + 𝑫𝑦 𝒗 + 𝑫𝑧 𝒘

4.4

Fig. 4 Staggered grids in the vertical. Full andhalf levels are represented by broken and solidlines, respectively. Prognostic variables anddiagnosed vertical velocity 𝜎̇ are placed on
full and half levels, respectively.

4.3

𝑘−

運動方程式の離散化

鉛直移流を除いた水平方向の運動方程式は3次元
デカルト座標系を用いて計算することができる(Swarztrauber et al., 1998) ．式（9）を用いるとデカルト座標系における運動方程式は以下のように計算
できる．

𝜅+1
𝜎 𝜅+1
1 −𝜎 1
𝑘+
𝑘−
1
2
2
𝜎𝑘𝜅 =
𝜅+1𝜎 1−𝜎 1
𝑘+
2

𝜕𝑢
𝜕𝑢
= −𝑑𝑢 + 𝜇𝑥 − 𝜎̇
𝜕𝑡
𝜕𝜎

2

球面の接平面に対する投影

デカルト座標系におけるベクトルは球面の接平面にあるとは限らないため，投影手法を適用する．
球面投影には投影行列𝑷をベクトルの左側からか𝜕𝑣
𝜕𝑣= −𝑑𝑣 + 𝜇𝑦 − 𝜎̇
𝜕𝑡
𝜕𝜎
𝜕𝑤
𝜕𝑤
= −𝑑𝑤 + 𝜇𝑧 − 𝜎̇𝜕𝑡𝜕𝜎
𝜇 = 𝑥 𝑑𝑢 + 𝑦 𝑑𝑣 + 𝑧 𝑑𝑤

ける．

1 − 𝑥2
𝑷 = [ −𝑥𝑦−𝑥𝑧−𝑥𝑦
1 − 𝑦2−𝑦𝑧

𝑑𝑢 = 𝑢 (𝑫x 𝐮) + 𝑣 (𝑫y 𝐮) + w(𝑫z 𝒖) + 𝑫x 𝜓

−𝑥𝑧
−𝑦𝑧 ]
1 − 𝑧2+RT (𝑫x ln (𝒑𝒔 )) + 𝑓(𝑦𝑤 − 𝑧𝑣) + 𝑰𝐮
𝑑𝑣 = 𝑢(𝑫x 𝒗) + 𝑣(𝑫y 𝒗) + 𝑤(𝑫z 𝒗) + 𝑫𝑦 𝜓

投影行列を用いると３次元デカルト座標系における勾配は以下の式で表される．

+RT (𝑫y ln (𝒑𝒔 )) + 𝑓(𝑢𝑧 − 𝑤𝑥) + 𝑰𝐯
𝑑𝑤 = 𝑢(𝑫x 𝑤) + 𝑣(𝑫y 𝑤) + 𝑤(𝑫z 𝒘) + 𝑫𝑧 𝜓

∇h = 𝑷∇c

(9)
(10)

(8)

― 154 ―+RT (𝑫z ln (𝒑𝒔 )) + 𝑓(𝑥𝑣 − 𝑢𝑦) + 𝑰𝑤

(11)

気温の鉛直移流は全ポテンシャルエンタルピーと𝑰に超粘性と場合によって拡散が適用される．
式（11）における−𝑑𝑢 + 𝜇𝑥の項は投影技法の未定乗数法である．ベクトルに対して位置ベクトルとの全エンストロフィーが保存するように以下の式で求める．

内積を取り，𝑥とかけることで𝑥軸の非接平面成分𝜇𝑥を計算でき，𝑑𝑥から引くことで接平面成分を計算で
きる．

4.5

連続の式の離散化𝜎̇𝜕𝑇
1= (𝜎̇ 1 (𝑇 1 − 𝑇𝑘 ) + 𝜎̇ 1 (𝑇𝑘 − 𝑇 1 )) ( )
𝑘+
𝑘+
𝑘−
𝑘−
𝜕𝜎
∆𝜎
2
2
2
2

気温のハーフレベルの値は以下の式で得られる．

連続の式（5）の離散化は式（10）を用いて以下の式で示される．

𝑇𝜕 ln 𝑝𝑠𝑁= −𝛴𝐾=1
𝐷𝑘 Δ𝜎𝑘𝜕𝑡

1 = 𝑎 1 𝑇𝑘 + 𝑏 1 𝑇𝑘+1
𝑘+
𝑘+
𝑘+
2
2
2

ここで𝑎𝑘 ，𝑏𝑘 は以下のように計算される．

𝐺k = 𝛁𝝈 ∙ 𝒖𝒌 + 𝒖𝒌 ∙ 𝛁ln(𝒑𝐬 )= 𝑫x 𝒖𝒌 + 𝑫y 𝒗𝒌 + 𝑫z 𝒘𝒌 + 𝑢𝑘 𝑫x ln 𝒑𝐬𝜎+𝑣𝑘 𝑫y ln 𝒑𝐬 + 𝑤𝑘 𝑫z ln 𝒑𝐬

4.6

𝑎𝑘 = [(1 𝜅𝑘+2

𝜎𝑘 𝜅 −1) − 1] [1 − (
) ]
𝜎𝑘
𝜎𝑘+1

熱力学方程式の離散化

𝜎

熱力学方程式(4)は，式(9)を用いて以下のように

𝑏𝑘+

離散化する．

4.10

𝜕𝑇= −𝑢(𝑫x 𝑻) − 𝑣(𝑫y 𝑻) − 𝑤(𝑫z 𝑻)𝜕𝑡

1 = [1 − (
2

1 𝜅
2

−1
𝜎𝑘+1 𝜅
) ] [() − 1]𝜎𝑘+1
𝜎𝑘𝑘+

ジオポテンシャルジオポテンシャルは以下の式で計算される．

𝜕𝑇 𝜅𝑇
+𝜎̇
+
𝝎
𝜕𝜎 𝜎𝑝𝑠̂𝜓 = 𝜓s + 𝑪𝑻
𝑻̂𝑖 = 𝑇𝑘 , 𝑘 = 1, ⋯ , 𝐾

4.7

𝜔の計算

ωは式（10）を用いて以下のようになる．

0, 𝑙 < 𝑘𝐶𝑝 𝐴𝑘+1/2 , 𝑙 = 𝑘

𝑪𝑘,𝑙 = {

𝐶𝑝 (𝐴

𝑙+

𝐶𝑝 𝑇𝑘
𝑅𝑇𝑘
𝜔 = 𝑅𝑇𝑘 𝑢𝑘 ∇ ln 𝒑𝒔 −+ 𝐵𝑘−1 𝑀𝑘−1 ][𝐴 1 𝑀𝑝𝑘 𝑘∆𝜎𝑘 𝑘+2 𝑘+1
2

1 + 𝐵 1) , 𝑙 > 𝑘
𝑙−
2
2

𝜎
𝐴𝑘+1/2 = (𝑘

𝑀𝑘 = ∑ 𝐺𝑖 ∆𝜎𝑖

𝑘+

1 𝜅
2

𝜎𝑘

𝑖=1

) −1

𝜎
𝐵𝑘−1/2 = 1 − (4.8 鉛直𝜎速度𝜎̇ は水平発散の項をデカルト座標系に変更して計

𝑘−

𝜎𝑘

1 𝜅
2

)

𝜓s は地表面のジオポテンシャルを表す．

算する．

5.
𝜎̇ = 𝜎𝐾 (𝜕𝜕𝑡ln 𝑝s ) − 𝑀𝑘

テストケースPolvani et al. (2004，以下P04) はプリミティブ方程式モデルのバグの発見とモデル間比較を容易に行え

4.9

鉛直移流項の離散化

るように設計されている．P04では12日間で成長する

運動量の鉛直移流はハーフレベルでの値の平均として，以下の式で計算する．

北半球中緯度ジェットにおける傾圧不安定波を計算する．積分時間で短期間で終わるためバグの修正が迅速に行え，時間・空間平均を用いないことでバグ(𝑢𝑘+1 − 𝑢𝑘 )(𝑢𝑘 − 𝑢𝑘−1 )𝜕𝑢𝜎̇
= 𝜎̇ 1+ 𝜎̇ 1𝑘+𝑘−
𝜕𝜎
2 Δσk+1 + Δσk
2 Δσk + Δσk−1

が平均に隠れることも防げる．バグを修正した後に，高精度なモデルで計算された数値的収束解と実験結― 155 ―

果を比較し，渦度場の再現性やモデルの特性を調べ

N10242

ることができるため，本研究ではP04を用いた．

5.1

初期値と摂動

P04において定義されている初期値について簡潔に述べる．風速は南半球全域で0とし，北半球では南北風は0，東西風𝑢𝑝 が以下の式で定義される．
N40962𝑢𝑝 (𝜃, 𝜎) = 𝑢0 sin3 (𝜋𝜇2 )𝐹(𝑧)1𝑧 − 22𝜋𝑧
𝐹(𝑧) = [1 − tan3 ()] sin ( )2530
𝜇 = sin(𝜃) , 𝑧 = −7.34 log(𝜎)𝜃は緯度である．気温場は以下の式で定義されている．

𝑇(𝜃, 𝑧) = ∫N163842𝜕𝑇(𝜃, 𝑧)𝑑𝜃 + 𝑇0 (𝑧)𝜕𝜃𝜕𝑇(𝜃, 𝑧)7.41
𝜕𝑢(𝑎𝑓 + 2𝑢 tan 𝜃)
=−
𝜕𝜃
𝑅𝑑
𝜕𝑧
𝑇0 (𝑧)はアメリカ標準大気，𝑎は地球平均半径である．
傾圧不安定を引き起こすための摂動は全層の気温場に与えられる．摂動は以下の式で計算される．

しかし，異なる離散化手法でも高解像度化によりおFig. 5 Surface vorticity field. The contour
interval is 1 × 10−5 s −1 . The horizontal andvertical axes are the longitude and latitude,
respectively. The contour almost converges at
N40962.

よそ等しい渦度場に収束することが報告されている．

いた．時間刻み幅は480,240,120[s]である．超粘性の渦度は最下層の値が用いられ，以下のノルムが計算

階数は6，任意係数𝛾は2 × 10−1である．ステンシルサ

される．

イズは51である．

𝑇′ (𝜆, 𝜃) = sech2 (3𝜆) sech2 {6(𝜃 − 𝜋/4)}5.2

数値的収束解本テストケースには厳密な解は示されていない．

ℓ2 =

1

√4𝜋

2

{∫ ∫ 𝜁 cos 𝜃 𝑑𝜃𝑑𝜑}

1
2

6.

実験結果

ℓ∞ = max|𝜁|
12日目における最下層の渦度場を節点数Nごとに|𝛻𝜁|∞:渦度の勾配ノルムの２ノルムの最大値

Fig. 5に示した．N10242においておよそ形を捉えてこれらの渦度ノルムの時間経過が高解像度化によりおよそ同じ時間変化に収束することを目指す．まいるがノイズがある．N40962において渦度場の等値
線の概形は収束している．

たノルムだけでは渦度場全体を把握できないので，Fig.6に渦度ノルムの時間発展を示す．ℓ2 は節点数決められた幅で描かれる等値線の概形が１つの形に10242の時点で高解像度のものと変わらない．ℓ∞ と収束することも収束判定に用いる．

|𝛻𝜁|∞ は40962において収束しており，12日目の値はわずかに異なる．次に12日目のノルム値をP04と比較5.3

する（Table.1）．ℓ∞ と|𝛻𝜁|∞ はおよそ等しい値に収モデルの実験設定

節点数は10242,40962,163842点の3通りを用いた．

束し，ℓ∞ がP04より小さい結果になった．

節点数に対して形状パラメーターは4.6,10,22を用― 156 ―

Table1 Vorticity norms at 12 days. The number ofnodes for ICOS and P04 are 163842 and 524288,
respectively.
ℓ2

7.

|𝛻𝜁|∞

ℓ∞
−6

P04

7.8 × 10ICOS

7.8 × 10−6

−5

3 × 10−10

7.3 × 10−5

3 × 10−10

7.4 × 10

スケーラビリティーの検証

Fig.7にプロセス数を変えたときの1日の計算時間を示す．用いた節点数は163842である．計算時間にはモデルの初期化からデータの読み書きまでの時間が含まれている．プロセス数を2倍にすると計算時間は2分の１になる．しかし，プロセス数およそ11.9倍の時にはおよそ10分の1の計算時間になった．

8.

考察

渦度場の再現性を12日目におけるノルムのP04とで比較すると，ℓ2 と|𝛻𝜁|∞ は同等の値で，ℓ∞ は低い結果であった．ℓ2 が計算できているため，大まかに渦度場は再現できている．しかし，ℓ∞ は低く，|𝛻𝜁|∞はおよそ同等であるため局所的な差がP04との間にある．
局所的な違いを離散化の差から検討する．Flyer etFig. 6 Vorticity norms vs time [d]. The verticaland horizontal axes are norm and day,respectively. N is the number of ICOS nodes.

al. (2012) の浅水波方程式モデルにおける実験からℓ2 , ℓ∞ and |𝛻𝜁|∞ are 2-norm, infinity norm and

道で風速が緯度方向に対して無限回微分不可になり，

infinitynormof

2-normof

gradients ,respectively.

スペクトル変換法とRBF-FDではGibbs現象による影響が異なると報告されている．P04の初期条件では赤Gibbs現象が起きうる．これによりGibbs現象の影響の差が離散化間で生じ，収束する場が異なっている可能性がある．
またRBF-FDでは安定化のために超粘性を適用しているため，全ての波が影響を受けている．これにより場の振幅がP04とずれて，最大値の差につながっている可能性がある．
次にスケーラビリティーを検証する．プロセス数の増加に従って，計算時間の減少率が鈍っていく原
因は，RBF-FD手法において並列化しても実行時間が頭打ちになる箇所があるためである．減少率は鈍る
が，並列化による高速化は十分に出来ているためモデルにはスケーラビリティーがある．

9.

Fig.7 Calculation time for 1 day for 21, 42, 83,
141, 166 and 249 processes.

まとめ

本研究では水平，鉛直方向の離散化にRBF-FD，と有限差分法を適用したプリミティブ方程式モデルを― 157 ―

開発した．PO4で提案された傾圧不安定波の実験を行った．実験の結果からP04のテストケースは先行研究とおよそ整合しており，局所的な差は離散化の特性の違いが原因と考えられる．また開発したモデルは並列化により高速化が可能であることが示された．

謝辞

本研究は京都大学防災研究所共同研究（令和 5 年
度萌芽的共同研究 2023NE-02）の支援を受けた．

参考文献
Arakawa, A., and Suarez M., 1983: Vertical
Differencing of the Primitive Equations in Sigma
Coordinates., Mon. Weather Rev., 111, pp. 34–45.
Bollig, E. F., Flyer N., and Erlebacher G., 2012:
Solution to PDEs using radial basis function finitedifferences (RBF-FD) on multiple GPUs., J.
Compt. Phys., 231, pp. 7133–7151.
Flyer, N., and Wright G. B., 2007: Transport
schemes on a sphere using radial basis functions.,
J. Comput. Phys., 226, pp. 1059–1084.
Flyer, N., and Wright G. B., 2009: A radial basis
function method for the shallow water equations
on a sphere., Proc. R. Soc. A: Math. Phys. Eng.,
465, pp. 1949–1976,.
Flyer, N., Lehto E., Blaise S., Wright G. B., and StCyr A., 2012: A guide to RBF-generated finite

differences for nonlinear transport: Shallow water
simulations on a sphere., J. Comput. Phys., 231, pp.
4078–4095.
Fornberg, B., and Lehto E., 2011: Stabilization of
RBF-generated finite difference methods forconvective PDEs., J. Comput. Phys, 230, pp. 22702285.
Gunderman, D., Flyer, N., and Fornberg, B., 2020:
Transport schemes in spherical geometries using
spline-based RBF-FD with polynomials., J.
Comput. Phys., pp. 408, 109256.
Polvani L. M., Scott R. K., and Thomas S. J., 2004:
Numerically Converged Solutions of the Global
Primitive Equations for Testing the DynamicalCore of Atmospheric GCMs., Mon. Weather Rev.,
132, pp. 2539–2552.
Swarztrauber, P. N., Williamson D. L., and Drake J.
B., 1998: The Cartesian method for solving partial
differential equations in spherical geometry. Dyn.
Atmos. Oceans., 27, pp. 679–706.
Tillenius, M., Larsson E., Lehto E., and Flyer, N.,
2015: A scalable RBF–FD method for atmosphericflow. J. Comput. Phys., 298, pp. 406–422.
Wright, G. B., and Fornberg, B., 2006: Scatterednode compact finite difference-type formulas
generated from radial basis functions., J. Comput.
Phys., 212, pp. 99–123.

― 158 ―

（論文受理日：2023年8月31日）