プログラムの高速化・並列化サービスの事例:チューニングレポート<要約>：天然磁化鉱物系の微視的磁気モデリングコードMicroMagの並列化

＊ここに掲載するのは、エジンバラ大学EPCCのChris M. Maynard博士によるHECToRレポート「MicroMagnetic modelling of naturally occurring magnetic mineral systems, Dr. Chris M. Maynard, EPCC, The University of Edinburgh, June 8, 2011」を要約したものです。

[2017年4月掲載]

概要

地球科学では、鉱物の磁気特性とそれらに含まれる古地磁気記録の安定性の正確な決定に依存する問題が多く存在します。例えば、地磁気の振舞いの研究は地表および地表面上の観測磁場の変動に大きく依存します。その直接の記録は高々200年前までしかなく、より詳細な解析は岩石中に保存された古い記録が頼りです。古代の場の方角的な記録の地質学的な解釈は、50年以上前に、地球上の海洋及び大陸プレートが移動し続けているという事実を導き出しました。磁気鉱物学のより完全な理解は、大陸ブロックの詳細な運動と回転だけでなく、岩石が出来た温度と熱時間履歴の詳細な分析を可能にします。さらに、磁気鉱物学は古細菌の変種を決定する手法としても頻繁に使用されています。

こうした研究の全ては、鉱物の微視的構造や化学、グレイン形状、グレイン間の磁気的相互作用と共に磁気的特性がどのように変化するかの理解に掛かっています。自然発生磁化鉱物の複雑性と多様性により、その数値的な微小磁気問題解析を、より理想的な人工的記録媒体よりもさらに複雑にしています。環境磁気研究の高度化は、天然磁化鉱物系で生じる高精度の磁気記録プロセスを基礎とします。磁化鉱物に局所的な地磁気方向と強さを記録して保存することを可能とする基本的なプロセスは、未だに完全に理解されているとは言えません。

平衡化した磁区構造は、Landau-Lifshitz-Gilbertの運動方程式を積分することで得られます。磁化領域形状はそれがモデル化された精度に敏感なため、こうした不均一な形状を扱うには有限要素法が適しています。コードMicroMagはプロジェクト開始時点では、Fortran90で記述されたシリアルコードでした。

アルゴリズム

微視的磁気数値モデルは、有限要素メッシュの頂点N個で配置された構造内の磁化を3Dカルテシアン単位ベクトルにより表現します。離散化誤差を小さくするために四面体要素を用います。計算セルのサイズは、交換長という微視的磁気量で決定され、これは最大のセルサイズに対する制約にもなります。

平衡磁区構造は系の全自由エネルギー最小化により求めます。これはLandau-Lifshitz-Gilbert $LLG$ 方程式を積分することで得られます。運動学的な振る舞いは、LLG方程式を各時間ステップで解くことで得られます。全有効場の寄与には、計算集約的な非局所的な場（静磁場）の計算と、局所交換場、異方性場、外場の計算も含まれます。これらの中で静磁場計算には工夫を要します。交換場と異方性場の計算には、有限体積法を適用します。これは本質的には、剛性マトリックスに磁化ベクトルを掛ける計算です。

静磁場はスカラーポテンシャルを用いて解きます。磁化の発散に対してはポアソン方程式、境界の頂点に対してはラプラス方程式を解きます。長距離相互作用の近似として境界要素法を用います。これらは、線形有限要素解のガラーキン法を用いて離散化します。その解法には標準的な共役勾配法を用います。ポアソン方程式は極めてスパースなため、効率的なスパース行列解法の選択が重要です。また、四面体有限要素基底関数は、残る3成分場を単純なスパース行列-ベクトル乗算で計算する効率的な計算手段を提供します。最初の有限要素基底関数の効率的な構築で生じるデータレイアウトから、専用のスパース行列管理が必要とされます。

一方、長距離相互作用に対する境界要素近似から生じる行列のサイズはN²です。ここでNは境界上の頂点数です。並列コードでは、この密行列は、個々の処理単位へ効率的に分散させることが可能です。LLG方程式である常微分方程式（ODE）の積分には、陰的ソルバーが必要であり、よって3N連立スパース線形方程式に対する近似逆行列も必要となります。時間ステップサイズが系の剛性で制約されるため、陽解法は有効ではありません。この陰解法は、可変係数を持つマルチステップ後退オイラー法をベースにした可変次数ソルバーです。

作業計画

当初の計画では、コードの完全な並列化作業のために、nAGとEPCCのスタッフ両方による16人月が要請されました。最終的には、以下の作業項目に対して6人月で決着しました。

·作業１：並列ライブラリ実装（3人月）

100000要素の相互作用無の単独の球、および100000要素の相互作用無の単独のフランボイドのテストケースによる、シリアルコードのプロファイリング
並列CVODEソルバー[1]呼出しの実装とプロファイリング
成果：並列CVODE実装MicroMagコード

