プログラムの高速化・並列化サービスの事例

＊ここに掲載するのは、インペリアルカレッジ・ロンドンのG. Marketos博士によるHECToRレポート「Opening up High Performance Computing to the Discrete Element Method User Community, G. Marketos, Civil and Environmental Engineering Department, Imperial College London」を要約したものです。

概要

個別要素法：Discrete Element Modelling (DEM)は分子動力学（MD)と多くの点で類似しています。DEMの応用として、土木および地盤工学における材料の粒状応答シミュレーションが有ります。しかしながら現在のシミュレーションは、そのHPCバージョンのアプリケーションで要求される機能の欠落のために制限されています。大規模原子分子超並列シミュレータLAMMPS(http://lammps.sandia.gov/)は、MDシミュレーションに広く利用されています。物理的かつアルゴリズム的なアナロジーから、LAMMPSはDEMシミュレーションの極めて良いプラットフォームの一つです。粒状物を対象にしたLAMMPSはそうしたコードの一つで、高並列化されていますがDEMに必要な機能が欠けています。

このプロジェクトの目的は"granular LAMMPS"の開発です。新たに追加するC++で実装するコードの機能は、2種の境界条件と2粒子間の結合に対する新しい接触モデルです。この開発により、HECToR上で高いスケール性を持つDEMシミュレーションを可能にします。

主要な目的と成果

現在英国の個別要素法（DEM）ユーザコミュニティでは、DEMシミュレーションはシリアルコードを用いる場合がほとんどです。これまでのEPSRCでの研究からGranular LAMMPSが大規模DEMシミュレーションにとってオープンソースの最適な並列コードであることが判っています。しかしながらHPC環境をフルに活用するには、Granular LAMMPSに境界条件を導入することが必須です。まず最初に固定壁の移動を許す、つまり境界の移動を制御する応力制御シミュレーションを可能にしました。次に、応力制御固定境界の実装を行いました。以上は最初の作業項目WK1で行われました。WK1の残りの作業で膜境界が実装されました。これはラテックス膜に覆われた材料実験シミュレーションをより現実的な物にします。この機能は他の如何なるオープンソースコードでも実現されておらず、必要となるボロノイ図形計算は高コストですが、DEMコミュニティでGranular LAMMPSとHPCの利用を促進することになると考えられます。

LAMMPSに開発が必要な別の機能として、DEMシミュレーションで用いられる粒子間相互作用モデルが有ります。分子動力学シミュレーションにおいては膨大な努力がこうした開発に充てられて来ました。しかしながらDEMシミュレーションではより少ない基本モデルのみが必要となります。それは線形および非線形弾性接触モデルで、各粒子間相互作用は個別の"Pair"クラスに実装されます。新しい非線形弾性接触バネモデルが、摩擦する粒子のDEMシミュレーションに最も適切なHertz-Mindlin弾性接触モデルとして作成されました。また、第2の作業項目の目的は、粒子間の結合相互作用を追加することです。これは多孔質岩や砂岩といった結合した粒状材料モデリングを可能にし、地質学のDEMユーザにHPC利用への道を切り開きます。これは工業的にも興味深く重要な様々な現象の研究にも役立ちます。例としては、石油工業でのドリルあるいは間隙流体抽出により生じる力学的負荷による砂岩の変形が有ります。

新しい境界条件の実装

DEMの開発を後押ししてきた分野の一つとして地盤工学が有ります(Cundall, 2001, Cundall and Strack, 1979)。最初のアルゴリズム(Cundall and Strack, 1979)が提唱されて以来DEMの利用は拡大し、様々な科学、工学分野に適用されてきました。しかしながら地盤の振る舞いの基礎的な理解のための地盤工学実験シミュレーションは、科学雑誌に掲載されるDEMシミュレーションの大きな割合を占めています(例えば、Cui et al, 2007, Marketos and Bolton, 2010, Shen, 2012)。こうしたシミュレーションは、如何に標本材料の応答が実験室で最も良く測定されるかという問題に重要な示唆を提供し、サンプル内の力の伝達やずれのパターンに境界が与える影響（これが粒状材料の実験の解釈を複雑にします）を解析するのに役立ちます。
作業開始前は、granular LAMMPSは固定壁を持つサンプルのみが可能でした。これは正弦曲線的な動きや、あるいは境界面に一致した（せん断）方向の動きのみが可能です。これは極めて強い制限で、DEMシミュレーションではより多くの境界制御が必要です。
そこで、LAMMPSのFixWallGranクラスに平坦境界を実装しました。新しい平坦境界利用には、LAMMPSの入力スクリプト"fix wall/gran"に追加された引数を指定して用いることが出来ます。

膜境界の実装

工業および研究における地盤力学の実験室での多くの物理実験では、サンプルは垂直に圧縮されるように柔軟で変形可能な膜で覆われます。この面外方向の自由度は上述の固定境界では表現不可能で、サンプル内のずれを過剰に制約します。例えば、実験で確認されている重要な多くの現象（例えば帯状せん断）は固定境界のシミュレーションでは観測されず、粒状材料のせん断特性に誤った予測を導きます。より正確に実験条件を模倣するには、Granular LAMMPSに擬似的な膜境界を実装する必要があります。

膜で覆われたサンプルのシミュレーションは多くなく、そうした機能を持つオープンソースや商用コードは存在しません。さらにこれが可能なコードは高コストのため、大量の粒子を用いた3次元シミュレーションはHPC環境でしか実行できません。このためGranular LAMMPSへの膜境界の実装は最重要項目であると判断されました。

擬似的な膜に対する多くの計算スキームが文献で提案されています(Cheung and O’Sullivan, 2008)。これらは皆、サンプルの境界から一定の距離内に存在する粒子の等価な重みあるいは領域サイズの計算が仮定されており、これらは境界の粒子に掛かる外力の計算に用いられます（通常、力はこの領域サイズと特定の圧力の積になります）。数値的な膜を実装する方法は概念的に3つの部分に分けられます：

• A:境界粒子の選択
• B:各境界粒子の等価領域サイズの計算
• C:膜力の計算

当初は膜アルゴリズムは2次元DEMシミュレーションに対して開発されました。この場合、3次元ほどは複雑でなく、境界粒子の特定は直截的です。Thomas and Bray (1999)およびCheung and O’Sullivan (2008)は境界粒子を繋ぐ線分長を用いて、これらを境界圧力値に乗算して膜力を計算しました。Kuhn (1995), O’Sullivan (2002), Cheung and O’Sullivan (2008)は、最外側粒子の重心を平面あるいは円筒形表面に射影することで、彼らのやり方を3次元シミュレーションへ拡張しました。次に、その平面上に射影された点の2次元ボロノイ図形を計算すれば、各粒子の重み(あるいは領域サイズ)が判ります。最後に膜力を、膜圧力値に各境界粒子の領域サイズを掛けて計算します。

上述は膜の一側面を表しますが、射影面に垂直な膜力のみに適用されるという点が欠点です。ここで、この条件を緩和する修正アルゴリズムを提案します。新しい膜の実装は、膜が曲がったり、面外に膨らんだりする場合や、力が膜射影面に垂直な方向から大きくずれる場合にも適用可能です。

このスキームでは、上記のAパートとBパートの出発点は粒子系の重み付きボロノイ図形(Laguerreあるいはradical Voronoi graphとも呼ばれます。Okabe et al, 2000を参照してください)の計算です。膜計算の全ては、LAMMPSのFixベースクラスから導かれるself-containedクラスFixMembraneGranとして実装されました。粒子に掛かる力は、接触バネにより粒子に掛かる力に対して粒子間接触モデルを適用した後に実行される関数"post_force()"を通して膜力分増加します。

この実装は膜力計算のためにボロノイ図形を利用します。DEMシミュレーションでは粒子は異なるサイズ（半径）を持つことが多く、半径で重みを付けたボロノイ図形の計算が最も適しています。こうしたグラフは以前にDEMで用いられていましたが(Chareyre et al, 2012, Rycroft et al, 2009等)著者の知識では、膜力計算に利用されるのは初めてです。Voro++はUC BerkeleyのChris Rycroft博士がメンテしており、LAMMPSへのリンクも容易に見えます。更に、Voro++はセルベース計算を基に、各粒子個別にボロノイセルを計算します。そのため、シミュレーション領域を空間分割して並列化を行うGranular LAMMPSに良く適合します。

LAMMPSユーザフォーラムでの議論からLAMMPS開発者によりVORONOIパッケージがリリースされました。これは粒子グループに対する通常のボロノイ図形計算コード（ComputeVoronoiAtomクラス）も含んでいます。このパッケージに膜を実装した新しいクラスFixMembraneGranが追加され、粒子グループに対する重み付きボロノイ計算をする新しいクラスComputeLaguerreAtomも実装されました。VORONOIパッケージはビルド時に利用を選択可能です。

膜計算コードの試験：
新しい膜計算コードを、6面全てを固定壁で覆った1600粒子で検証しました。サンプルを平衡化した後に、6面全てを膜に置き換えて、膜圧力は固定壁に適用した値と同じにしました。LAMMPSの繰返し後、サンプルは応力を維持し続けている様子が見られました。

Granular LAMMPSの接触モデルの拡張

インペリアルカレッジのLAMMPS開発チームメンバのKevin Hanley博士とXin Huang氏は、クラスPairGranHertzHistory内の非線形弾性相互作用は通常DEMシミュレーションで用いられるHertz-Mindlin相互作用(Mindlin and Deresiewicz, 1953、あるいはItasca, 2003を参照してください)とは異なることを指摘しました。法線力(即ち、2つの接触粒子の中心を結ぶベクトル方向の力)の計算はHertzian理論(Johnson, 1985)で可能ですが、それに垂直な面内の相互作用(即ち、せん断あるいは接線方向の相互作用)は異なります。そこでHertz-Mindlin相互作用を実装した新しいペアクラス"PairGranShmHistory"の実装が好ましいと判断しました。既存の修正よりも新規開発の方が後方互換性に優れると考えられます。

新しい結合モデルの実装

多くの結合モデルがDEMシミュレーションで提案されてきました(Potyondy and Cundall, 2004, Utili and Nova, 2008, Weatherley, 2009)。これらの内でPotyondy and Cundallによるモデルが最も頻繁に利用されています。インペリアルカレッジのLAMMPSユーザグループメンバはこれまで多くのシミュレーションでこのモデルを用いており(Marketos and Bolton, 2009, Hanley et al, 2010, Cheung and O’Sullivan, 2008)、Granular LAMMPSへ追加するモデルとしてはこのモデルが最も適切です。

Granular LAMMPSへ新しく結合モデルを実装するには2つのやり方が有ります。それはbondスタイルかpairスタイルのどちらかです。bondスタイルは、結合粒子間の相互作用力計算に自己完結型スキームの利点を生かせます。しかしながらLAMMPSにおいてはこのようなbondスタイルは分子動力学でしか使えません。LAMMPSに個別要素法用のbondを適用する際の問題は、結合状態に対して結合毎の追加メモリーが必要となることです。このメモリーは、粒子がプロセッサを跨ぐ際にコピーしなくてはならず、結合リストから結合が追加あるいは削除される際には注意深い取り扱いが必要です。

よって、LAMMPSには"Pair"クラスとして結合モデルを実装する選択をしました。これにより、結合状態情報を保管するのに接触リスト(LAMMPSでは"neighbor list"と呼ばれる)を利用する利点が生じます。粒子の"neighbor list"のやり方はLAMMPSにおいてよくテストされ、過去何年にも渡り多くのユーザによりデバッグされて来ました。よってこれは開発時間を考えるとこれが結合モデルを実装する安全な方法と言えます。

結合相互作用モデルの検証：
Cheung (2010)はその文献において、Potyondy and Cundall (2004)のモデルを実装した商用コードPFC3Dにおける2つの結合モデルの振る舞いの調査結果が詳細に述べられています。ここで、結合モデルの実装検証のために、Cheungの文献で記載されたものと同等の2粒子シミュレーションを実行しました。これはLAMMPSのシリアルバージョンの実行です。このシミュレーションにおいて、丁度接触する一で半径0.5mの2つの粒子を作成しました。粒子Aの測度と角速度はゼロに設定し、粒子Bを接触バネと接触位置に設定された平行結合の両方で動くようにしました。 結果はCheung (2010)のものと非常によく一致しています。

新しい機能を用いた大規模シミュレーション

granular LAMMPSの性能確認のため、HPC環境で大規模シミュレーションを実行しました。1.95百万粒子を重力下で降下させ、立方体容器内で平衡化させました。用いた粒子特性(例えば、粒子サイズ、密度、剛性)は、EPSRCプロジェクトにおいてブリストル大学の実験研究者が用いた細粒ガラスビーズのサンプルです。サンプル容器は実験容器と同じに設定しました。容器上部の余分な約500,000粒子を除いて平坦にし、全粒子数1,428,532でシミュレーションを行いました。時間ステップは既存のガイドライン(Itasca, 2003)に従い、応力シミュレーションを、6つの側面に対してレート200kPa/secを掛けて行いました。

次に100kPa下のサンプルを、ベンダーエレメント実験のアナログである点源からの波伝播シミュレーションに用いました。この実験は地盤のせん断剛性の決定に通常良く用いられる手法ですが、生成する信号の理解が課題となっています。 結果として、HPC環境において以前より大規模なシミュレーションが可能になりました。

成果

上述のシミュレーションはこの分野のマイルストーンを表しています。以前に比べ30倍の粒子を扱っています。全シミュレーション時間は、サンプル調整を含めて6週間を超えない程度です。サンプル調整にはImperial College Londonのクラスターマシン（72 cores）を用いて、約1か月かかりました。次の応力シミュレーションには、1MPaに達するのに約48CPU時間掛かりました。この実シミュレーション時間は5秒で、20百万時間ステップ以上掛かりました。これまでに地盤力学分野では百万粒子を超える解析は他になく、本計算はこの分野の大きな成果となりました。

謝辞

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

文献

[1]	Chareyre B, Cortis A, Catalano E & Barthelemy E (2012). Pore-Scale Modeling of Viscous Flow and Induced Forces in Dense Sphere Packings. Transport in Porous Media, Vol. 92, No. 2, pp 473-493.
[2]	Cheung G & O'Sullivan C (2008). Effective simulation of flexible lateral boundaries in two- and three-dimensional DEM simulations, Particuology, Vol. 6, pp. 483-500.
[3]	Cheung, G. (2010). Micromechanics of sand production in oil wells. Ph.D. thesis, Imperial College London.
[4]	Cundall, P. (2001).A discontinuous future for numerical modelling in geomechanics? Geotechnical Engineering Proceedings of the Institution of Civil Engineers 149 (1), 41.47.
[5]	Cundall, P. and O. Strack (1979a).A discrete numerical model for granular assemblies.Geotechnique 29 (1), 47.65.
[6]	Cui L, O’Sullivan C & O’Neill S (2007). An analysis of the triaxial apparatus using a mixed boundary three-dimensional discrete element model. Geotechnique, Vol. 57, No. 10, pp. 831-844.
[7]	Hanley (2013) personal communication.
[8]	Hanley, K., O’Sullivan, C., Byrne, E. P. , Cronin, K. (2012) Discrete element modelling of individual agglomerates of infant formula. Particuology, 10(5), pp 523-531.
[9]	Itasca Consulting Group Inc. (2003). PFC3D: Particle Flow Code in 3 Dimensions, Version 3.0, Minneapolis, USA.
[10]	Johnson, K. (1985). Contact Mechanics. Cambridge University Press.
[11]	Kuhn MR (1995). A flexible boundary for three-dimensional DEM particle assemblies. Engineering Computations, Vol, 12(2), pp. 175-183.
[12]	Marketos G & Bolton MD (2009). Compaction bands simulated in Discrete Element Models, Journal of Structural Geology, Vol. 31, pp. 479-490.
[13]	Marketos G & Bolton MD (2010). Flat boundaries and their effect on sand testing, International Journal for Numerical and Analytical Methods in Geomechanics, Vol. 34, pp. 821-837.
[14]	Mindlin RD & Deresiewicz H (1953). Elastic spheres in contact under varying oblique forces. ASME Journal of Applied Mechanics Vol. 20, pp. 327-344.
[15]	O'Donovan J, O'Sullivan C & Marketos G (2012). Two-dimensional discrete element modelling of bender element tests on an idealised granular material. Granular Matter, Vol. 14, pp. 733-747.
[16]	O’Sullivan, C (2002). The application of discrete element modelling to finite deformation problems in geomechanics. Doctoral Dissertation, University of California, Berkeley.
[17]	O’Sullivan C (2011). Particulate Discrete Element Modelling: A Geomechanics Perspective. Taylor and Francis, London.
[18]	Okabe A, Boots B, Sugihara K & Chiu SN (2000) Spatial tessellations: concepts and applications of Voronoi diagrams. 2nd Edition, Wiley.
[19]	Potyondy, DO and PA Cundall (2004). A bonded-particle model for rock. International Journal of Rock Mechanics and Mining Sciences, Vol. 41(8), pp. 1329-1364.
[20]	Rycroft CH (2009). Voro++: A three-dimensional Voronoi cell library in C++, Chaos Vol.19, 041111
[21]	Shen C-K (2012) PhD Thesis, Imperial College London
[22]	Thomas P & Bray J (1999). Capturing nonspherical shape of granular media with disk clusters. Journal of Geotechnical and Geoenvironmental Engineering, Vol. 125(3), pp.169-178.
[23]	Utili, S. and R. Nova (2008). DEM analysis of bonded granular geomaterials. International Journal for Numerical and Analytical Methods in Geomechanics 32 (17), 1997.2031.
[24]	Weatherley, D. (2009). Esys-particle v2.0 users guide. Technical report, Earth Systems Science Computational Centre, University of Queensland.