＊ここに掲載するのは、エジンバラ大学 EPCCのChris M. Maynard博士らによるHECToRレポート「MicroMagnetic modelling of naturally occurring magnetic mineral systems: II, Dr. Chris M. Maynard, Mr. Paul J. Graham, EPCC, The University of Edinburgh, 21/02/2013」を要約したものです。
概要
鉱物の磁気特性と古磁気記録の安定性の正確な決定に関して、多くの地球科学的問題が存在します。例えば、地磁気の研究には地表面上とその上部で観測された場の変化に大きく依存します。直接観察の記録は200年前に遡りますが、より詳細な分析は岩石内に形成された古代の記録に依存します。古代の場の直接的な記録の地質学的解釈からの発見が50年以上前にありました。それは海洋と大陸のプレートが移動し続けていることです。今日の磁気鉱物学の進んだ理解は、大陸ブロックの運動と回転の詳細スケールの解析だけでなく、岩石の定値温度と熱履歴の解析も可能にします。さらに磁気鉱物学は古気候変化を決定する際の代替手法として利用され始めました。
こうしたアプリケーションは、鉱物の微小構造や化学的性質、グレイン形状およびグレイン間の磁気相互作用によって、その磁気的性質がどのように変化するかの理解に依存します。その数値的な微小磁気的問題の検討は、理想的な人工的記録媒体を考える場合よりも自発的磁化鉱物の複雑性と多様性により困難さが増します。環境磁場研究の発展は、天然磁化鉱物系で生じる高い忠実性を持つ磁化記録プロセスが基礎です。通常の磁化鉱物が局所的な地磁気方向と強さを記録し、地質学的時間スケールでそれを保持することを可能とする基礎的なプロセスに関しては、完全な理解に程遠い状態です。天然鉱物内の記録の信頼度を測るには、その記録プロセスに対するより深い理解が必要です。
平衡磁区構造はLandau-Lifshitz-Gilbert運動方程式で決定されます。不規則な構造を扱う場合に有用なのが有限要素法です。それは磁区構造がモデル化された精度に敏感なためです。このプロジェクト以前のmicroMagはFortran90で記述されたシリアルコードでした。
アルゴリズムの詳細
数値的な微小磁化モデルは、有限要素（FE)メッシュのN個の頂点に置かれた3Dカルテシアンベクトルで磁化を表現します。離散化誤差を最小にするために四面体要素を用います。計算セルのサイズは交換長と呼ばれる微小磁化的物理量で支配され、これは最大セルサイズを制限します。
平衡磁区構造は系の全磁化自由エネルギーを最小化して得られますが、これはLandau-Lifshitz-Gilbert（LLG)運動方程式を積分することで達成されます。そのダイナミクスは各時間ステップでLLG運動方程式を解くことで得られます。有効場の計算には、計算集約的な非局所的場静磁場の計算や、局所交換場、異方性場および外場の計算が含まれます。また、全有効場の計算も演算集約的です。交換場と異方性場の計算には有限体積法が用いられますが、これは本質的には全磁化ベクトルに剛性マトリックスを乗算することで計算されます。
静磁場はスカラーポテンシャルから計算されます。ここでは、磁化の発散に対するポアソン方程式と境界の頂点に対するラプラス方程式の求解が必要です。長距離相互作用の近似には境界要素法を用います。これら方程式は、線形有限要素解に対するガラーキン法を用いて離散化します。求解には標準的な共役勾配（BG）法を用います。ポアソン方程式はスパース行列なので、その手法の効率性が重要となります。また四面体有限要素基底関数は、残る3成分場を簡潔な行列ベクトル積で計算する効率的な手法を提供します。初期有限要素基底関数の効率的な構築のためには、専用のスパース行列構築が必要です。
長距離相互作用に対する境界要素近似で生じる行列は、対象物質の境界上の頂点数Nの2乗サイズです。実際のコードにおいては、この密行列は各プロセッサへ効率的に分散可能です。剛性システムのLLG運動方程式による常微分方程式（ODE）の積分には陰解法ソルバーを用いるため、スパースな3N線形系に対する近似逆行列計算が必要になります。時間ステップサイズが剛性で制御されるので陽解法では非効率です。この陰解法は可変係数を持つ多段後退オイラー法をベースにした可変次数ソルバーです。
作業計画
当初の作業計画は以下の作業項目からなる6人月のプロジェクトです：
作業項目1：効率的かつ高精度の剛性マトリックスコードの構築[1人月]
作業項目2：行列の並列生成[2人月]
a：剛性行列のスパース性の決定
b：PETSc行列データオブジェクトの生成
c：行列の計算とPETSc行列オブジェクトの配置
作業項目3：メッシュ生成用にParMETISを実装[2.5人月]
a：ParMETIS呼び出しの実装
b：他のメッシュ分割サポート実装
c：ベンチマーク測定
ただし、作業項目2において、行列オブジェクトの初期化は計算全体に掛かる時間の内ごく僅かでした。
またプロファイル測定により、ボトルネックの一つであるデータ変換処理が見つかりました。このデータ変換ルーチンの最適化も作業項目に追加しました。
実装
既存コードはMPIで並列化されていますが、初期化フェーズがほぼシリアル実行です。並列性能は非常に悪く、結果の正確性にも問題がありました。
·結果の検証とプロファイリング
シリアルコードと並列コードの結果の一致が不十分でした。シリアルコードは不完全コレスキー分解をプレコンディショナにしたCGソルバーを用いています。プレコンディショナを用いないと結果が異なります。特にプロコンディショナを用いたCGソルバーの結果ベクトルは、並列コードに比べて有効数字で3桁異なり1％の相対誤差を持ちます。プレコンディショナを用いないと12桁有効数字が異なります。プレコンディショナの有無で結果は異なるべきではないですが、こうした系では丸め誤差が重要になります。これは条件数の多い線形系や、CGソルバーでなくてプレコンディショナ自身の計算で丸め誤差が大きくなる場合に、よく見られる現象です。多くの検証作業の結果、この事象は数値的丸め誤差であることが判りました。
PETScの自己プロファイル機能を用いたところ、行列ベクトル積のコストが大きいことが判りました。さらにCrayPATでプロファイルを調査したところ、MPI通信ルーチンが大きなコストを持つことが判りました。これは剛性マトリックス計算で大きな通信負荷があることが判りました。
PETScはデフォルトで、行列をプロセッサランクに渡って列で分散します。スパース行列は、プロセッサに配置される対角ブロックの列セットと非対角要素という行列要素の区別があります。対角ブロック上の行列ベクトル積は完全にローカルです。非対角部については通信が発生します。既存コードは剛性マトリックスをFrotran配列としてシリアル実行で構築していましたが、PETScはデータオブジェクトを分散配置して構築します。これは極めて非効率的な処理ですが、それ自体は性能劣化の原因ではありません。
性能劣化の原因として考えられたのはメッシュ分割です。メッシュはどのようなジオメトリーに対してもサードパーティライブラリのCubitで生成されています。これはシリアル実行コードです。その出力は、頂点座標と四面体結合情報、つまり四面体頂点番号と境界上の要素とそれらが結合する面数のリストです。microMagは並列実行においてこのメッシュの分割について何も意識しません。単純にファイルを読み込み、メッシュファイルの値の剛性マトリックスを作成します。PETScはこの行列を列で分割し、結果的にメッシュの分割が行われます。これが意味するのは、互いに結合する四面体が多くのプロセッサに渡って分散的に分割されることです。つまり行列の対角部より非対角部が多く存在することになり、結果として通信オーバーヘッドが大きくなります。通信オーバーヘッドを削減するには、より効果的な分割が必要です。
これがParMETISによる分割を導入する理由です。これは新たなプログラムを独立に実装することで行いました。こうしてmicroMagは多くの修正なしに新たなメッシュファイルを読み込み、そのテストも容易になります。その後にこの領域分割コードをmicroMagへ組み込みました。
·ParMETIS、PETScと領域分割
PETScは、サードパーティーの領域分割ライブラリとの簡便なインターフェイス機能を持ちます。ParMETISはParMATISについては、デフォルト分割ライブラリとしてPETScインターフェイスを用い、後続の開発でのコード修正量を削減します。分割は頂点でなく要素で行います。
この工夫はメッシュを都合の良い順序で分割し、結果として行列の並びも適した形になり、通信は最小化されるというものです。ここで下記の2つの新しいサブルーチンを作成しました：
repartVert：これは独立したプログラムで、モデルデータを読み込み、サブルーチンreorderを呼び出して、再分割したデータを書き出します。
reorder：このルーチンでは、最初に頂点の結合情報データを生成し、PETSc/ParMETISで必要な隣接リストの作成に用いられます。次にPETScのデータ構造が生成され、領域分割処理が呼び出されます。次に新しい分割データ（これは頂点リストとそのプロセッサ上の新しい配置です）が生成され、ランク0のプロセスへ集約され、分割処理が適用されて頂点と四面体要素が適切に順序付けられます。最後にこれら頂点と四面体要素は適切なプロセッサへ配置されます。
現在MicroMagには論理変数の分割処理オプションが装備され、これがtrueであればreorderルーチンが呼び出されて、さらに新しい分割における境界ノードの再定義を行います。
·データ変換処理の最適化
残念ながらMicroMagの並列コードは高速化を示さず、場合によっては遅くなっていました。詳細に調査したところ、この原因はデータ変換ルーチンにあることが判りました。MicroMagは磁場を2つの異なるデータ構造で保持します。その一つはサイズ（NMAX,3）の2次元配列で、NMAXは頂点数、3はカルテシアン空間の物理次元を表します。第2の構造は、1次元化された配列3∗NMAXです。コード内では幾度もこの2つのデータ構造が入れ替わる中で、Sundialsソルバーは1次元配列を用います。この変換は極めて高コストで、既存のコードでは複数のPETSc VecScatter呼び出しを要求しており、呼ばれた回数分の大きな通信オーバーヘッドが生じます。さらにこのオーバーヘッドはコア数が増えるにつれて増加し、他のコード部分の高速化による利得を打ち消してしまいます。
コード内の全てを1次元配列に書き直すには、根本的な書き直しが必要なため、別の方策を用いました。調査の結果、この2つの配列表現の変換は、コア間のデータ受け渡し無しに各プロセッサ上で局所的に実行できることが判り、新しい変換ルーチンを作成しました。
測定の結果、10955要素のテストデータでは8コアで2倍の並列による高速化が示されました。これ以上にコア数を増やすと効率は減少していきます。
さらに大規模な100245要素ケースでは、初期化処理がシリアルケースでは約3％を占めます。時間の制約上、収束まで実行せずに1,000時間ステップを測定しました（シングルコアでは収束に10時間かかります）。初期化を省いた実行時間では、16コア上での高速化は5倍に達しました。
謝辞
このプロジェクトは、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/ |