·作業２：完全な並列化(2人月）

Cubitの四面体メッシュ分割をMetisライブラリー実装により自動化
有限要素番号付けシステムのためのローカルからグローバルへのマッピングを実装
完全並列化MicroMagのテスト
成果：MicroMag並列バージョンおよび、HECToRの256コアを用いたデモンストレーション

既存コード概要

MicroMagは長い間、他の科学研究コードと同様に、複数の人物により、明確な開発計画なしに多くの機能が開発されているため、旧式の特性を持ち保守が面倒になっていました。

最初に、既存のODEソルバーを現代的なライブラリCVODEへ置き換えました。これだけでもシリアル性能が高速化しました。更にCVODEは並列実行も可能です。この後にnAG CSEチームがシリアルコードをリファクタリングしました。これにより、不要な機能と歴史的な由来を持つ部分が省かれて、簡潔な構造になりました。

並列ライブラリ実装

CVODE実装の結果、XT4上で20K頂点、100K要素の球ケースの外場計算では、既存コードが3098.86秒だったのに対し、CVODE版は1890.35秒に減少し1.64倍高速化しました。

この計算は2つの部分、陰解法ODEソルバーおよびその時間ステップ更新に必要となる右辺項（RHS）の計算から成ります。RHSはエネルギーを求めるもので、最も計算集約的なコード部分です。ここでは、RHS計算はシリアル実行してその結果のエネルギーを並列要素へ分散します。CVODEソルバーが並列に実行されて、結果が集約されて次のRHSシリアル計算に用いられます。

PETScによる並列実装

nAG CSEチームとの協議の結果、PETScライブラリ[3]を利用することとしました。PETScはMicroMagで必要となる、ベクトルおよびスパース行列のデータ型をサポートします。重要なのはSUNDIALS CVODE library [1]をPETScから呼べることです。RHS計算に必要なクリロフ部分空間ソルバーのサポート可能です。さらに、有限要素分割ツールMetisとその並列版ParMetisもPETScから呼び出すことが可能です。

·行列/ベクトルの並列構築

作業のほとんどが、PETScを効率的に用いるための並列データ構造の構築に費やされました。MicroMagは、Metis領域分割メッシュを読み込み、鉱物グレインに関する情報とともにこれを用いて一連のスパース行列を作成し、これをベクトルである磁場の更新に使用します。MicroMagは結果を、SLAPカラム形式で行列フォーマットに保存します。
MicroMagの扱う行列が対称行列であることを利用すれば、PETScの列方向順序への対応は単純な変換のみで済みます。

·時間発展に対する並列CVODE

データ構造の並列構築の後にやるべきことは、CVODEソルバー並列実行とPETScのクリロフ部分空間ソルバーによるRHS計算の並列化です。しかしながら、ここで発見されたインターフェイス不具合のPETSc開発者からの解決を待つ状況の中、CVODEの代わりに単純なオイラー法で代用せざるを得ませんでした。

·RHS並列計算

PETScの行列ベクトル積を用いて、RHS計算を実装します。MicroMagは、磁場を2種のデータ構造として保持しています。一つは、有限要素数NMAX×物理次元数３の2次元行列m $NMAX,3$ です。もう一つはこれの一次元版Y $3*NMAX$ です。これをPETScデータ形式へ適用しました。この他に行列ベクトル積等へのPETSc適用を実装しました。この時点で、作業時間の制約により一旦プロジェクトは終結しました。

後続プロジェクトでの結果

2011年第一四半期に新たにプロジェクトが組まれ[4]、PETSc開発者から入手したCVODEインターフェイスのバグ修正パッチが適用されました。その後、XE6上の16ノード $384コア$ によるベンチマークが行われましたが、残存バグが残り、性能の改善が見られない状態で作業は一旦終了しています。今後の更なる性能調査が必要です。

謝辞

このプロジェクトは、nAG Ltd.が運営するHECToRの分散計算科学および工学 $CSE$ サービスの基に実行されました。英国の国立スーパーコンピューティング・サービスである、HECToR：英国リサーチ・カウンシル・ハイエンド計算サービスは、リサーチ・カウンシルを代行するEPSRCが管理しています。そのミッションは英国学術界の科学および工学の研究支援です。HECToRスーパーコンピューターは、UoE HPCx Ltd.およびnAG Ltd.のCSEサポートサービスにより管理運営されています。

文献

[1]	The SUNDIALs library https://computation.llnl.gov/casc/sundials/main.html
[2]	The Metis library http://glaros.dtc.umn.edu/gkhome/views/metis
[3]	The PETSc library http://www.mcs.anl.gov/petsc/petsc-as/
[4]	The Edikt project http://www.edikt.org.uk/edikt2/