（続）スーパーコンピュータ「京」の利用：5．1km格子間隔を切る大気大循環シミュレーションへの道のり

全文

(1)特集. 続. スーパーコンピュータ「京」の利用. 5. 基応専般. 1km 格子間隔を切る大気大循環シミュレーションへの道のり. 八代尚. 1,2. 富田浩文. 1,2. 宮本佳明. 1. 1 理化学研究所計算科学研究機構 2 海洋研究開発機構. 全球雲解像モデルによる大気シミュレーション. 持っている．凝結した水は潜熱を開放し，浮力を持. は，気象・気候学分野では国内外問わず多くの反. って上昇する．より高い高度ではさらに氷や雪に. 響とともに，解像度不足であるとの批判も受けて. 変化する．こうして水の相変化によって引き起こ. きた．積乱雲の水平スケールは 1km のオーダであ. された対流運動は時に高度 10km 以上にも達する. り，2004 年に達成した 3.5 ㎞という解像度で雲を. 雄大な積乱雲を作り出し，大気の大循環を駆動す. 陽に表現することでは雲解像とは言えないとの批. る源にもなっている．このミクロスケールから惑. 判である．この批判に対する 1 つの解答として，我々. 星スケールまでにおよぶ一連の過程をモデルシミ. は「京」を用い 1㎞を切る解像度での全球シミュレ. ュレーションで現実的に表現するには，一体どの. ーションを行い，積乱雲の統計的性質が 2km あた. 程度の格子（メッシュ）間隔が必要なのであろうか．. りを境に変化し，積乱雲ひとつひとつをより正確. 今から 20 年前の計算機の能力では，そもそも粗い. に表現する足掛かりを得た．本稿では超高解像度. メッシュ間隔（数百 km）でなければ全球シミュレ. 大気実験の超並列実行に至るさまざまな問題，克. ーションを行うことができなかった．積乱雲の水. 服点を紹介する．. 平スケールは 1km のオーダであるため，数百 km のメッシュの中に存在するはずの積乱雲群を，エ. 全球気象モデルは大気現象を十分に表現しているか. ネルギー・質量収支の観点からモデル化し，鉛直. 気象の移り変わりを第一原理に基づいて表現する. れ，今も多くの改良を経て広く用いられている．し. としたら，大気は分子にまで還元される．全大気の. かし積乱雲群を統計的に扱っているため，大きな不. 分子運動をそれぞれ解くような計算ができる未来は. 確定性が残る．日本では海洋研究開発機構におい. まだまだ先であろう．そこで大気を秒速数 m から. て，熱帯の積乱雲を表現するのに十分な空間解像度. 数百 m で移動する流体として扱うのだが，気象現. でのシミュレーションを可能にする「全球雲解像モ. 象はこの大気の大きな運動だけにはとどまらない．. デル」の開発が 2000 年より開始された．全球非静. たとえば水蒸気は数十～数百ナノメートルの粒子を. 力学モデル NICAM（Nonhydrostatic Icosahedral. 種にして液体の雲となり，それらがぶつかり合って. Atmospheric Model）と名付けられたこの大気モデ. やがて雨となる．また，太陽から降り注ぐ光は，ま. ルは，地球シミュレータの性能を十分に活用し，世. さに光の速さで到達して大気や地表を暖める．我々. 界で初めて水平メッシュ 3.5km での全球シミュレ. はこれら時空間スケールの異なる多くの過程をモデ. ーションを達成した（Tomita et al., 2005 など）．. ル化・粗視化して，コンピュータ上で気象シミュレ. それまで地球上の一部分の領域のみを区切って実行. ーションを行っている．その中でも雲の生成消滅は，. していた高解像度の実験が，全球をすっぽり覆う形. シミュレーション結果に対して大きなインパクトを. で実現したのである．成果は国内外の気象・気候コ. 方向に再分配することで積雲対流の効果を表現した．この手法は積雲パラメタリゼーションと呼ば. 1）. 情報処理 Vol.55 No.8 Aug. 2014. 811.

