（様式第 20）成果報告書

共同利用（産業利用トライアルユース：先端研究施設共用促進事業 『みんなのスパコン』TSUBAME によるペタスケールへの飛翔） 成果報告書 平成２３年度 課題種別

利用課題名 個別要素法を用いた粉末充填シミュレーションプログラムの並列化とその評価
英文：Parallel benchmark for DEM code in the process of compacting metal powder
大塚 順
Jun Otsuka
住友電気工業株式会社 解析技術研究センター
Sumitomo Electric Industries, Ltd.
http://www.sei.co.jp/index.ja.html
邦文抄録（300 字程度）
個別要素法のシミュレーションプログラムに OpenMP 並列化を適用した。接触力の計算と粒子の座標・速度の更新に対して並列化を実施した。ベンチマークの結果、90 万粒子系で単体計算の 2.92 倍の高速化を達成した。本プロジェクトの取り組みにより、実プロセス規模の解析が可能になった。
英文抄録（100 words 程度）
We adapted the OpenMP parallel method to simulate the Discreet Element Method (DEM)
simulation program. We did parallel calculations for contact force calculation and updating thecoordination and velocity of the particles. We achieved a 2.92 times speed increase for parallel
calculations, compared with single calculations, in the case of doing a 900,000 particle simulation.
An analysis of the real process scale became possible through this simulation project.

Keywords: 個別要素法, DEM, 並列化, 高速化, 粉末冶金

背景と目的を目的とし、計算プログラムに OpenMP 並列化を適用

原料に粉末を用い、これを添加物と混合、成形してした。その結果、単体計算の 2.92 倍の高速化を達成で最後に焼結することで製造される粉末冶金製品は、各き、大規模モデルをシミュレーションできるようになっ種の機械装置に部品として組み込まれており、現代のた。

我々の生活において必要不可欠である。粉末を金型に充填する給粉工程におけるプロセス技術が製品の性概要能に大きく影響するが、これまで暗黙知に依存する部個別要素法を高速化するため、OpenMP 並列化を分が大きかった。新たな製品開発やさらなる生産性向適用し、その評価を行った。単体計算で 9.5 割の計算上を図るにはシミュレーション技術の開発が重要な課時間を占める接触力計算と粒子の座標・速度の更新に題である。そこで、有力なシミュレーション手法として、対して並列化を実施した。ベンチマークのため粉粒体をお互いに相互作用する個々の粒子の運動をコンピュー貯蔵装置から排出するときの粉体挙動をシミュレーショタ上で逐次追跡する個別要素法によるプロセス設計がンし、90 万粒子系で単体計算の 2.92 倍の高速化を確検討されている。

認した。この高速化によって、従来 1s 間のシミュレーシ個別要素法では粒子の挙動をニュートンの運動方程式を数値積分することで求めるが、その時間刻みは非ョンを実施するために 24 日かかっていた計算が 8 日間で完了できるようになった。

常に小さいことに加えて、シミュレーションモデルを現実に近づけるためには扱う粒子数を増やす必要があるた
め、膨大な計算時間が必要になる。

結果および考察
個別要素法は、粒子に働く力とニュートンの運動方

本プロジェクトでは、個別要素法で現実を反映した膨程式を陽解法で繰り返し計算する手法である。フローチ

大な粒子数を TSUBAME で計算できるようにすることャートを図 1 に示す。膨大な計算時間を短縮するため、

（様式第 20）成果報告書
OpenMP 並列の適用により高速化を実施した。並列化

ここまでの変更を施したプログラムで 8 並列計算したを適用した単体計算プログラムは、参考文献[1]に紹介と き の プ ロ フ ァ イ ル を 図 2 に 示 す 。 図 は 、 VTune

されているプログラムの様に計算フローを忠実にプログ

Amplifier の Lock & Wait 機能の解析結果であり、スラミングしている。

レッド間通信の並列性の情報を時系列で示している。

まず始めに、単体計算のプロファイルを取得し、ボト粒子－壁間の接触力計算以外では、1 step の処理全ルネックを特定した。プロファイルの取得には Intel体にわたって並列化できていることがわかる。粒子間VTune Amplifier を使用した。計算条件は、粒子数はの接触力計算の前半と後半において並列計算と単体300 とし、初期状態は密な配置とした。解析時間は計算が 8 回交互に現れているが、これはスレッド毎に1.0(s)、時間ステップは 2.5e-7(s)。各サブルーチンの計用意した接触力の計算結果を格納するための配列に算時間が全体に占める割合を調べたところ、接触力の対してゼロクリアおよびスレッド間の総和を計算すると計算が 0.8、粒子座標と粒子速度の更新が 0.15、接触ころで同期による待ちが発生していることを意味してい候補の判定が 0.05 であった。力の計算では粉末を構る。

