(57)【特許請求の範囲】
【請求項１】
複数のノードを有する並列計算機のための複数種類の粒子が分布する領域に対する粒子シミュレーションプログラムであって、
前記領域を分割した複数の部分領域各々を各ノードに割り当てて、粒子シミュレーションのための複数のプログラムを各ノード毎に実行し、
前記複数のプログラムそれぞれの実行時間を、前記複数のプログラムそれぞれの処理対象粒子と同じ種類の粒子の個数と、前記複数のプログラムそれぞれの比例係数との積で表した計算式と、前記複数のプログラムそれぞれの実行時間の総和を、各種類の粒子１個当たりの計算コストと、各種類の粒子の個数との積の総和で表した計算式とを用いて、前記複数のプログラムそれぞれの実行時間と各種類の粒子の個数とから、各種類の粒子１個当たりの計算コストを求め、
各部分領域における計算量を、各種類の粒子１個当たりの計算コストと各部分領域内における各種類の粒子の個数とを乗算した値の総和により求め、
前記複数の部分領域それぞれにおける計算量が、前記複数の部分領域それぞれを割り当てられたノードが有する計算機資源に比例するように、いずれかの部分領域の領域境界の位置を変更し、
変更された領域境界によって決定された各部分領域を各ノードに割り当てて、前記粒子シミュレーションのための前記複数のプログラムを各ノード毎に実行する、
処理を前記並列計算機に実行させる粒子シミュレーションプログラム。
【請求項２】
前記並列計算機は、各ノードで実行された前記複数のプログラムそれぞれの前記実行時間のばらつきが所定値以内の場合には、前記領域境界の位置の変更は行わないことを特徴とする請求項１記載の粒子シミュレーションプログラム。
【請求項３】
複数のノードを有する粒子シミュレーション装置であって、
前記粒子シミュレーション装置は、複数種類の粒子が分布する領域を分割した複数の部分領域各々を各ノードに割り当てて、粒子シミュレーションのための複数のプログラムを各ノード毎に実行し、
前記粒子シミュレーション装置は、
前記複数のプログラムそれぞれの実行時間を、前記複数のプログラムそれぞれの処理対象粒子と同じ種類の粒子の個数と、前記複数のプログラムそれぞれの比例係数との積で表した計算式と、前記複数のプログラムそれぞれの実行時間の総和を、各種類の粒子１個当たりの計算コストと、各種類の粒子の個数との積の総和で表した計算式とを用いて、前記複数のプログラムそれぞれの実行時間と各種類の粒子の個数とから、各種類の粒子１個当たりの計算コストを求める計算部と、
各部分領域における計算量を、各種類の粒子１個当たりの計算コストと各部分領域内における各種類の粒子の個数とを乗算した値の総和により求め、前記複数の部分領域それぞれにおける計算量が、前記複数の部分領域それぞれを割り当てられたノードが有する計算機資源に比例するように、いずれかの部分領域の領域境界の位置を変更する変更部と、
を備え、
変更された領域境界によって決定された各部分領域を各ノードに割り当てて、前記粒子シミュレーションのための前記複数のプログラムを各ノード毎に実行することを特徴とする粒子シミュレーション装置。
【請求項４】
複数のノードを有する並列計算機が、
複数種類の粒子が分布する領域を分割した複数の部分領域各々を各ノードに割り当てて、粒子シミュレーションのための複数のプログラムを各ノード毎に実行し、
前記複数のプログラムそれぞれの実行時間を、前記複数のプログラムそれぞれの処理対象粒子と同じ種類の粒子の個数と、前記複数のプログラムそれぞれの比例係数との積で表した計算式と、前記複数のプログラムそれぞれの実行時間の総和を、各種類の粒子１個当たりの計算コストと、各種類の粒子の個数との積の総和で表した計算式とを用いて、前記複数のプログラムそれぞれの実行時間と各種類の粒子の個数とから、各種類の粒子１個当たりの計算コストを求め、
各部分領域における計算量を、各種類の粒子１個当たりの計算コストと各部分領域内における各種類の粒子の個数とを乗算した値の総和により求め、
前記複数の部分領域それぞれにおける計算量が、前記複数の部分領域それぞれを割り当てられたノードが有する計算機資源に比例するように、いずれかの部分領域の領域境界の位置を変更し、
変更された領域境界によって決定された各部分領域を各ノードに割り当てて、前記粒子シミュレーションのための前記複数のプログラムを各ノード毎に実行する、
ことを特徴とする粒子シミュレーションにおける計算機資源配分方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、粒子シミュレーションプログラム、粒子シミュレーション装置、及び計算機資源配分方法に関する。
【背景技術】
【０００２】
数値計算を用いて水、空気等の流れを調べる流体解析、又は圧縮されたゴム等の振る舞いを調べる弾性体解析において、流体又は弾性体のような連続体を粒子の分布で表現する粒子法が用いられることがある。粒子法の具体例としては、Smoothed Particles Hydrodynamics（ＳＰＨ）法、Moving Particles Semi-implicit（ＭＰＳ）法等が知られている。
【０００３】
分散メモリ型並列計算により粒子の状態を計算する粒子状態計算装置も知られている（例えば、特許文献１を参照）。シミュレーション空間内に分布する粒子に対して、同数の粒子を含むようにシミュレーション空間を分割する分子動力学計算装置も知られている（例えば、特許文献２を参照）。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特開２０１２−１２８７９３号公報
【特許文献２】特開平５−２７４２７７号公報
【発明の概要】
【発明が解決しようとする課題】
【０００５】
複数のノードを有する並列計算機を用いて、計算量が異なる複数種類の粒子が混在する空間に対する粒子法のシミュレーションを行う際、それぞれの粒子に対して計算機資源を効率よく配分する技術は知られていない。
【０００６】
なお、かかる問題は、流体解析又は弾性体解析に限らず、他の連続体を含む粒子シミュレーションにおいても生ずるものである。
【０００７】
１つの側面では、本発明は、並列計算機による複数種類の粒子に対する粒子シミュレーションにおいて、計算機資源を効率よく配分することを目的とする。
【課題を解決するための手段】
【０００８】
１つの案では、粒子シミュレーションプログラムは、複数のノードを有する並列計算機に以下の処理を実行させる。
（１）並列計算機は、複数種類の粒子が分布する領域を分割した部分領域を各ノードに割り当てて、粒子シミュレーションのための複数のプログラムを各ノード毎に実行する。
（２）並列計算機は、複数のプログラムそれぞれの処理対象粒子の種類と、複数のプログラムそれぞれの実行時間とに基づいて、各処理対象粒子に対する演算コストを求める。
（３）並列計算機は、演算コストと各処理対象粒子の粒子数を基に部分領域の領域境界の位置を変更する。
【発明の効果】
【０００９】
実施形態によれば、並列計算機による複数種類の粒子に対する粒子シミュレーションにおいて、計算機資源を効率よく配分することができる。
【図面の簡単な説明】
【００１０】
【図１】影響領域を示す図である。
【図２】領域分割を示す図である。
【図３】袖粒子を示す図である。
【図４】部分領域間における粒子の移動を示す図である。
【図５】複数種類の粒子が混在する計算対象領域を示す図である。
【図６】粒子シミュレーション装置の機能的構成図である。
【図７】計算機資源配分処理のフローチャートである。
【図８】並列計算機の構成図である。
【図９】ノードの構成図である。
【図１０】粒子シミュレーションのフローチャートである。
【図１１】サブルーチンデータを示す図である。
【図１２】サブルーチンの実行時間を示す図である。
【図１３】津波シミュレーションにおける粒子の分布を示す図である。
【図１４】コスト係数に基づく部分領域を示す図である。
【図１５】粒子の個数に基づく部分領域を示す図である。
【図１６】４個のノードにおけるサブルーチンの実行時間を示す図である。
【図１７】津波シミュレーションが進行した後の部分領域を示す図である。
【発明を実施するための形態】
【００１１】
以下、図面を参照しながら、実施形態を詳細に説明する。
実施形態において、粒子が分布する計算対象領域は、３次元空間であってもよく、２次元空間であってもよい。粒子法では、ある粒子に対してあらかじめ設定された影響領域内に存在する近傍粒子のみからの相互作用を計算することで、粒子の運動が解析される。
【００１２】
図１は、このような影響領域の例を示している。粒子１０１に対する影響領域１１１は、粒子１０１を中心とする半径ｈの円又球で表され、粒子１０１の運動は、影響領域１１１内の複数の近傍粒子１０２からの相互作用に基づいて計算される。
【００１３】
粒子法の計算では、計算対象領域を広くとりながら所望の解像度を実現するために、多数の粒子を扱うことが多い。例えば、津波の浸水計算において５ｋｍ四方の領域を１ｍの解像度で扱う場合、数億個の粒子が用いられる。このような大規模な計算を行うために、多数のノードを有する並列計算機が広く用いられている。特許文献１等で用いられている領域分割法では、計算対象領域を複数の部分領域に分割して、それぞれの部分領域を複数のノードに割り当てることで、並列計算が行われる。
【００１４】
図２は、このような領域分割の例を示している。図２の領域分割では、計算対象領域が部分領域２０１〜部分領域２０４に分割され、４個のノードがそれぞれの部分領域を担当し、担当部分領域に含まれる粒子の計算を行う。各部分領域内のメッシュは、図１の影響領域１１１の半径ｈを縦横のサイズとする正方形又は立方体を表す。
【００１５】
各粒子の運動は影響領域内に存在する近傍粒子からの影響を受けるため、部分領域間の境界から距離ｈまでの範囲に存在する粒子（袖粒子）は、隣接する部分領域に含まれる粒子に対して影響を及ぼすことになる。そこで、各ノードは、担当部分領域に含まれる袖粒子を特定し、その袖粒子の情報を隣接する部分領域を担当するノードへ送信する。
【００１６】
図３は、このような袖粒子の例を示している。部分領域２０２を担当するノードは、境界３０１から距離ｈまでの範囲３０２に存在する袖粒子を特定し、その袖粒子の情報を、部分領域２０１を担当するノードへ送信する。袖粒子の情報を受信したノードは、その情報を用いて部分領域２０１内の粒子の運動を解析する。このとき、部分領域２０１を担当するノードは、例えば、連続方程式を解いて密度を更新したり、運動方程式を解いて得られた加速度を用いて、速度及び位置座標を更新したりする。
【００１７】
図４は、部分領域間における粒子の移動の例を示している。並列計算では、粒子シミュレーションの時間ステップ毎に粒子の運動が計算される。ある時間ステップにおいて部分領域２０１内に存在していた粒子４０１が、次の時間ステップにおいて、境界３０１を超えて部分領域２０２へ移動した場合、部分領域２０１を担当するノードは、粒子４０１の情報を部分領域２０２を担当するノードへ送信する。これにより、部分領域２０２を担当するノードが粒子４０１の運動を解析することが可能になる。
【００１８】
ある時間ステップの並列計算においてノード毎に計算が終了するタイミングが異なると、無駄な待ち時間が発生してしまうため、並列計算の計算効率（並列化効率）が低下する原因となる。そこで、特許文献２等のように、複数のノードの計算時間を均等にするための技術がいくつか提案されている。
【００１９】
しかしながら、従来の領域分割法では、粒子の個数を基準として部分領域間の境界を決定しているため、次のような場合には、必ずしも最適な並列化効率が得られるとは限らない。
（ａ）流体を表す粒子、構造物を表す粒子、弾性体を表す粒子等のように、複数種類の粒子が混在する。
（ｂ）粒子の種類毎に計算量（計算コスト）が異なる。
（ｃ）粒子の種類毎に粒子の分布の様子が異なる。
【００２０】
図５は、複数種類の粒子が混在する計算対象領域の例を示している。図５の計算対象領域内には、流体を表す種類Ａの粒子５０１と、構造物を表す種類Ｂの粒子５０２とが混在している。流体を表す粒子５０１は時間ステップ毎に移動するのに対して、構造物を表す粒子５０２の位置は変化しないため、粒子５０１の計算量は粒子５０２の計算量よりも大きくなると考えられる。
【００２１】
この計算対象領域を、粒子数が均等になるように部分領域５１１〜部分領域５１４に分割した場合、部分領域５１４の計算量が部分領域５１１〜部分領域５１３の計算量よりも極めて大きくなる。このため、ノード間のロードバランスが悪くなって、並列化効率が低下する。
【００２２】
図６は、実施形態の粒子シミュレーション装置の機能的構成例を示している。図６の粒子シミュレーション装置６０１は、例えば、並列計算機に対応し、ノード６１１−１〜ノード６１１−Ｎ（Ｎは２以上の整数）のＮ個のノードを有する。粒子シミュレーション装置６０１は、複数種類の粒子が分布する領域（計算対象領域）を分割した部分領域を各ノードに割り当てて、粒子シミュレーションのための複数のプログラムを各ノード毎に実行する。ノード６１１−１は、計算部６２１及び変更部６２２を含み、粒子シミュレーションにおける計算機資源配分処理を行う。
【００２３】
図７は、図６のノード６１１−１が行う計算機資源配分処理の例を示すフローチャートである。まず、計算部６２１は、複数のプログラムそれぞれの処理対象粒子の種類と、複数のプログラムそれぞれの実行時間とに基づいて、各処理対象粒子に対する演算コストを求める（ステップ７０１）。次に、変更部６２２は、演算コストと各処理対象粒子の粒子数を基に、部分領域の領域境界の位置を変更する（ステップ７０２）。
【００２４】
図６の粒子シミュレーション装置６０１によれば、並列計算機による複数種類の粒子に対する粒子シミュレーションにおいて、計算機資源を効率よく配分することができる。
【００２５】
図８は、並列計算機の構成例を示している。図８の並列計算機８０１は、図６の粒子シミュレーション装置６０１に対応し、ノード８１１−１〜ノード８１１−Ｎ、入力装置８１２、及び出力装置８１３を含む。ノード８１１−１〜ノード８１１−Ｎ、入力装置８１２、及び出力装置８１３は、通信ネットワーク８１４により互いに接続されている。
【００２６】
ノード８１１−１〜ノード８１１−Ｎは、図６のノード６１１−１〜ノード６１１−Ｎにそれぞれ対応し、粒子シミュレーションのための複数のプログラムに対応する複数のサブルーチンを実行する。
【００２７】
入力装置８１２は、例えば、キーボード、ポインティングデバイス等であり、オペレータ又はユーザからの指示又は情報の入力に用いられる。出力装置８１３は、例えば、表示装置、プリンタ、スピーカ等であり、オペレータ又はユーザへの問い合わせ又は処理結果の出力に用いられる。処理結果は、粒子シミュレーションのシミュレーション結果であってもよい。
【００２８】
なお、並列計算機８０１が図８のすべての構成要素を含む必要はなく、用途又は条件に応じて一部の構成要素を省略することも可能である。例えば、オペレータ又はユーザからの指示又は情報を入力する必要がない場合は、入力装置８１２を省略してもよい。オペレータ又はユーザへの問い合わせ又は処理結果を出力する必要がない場合は、出力装置８１３を省略してもよい。
【００２９】
図９は、図８のノード８１１−ｉ（ｉ＝１〜Ｎ）に対応する情報処理装置の構成例を示している。図９のノード８１１−ｉは、Central Processing Unit（ＣＰＵ）９０１、メモリ９０２、補助記憶装置９０３、ネットワーク接続装置９０４、及び媒体駆動装置９０５を含む。ＣＰＵ９０１、メモリ９０２、補助記憶装置９０３、ネットワーク接続装置９０４、及び媒体駆動装置９０５は、バス９０６により互いに接続されている。
【００３０】
メモリ９０２は、例えば、Read Only Memory（ＲＯＭ）、Random Access Memory（ＲＡＭ）等の半導体メモリであり、粒子シミュレーションに用いられるプログラム及びデータを記憶する。
【００３１】
ＣＰＵ９０１（プロセッサ）は、例えば、メモリ９０２を利用してプログラムを実行することにより、図６の計算部６２１及び変更部６２２として動作する。ＣＰＵ９０１は、粒子シミュレーションのための複数のサブルーチンも実行する。
【００３２】
補助記憶装置９０３は、例えば、磁気ディスク装置、光ディスク装置、光磁気ディスク装置等である。補助記憶装置９０３は、ハードディスクドライブであってもよい。ノード８１１−ｉは、補助記憶装置９０３にプログラム及びデータを格納しておき、それらをメモリ９０２にロードして使用することができる。
【００３３】
ネットワーク接続装置９０４は、図８の通信ネットワーク８１４に接続され、通信に伴うデータ変換を行う通信インタフェースである。ノード８１１−ｉは、他のノード又は並列計算機８０１の外部に設けられた装置から、通信ネットワーク８１４を介してプログラム及びデータを受信し、それらをメモリ９０２にロードして使用することができる。
【００３４】
媒体駆動装置９０５は、可搬型記録媒体９０７を駆動し、その記録内容にアクセスする。可搬型記録媒体９０７は、メモリデバイス、フレキシブルディスク、光ディスク、光磁気ディスク等である。可搬型記録媒体９０７は、Compact Disk Read Only Memory（ＣＤ−ＲＯＭ）、Digital Versatile Disk（ＤＶＤ）、Universal Serial Bus（ＵＳＢ）メモリ等であってもよい。オペレータ又はユーザは、この可搬型記録媒体９０７にプログラム及びデータを格納しておき、それらをメモリ９０２にロードして使用することができる。
【００３５】
このように、粒子シミュレーションに用いられるプログラム及びデータを格納するコンピュータ読み取り可能な記録媒体は、メモリ９０２、補助記憶装置９０３、又は可搬型記録媒体９０７のような、物理的な（非一時的な）記録媒体である。
【００３６】
なお、ノード８１１−ｉが図９のすべての構成要素を含む必要はなく、用途又は条件に応じて一部の構成要素を省略することも可能である。例えば、通信ネットワーク８１４を介してプログラム及びデータを受信する場合は、補助記憶装置９０３を省略してもよく、可搬型記録媒体９０７を使用しない場合は、媒体駆動装置９０５を省略してもよい。
【００３７】
並列計算機８０１のノード８１１−１〜ノード８１１−Ｎは、粒子シミュレーションの時間ステップ（シミュレーションステップ）毎に複数のサブルーチンを実行し、複数種類の粒子の運動を計算する。このとき、ノード８１１−１は、複数のサブルーチンそれぞれの処理対象粒子の種類を示す情報と、ある時間ステップにおける複数のサブルーチンそれぞれの実行時間とを用いて、粒子の種類毎の重みを表すコスト係数を求める。コスト係数は、粒子１個当たりの計算量を表す指標である。
【００３８】
そして、ノード８１１−１は、それぞれの種類の粒子の計算量に基づいて、計算対象領域をＮ個の部分領域に分割する。ある種類の粒子の計算量は、その種類の粒子の個数をコスト係数により重み付けすることで、求めることができる。次の時間ステップにおいて、ノード８１１−１〜ノード８１１−Ｎは、分割されたＮ個の部分領域それぞれに含まれる粒子を対象として、複数のサブルーチンを実行する。これにより、粒子の分布の変化に応じて、計算対象領域の分割方法を変更することが可能になる。
【００３９】
並列計算機８０１は、粒子法を用いた数値計算において広く利用可能である。例えば、防災計画において河川又は海の水に対する流体解析を行ったり、金属製品の設計において鋳造過程を解析したりすることができる。また、製品設計の際に弾性体解析を行うことで、封止ゲルの形状等を適切に決定することもできる。
【００４０】
図１０は、図８の並列計算機８０１が行う粒子シミュレーションの例を示すフローチャートである。まず、ノード８１１−１のＣＰＵ９０１は、入力装置８１２又は並列計算機８０１の外部の装置から入力される粒子データを取得し、メモリ９０２に格納する（ステップ１００１）。
【００４１】
粒子データは、計算対象領域内の流体、弾性体、構造物等の連続体を表す複数の粒子それぞれの位置、速度等を表すデータである。例えば、数値計算の直接の対象である、位置の変化を伴う移動粒子は、流体、弾性体等をモデル化した粒子である。各移動粒子の粒子データは、中心座標、速度、影響領域の半径、密度、質量、変形勾配テンソル、材料物性、温度等の物理量を含んでいてもよい。計算対象領域内には、移動速度の異なる複数種類の移動粒子が含まれることもある。
【００４２】
一方、位置の変化を伴わず、移動粒子の運動の境界条件を規定するに過ぎない構造物粒子は、流体を格納する容器の壁等をモデル化した粒子である。各構造物粒子の粒子データは、中心座標、面要素の法線ベクトル、面要素の面積等の物理量を含んでいてもよい。
【００４３】
次に、ノード８１１−１のＣＰＵ９０１は、入力装置８１２又は外部の装置から入力されるサブルーチンデータを取得し、メモリ９０２に格納する（ステップ１００２）。サブルーチンデータは、複数のサブルーチンそれぞれの処理対象粒子の種類を表す。
【００４４】
図１１は、サブルーチンデータの例を示している。図１１のサブルーチンデータは、実行順序、サブルーチンＩＤ、及び粒子の種類を含む。実行順序は、各サブルーチンの実行順序を表し、サブルーチンＩＤは、各サブルーチンの識別情報である。粒子の種類は、各サブルーチンの処理対象粒子の種類を表す。各サブルーチンは、処理対象粒子の種類に応じて異なる方程式をモデル化することで実装される。例えば、１番目のサブルーチンＲ１を実行すると、種類Ａ及び種類Ｂの２種類の粒子を対象とする計算が行われる。
【００４５】
次に、ノード８１１−１〜ノード８１１−Ｎは、それぞれが担当する部分領域内の移動粒子に対して設定された物理モデルに従って、１時間ステップに対する移動粒子の運動を計算する（ステップ１００３）。
【００４６】
粒子シミュレーションの最初の時間ステップにおいては、まだコスト係数が求められていないため、ノード８１１−１のＣＰＵ９０１は、所定の分割方法に従って、計算対象領域をＮ個の部分領域に分割する。そして、ノード８１１−１は、ノード８１１−２〜ノード８１１−Ｎに対して、サブルーチンデータとそれぞれの担当部分領域の粒子データとを送信する。
【００４７】
各ノードのＣＰＵ９０１は、各担当部分領域の粒子データを用いて、サブルーチンデータが示す複数のサブルーチンを実行順序に従って実行し、各サブルーチンの実行時間を計測してメモリ９０２に格納する。このとき、ノード８１１−２〜ノード８１１−Ｎは、各サブルーチンの実行時間をノード８１１−１へ送信する。そして、各ノードのＣＰＵ９０１は、計算した加速度を含む時間微分項を用いて、各粒子の物理量を更新する（ステップ１００４）。
【００４８】
次に、ノード８１１−１のＣＰＵ９０１は、領域境界の更新を行うか否かを判定する（ステップ１００５）。領域境界の更新は、事前に設定した所定数の時間ステップ毎に１回行ってもよい。あるいは、ノード８１１−１〜ノード８１１−Ｎにおける各サブルーチンの実行時間のばらつきが所定値以上だった場合に、領域境界の更新を行う。
【００４９】
領域境界の更新を行う場合（ステップ１００５，ＹＥＳ）、ノード８１１−１のＣＰＵ９０１は、粒子の種類毎のコスト係数を求める（ステップ１００６）。コスト係数は、ノード８１１−１〜ノード８１１−Ｎにおける各サブルーチンの実行時間と、サブルーチンデータが示す各サブルーチンの粒子の種類とに基づいて求められる。
【００５０】
次に、ノード８１１−１のＣＰＵ９０１は、それぞれの種類の粒子の計算量に基づいて、部分領域間の境界の位置を決定する（ステップ１００７）。ある種類の粒子の計算量は、その種類の粒子の個数にコスト係数を乗算することで、求めることができる。例えば、ノード８１１−１〜ノード８１１−Ｎが有する計算機資源が同程度である場合、それぞれの部分領域内における計算量が均等になるように、境界の位置が決定される。各部分領域内における計算量は、例えば、種類毎の粒子の個数とコスト係数とを乗算した値の総和で表される。
【００５１】
次に、各ノードは、領域境界の更新により担当ノードが変更された領域内の各粒子について、ステップ１００４で更新した物理量を表す粒子データを、更新後の担当ノードへ送信する（ステップ１００８）。
【００５２】
次に、ノード８１１−２〜ノード８１１−Ｎは、担当するすべての粒子の粒子データをノード８１１−１へ送信する（ステップ１００９）。そして、ノード８１１−１は、ノード８１１−１〜ノード８１１−Ｎの粒子データを含む計算結果ファイルを生成して、補助記憶装置９０３に格納する。出力装置８１３は、計算結果ファイルに含まれる粒子データを、数値データ、グラフ、画像等の形式で出力することができる。
【００５３】
一方、ノード８１１−１〜ノード８１１−Ｎにおける各サブルーチンの実行時間のばらつきが所定値以内の場合は、領域境界の更新を行わない。領域境界の更新を行わない場合（ステップ１００５，ＮＯ）、並列計算機８０１は、ステップ１００９以降の処理を行う。
【００５４】
次に、ノード８１１−１のＣＰＵ９０１は、粒子シミュレーションを終了するか否かを判定する（ステップ１０１０）。並列計算機８０１は、時間ステップの数が上限値に達したときに粒子シミュレーションを終了してもよく、粒子の物理量の変化が閾値未満になったときに粒子シミュレーションを終了してもよい。
【００５５】
粒子シミュレーションを終了しない場合（ステップ１０１０，ＮＯ）、並列計算機８０１は、次の時間ステップについて、ステップ１００３以降の処理を繰り返す。次の時間ステップの計算では、更新後の境界で分割されるＮ個の部分領域が用いられる。粒子シミュレーションを終了する場合（ステップ１０１０，ＹＥＳ）、並列計算機８０１は、処理を終了する。
【００５６】
図１０の粒子シミュレーションによれば、計算量の異なる複数種類の粒子が混在する計算対象領域において適切な領域境界を設定することができ、並列計算機の計算機資源を複数の粒子に対して効率よく配分することができる。これにより、並列化効率が向上するため、各ノードに搭載する計算機資源を削減することが可能になる。
【００５７】
例えば、津波シミュレーションにおいて、海水を表す種類Ａの流体粒子と、地形及び建物を表す種類Ｂの構造物粒子とが、計算対象領域内に分布しており、流体粒子の個数がＮＡ個であり、構造物粒子の個数がＮＢ個であるものとする。また、図１１のサブルーチンデータが入力されたものとする。
【００５８】
この場合、サブルーチンＲ１は、例えば、ある粒子の影響領域内に存在する近傍粒子を探索するためのデータ構造を生成する前処理のサブルーチンであり、種類Ａ及び種類Ｂの粒子を対象とする。サブルーチンＲ２は、例えば、流体粒子の連続方程式を解いて密度を計算するサブルーチンであり、種類Ａの粒子を対象とする。
【００５９】
サブルーチンＲ３は、例えば、構造物粒子から流体粒子への反発力を計算するサブルーチンであり、種類Ｂの粒子を対象とする。サブルーチンＲ４は、例えば、流体粒子の運動方程式を解いて加速度を計算するサブルーチンであり、種類Ａの粒子を対象とする。流体粒子の運動方程式として、ナヴィエ・ストークス方程式を用いてもよい。
【００６０】
図１２は、ある時間ステップにおけるこれらのサブルーチンの実行時間の例を示している。サブルーチンＲ１〜サブルーチンＲ４の実行時間は、それぞれ、Ｔ１〜Ｔ４である。実行時間Ｔｊ（ｊ＝１〜４）は、ノード８１１−１〜ノード８１１−ＮにおけるサブルーチンＲｊの実行時間の総和を表す。この場合、ある時間ステップにおけるサブルーチンＲ１〜サブルーチンＲ４の総実行時間Ｔ０は、次式により表すことができる。
Ｔ０＝Ｔ１＋Ｔ２＋Ｔ３＋Ｔ４      （１）
【００６１】
実行時間Ｔｊ（ｊ＝１〜４）は、サブルーチンＲｊの計算量に比例し、総実行時間Ｔ０は、１時間ステップにおけるサブルーチンＲ１〜サブルーチンＲ４の総計算量に比例すると考えられる。ここで、各サブルーチンの実行時間が処理対象粒子の個数に比例すると仮定すると、コスト係数は、粒子の種類毎の実行時間と個数との比例係数に対応する。このとき、総実行時間Ｔ０は、種類Ａの粒子のコスト係数ＣＡと、種類Ｂの粒子のコスト係数ＣＢとを用いて、次式により表すことができる。
Ｔ０＝ＣＡ＊ＮＡ＋ＣＢ＊ＮＢ      （２）
【００６２】
また、各サブルーチンの実行時間Ｔ１〜実行時間Ｔ４は、次式により表すことができる。
Ｔ１＝Ｇ１＊（ＮＡ＋ＮＢ）      （３）
Ｔ２＝Ｇ２＊ＮＡ                （４）
Ｔ３＝Ｇ３＊ＮＢ                （５）
Ｔ４＝Ｇ４＊ＮＡ                （６）
【００６３】
Ｇ１〜Ｇ４は、各サブルーチンの比例係数を表す。式（１）〜式（６）より、次式が得られる。
Ｔ０＝ＣＡ＊ＮＡ＋ＣＢ＊ＮＢ
＝Ｔ１＋Ｔ２＋Ｔ３＋Ｔ４
＝Ｇ１＊（ＮＡ＋ＮＢ）＋Ｇ２＊ＮＡ＋Ｇ３＊ＮＢ＋Ｇ４＊ＮＡ
＝（Ｇ１＋Ｇ２＋Ｇ４）＊ＮＡ＋（Ｇ１＋Ｇ３）＊ＮＢ          （７）
【００６４】
式（７）より、次式が得られる。
ＣＡ＝Ｇ１＋Ｇ２＋Ｇ４      （８）
ＣＢ＝Ｇ１＋Ｇ３            （９）
【００６５】
式（３）〜式（６）を変形すると、次式が得られる。
Ｇ１＝Ｔ１／（ＮＡ＋ＮＢ）      （１０）
Ｇ２＝Ｔ２／ＮＡ                （１１）
Ｇ３＝Ｔ３／ＮＢ                （１２）
Ｇ４＝Ｔ４／ＮＡ                （１３）
【００６６】
式（８）及び式（９）のＧ１〜Ｇ４に対して、式（１０）〜式（１３）を適用することで、次式が得られる。
ＣＡ＝Ｔ１／（ＮＡ＋ＮＢ）＋Ｔ２／ＮＡ＋Ｔ４／ＮＡ      （１４）
ＣＢ＝Ｔ１／（ＮＡ＋ＮＢ）＋Ｔ３／ＮＢ                  （１５）
【００６７】
ノード８１１−１のＣＰＵ９０１は、式（１４）及び式（１５）を用いて、各サブルーチンの実行時間Ｔ１〜実行時間Ｔ４から、コスト係数ＣＡ及びコスト係数ＣＢを求めることができる。サブルーチンの実行時間の測定値に基づいてコスト係数を決定することで、粒子シミュレーションの実行環境に合わせてコスト係数を調整することが可能になる。
【００６８】
図１３は、津波シミュレーションの初期における粒子の分布の例を示している。津波シミュレーションの初期においては、陸上に海水が存在しないため、流体粒子は計算対象領域の一部の領域１３０１に局在するのに対して、地形及び建物を表す構造物粒子は計算対象領域全体にほぼ一様に分布する。
【００６９】
ここで、ＮＡ＝１００００、ＮＢ＝４０００００とすると、計算対象領域の１／４に相当する領域１３０１に種類Ａの流体粒子が１万個存在し、計算対象領域全体に種類Ｂの構造物粒子が４０万個存在する。Ｔ１＝１秒、Ｔ２＝２秒、Ｔ３＝０．５秒、及びＴ４＝２秒とすると、式（１４）及び式（１５）より、ＣＡ＝４．０２ｅ−４及びＣＢ＝３．６９ｅ−６となる。この場合、流体粒子のコスト係数ＣＡは、構造物粒子のコスト係数ＣＢに対して百倍程度大きくなる。
【００７０】
図１４は、図１３の計算対象領域をコスト係数に基づいて分割した部分領域の例を示している。この例では、並列計算機８０１のノードの個数が４個であり、境界１４１１及び境界１４１２によって、計算対象領域が部分領域１４０１〜部分領域１４０４に分割されている。それぞれの部分領域内において、流体粒子の個数とコスト係数ＣＡとを乗算した値と、構造物粒子の個数とコスト係数ＣＢとを乗算した値との和が均等になるように、境界１４１１及び境界１４１２の位置が決定される。
【００７１】
この場合、ノード８１１−１〜ノード８１１−４は、それぞれ、部分領域１４０１〜部分領域１４０４の計算を担当し、ノード８１１−１〜ノード８１１−４の計算量が均等化される。このため、ノード間のロードバランスは良好であり、並列化効率が向上する。
【００７２】
一方、図１５は、図１３の計算対象領域を粒子の個数のみに基づいて分割した部分領域の例を示している。この例では、境界１５１１及び境界１５１２によって、計算対象領域が部分領域１５０１〜部分領域１５０４に分割されている。それぞれの部分領域内において、流体粒子の個数と構造物粒子の個数との和が均等になるように、境界１５１１及び境界１５１２の位置が決定される。
【００７３】
この場合、ノード８１１−１〜ノード８１１−４は、それぞれ、部分領域１５０１〜部分領域１５０４の計算を担当し、ノード８１１−４の計算量がノード８１１−１〜ノード８１１−３の計算量よりも極めて大きくなる。このため、ノード間のロードバランスが悪く、並列化効率が低下する。
【００７４】
図１６は、ノード８１１−１〜ノード８１１−４におけるサブルーチンの実行時間の例を示している。プロセスＰ１〜プロセスＰ４は、それぞれ、ノード８１１−１〜ノード８１１−４内で動作し、ある時間ステップにおいてサブルーチンＲ１〜サブルーチンＲ４を実行する。
【００７５】
このとき、各サブルーチンの実行時間は、プロセス毎に異なることがある。例えば、プロセスＰ１によるサブルーチンＲ１の実行時間は、プロセスＰ２〜プロセスＰ４によるサブルーチンＲ１の実行時間よりも長い。また、プロセスＰ１によるサブルーチンＲ２の実行時間は、プロセスＰ２〜プロセスＰ４によるサブルーチンＲ２の実行時間よりも短い。
【００７６】
そこで、ノード８１１−１のＣＰＵ９０１は、ノード８１１−１〜ノード８１１−４におけるサブルーチンＲｊ（ｊ＝１〜４）の実行時間の総和を、式（１）の実行時間Ｔｊとして用いる。これにより、ノード間における計算量の偏りが均等化されるように、コスト係数ＣＡ及びコスト係数ＣＢを決定することができる。
【００７７】
図１７は、図１４の計算対象領域に対する津波シミュレーションが進行した後の部分領域の例を示している。計算が進行して津波による陸上への浸水が進むと、流体粒子は広い領域１７２１に分布するようになる。
【００７８】
この場合、コスト係数に基づいて決定された境界１７１１及び境界１７１２によって、計算対象領域が部分領域１７０１〜部分領域１７０４に分割される。ノード８１１−１〜ノード８１１−４は、それぞれ、部分領域１７０１〜部分領域１７０４の計算を担当する。このように、粒子の分布が変化すると、それぞれの部分領域の計算量が均等になるように、領域境界の位置も変化する。
【００７９】
計算量の異なる３種類以上の粒子が混在する場合でも、図１０の粒子シミュレーションによって、それぞれの部分領域内における計算量が均等になるように領域境界の位置を決定することが可能である。
【００８０】
例えば、計算対象領域内にＬ種類の粒子が分布しており、種類ｋ（ｋ＝１〜Ｌ）の粒子の個数がＮｋ個であり、ノード８１１−１〜ノード８１１−ＮがサブルーチンＲ１〜サブルーチンＲＭを実行する場合について考察する。この場合、ある時間ステップにおけるサブルーチンＲ１〜サブルーチンＲＭの総実行時間Ｔ０は、次式により表すことができる。
Ｔ０＝Ｔ１＋Ｔ２＋・・・＋ＴＭ      （２１）
【００８１】
式（２１）の実行時間Ｔｊ（ｊ＝１〜Ｍ）は、サブルーチンＲｊの実行時間を表す。このとき、総実行時間Ｔ０は、種類ｋ（ｋ＝１〜Ｌ）の粒子のコスト係数Ｃｋを用いて、次式により表すことができる。
Ｔ０＝Ｃ１＊Ｎ１＋Ｃ２＊Ｎ２＋・・・＋ＣＬ＊ＮＬ      （２２）
【００８２】
また、各サブルーチンの実行時間Ｔｊは、次式により表すことができる。
Ｔｊ＝Ｇｊ＊（Ｂｊ１＊Ｎ１＋Ｂｊ２＊Ｎ２＋・・・＋ＢｊＬ＊ＮＬ）      （２３）
【００８３】
式（２３）のＢｊｋは、０又は１のいずれか一方の値を表し、種類ｋの粒子がサブルーチンＲｊの処理対象粒子である場合、Ｂｊｋ＝１となり、種類ｋの粒子がサブルーチンＲｊの処理対象粒子ではない場合、Ｂｊｋ＝０となる。Ｇｊは、サブルーチンＲｊの比例係数を表す。式（２１）〜式（２３）より、次式が得られる。
Ｔ０＝Ｃ１＊Ｎ１＋Ｃ２＊Ｎ２＋・・・＋ＣＬ＊ＮＬ
＝Ｔ１＋Ｔ２＋・・・＋ＴＭ
＝Ｇ１＊（Ｂ１１＊Ｎ１＋Ｂ１２＊Ｎ２＋・・・＋Ｂ１Ｌ＊ＮＬ）
＋Ｇ２＊（Ｂ２１＊Ｎ１＋Ｂ２２＊Ｎ２＋・・・＋Ｂ２Ｌ＊ＮＬ）
＋・・・＋ＧＭ（ＢＭ１＊Ｎ１＋ＢＭ２＊Ｎ２＋・・・＋ＢＭＬ＊ＮＬ）
＝（Ｇ１＊Ｂ１１＋Ｇ２＊Ｂ２１＋・・・＋ＧＭ＊ＢＭ１）＊Ｎ１
＋（Ｇ１＊Ｂ１２＋Ｇ２＊Ｂ２２＋・・・＋ＧＭ＊ＢＭ２）＊Ｎ２
＋・・・＋（Ｇ１＊Ｂ１Ｌ＋Ｇ２＊Ｂ２Ｌ＋・・・＋ＧＭ＊ＢＭＬ）＊ＮＬ
（２４）
【００８４】
式（２４）より、次式が得られる。
Ｃｋ＝Ｇ１＊Ｂ１ｋ＋Ｇ２＊Ｂ２ｋ＋・・・＋ＧＭ＊ＢＭｋ      （２５）
【００８５】
式（２３）を変形すると、次式が得られる。
Ｇｊ＝Ｔｊ／（Ｂｊ１＊Ｎ１＋Ｂｊ２＊Ｎ２＋・・・＋ＢｊＬ＊ＮＬ）      （２６）
【００８６】
式（２５）のＧ１〜ＧＭに対して、式（２６）を適用することで、次式が得られる。
Ｃｋ＝Ｔ１＊Ｂ１ｋ／（Ｂ１１＊Ｎ１＋Ｂ１２＊Ｎ２＋・・・＋Ｂ１Ｌ＊ＮＬ）
＋Ｔ２＊Ｂ２ｋ／（Ｂ２１＊Ｎ１＋Ｂ２２＊Ｎ２＋・・・＋Ｂ２Ｌ＊ＮＬ）
＋・・・
＋ＴＭ＊ＢＭｋ／（ＢＭ１＊Ｎ１＋ＢＭ２＊Ｎ２＋・・・＋ＢＭＬ＊ＮＬ）
（２７）
【００８７】
ノード８１１−１のＣＰＵ９０１は、式（２７）を用いて、各サブルーチンの実行時間Ｔ１〜実行時間ＴＭから、コスト係数Ｃ１〜コスト係数ＣＬを求めることができる。
【００８８】
鋳造過程のシミュレーションでは、例えば、溶融金属を表す種類Ａの流体粒子と、金型を表す種類Ｂの構造物粒子とが分布する計算対象領域を用いることができる。この場合、例えば、津波シミュレーションで用いたサブルーチンＲ１〜サブルーチンＲ４に続いて、サブルーチンＲ５が実行される。サブルーチンＲ５は、例えば、エネルギー方程式を解いて熱伝導による温度変化を計算するサブルーチンであり、種類Ａ及び種類Ｂの粒子を対象とする。
【００８９】
弾性体のシミュレーションでは、例えば、弾性体を表す移動粒子を対象として変形勾配テンソルを計算するサブルーチンと、変形勾配テンソルからコーシー応力を求めるサブルーチンとを用いることができる。
【００９０】
図６の粒子シミュレーション装置６０１、図８の並列計算機８０１、及び図９のノード８１１−ｉの構成は一例に過ぎず、粒子シミュレーション装置６０１及び並列計算機８０１の用途又は条件に応じて一部の構成要素を省略又は変更してもよい。例えば、図６の粒子シミュレーション装置６０１において、ノード６１１−１は複数のプログラムを実行せず、ノード６１１−２〜ノード６１１−Ｎが複数のプログラムを実行してもよい。同様に、図８の並列計算機８０１において、ノード８１１−１は複数のサブルーチンを実行せず、ノード８１１−２〜ノード８１１−Ｎが複数のサブルーチンを実行してもよい。
【００９１】
図７及び図１０のフローチャートは一例に過ぎず、粒子シミュレーション装置６０１及び並列計算機８０１の構成又は条件に応じて一部の処理を省略又は変更してもよい。例えば、図１０の粒子シミュレーションにおいて、１時間ステップ毎に領域境界の更新を行う場合は、ステップ１００５の処理を省略することができる。
【００９２】
ノード８１１−１〜ノード８１１−Ｎが有する計算機資源に差がある場合、図１０のステップ１００７において、ノード８１１−１は、それぞれの部分領域内における計算量がそれぞれのノードの計算機資源に比例するように、境界の位置を決定してもよい。
【００９３】
図１１のサブルーチンデータは一例に過ぎず、別のサブルーチンデータを用いてもよい。例えば、サブルーチンの個数は、３個以下であってもよく、５個以上であってもよい。サブルーチンの処理対象粒子の種類は、３種類以上であってもよい。図１２及び図１６の実行時間は一例に過ぎず、実行時間は使用するサブルーチンによって変化し、時間ステップ毎にも変化する。
【００９４】
図１〜図５、図１３〜図１５、及び図１７の粒子、計算対象領域、及び部分領域は一例に過ぎず、別のモデルを用いて粒子シミュレーションを行ってもよい。津波又は鋳造過程以外の流体シミュレーション、弾性体シミュレーション、又は他の連続体を含む粒子シミュレーションを行うことも可能である。
【００９５】
開示の実施形態とその利点について詳しく説明したが、当業者は、特許請求の範囲に明確に記載した本発明の範囲から逸脱することなく、様々な変更、追加、省略をすることができるであろう。
【００９６】
図１乃至図１７を参照しながら説明した実施形態に関し、さらに以下の付記を開示する。
（付記１）
複数のノードを有する並列計算機のための複数種類の粒子が分布する領域に対する粒子シミュレーションプログラムであって、
前記領域を分割した部分領域を各ノードに割り当てて、粒子シミュレーションのための複数のプログラムを各ノード毎に実行し、
前記複数のプログラムそれぞれの処理対象粒子の種類と、前記複数のプログラムそれぞれの実行時間とに基づいて前記各処理対象粒子に対する演算コストを求め、
前記演算コストと前記各処理対象粒子の粒子数を基に前記部分領域の領域境界の位置を変更する、
処理を前記並列計算機に実行させる粒子シミュレーションプログラム。
（付記２）
前記並列計算機は、変更された前記領域境界によって決定された前記部分領域を前記各ノードに割り当てて、前記粒子シミュレーションのための前記複数のプログラムを前記各ノード毎に実行することを特徴とする付記１記載の粒子シミュレーションプログラム。
（付記３）
前記並列計算機は、各ノードで実行された前記複数のプログラムそれぞれの前記実行時間のばらつきが所定値以内の場合には、前記領域境界の位置の変更は行わないことを特徴とする付記１又は２記載の粒子シミュレーションプログラム。
（付記４）
複数のノードを有する粒子シミュレーション装置であって、
前記粒子シミュレーション装置は、複数種類の粒子が分布する領域を分割した部分領域を各ノードに割り当てて、粒子シミュレーションのための複数のプログラムを各ノード毎に実行し、
前記粒子シミュレーション装置は、
前記複数のプログラムそれぞれの処理対象粒子の種類と、前記複数のプログラムそれぞれの実行時間とに基づいて前記各処理対象粒子に対する演算コストを求める計算部と、
前記演算コストと前記各処理対象粒子の粒子数を基に前記部分領域の領域境界の位置を変更する変更部と、
を備えることを特徴とする粒子シミュレーション装置。
（付記５）
前記粒子シミュレーション装置は、変更された前記領域境界によって決定された前記部分領域を前記各ノードに割り当てて、前記粒子シミュレーションのための前記複数のプログラムを前記各ノード毎に実行することを特徴とする付記４記載の粒子シミュレーション装置。
（付記６）
前記変更部は、各ノードで実行された前記複数のプログラムそれぞれの前記実行時間のばらつきが所定値以内の場合には、前記領域境界の位置の変更は行わないことを特徴とする付記４又は５記載の粒子シミュレーション装置。
（付記７）
複数のノードを有する並列計算機が、
複数種類の粒子が分布する領域を分割した部分領域を各ノードに割り当てて、粒子シミュレーションのための複数のプログラムを各ノード毎に実行し、
前記複数のプログラムそれぞれの処理対象粒子の種類と、前記複数のプログラムそれぞれの実行時間とに基づいて前記各処理対象粒子に対する演算コストを求め、
前記演算コストと前記各処理対象粒子の粒子数を基に前記部分領域の領域境界の位置を変更する、
ことを特徴とする粒子シミュレーションにおける計算機資源配分方法。
（付記８）
前記並列計算機は、変更された前記領域境界によって決定された前記部分領域を前記各ノードに割り当てて、前記粒子シミュレーションのための前記複数のプログラムを前記各ノード毎に実行することを特徴とする付記７記載の計算機資源配分方法。
（付記９）
前記並列計算機は、各ノードで実行された前記複数のプログラムそれぞれの前記実行時間のばらつきが所定値以内の場合には、前記領域境界の位置の変更は行わないことを特徴とする付記７又は８記載の計算機資源配分方法。
【符号の説明】
【００９７】
１０１、４０１、５０１、５０２  粒子
１０２  近傍粒子
１１１  影響領域
２０１〜２０４、５１１〜５１４、１４０１〜１４０４、１５０１〜１５０４、１７０１〜１７０４  部分領域
３０１、１４１１、１４１２、１５１１、１５１２、１７１１、１７１２  境界
３０２  範囲
６０１  粒子シミュレーション装置
６１１−１〜６１１−Ｎ、８１１−１〜８１１−Ｎ  ノード
６２１  計算部
６２２  変更部
８０１  並列計算機
８１２  入力装置
８１３  出力装置
８１４  通信ネットワーク
９０１  ＣＰＵ
９０２  メモリ
９０３  補助記憶装置
９０４  ネットワーク接続装置
９０５  媒体駆動装置
９０６  バス
９０７  可搬型記録媒体
１３０１、１７２１  領域
【図１】
0006543557.tif 000002
【図２】
0006543557.tif 000003
【図３】
0006543557.tif 000004
【図４】
0006543557.tif 000005
【図５】
0006543557.tif 000006
【図６】
0006543557.tif 000007
【図７】
0006543557.tif 000008
【図８】
0006543557.tif 000009
【図９】
0006543557.tif 000010
【図１０】
0006543557.tif 000011
【図１１】
0006543557.tif 000012
【図１２】
0006543557.tif 000013
【図１３】
0006543557.tif 000014
【図１４】
0006543557.tif 000015
【図１５】
0006543557.tif 000016
【図１６】
0006543557.tif 000017
【図１７】
0006543557.tif 000018