(2) 特集. 続. スーパーコンピュータ「京」の利用ミュニティに大きな驚きをもって迎えられた．それ. し 10 ％を超える実行効率を得た．しかし，アプリ. と同時に，3.5km ではまだ積乱雲 1 つ 1 つを表現し. ケーション全体に戻ると効率はピーク性能比で 1%. ているとはいえないという批判も展開された．水平. 程度しか向上しなかった．実行速度でみると 2 割. メッシュ数百 m での理想実験においても，積乱雲. 程度しか高速化しなかったということである．これ. の構造の表現はいまだ不十分であるとの研究結果も. はアムダールの法則に従うもので，高速化した部分. 報告されている．このような批判をふまえて，我々. の実行時間に占める割合が相対的に小さくなったた. は全球モデルの解像度のさらなる向上によって積乱. めであった．NICAM のソースコードは数十万行あ. 雲の表現がどう変化していくのかを知るべく，新た. るが，その中で総計算時間の大部分を占めるような. なスーパーコンピュータである「京」の性能を活か. 計算カーネルはなく，演算量の少ないカーネルがい. した超高解像度・大規模実験に挑戦することにした．. くつも存在している．そのため一般的な性能プロファイラによるコストランキングでは，性能向上のた. 「京」における NICAM の性能. めの糸口を捉えることが難しかった．そこで我々は，. NICAM は 2006 年より始まった次世代スーパー. み，浮動小数点演算と経過時間の計測を徹底的に行. コンピュータ開発プロジェクトにおいて，初期の段. った．その結果，これまでベクトル計算機での最適. 階からターゲットアプリケーションの 1 つとして選. 化では気にしていなかった中間配列の多用や配列の. 定され，開発チームと計算科学研究機構・運用技術. ゼロ埋め等が，メモリバンド幅の小さいスカラー機. 部門との共同で「京」での性能最適化が進められた．. では無視できないほどの時間ロスと実行効率の低下. これまで NICAM はベクトル機である地球シミュ. の原因になっていることが分かった．これら演算密. レータをメインマシンとして開発されてきた．その. 度の低い部分を根気よく排除することで，実行時間. ため，「京」への移植ではメモリバンド幅の小さい. の短縮と高効率化の両方を得ることができた．加え. スカラー機でいかに性能を向上させるかという点が. て，さらなる高効率化のために，コードのリファク. 鍵となった．大気モデルのメイン計算部分は流体力. タリングを通して条件分岐文の挿入位置の改善やサ. 学の方程式系を時間発展させて解く部分（力学過. ブルーチン階層の見直しを行った．その結果，京コ. 程）と，力学過程以外の諸過程をモデル化した部分. ンピュータのピーク性能に対し，アプリケーション. （物理過程）に大別される．物理過程には，雲の相. 全体で最大 10％（倍精度）の実行効率を達成した．. 変化や雨の落下，地表面を通した熱・水蒸気の流出. 小さな時間ロスの積み重ねにより，全体性能が容易. 入，太陽・赤外放射による大気の加熱・冷却，そし. に低下するということは，今回の最適化の経験で得. て格子幅よりも小さいスケールの流体（乱流）運動. られた大きな教訓である．一般的に気象・気候モデ. の効果等が含まれる．力学過程が主にステンシル計. ルは学際的な連携によって開発・統合されたアプリ. 算（隣接する格子の値を使った計算）を行うのに対. ケーションであり，今後も新しいコンポーネントや. し，物理過程の各コンポーネントは水平方向には独. スキームが追加されていく．継続的に高い演算性能. 立して計算され，ループの内部にたくさんの演算と. を維持するには，どのコンポーネントにおいても高. 条件分岐が記述される傾向にある．NICAM では力. 効率性に留意し，それらを組み合わせていく努力が. 学過程と物理過程の所要時間がおおむね 1:1 であり，. 必要である．. それぞれに異なる最適化方針がとられた．詳細につ. 「京」での単体性能，スケーリング性能の検証を. 2）. いては南ら（2012）. 812. 百カ所を超えるカウンタを手動でコード内に埋め込. を参照されたい．最適化に. 行うにあたり，我々はベンチマークテストをできる. よって，発散項や勾配項を計算する主要なステンシ. だけ現実の気象場を想定した設定で行うことにし. ル計算カーネルは京コンピュータのピーク性能に対. た．これは大気モデルが内部にたくさんの条件分岐. 情報処理 Vol.55 No.8 Aug. 2014.