成する全ての粒子に対して粘弾性モデルで表現した接今回実施した並列化の効果を評価するため、粉粒体触力を計算するため、計算フローの中で最も計算コストを貯蔵装置から排出するときの粉体挙動をシミュレーが高い。並列化によるコンピュータの性能効果についションした。シミュレーションモデルを図 3 に示す。粒子て示したアムダールの法則より、並列効率の理論値は数は 90 万であり、粉末パラメータは測定値に基づき設

接触力計算のみを並列化した場合が 8 並列で 3.3 倍、定した。並列数を変更した時のスピードアップを図 4 に接触力計算と位置・速度の更新を並列化した場合が示す。スピードアップは、単体計算の時間÷並列計算5.9 倍となる。の時間とした。PGI コンパイラを使用した場合、8 並列粒子間の接触力計算に OpenMP 並列化を適用した。

計算でスピードアップが最大となり、単体計算の 2.92

接触力計算は、粒子番号と接触候補番号の２重ループ倍となった。一方、Intel コンパイラを使用した場合、並である。粒子番号のループに対して、OpenMP 並列し列数を増やした時の性能低下が顕著だった。今回達成た時に各スレッドが分担する接触力計算の負荷を均すした 2.92 倍に対して、アムダールの法則による理論値ため、schedule(dynamic)オプションを使用した。このが 3.3 倍なので、理論性能に近い並列効率を達成できことで、粗粒を多く担当するスレッドで、細粒よりも接触ている。本プロジェクトの取り組みによって、並列化適数が多いことが原因で計算負荷が大きくなっていた問用前には 1s 間のシミュレーションを実施するために単題を回避できた。また、各スレッドで計算した接触力を体計算で 24 日かかっていたのが、8 並列計算で 8 日間足し合わせるとき、メモリアクセスが競合するのを避けに短縮できた。

るため、スレッド毎に異なるメモリ領域を確保するように個別要素法によるシミュレーション事例を示す。小型した。この計算手順は、①接触力の計算結果を格納すの粉末貯蔵装置(65000 粒子系)でシミュレーションした。

るための配列をスレッド数分用意する、②スレッド毎に粉末断面のスナップショットを時系列で図 5 に示す。実異なる領域に加算する、③スレッド間で総和する。一方、験と同じように”あり地獄”が形成されることを確認でき粒子－壁間の接触力計算の並列化についても検討した。

たものの、粒子－粒子の接触に比べて接触数が少なく正味の負荷が小さいため、計算時間の短縮が並列化のオーバーヘッドを下回る。そのため、粒子－壁間の接触力計算には並列化を適用しないことにした。
並列化の効果を大きくするため、粒子の座標・速度の更新についても並列化を適用した。ここでは、不要な同期を除くため、nowait 指示句を挿入した。

まとめ、今後の課題
簡単に実施できるプログラム変更の範囲では、一定レベルの並列効率を達成でき、実際の計算に役立つレベルに達することができた。
今後の課題として、高速化については、単体計算性能の改善がさらなる高速化のために必要となる。現行

（様式第 20）成果報告書のプログラムは、接触力計算のループが長すぎるためコンパイラの解析に不都合であること、条件分岐が残っているため命令のスケジューリングが非最適であることが改善余地と考えている。プリフェッチやパイプライン処理の有効利用で高速化できそうと見込んでいる。具体的な解析事例への取り組みについては、実際のプロセスの改良に役立つ知見を得ることを目指す。
TSUBAME で多数粒子のシミュレーションを実施することで、より現実に近いモデルで紛体挙動を解析できると考えられる。解析により得られる知見は個別要素法

図 3 大規模計算のベンチマークモデル

で紛体貯蔵装置を解析するときに重要となる学術的にも価値のある成果になると期待できる。

本研究は、神戸大学と兵庫県の「企業を牽引する計算科学高度技術者の要請」プロジェクトの支援を受けて実施された。ここに記して、感謝の意を表したい。

参考文献[1] 粉体工学会 編(1998), 『粉体シミュレーション入門―コンピュータで粉体技術を創造する』, 産業図書
初期設定

図 4 並列計算によるスピードアップ

2粒子衝突モデル

接触候補の判定

垂直力

接触力の計算
no

せん断力

全粒子終了
yes
粒子位置の更新

垂直方向成分せん断方向



スライダ（摩擦）

粒子速度の更新
noダッシュポット（粘性）

ステップ終了yes
終了

バネ(弾性)I    r   kr   022

ダッシュポット

バネmu u  ku  0スライダ ダッシュ バネ
ポット

図 1 個別要素法の計算フロー

図 2 8 並列計算実行時のプロファイル

図 5 小型粉末貯蔵装置のシミュレーション