(3) 間隔 0.87km，鉛直 94 層とした．このときの 3 次元. ／低い，高度に雲がある／ないかによって，演算量. 格子数は 630 億，時間ステップの刻みは 2 秒である．. が大きく変わってくる．図 -1 にウィークスケーリ. 「京」全体のおよそ 1/4 である 20,480 ノード（163,840. ングの結果を示す．この実験では，ノードあたり. コア）を用い，1 回のジョブ単位で 2 時間の時間発. の問題サイズが一定になるように，水平メッシュ. 展を計算した．このときの演算性能は 230TFLOPS. 間隔を狭めて総格子数を増やしつつノード数を増. （ピーク性能比約 8%），メモリ使用量は 250TB で. 加させている．結果は良好であり，81,920 ノード. あった．加えてこれまで行ってきた解像度を含む，. （655,360 コア）まで十分な実行効率を保って実行. メッシュ間隔 14km から 1.7km での比較実験も並. されることが確かめられた．ノード数が増加する. 行して行った．図 -2 に 0.87km メッシュ実験によ. につれて，実行時間が増加し効率が若干下がる傾向. って得られた全球の雲分布を示す．図をみると，惑. にあるが，これはノード間通信時間やノードごとの. 星スケールの現象から非常に細かい雲の分布までの. 演算量のばらつきがより顕著になるためであると考. 再現に成功したことが分かる．特に雲の水平・鉛直. えられる．演算量は前述のとおり，各ノードが担当. 構造は，これまでに得られたどの全球大気シミュレ. する格子の中にどれだけ雲が存在しているかの違い. ーションの結果よりも精細なものと判断できる．さ. で大きなインバランスが発生し得る．将来のスパコ. らに各領域に注目すると，熱帯と中高緯度帯では異. ンでさらにプロセス分割数が増加したり，より精緻. なるタイプの雲の集団があることが見てとれる．ま. で演算量の変動が激しいスキームを採用したりした. た日本付近には大型で非常に強い台風（平成 24 年. 場合，インバランスはより顕著になると予想される．. 台風第 15 号）が接近しており，その南西にはもう. 5 1km 格子間隔を切る大気大循環シミュレーションへの道のり. を持っているためである．たとえば，零度より高い. ひとつの台風（平成 24 年台風第 14 号）も存在し. 格子間隔が 1km を切る全球大気シミュレーション. ている．図 -3 に台風の付近を拡大し，より粗い解像度（3.5km メッシュ）での結果を並べたものを示す．解像度の向上により，荒い解像度では捉えきれ. る超高解像度実験が進められた．計算資源や実行に. に再現されている．この台風だけでなく，異なる雲. 要する時間を考慮し，最大実験規模は水平メッシュ. の大規模システムの様相を全球同時に得ることがで. GL09RL02Z94 14km. 13. GL13RL06Z94 870m. GL11RL04Z94 3.5km. 12. 11. 10. 9. 8. 7. 6. 5. 4. 3. 2. 1. 40. 160. 640. 2560. 10240. Number of cores. 40960. 163840. 655360 実行時間. Performance eﬃcinecy［%］. なかった台風の内部構造が，よりリアリスティック. Elapsed time［sec/step］. こうして得られた演算性能をもとに，「京」におけ. 0 実行効率. 図 -1 「京」における NICAM のウィークスケーリング性能．横軸は使用コア数（使用ノード数の 8 倍）を表す．棒グラフが 1 ステップあたりの経過時間，折れ線が実行効率を表す．. 情報処理 Vol.55 No.8 Aug. 2014. 813.

(4) 特集. 続. スーパーコンピュータ「京」の利用. 図 -2 全球非静力学モデル NICAM による水平メッシュ 0.87km シミュレーション： 2012 年８月 25 日における雲分布（理化学研究所吉田龍二氏提供）. きたのもこの実験の大きな特徴である．我々はこの利点を活かし，それぞれの雲集団からさらに，その中に含まれる積乱雲 1 つ 1 つを抽出することを試みた．大きく発達した積乱雲が存在する格子を客観. 0.87km メッシュ実験における入力データサイズ. 的に抽出し，それらの点を中心として周囲の場がど. は，境界条件が 1TB, リスタートファイルが 7TB. う表現されているかについて，全球での統計的性質. であった．NICAM ではこれらの入力ファイルを計. 3）. を解析した．図 -4 は Miyamoto et al.（2013）. 814. 膨大なデータ量のシミュレーション結果とその解析. に. 算する際のプロセス数で分割して保持する形をとっ. よって示された，積乱雲とその周辺域の上昇流を全. ている．たとえば「京」では 1 ノードあたり 1 プロ. 球平均した結果である．図より，3.5km よりも粗い. セスで 20,480 ノードを利用する際に 20,480 ファイ. 格子間隔では上昇流の強い領域が積乱雲中心の格子. ル準備し，各計算ノードのローカルファイルシステ. 点だけで強く現れるのに対して，1.7km よりも細か. ムへステージインを行う．出力においても同様に分. い格子間隔では複数の格子点で表現されていること. 散ファイルでのアウトプットを行う．今回の実験で. が示されている．すなわち，NICAM は水平格子間. はジョブあたり 2 時間のシミュレーションを行う中. 隔 2km を切るあたりから積乱雲をより正確に表現. で，30 分ごとに 3 次元の風速，気温，気圧，各相. し始めていることを発見したのである．また，積乱. での水物質の混合比や 2 次元の降水量，海面更正気. 雲同士の最短間隔についても解析を行い，1.7km と. 圧などの変数出力を行った．出力タイミング１回に. 0.87km では粗い解像度の時よりも，積乱雲が互い. つき書き出されるデータ総量は 3TB であり，1 ジ. に格子点間隔を広くとる傾向にあることが分かった．. ョブでは 13TB に達した．入出力ファイルをどの. このことは，今までのシミュレーションでは積乱雲. ような形で扱うかについては，いくつかの選択肢が. 同士の距離が粗いメッシュ間隔によって規定されて. ある．国内外の気象・気候モデルでは，シミュレー. いたのに対し，解像度の向上によってより物理法則. ション実行時に 1 ファイルへ集約して出力する方式. に適った距離を再現することが可能になったことを. をとることが多い．出力後のデータ解析や可視化で. 示唆する結果である．. 扱いやすいからである．しかしこのやり方は，高解. 情報処理 Vol.55 No.8 Aug. 2014.

(5) 5 1km 格子間隔を切る大気大循環シミュレーションへの道のり. 図 -3 NICAM によって再現された平成 24 年台風第 15 号（アジア名：ボラヴェン）．左が 3.5km メッシュ，右が 0.87km メッシュのシミュレーション結果を表す．（理化学研究所吉田龍二氏提供）. 像度化に伴いデータ量が膨大になると，計算実行時. 積乱雲間の距離計算をはじめとする一連の解析が行. にファイル I/O スループットが全体の律速となっ. われた．1 つのタイムステップの出力からの計算に. てしまう．これを回避するためには MPI-IO のよう. 32PE，メモリ 1TB 搭載の Linux マシンで約 300. な並列 I/O を採用する必要が出てくる．もう 1 つ，. 時間かかっていたものを，「京」を用いて約 5 分で. 膨大な入出力データが抱える問題として，データス. 行うことが可能になり，積乱雲抽出の手法や，しき. トレージの圧迫がある．ディスクスペースの削減の. い値を変えた解析を容易に行うことができた．現在. ために今後は圧縮フォーマットが主流になると考え. はさらに解析ツールを拡充し，各格子点での地方時. られるが，たとえば気象・気候モデルで広く使われ. に沿った時系列データの並べ替えと全球統計を組み. る NetCDF は MPI-IO での出力と圧縮を同時に行. 合わせた解析を行っている．この解析では，1 回の. うことができない．そこで，今後は圧縮と並列 I/O. ジョブ計算のためにシミュレーション 1 日分にあた. の両方が同時に可能なライブラリの普及が望まれる．. る，延べ 160TB の出力データを一度にステージイ. NICAM はデータが正 20 面体格子に沿う形で格納. ンしている．. されていて，気象解析で一般的に用いられる解析ソフトウェアを活用するには緯度経度格子に変換する必要がある．性能と利便性を考慮した結果，シミュ. まとめ. レーション時にはファイル I/O スループットを重. 全球非静力学モデル NICAM は「京」での高い. 視して，ノードごとに独立した入出力を分散ファイ. 演算性能，スケーリング性能を達成し，水平メッシ. ルに対して行う方式を採用し，一方ポストプロセ. ュ間隔 1km を切る全球大気計算により気象モデル. スで格子系の変換と 1 ファイル化を行うことにし. シミュレーション研究を新たなステージに押し上げ. ている．. ようとしている．今後も高効率・高並列性を意識し. 0.87km 実験で得られた出力データは，日々巨大. たモデルコンポーネントおよび解析ツールの継続的. なデータサイズに慣れている我々からしても簡単に. な開発によって，気象・気候研究のみならず，防災・. 取り扱うのが難しい量であり，通常の感覚でポスト. 減災を支える強力なアプリケーションソフトウェア. プロセスや解析作業を行おうとしたところ，並の. として発展していくことが期待される．我々は将来. Linux クラスタマシンでは現実的な時間で完了しな. の計算機環境においてもその性能を最大限に活かし，. いという事態に陥ってしまった．これを解決するた. さらなる超高解像度シミュレーションを追求するつ. めに，「京」を用いた大規模並列でのポストプロセ. もりである．また今回の 0.87km 実験は「京」を用. スと解析を行うためのツールプログラムを新たに用. いて数日分を計算するのがやっとであった．台風や. 意した．開発した解析プログラムを用いて，前述の. 熱帯域での変動現象の研究のためには長期間のシミ. 情報処理 Vol.55 No.8 Aug. 2014. 815.

(6) 特集. 続. スーパーコンピュータ「京」の利用. 16. (a) 14. 16. (b) 7.0. 16. (c) 3.5. 16. (d) 1.7. 16. 12. 12. 12. 12. 12. 10. 10. 10. 10. 10. 8. 8. 8. 8. z (km). 14. z (km). 14. z (km). 14. z (km). 14. z (km). 14. (e) 0.87. 8. w(ms-1) 0.5 0.4 0.3. 6. 6. 6. 6. 6. 0.2. 4. 4. 4. 4. 4. 0.1. 2. 2. 2. 2. 2. 0 -6 -4 -2 0 2 4 6 x (grid). 0 -6 -4 -2 0 2 4 6 x (grid). 0 -6 -4 -2 0 2 4 6 x (grid). 0 -6 -4 -2 0 2 4 6 x (grid). 0 -6 -4 -2 0 2 4 6 x (grid). 0. 図 -4 全球に存在する積乱雲をすべて平均した鉛直速度の高度分布，左からメッシュ間隔 14km，7km，3.5km，1.7km，0.87km の結果を表す．横軸は積乱雲中心からの格子距離を示す， 3）（Miyamoto et al., 2013 より引用）. ュレーションを必要とする．今後は１ステップの計算時間をできるだけ短縮するための最適化，改良も進めていく．. 参考文献 1） Tomita, H., Miura, H., Iga, S., Nasuno, T. and Satoh, M. : A Global Cloud-Resolving Simulation : Preliminary Results from an Aqua Planet Experiment, Geophys. Res. Lett. Vol.32, L08805, doi : 10.1029/2005GL022459 (2005). 2）南一生，井上俊介，熊畑清，黒田明義，寺井優晃，長谷川幸弘 : スーパーコンピュータ「京」: 6. アプリケーションの性能と最適化─アプリケーションによるペタ FLOPS 性能の実証─，情報処理 ,Vol.53, No.8, pp.787-794 (2012). 3） Miyamoto, Y., Kajikawa, Y., Yoshida, R., Yamaura, T., Yashiro, H. and Tomita, H. : Deep Moist Atmospheric Convection in a Sub-Kilometer Global Simulation, Geophys. Res. Lett. Vol.40, doi : 10.1002/grl.50944 (2013). （2014 年 3 月 1 日受付）. 816. 情報処理 Vol.55 No.8 Aug. 2014. 八代尚 [email protected] 2007 年東北大学大学院理学研究科地球物理学専攻修了．博士（理学）．現在，（独）理化学研究所計算科学研究機構複合系気候科学研究チーム研究員．高解像度気候・気象モデルを用いた研究に従事．富田浩文 [email protected] 1999 年東京大学大学院工学系研究科航空宇宙工学科修了．博士（工学）．現在，（独）理化学研究所計算科学研究機構複合系気候科学研究チームリーダー．専門は気象・気候学．宮本佳明 [email protected] 2011 年京都大学大学院理学研究科地球惑星科学専攻修了．博士（理学）．現在 ,（独）理化学研究所計算科学研究機構複合系気候科学研究チーム基礎科学特別研究員．専門は気象学．.

(7)