粒子状態計算装置及び粒子状態計算方法
【課題】 効率的に分散メモリ型並列計算により粒子状態を計算する。
【解決手段】 複数のノード１０は、粒子状態計算システム１を構成し、複数に分割された計算領域である自領域に配置される粒子の状態を時間ステップ毎に計算する装置である。ノード１０は、自領域の粒子の状態を計算する粒子状態計算部１１と、時間ステップが所定の周期の時間ステップか否かを判断する周期判断部１２と、所定の周期の時間ステップに別のノードで用いられる自領域仮想粒子を示す情報を出力する自仮想粒子情報出力部１３と、自領域において計算に用いられる他領域仮想粒子情報を入力する他仮想粒子情報入力部１４と、自領域仮想粒子の状態を出力する自仮想粒子状態出力部１５と、他領域仮想粒子の状態を入力する他仮想粒子状態入力部１６とを備える。粒子状態計算部１１は他領域仮想粒子の情報を用いて粒子の状態を計算する。
【解決手段】 複数のノード１０は、粒子状態計算システム１を構成し、複数に分割された計算領域である自領域に配置される粒子の状態を時間ステップ毎に計算する装置である。ノード１０は、自領域の粒子の状態を計算する粒子状態計算部１１と、時間ステップが所定の周期の時間ステップか否かを判断する周期判断部１２と、所定の周期の時間ステップに別のノードで用いられる自領域仮想粒子を示す情報を出力する自仮想粒子情報出力部１３と、自領域において計算に用いられる他領域仮想粒子情報を入力する他仮想粒子情報入力部１４と、自領域仮想粒子の状態を出力する自仮想粒子状態出力部１５と、他領域仮想粒子の状態を入力する他仮想粒子状態入力部１６とを備える。粒子状態計算部１１は他領域仮想粒子の情報を用いて粒子の状態を計算する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、分散メモリ型並列計算により粒子の状態を計算する粒子状態計算装置及び粒子状態計算方法に関する。
【背景技術】
【０００２】
雲の発達初期において雲粒子（水粒子）の粒子系分布が急激に広がり（粒子系分布の分散が急激に大きくなり）、速やかな降雨開始をもたらす現象が長い間議論されている。その原因として、雲乱流による雲粒子の衝突促進効果、乱流エントレインメント、巨大凝結核、雲粒子の乱流拡散等が挙げられている。中でも、初めに挙げた雲乱流による衝突促進効果が最も長い間議論されてきた。それに伴い、乱流中での粒子の衝突頻度に関する研究が盛んに行われてきた。
【０００３】
これまでに粒子の乱流衝突頻度を予測するモデルがいくつか提案されてきた。例えば、非特許文献１に記載されたような、直接数値計算（Direct Numerical Simulation，ＤＮＳ）によって求められた衝突頻度因子データを基にしたものが提案されている。
【先行技術文献】
【非特許文献】
【０００４】
【非特許文献１】Woittiez, E.J.P., Jonker, H.J.J. and Portela, L.M., “On the Combined Effects of Turbulence and Gravity onDroplet Collisions in Clouds: a Numerical Study”, Journal of the AtmosphericSciences, Vol. 66 (2009), pp.1926-1943.
【発明の概要】
【発明が解決しようとする課題】
【０００５】
しかしながら、上記の従来の方法で得られる予測結果は、以下のように高レイノルズ数での信頼性が確かめられていない。従来の研究で得られたデータは、非特許文献１に示されるように最大でも８５というテイラーマイクロスケール基準レイノルズ数Ｒｅλ（＝ｕ´ｌλ／ν、ここでｕ´は速度変動強度、ｌλはテイラーマイクロスケール、νは動粘度係数である）で得られたデータである。そのレイノルズ数は、実際の雲乱流のＲｅλ＝１０３〜４という高レイノルズ数とはほど遠い。そのため、従来の乱流衝突モデルをそのまま雲乱流のシミュレーションに適用することの妥当性には疑問が残る。
【０００６】
これを解決するためには、できるだけ高いレイノルズ数における衝突頻度因子データが必要である。そして、そのためには大規模な並列計算が必須となる。並列計算には、各プロセスが同じ物理メモリを共有する共有メモリ型並列計算（ｏｐｅｎＭＰや自動並列化ライブラリなど）と、各プロセスがそれぞれ別個の物理メモリを使用する分散メモリ型並列計算（Message Passing Interface（ＭＰＩ）ライブラリが一般的）の２種類に分けられる。
【０００７】
前者は、せいぜい数十プロセスを使った並列計算しか実行できず、また使用できるメモリも限られる。従って、大規模な並列計算を実行できないという短所がある。これまでに、乱流中での粒子衝突に対する共有メモリ型（ｏｐｅｎＭＰ等）の並列コードは開発されてきた。しかし、上記の理由により、共有メモリ型の並列計算では、高レイノルズ数における粒子衝突データを取得することは難しい。
【０００８】
一方、後者の分散メモリ型並列計算では、ハードウェアを寄せ集めれば、多くのプロセスを使用した並列計算が可能である。分散メモリ型並列計算では、それぞれのプロセスが（分散した）メモリを使用して、必要な時にデータ通信を行いながら計算を進める。データ通信速度はメモリアクセス速度に比べて極端に遅いため、効率的なデータ通信が求められる。
【０００９】
しかし、粒子衝突に対する分散メモリ型並列計算コードは開発されてこなかった。その一つの原因として、気相乱流計算に対して主にスペクトル法を用いたＤＮＳが行われてきたことが挙げられる。分散メモリ型並列環境でスペクトル計算する場合には、全領域データの（全対全の）通信が発生する。膨大なデータ通信を効率的に処理する技術が必要となるため、現状では３次元領域分割に対応したスペクトル法は開発されていない。２次元領域分割による大規模並列計算では、一つ一つの領域は細長い、いわゆるペンシル型の形状になり、（表面積）／（体積）が大きくなる。粒子計算のような近傍プロセス間で情報の受け渡しが必要となる並列計算の場合、通信量は領域の表面積に比例するため、（表面積）／（体積）の増大は（通信量）／（計算量）の増大につながり、結果として並列効率の低下を引き起こす。そのため２次元分割では、粒子運動に対する大規模並列計算に向かない。これらのことから、２次元分割にしか対応していないスペクトル法を用いて乱流を計算し、高レイノルズ数における粒子衝突データを取得するには大きな困難があると考えられる。
【００１０】
本発明は、上記を鑑みてなされたものであり、効率的に分散メモリ型並列計算により粒子状態を計算し、その結果より大規模なモデルに対する計算を可能とする粒子状態計算装置及び粒子状態計算方法を提供することを目的とする。
【課題を解決するための手段】
【００１１】
上記目的を達成するために、本発明に係る粒子状態計算装置は、複数の分割された計算領域の一つである自領域に配置される粒子の状態を時間ステップ毎に繰り返し計算する粒子状態計算装置であって、各時間ステップにおける自領域の粒子の状態を計算する粒子状態計算手段と、粒子状態計算手段における対象の時間ステップが、予め設定された周期の時間ステップか否かを判断する周期判断手段と、周期判断手段によって予め設定された周期の時間ステップであると判断された時間ステップに、粒子状態計算手段によって計算された粒子の位置に基づいて、当該粒子から、自領域に隣接する分割された計算領域である隣接領域において計算に用いられる自領域仮想粒子を特定して、特定した自領域仮想粒子を示す情報を当該隣接領域の粒子の状態を計算する別の粒子状態計算装置に出力する自仮想粒子情報出力手段と、別の粒子状態計算装置から、自領域において計算に用いられる他領域仮想粒子を特定する情報を入力する他仮想粒子情報入力手段と、各時間ステップにおいて粒子状態計算手段によって算出された自領域仮想粒子の状態を示す情報を別の粒子状態計算装置に出力する自仮想粒子状態出力手段と、各時間ステップにおいて別の粒子状態計算装置から他領域仮想粒子の状態を示す情報を入力する他仮想粒子状態入力手段と、粒子状態計算手段による計算に応じた情報を出力する出力手段と、を備え、粒子状態計算手段は、他仮想粒子情報入力手段によって入力される他領域仮想粒子を特定する情報、及び他仮想粒子状態入力手段によって入力された他領域仮想粒子の状態を示す情報を用いて粒子の状態を計算することを特徴とする。
【００１２】
本発明に係る粒子状態計算装置では、予め設定された周期の時間ステップにおいて自領域における自領域仮想粒子が特定されて、特定された自領域仮想粒子を示す情報が隣接領域の状態を計算する別の粒子状態計算装置に出力される。その一方で、別の粒子状態計算装置から、自領域において計算に用いられる他領域仮想粒子を特定する情報が入力される。そして、各時間ステップにおいて、自領域仮想粒子の状態を示す情報が出力されると共に、他領域仮想粒子の状態を示す情報が入力される。入力された他領域仮想粒子を特定する情報、及び当該粒子の状態を示す情報を用いて自領域の粒子の状態が計算される。即ち、粒子状態計算装置間でやり取りする情報に係る粒子を限定された仮想粒子とすることで、各粒子状態計算装置で必要な情報を確実に入力させると共に粒子状態計算装置間の通信を効率的にする。また、仮想粒子を特定する情報のやりとりを、予め設定された周期の時間ステップとすることで粒子状態計算装置間の通信が更に効率的になる。これによって、本発明に係る粒子状態計算装置によれば、効率的に分散メモリ型並列計算により粒子状態を計算し、その結果より大規模なモデルに対する計算が可能となる。
【００１３】
粒子は、流体中の粒子であり、粒子状態計算手段は、有限差分法によって流体の状態を算出して、当該流体の影響を考慮して粒子の状態を計算する、ことが望ましい。この構成によれば、流体中の粒子に係る大規模なモデル計算が可能となる。
【００１４】
粒子状態計算手段は、他仮想粒子情報入力手段によって入力される他領域仮想粒子を特定する情報に基づいて、粒子間の近傍関係を特定して、当該近傍関係に基づいて粒子の状態を計算することが望ましい。この構成によれば、粒子間の衝突等の粒子間の相互作用が生じるモデルの計算を確実に行うことができる。
【００１５】
ところで、本発明は、上記のように粒子状態計算装置の発明として記述できる他に、以下のように粒子状態計算方法の発明としても記述することができる。これはカテゴリが異なるだけで、実質的に同一の発明であり、同様の作用及び効果を奏する。
【００１６】
即ち、本発明に係る粒子状態計算方法は、複数の分割された計算領域の一つである自領域に配置される粒子の状態を時間ステップ毎に繰り返し計算する粒子状態計算装置による粒子状態計算方法であって、各時間ステップにおける自領域の粒子の状態を計算する粒子状態計算ステップと、粒子状態計算ステップにおける対象の時間ステップが、予め設定された周期の時間ステップか否かを判断する周期判断ステップと、周期判断ステップにおいて予め設定された周期の時間ステップであると判断された時間ステップに、粒子状態計算ステップにおいて計算された粒子の位置に基づいて、当該粒子から、自領域に隣接する分割された計算領域である隣接領域において計算に用いられる自領域仮想粒子を特定して、特定した自領域仮想粒子を示す情報を当該隣接領域の粒子の状態を計算する別の粒子状態計算装置に出力する自仮想粒子情報出力ステップと、別の粒子状態計算装置から、自領域において計算に用いられる他領域仮想粒子を特定する情報を入力する他仮想粒子情報入力ステップと、各時間ステップにおいて粒子状態計算ステップにおいて算出された自領域仮想粒子の状態を示す情報を別の粒子状態計算装置に出力する自仮想粒子状態出力ステップと、各時間ステップにおいて別の粒子状態計算装置から他領域仮想粒子の状態を示す情報を入力する他仮想粒子状態入力ステップと、粒子状態計算ステップにおける計算に応じた情報を出力する出力ステップと、を含み、粒子状態計算ステップにおいて、他仮想粒子情報入力ステップにおいて入力される他領域仮想粒子を特定する情報、及び他仮想粒子状態入力ステップにおいて入力された他領域仮想粒子の状態を示す情報を用いて粒子の状態を計算することを特徴とする。
【発明の効果】
【００１７】
本発明によれば、粒子状態計算装置間でやり取りする情報に係る粒子を限定された仮想粒子とすることで、各粒子状態計算装置で必要な情報を確実に入力させると共に粒子状態計算装置間の通信を効率的にする。また、仮想粒子を特定する情報のやりとりを、予め設定された周期の時間ステップとすることで粒子状態計算装置間の通信が更に効率的になる。これによって、本発明によれば、効率的に分散メモリ型並列計算により粒子状態を計算し、その結果より大規模なモデルに対する計算が可能となる。
【図面の簡単な説明】
【００１８】
【図１】本発明の実施形態に係る粒子状態計算装置であるノードの機能構成を示す図である。
【図２】流体及び粒子の計算領域及び（流体計算における）計算領域の分割を示す図である。
【図３】計算領域のセルへの分割を示す図である。
【図４】複数のノードにより構成される粒子状態計算システム全体で実行される処理を示すフローチャートである。
【図５】本発明の実施形態に係る粒子状態計算装置である各ノードで実行される処理（粒子状態計算方法）を示すフローチャートである。
【図６】図５のフローチャートの各処理に応じた計算領域における粒子を示す図であり、説明のために平面図としている。
【図７】本実施形態に係る粒子状態計算装置及び粒子状態計算方法により得られる粒子の衝突頻度因子の結果を示すグラフである。
【図８】本実施形態に係る粒子状態計算装置及び粒子状態計算方法により得られるストークス数が１の粒子の衝突頻度因子の結果を示すグラフである。
【図９】Ｎ２５６計算に対して演算プロセス数を変化させた場合の（ａ）１時間ステップに要した演算時間、（ｂ）演算速度を示すグラフである。
【発明を実施するための形態】
【００１９】
以下、図面と共に本発明に係る粒子状態計算装置及び粒子状態計算方法の好適な実施形態について詳細に説明する。なお、図面の説明においては同一要素には同一符号を付し、重複する説明を省略する。
【００２０】
図１に本実施形態に係る粒子状態計算装置であるノード１０を示す。ノード１０は、別のノード１０と共に粒子状態計算システム１を構成する。粒子状態計算システム１は、例えば、雲を構成する水粒子等の粒子の状態を計算（シミュレーション）するシステムである。各ノード１０は、ＣＰＵ（Central Processing Unit）、メモリ、ハードディスク等のハードウェアを備えるコンピュータである。各ノード１０はイントラネットやインターネット等のネットワークで接続されており、各ノード１０間で情報の送受信が可能になっている。粒子の状態の計算は、各ノード１０が別個の物理メモリを用いて分散して計算を行う分散メモリ型並列計算によって行われる。ここで各ノード１０における、計算の処理を（演算）プロセスと呼ぶ。
【００２１】
計算対象の粒子は、予め大きさが設定される計算領域に配置される。計算領域は複数の領域に分割される。各ノード１０は、分割された領域の一つが計算対象として割り当てられて、割り当てられた領域である自領域に配置される粒子の状態を時間ステップ毎に繰り返し計算する。計算される粒子の状態は、具体的には粒子の位置及び速度である。また、計算領域には流体（より具体的には、非圧縮流体）が流れており、粒子は流体中にあり流体の影響を受けることとしてもよい（即ち、流体についての計算も行う）。流体は、例えば、水粒子の状態に影響を与える気流である。計算される流体の状態は、具体的には流体の速度及び圧力である。
【００２２】
計算対象の領域が隣接しているノード１０同士は、粒子の移動等に応じて粒子の情報を送受信して、粒子の計算に利用する。
【００２３】
なお、粒子状態計算システム１には、分散メモリ型並列計算を統括的に制御する装置も含まれている。その装置が計算のためのパラメータの入力、計算領域の分割、及び各ノード１０への計算の指示等の分散メモリ型並列計算を統括的に制御する処理を行う。当該装置が複数のノード１０のうちの何れかであってもよい。
【００２４】
ここで、粒子及び流体の計算について説明する。まず、流体の支配方程式と数値計算法について説明する。非圧縮流体の方程式は、次の連続の式と、ナビエ−ストークス方程式である。
【数１】
ここで、ｕｉは流体の速度ベクトル、ｐは流体の圧力を示す変数である。ｔは時間ステップを示す変数である。添え字ｉ，ｊは、ｘ，ｙ，ｚ何れかの方向を表しており、式（１）及び（２）は総和規約に従って表記されている。全ての変数は、代表長さＬ０と代表速度Ｕ０とによって無次元化されている。また、Ｒｅは動粘性係数νを用いてＲｅ＝Ｕ０Ｌ０／νと定義されるパラメータとして与えられる値である。例えば、νを１気圧２９８Ｋでの値１．５×１０−５ｍ２／ｓに設定される。式（２）の右辺最終項ｆｉは外力項である。
【００２５】
対流項の差分の計算には、Morinishi, Y., Lund, T.,Vasilyev, O. and Moin, P., “Fully Conservative Higher Order Finite DifferenceSchemes for Incompressible Flow”, Journal of Computational Physics, Vol. 143(1998), pp.90-124. (4th-order scheme for Non-linear term)に示される保存型４次中心差分スキーム、粘性項の差分の計算には４次中心差分法を用いることができる。時間進行（時間発展）の計算には、２次のルンゲ−クッタ法を用い、時間発展後の速度が連続の式を満たすように速度と圧力とを同時緩和させながら反復計算するＨＳＭＡＣ法（Hirt, C.W. and Cook, J.L., “Calculating Three-Dimensional Flowsaround Structures and over Rough Terrain”, Journal of Computational Physics,Vol. 10 (1972), pp.324-340）を用いることができる。その反復計算の中での空間差分法には２次中心差分法を用いることができる。反復計算の収束判定条件は、例えば、∂ｉｕｉ／（Ｕ０／Δ）＜１０−３（Δは計算格子幅）とすることができる。
【００２６】
外力項ｆｉには、Reduced-CommunicationForcing（ＲＣＦ）法を用いることができる（Onishi, R, Baba, Y. andTakahashi, K., “Efficient Large-Scale Forcing in Finite-Difference Simulationsof Steady Isotropic Turbulence”, Topics on Chaotic Systems -selected papersfrom CHAOS 2010 international conference, Skiadas, C.H. et al. eds.,World Scientific Publishing Co. (in press)）。これは、大スケール運動にのみエネルギーを注入する大スケール強制法の一つであり、有限差分法において高並列効率を達成可能な強制法である。例えば、波数ｋの大きさが｜ｋ｜＜２．５の大スケール運動のエネルギーが保たれるように設定される。
【００２７】
続いて、流体の計算領域と並列計算法について説明する。図２に示すように、一辺２πＬ０の立方体計算領域に対して、Ｎ３個のスタッガード格子を配する。つまり、計算格子幅Δ＝２πＬ０／Ｎとなる。ここで、Ｎは予め設定あるいは入力されるパラメータである。具体的な計算コードは、ＭＰＩライブラリを用いたＦｏｒｔｒａｎ９０コードで記述することができる。これは、３次元領域分割に対応しており、立方体計算領域の座標軸であるｘ、ｙ及びｚの各方向を等間隔に分割することができる。各方向の分割数Ｍｘ、Ｍｙ及びＭｚは、ＭＰＩプロセス総数ＭｐｒｏｃとＭｐｒｏｃ＝ＭｘＭｙＭｚの関係にある。このとき、各プロセスはｎｘ×ｎｙ×ｎｚ（ｎｘ＝Ｎ／Ｍｘ、ｎｙ＝Ｎ／Ｍｙ及びｎｚ＝Ｎ／Ｍｚ）の大きさを持った“ＭＰＩ分割領域”内の流体及び粒子運動の計算を担当する。Ｍｘ、Ｍｙ及びＭｚは、任意のＮの約数の値に予め設定あるいは入力されるパラメータである。なお、上記の例では、Ｎ３に等間隔に分割することとしているが、Ｎｘ，Ｎｙ，Ｎｚと各方向に分割数を変えることもできる。
【００２８】
例えば、以下の表に示す計算条件で気相乱流場の計算が行われる。全てのケースで、クーラン数（Ｕ０Δｔ／Δ）が０．２以下になるように時間刻みΔｔを設定した。ここで、時間刻みΔｔは繰り返し計算の時間ステップの間隔である。但し、粒子の大きさによっては、口述するようにΔｔを更に小さく設定しなければならない場合もある。
【表１】
【００２９】
引き続いて、粒子の状態の計算について説明する。粒子の運動は、例えば、ラグラジアン粒子追跡法を用いて計算することができる。線形ストークス抗力を仮定すると支配方程式は次式となる。
【数２】
ここで、ｖｐは粒子速度、ｕｉ（アンダーバーの記載は省略する）は粒子位置における流体速度を示す変数である。ここで、添え字ｉはｘ，ｙ，ｚ方向を表す。τｐは粒子緩和時間である。式（１）と同様に全ての変数は、代表長さＬ０と代表速度Ｕ０とによって無次元化されている。また、粒子に働く重力は無視する。粒子緩和時間τｐは、次式で計算される。
【数３】
ここで、ｒは粒子半径、ρｐ及びρｆはそれぞれ粒子及び流体の密度であり、予め設定あるいは入力されるパラメータである。例えば、雲粒を想定して１気圧２９８Ｋでの水及び空気の密度の値から、ρｐ／ρｆ＝８４０とする。線形ストークス抗力モデルは粒子レイノルズ数Ｒｅｐ（＝２ｒ｜ｖｐ−ｕ｜／ν）が１に比べて十分小さい場合に適用可能である。雲粒を想定した場合、およそｒ＝４０μｍのときＲｅｐ＝１である。それ以上大きな粒子の場合には非線形抗力モデルを用いるべきであるが、簡単のために全ての場合で線形ストークス抗力モデルを用いてもよい。また、現実の雲では水滴の空気に対する質量割合最大でも１０−３程度であるので、粒子は希薄であり流体運動に影響を与えないとしてもよい。粒子位置における流体速度ｕは、粒子近傍の８点の速度データから線形補間により算出される。時間進行の計算には、流体計算と同じく、２次のルンゲ−クッタ法を用いる。
【００３０】
続いて、粒子の衝突の計算について説明する。衝突は２粒子間衝突のみを考慮する。２粒子がΔｔ間に移動する軌跡を求めて、それらの軌跡の最小距離が粒子直径以下であるか否かを判断して、粒子直径以下であれば衝突とみなす。衝突後には一方の粒子を消滅させる。従って、一回衝突が起こると粒子の総数が１個減少する。これを利用して、粒子総数の減少率から衝突頻度Ｎｃ及び衝突頻度因子Ｋｃ（＝２Ｎｃ／ｎｐ２、ｎｐは粒子数密度）を算出する。ある時間ステップｎにおける粒子総数Ｎｐｎ、１ステップ前の粒子総数Ｎｐｎ−１及び計算領域の体積Ｖｄ（＝（２πＬ０）３）から、ｎ−１／２ステップにおける衝突頻度Ｎｃｎ−１／２は以下のように算出される。
【数４】
この衝突頻度とｎ−１／２ステップにおける粒子数密度ｎｐｎ−１／２（＝０．５（Ｎｐｎ−１＋Ｎｐｎ）／Ｖｄ）から、ｎ−１／２ステップにおける衝突頻度因子を次の式により算出する。
【数５】
この衝突頻度因子を時間平均することによって平均衝突因子を算出する。
【００３１】
計算領域内にＮｐ個の粒子が存在する場合、特殊な方法を使わなければ、Ｎｐ（Ｎｐ−１）／２〜０．５Ｎｐ２回の衝突判定計算をしなければならない。衝突判定回数を削減する方法として分子動力学法で使用されるセルインデックス法（Allen, M.P., and Tildesley, D.J., “Computer Simulation of Liquids”,Oxford University Press, (1987), p.408）がある。これは、図３に示すように衝突頻度領域をセルに小分けし、ある粒子を含むセルとそのセルとを囲むｘ，ｙ，ｚ各方向の前後、合計２７個のセルの中の粒子との衝突を判定する。但し、逐次的に衝突を判定していく場合には、２７個ではなく、８個のセル（例えば、対象粒子を含むセルとｘ，ｙ，ｚ各方向前方側のセル）を判定の対象とすれば、重複判定がなく効率がよい。このセルインデックス法では、粒子一つ一つに番号を付けた上で、何番と何番の粒子がどのセル内にあるかを管理するリスト（セルインデックス）を作成する。
【００３２】
図３に示すように計算領域をＭｃｅｌｌ３個のセルに分割した場合、各セル内の平均粒子数Ｎｐ，ｃｅｌｌはＮｐ，ｃｅｌｌ＝Ｎｐ／Ｍｃｅｌｌ３となる。ある粒子一つに対する衝突判定回数は、粒子を含むセル内の他粒子との衝突判定回数（Ｎｐ，ｃｅｌｌ−１）／２と、ｘ，ｙ，ｚ各方向後方側のセル内粒子との衝突判定回数７Ｎｐ，ｃｅｌｌの和、つまり（Ｎｐ，ｃｅｌｌ−１）／２＋７Ｎｐ，ｃｅｌｌと見積もられる。結局、衝突判定回数の総数は、Ｎｐ×（（Ｎｐ，ｃｅｌｌ−１）／２＋７Ｎｐ，ｃｅｌｌ）〜７．５Ｎｐ２／Ｍｃｅｌｌ３である。つまり、Ｍｃｅｌｌ＞２．５に設定すれば衝突判定回数の回数を削減できる。この見積もりでは、Ｍｃｅｌｌを大きくすればするほど衝突判定回数の回数は小さくなる。しかし、セルの中に粒子がほとんどないような状況では、セルの中に粒子があるかないかを調べる処理が無駄になる。そのため、Ｎｐ，ｃｅｌｌが１以上になるようにＭｃｅｌｌを設定するのが望ましい。
【００３３】
本実施形態では、有限差分法を用いて気相乱流を計算し、乱流中で運動する粒子同士の衝突を検出する。ＭＰＩライブラリを用いた分散メモリ型並列計算に対応させるためには、下記の３つの対応するＭＰＩ通信の機能を実現しなければならない。即ち、当該機能に対応するＭＰＩ通信プログラムを記述しなければならない。
（Ａ）差分計算に必要な袖領域の気相データの取得
（Ｂ）ＭＰＩ分割領域（個々のノード１０に割り当てられた計算領域）からはみ出した粒子を隣のプロセスに移動させる処理
（Ｃ）ＭＰＩ分割領域の境界近傍にある粒子情報の共有
ある粒子ペアがＭＰＩ分割領域の境界をまたいで近接する場合、その情報を両側のプロセスが共有していなければ衝突を検出できない可能性がある。この問題に対処するため上記の（Ｃ）の操作が必要になる。本発明は、セルインデックスを活用して、（Ｂ）（Ｃ）を効率的に処理するためものである。
【００３４】
引き続いて、ノード１０の本発明の係る機能を詳細に説明する。図１に示すようにノード１０は、粒子状態計算部１１と、周期判断部１２と、自仮想粒子情報出力部１３と、他仮想粒子情報入力部１４と、自仮想粒子状態出力部１５と、他仮想粒子状態入力部１６と、出力部１７とを備えて構成される。なお、各構成要素の更に具体的な機能、各構成要素によって行われる更に具体的な処理については、後述する処理フローの説明で説明する。これらの機能は、例えばノード１０においてプログラムが実行されることによって実現される。
【００３５】
粒子状態計算部１１は、各時間ステップにおける自領域の粒子の状態を計算する粒子状態計算手段である。自領域とは、自ノード１０に割り当てられた計算領域である。粒子状態計算部１１は、具体的には上述した方法によって粒子の状態を計算する。なお、計算に必要なパラメータについては、上述した分散メモリ型並列計算を統括的に制御する装置から入力される。また、計算に必要な計算式等は予め粒子状態計算部１１に記憶されており、それらが読み出されて計算が行われる（他の構成要素についても同様）。
【００３６】
粒子状態計算部１１は、粒子間の相互作用（本実施形態では、粒子同士の衝突）の判断、計算のために上述したようにセルインデックス法に応じた処理を行う。即ち、各粒子に粒子を特定する情報である番号を付与する。粒子状態計算部１１は、自領域にある粒子の番号を、当該粒子が含まれるセルを特定する情報に対応付けて記憶、管理する。セルは、上述したように自領域を更に分割した領域である。
【００３７】
粒子状態計算部１１は、後述するように、他仮想粒子情報入力部１４によって入力される他領域仮想粒子を特定する情報、及び他仮想粒子状態入力部１６によって入力された他領域仮想粒子の状態を示す情報を用いて粒子の状態を計算する。仮想粒子とは、自領域（自ノード１０）における状態の計算対象の粒子ではあるが、自領域に隣接する分割された計算領域である隣接領域（別のノード１０）における粒子の状態の計算に必要な粒子である。各ノード１０間で、仮想粒子の情報を送受信して、その情報を粒子の計算に利用する。自領域にある粒子のうち、隣接領域において計算に用いられる仮想粒子を自領域仮想粒子と呼ぶ。隣接領域にある粒子のうち、自領域において計算に用いられる仮想粒子を他領域仮想粒子と呼ぶ。
【００３８】
粒子状態計算部１１は、隣接領域（別のノード１０）との間で、自領域からはみ出して隣接領域に入った粒子の情報及び隣接領域からはみ出して自領域に入った粒子の情報の送受信も行う。
【００３９】
また、粒子状態計算部１１は、上述したように、有限差分法によって流体の状態を算出して、当該流体の影響を考慮して粒子の状態を計算する。
【００４０】
粒子状態計算部１１は、他仮想粒子情報入力部１４によって入力される他領域仮想粒子を特定する情報に基づいて、粒子間の近傍関係を特定して、当該近傍関係に基づいて粒子の状態を計算する。本実施形態では近傍関係とは衝突が判断される粒子のペアである。具体的には、後述するようにセルインデックスを用いる。
【００４１】
周期判断部１２は、粒子状態計算部１１における計算の時間ステップが、予め設定された周期の時間ステップか否かを判断する周期判断手段である。この周期を示す値は２以上の自然数であり、予め周期判断部１２に記憶されている。周期判断部１２は、例えば、粒子状態計算部１１によって管理される現時点での時間ステップの情報を参照することによって上記の判断を行う。周期判断部１２は、計算の時間ステップが予め設定された周期の時間ステップであると判断するとその旨を粒子状態計算部１１及び自仮想粒子情報出力部１３に通知する。
【００４２】
自仮想粒子情報出力部１３は、周期判断部１２によって予め設定された周期の時間ステップであると判断された時間ステップに、自領域仮想粒子を特定して、特定した自領域仮想粒子を示す情報を隣接領域の粒子の状態を計算する別のノード１０に出力する自仮想粒子情報出力手段である。自領域仮想粒子の特定は、粒子状態計算部１１によって計算された粒子の位置に基づいて行われる。自領域中のセルのうち、隣接領域に隣接するセルに位置する粒子が自領域仮想粒子とされる。自仮想粒子情報出力部１３は、隣接領域毎に自領域仮想粒子の特定し、自領域仮想粒子の番号を、隣接領域の粒子の状態を計算する別のノード１０に出力する。また、自仮想粒子情報出力部１３は、自領域仮想粒子の番号を自仮想粒子状態出力部１５に通知する。
【００４３】
他仮想粒子情報入力部１４は、隣接領域の粒子の状態を計算する別のノード１０から、他領域仮想粒子を特定する情報である他領域仮想粒子の番号を入力する他仮想粒子情報入力手段である。他仮想粒子情報入力部１４は、入力した他領域仮想粒子の番号を粒子状態計算部１１に通知する。
【００４４】
自仮想粒子状態出力部１５は、各時間ステップにおいて粒子状態計算部１１によって算出された自領域仮想粒子の状態を示す情報を別のノード１０に出力する自仮想粒子状態出力手段である。自仮想粒子状態出力部１５は、自仮想粒子情報出力部１３から通知された自領域仮想粒子の番号によって示される粒子の状態の情報を粒子状態計算部１１から取得して別のノード１０に出力する。出力対象のノード１０は、自仮想粒子情報出力部１３が番号の情報を出力したノード１０である。出力する粒子の状態の情報は、粒子の位置（座標）を示す情報と、粒子の速度（方向、大きさ）を示す情報である。
【００４５】
他仮想粒子状態入力部１６は、各時間ステップにおいて別のノード１０から他領域仮想粒子の状態を示す情報を入力する他仮想粒子状態入力手段である。自仮想粒子状態出力部１５は、入力した他領域仮想粒子の状態を示す情報を粒子状態計算部１１に出力する。粒子の情報を入力する他領域仮想粒子は、他仮想粒子情報入力部１４によって特定する情報が入力された他領域仮想粒子である。他仮想粒子状態入力部１６は、入力した他領域仮想粒子の状態を示す情報を粒子状態計算部１１に出力する。
【００４６】
出力部１７は、粒子状態計算部１１による計算に応じた情報を出力する出力手段である。具体的には例えば、出力部１７は、粒子状態計算部１１による計算により得られる、上述した平均衝突頻度因子を計算するための情報を、分散メモリ型並列計算を統括的に制御する装置に出力する。
【００４７】
引き続いて、図４及び図５のフローチャートを用いて、本実施形態に係るノード１０で実行される処理（粒子状態計算方法）を説明する。まず、図４のフローチャートを用いて、粒子状態計算システム１全体の処理の概要について説明し、続いて、図５のフローチャートを用いて各ノード１０による本発明に係る処理について説明する。
【００４８】
粒子状態計算システム１では、まず、分散メモリ型並列計算を統括的に制御する装置に計算用のパラメータが入力される（Ｓ０１）。このパラメータは、計算領域の大きさ、流体計算のための格子数、粒子計算のためのセル数、粒子数、各粒子の大きさ、密度、速度、分散メモリ型並列計算の分割数及び時間刻み（時間ステップの間隔）等の情報である。この入力は、例えば、ユーザによる当該装置の操作により行われる。
【００４９】
分散メモリ型並列計算を統括的に制御する装置では、入力されたパラメータに基づいて、各ノード１０が担当する計算領域を決定して、各ノードに対して計算に必要な情報を出力する（Ｓ０２）。各ノード１０では計算が行われる（Ｓ０３）。続いて、各ノード１０の計算結果に応じた情報（例えば、平均衝突頻度因子）が出力される（Ｓ０４）。結果の出力は、例えば、各ノード１０による計算結果が、分散メモリ型並列計算を統括的に制御する装置に集められて、当該装置によって最終結果が生成されて出力される。
【００５０】
引き続いて、図５のフローチャートを用いて、各ノード１０における処理（図４のフローチャートにおけるＳ０３に相当）について説明する。また、図６に各処理に対応する計算領域を示す。図６に示す図では、計算領域は２次元領域としており、あるプロセスが受け持つＭｃｅｌｌ，ｘ×Ｍｃｅｌｌ，ｙ個のセル（Ｍｃｅｌｌ，ｘ＝Ｍｃｅｌｌ／Ｍｘ及びＭｃｅｌｌ，ｙ＝Ｍｃｅｌｌ／Ｍｙ）が実線で描かれている。実線で描かれた領域は対象プロセスが担当するＭＰＩ分割領域（自領域）でもあり、点線部分は別のプロセス（近傍プロセス）が計算を担当している領域を示している。以下では、セルインデックス作成時にＭＰＩ分割領域内にある粒子を実粒子と呼ぶ。衝突判定のためには、Ｍｃｅｌｌ，ｘ＋１やＭｃｅｌｌ，ｙ＋１の領域にある粒子も考慮する必要がある。それらの粒子が、仮想粒子である。仮想粒子は衝突検出のためだけに用いられて、その運動は計算されない。図６において黒塗りの点は実粒子であり、白抜きの点は（他）仮想粒子である。
【００５１】
本処理では、まず、粒子状態計算部１１によって所定の時間ステップにおける流体の状態の計算が行われる（Ｓ１１、粒子状態計算ステップ）。計算される流体の状態は流体の速度及び圧力である。また、粒子状態計算部１１によって所定の時間ステップにおける実粒子の状態の計算が行われる（Ｓ１２、粒子状態計算ステップ）。計算される実粒子の状態は実粒子の位置及び速度である。これらの計算は、前の時間ステップ又は初期状態における流体及び実粒子の状態、並びに計算後の流体の状態に基づいて、上述したようにルンゲ−クッタ法により行われる。なお、流体及び実粒子の状態の計算は、必ずしも順番に行われるわけではなく、流体及び実粒子それぞれについてルンゲ−クッタ法の１次ステップ及び２次ステップ等のように交互あるいは同時に計算される。図６（ａ）に示す状態が前の時間ステップ又は初期状態（即ち、Ｓ１１の直前）であり、図６（ｂ）に示す状態が計算後の状態（即ち、Ｓ１２の直後）である。
【００５２】
Ｓ１２の後、図６（ｂ）に示すようにＭＰＩ分割計算領域からはみ出す（実）粒子が現れる。はみ出した粒子の処理（上述した（Ｂ）の処理）を毎時間ステップ実行すると、多大なノード１０間の通信が発生して、並列性能を低下させるおそれがある。そこで、以下のように予め設定された周期のみにはみ出した粒子の処理を行わせる。はみ出した実粒子の計算を続けるためにＭＰＩ分割計算領域外の流体速度データが必要になる。各プロセス（粒子状態計算部１１）は、上述（Ａ）の処理に伴って、４次中心差分計算のために４格子分の袖データを持つ。そのため、粒子が３．５Δ分（スタッガード格子を採用しているため４Δではなく３．５Δ）以上はみ出さない限りは、余分なＭＰＩ通信（ノード１０間の通信）をすることなくｕを算出することができる。つまり、流体計算を有限差分法によって計算するからこそ可能な効率化手法である。また、以下に示すようにセルインデックス作成も予め設定された周期毎に行う。
【００５３】
Ｓ１２の後、周期判断部１２によって、現時点の時間ステップが予め設定された周期の時間ステップか否かが判断される（Ｓ１３、周期判断ステップ）。現時点の時間ステップが予め設定された周期の時間ステップであると判断されると（Ｓ１３のＹｅｓ）、周期判断部１２から粒子状態計算部１１及び自仮想粒子情報出力部１３にその旨が通知され、その場合の処理（Ｓ１４〜Ｓ１９）が行われる。
【００５４】
まず、図６（ｃ）に示すように、粒子状態計算部１１において記憶されている他仮想粒子に関する情報が除去される（Ｓ１４）。続いて、粒子状態計算部１１によって自領域から出た実粒子の情報が入出力される（Ｓ１５、粒子状態計算ステップ）。具体的には、自領域から出た実粒子の情報が、当該実粒子が位置する領域の計算を担当するノードに出力され、自領域における実粒子から削除される。この実粒子は、他領域における実粒子となる。一方、他領域から出て自領域に入った粒子の情報が、他領域の計算を担当するノード１０から受け取られ、自領域の実粒子とされる（図６（ｄ）の矢印で示される粒子である）。即ち、実粒子のノード（領域）間での引き渡しが行われる。ここでノード１０間の通信（ＭＰＩ通信）が発生する。
【００５５】
続いて、粒子状態計算部１１によってセルインデックスが生成（更新）される。セルインデックスは、セルに含まれる（セル上に位置する）実粒子の番号をそのセルを特定する情報に対応付けた情報である（Ｓ１６、粒子状態計算ステップ）。図６（ｅ）に示す（２，２）のセルの実粒子のリストは、［４，１０２，６０５，６５１］である。セルインデックスは、粒子状態計算部１１によって記憶され管理される。なお、（２，２）のセルとは、ｘ軸において２番目、ｙ軸において２番目のセルであることを示す（以下同様）。
【００５６】
続いて、自仮想粒子情報出力部１３によって、自領域仮想粒子が特定されて、自領域仮想粒子を示す情報が隣接領域の粒子の状態を計算する別のノード１０に出力される（Ｓ１７、自仮想粒子情報出力ステップ）。具体的には、粒子状態計算部１１によってセルインデックスにおいて、そのセルに含むセルが自領域仮想粒子であるとされるセルに対応付けられた実粒子が自領域仮想粒子とされる。そのようなセルは、領域間の位置関係によって決定される。当該セルを特定する情報と自領域仮想粒子の番号とが自仮想粒子情報出力部１３から別のノード１０に出力される。また、自領域仮想粒子の番号は、自仮想粒子情報出力部１３から自仮想粒子状態出力部１５に通知される。
【００５７】
その一方で、別のノード１０から送信される他領域仮想粒子の番号と当該粒子を含むセルとが対応付けられた情報が、他仮想粒子情報入力部１４によって入力される（Ｓ１８、他仮想粒子情報入力ステップ）。この情報は、別のノード１０におけるＳ１７に相当する処理によって送信されたものである。上記の情報は、他仮想粒子情報入力部１４から粒子状態計算部１１に入力される。これにより、図６（ｆ）に示すように、粒子状態計算部１１によって自領域に隣接するセルに位置する他仮想粒子の番号が把握される。Ｓ１７及びＳ１８ではノード１０間の通信（ＭＰＩ通信）が発生する。
【００５８】
続いて、粒子状態計算部１１によって、自領域の各粒子について、衝突の有無をチェックする粒子のペアが生成される（Ｓ１９、粒子状態計算ステップ）。この粒子のペアは、対象粒子と当該対象粒子を含むセル及び当該セルに隣接しているセルに位置する粒子とのペアである。このペアの生成の際には、他仮想粒子も考慮される（実粒子と他仮想粒子のペア）。
【００５９】
（ｉ，ｊ）番目のセル（（ｉ，ｊ）セル）は、自分自身を含めて、（ｉ，ｊ）、（ｉ，ｊ±１）、（ｉ±１，ｊ）、（ｉ±１，ｊ±１）の計９個のセルに囲まれている（３次元の場合は２７個）。しかし、セル番号を順に増やして衝突をチェックする場合には、（ｉ，ｊ）、（ｉ，ｊ＋１）、（ｉ＋１，ｊ）、（ｉ＋１，ｊ＋１）の計４個（３次元の場合は８個）のセルを調べていけば、重複することなく全ての組み合わせをチェックできる。そのため粒子のペアのリストを作る際には、ｉ＝０やｊ＝０のセルを含める必要はない。図６（ｆ）に示す番号が１３０のペアとなるリストは、上記の４つのセルに含まれる粒子である［３３２，６９，８７１，７７７，１０９，４，６０５，１０２，６５１］である。
【００６０】
Ｓ１９の後、あるいはＳ１２の後、周期判断部１２によって現時点の時間ステップが予め設定された周期の時間ステップでないと判断された場合（Ｓ１３のＮｏ）、以下の処理が行われる。
【００６１】
まず、自仮想粒子状態出力部１５によって、自仮想粒子情報出力部１３から通知された自領域仮想粒子のその時点の時間ステップの状態を示す情報が取得される。取得される情報は、Ｓ１２において粒子状態計算部１１よって計算されたものである。取得された情報は、自仮想粒子状態出力部１５から上記の別のノード１０に出力される（Ｓ２０、自仮想粒子状態出力ステップ）。
【００６２】
その一方で、別のノード１０から送信される他領域仮想粒子の状態を示す情報が、他仮想粒子状態入力部１６によって入力される（Ｓ２１、他仮想粒子状態入力ステップ）。この情報は、別のノード１０におけるＳ２０に相当する処理によって送信されたものである。上記の情報は、他仮想粒子状態入力部１６から粒子状態計算部１１に入力される。これにより、図６（ｇ）に示すように、粒子状態計算部１１によって自領域に隣接するセルに位置する他仮想粒子の状態である位置及び速度が把握される。Ｓ２１及びＳ２２ではノード１０間の通信（ＭＰＩ通信）が発生する。
【００６３】
続いて、粒子状態計算部１１では、図６（ｈ）に示すように、上述したように粒子間の衝突の判断が行われて、その判断に基づく計算が行われる（Ｓ２２、粒子状態計算ステップ）。この判断は、Ｓ１９で生成された全てのペアの間で行われる。実粒子間の衝突の判断は、それぞれＳ１２で計算された粒子の状態に基づいて行われる。実粒子と他領域仮想粒子との衝突の判断は、実粒子についてはＳ１２で計算された粒子の状態に基づいて、他領域仮想粒子についてはＳ２１で入力された粒子の状態に基づいて行われる。
【００６４】
続いて、粒子状態計算部１１では、計算の終了の判断が行われる（Ｓ２３、粒子状態計算ステップ）。この判断は、粒子状態計算部１１に予め設定あるいは入力された判断基準に基づいて行われる。計算が終了しないと判断された場合（Ｓ２３のＮｏ）には、次の時間ステップに進められてＳ１１からの処理が繰り返される（Ｓ２４、粒子状態計算ステップ）。
【００６５】
計算が終了する判断された場合（Ｓ２３のＹｅｓ）には、出力部１７によって、粒子状態計算部１１による計算に応じた情報が出力される（Ｓ２５、出力ステップ）。以上で、ノード１０における処理が終了する。なお、上記の処理における時間ステップの更新（Ｓ２４）の後、毎ステップあるいは数ステップに一度、例えば、粒子の衝突、粒子の速度及び粒子のエネルギー等の統計量を算出して出力することとしてもよい。
【００６６】
上述した処理では、一度セルインデックスが生成された後、予め設定された周期の時間ステップを経過するまで更新されないので、その間に粒子は当初のセルから移動する可能性がある。しかし、上記の周期の時間ステップの間の粒子の移動距離がセル幅Δｃｅｌｌ（＝２πＬ０／Ｍｃｅｌｌ）よりも小さければ、衝突する可能性のある粒子ペアを漏らす可能性を無視できる。本実施形態では、例えば、Δｃｅｌｌ＝４Δ（つまり、Ｍｃｅｌｌ＝Ｎ／４）のとき、上記の周期を８に設定し、クーラン数を０．２以下に設定すると、粒子は流体速度程度で移動すると考えると上記の周期の時間ステップの間に粒子はせいぜい１．６Δ程度しか移動しない。つまり、この設定では、衝突検出の漏れは無視できるほど小さい。実際、いくつかのケースで上記の周期を１に設定した場合と、８に設定した場合とで比較して、衝突ペアに変化が無いことを確認した。
【００６７】
上述したように本実施形態では、ノード１０間でやり取りする情報に係る粒子を限定された仮想粒子とすることで、各ノード１０で必要な情報を確実に入力させると共にノード１０間の通信を効率的にする。また、仮想粒子を特定する情報のやりとりを、予め設定された周期の時間ステップとすることでノード１０間の通信が更に効率的になる。これによって、本発明に係る本実施形態によれば、効率的に分散メモリ型並列計算により粒子状態を計算し、その結果より大規模なモデルに対する計算が可能となる。
【００６８】
また、本実施形態のように流体中の粒子を計算することとすれば、流体中の粒子に係る大規模なモデル計算が可能となる。但し、必ずしも流体中の粒子の計算でなくてもよく、分散メモリ型並列計算により粒子状態を計算するものに対して適用可能である。
【００６９】
また、本実施形態のように、セルインデックスを用いて粒子間のペアの相互作用を判断することとすれば、粒子間の衝突等の粒子間の相互作用が生じるモデルの計算を確実に行うことができる。
【００７０】
なお、上述した実施形態における粒子の計算については、定常な気相乱流場が形成されてから、粒子を混入させて行うこととしてもよい。また、衝突の判定については、粒子の混入から一定時間経過後（例えば、Ｔ０（＝Ｌ０／Ｕ０）以上）としてもよい。
【００７１】
なお、本実施形態においては雲をシミュレーションする例を示したが、必ずしも雲を対象とした計算に限られず、粒子の状態を計算するものであれば適用が可能である。例えば、小惑星の生成シミュレーション等にも適用可能である。
【００７２】
以下に、上述した実施形態による具体的な計算結果を示す。なお、定常な気相乱流相場が形成されたことを確認後、粒子を混入させた。粒子運動を計算する場合、時間刻みΔｔを０．５τｐ以下に設定した。この条件と流体計算のクーラン数による条件を同時に満たすΔｔを流体計算と粒子計算との両方に用いた。また、粒子混入からＴ０（＝Ｌ０／Ｕ０）以上経過後から衝突判定を開始した。
【００７３】
以下の表に粒子の計算条件を示す。最左列に示されているケース名は、上記の表と対応している。全てのケースで粒子数密度ｎｐが雲の代表値である１０８オーダーになるように設定した。セルに含まれる平均粒子数はＮ６４及びＮ２５６の場合は１、Ｎ５１２の場合は４．８であった。
【表２】
【００７４】
以下の表に、統計的に定常状態に達した気相等方性乱流場の乱流特性量を示す。表中のｕ´は速度変動強度、Ｒｅλ（＝ｕ´ｌλ／ν、ｌλはテイラーマイクロスケール）はテイラーマイクロスケール基準乱流レイノルズ数、ｋｍａｘ（＝Ｎ／２）は最大有効波数、ｌη（＝（ν３／ε）１／４、εはエネルギー離散率）はコルモゴロフスケールである。ｋｍａｘｌηは計算解像度の指標であり、運動エネルギー等の２次元量を議論する場合には１程度、スキューネスやフラットネス等の高次元を議論する場合には２程度必要だと言われている。従来の方法においては粒子衝突の計算では１．３〜１．４程度に設定されることが多いが、本例では通常よりも高い解像度である２程度になるように設定した。
【表３】
【００７５】
図７にＮ６４の場合（Ｒｅλ＝５２．３）の衝突頻度因子の結果を示す。縦軸は局所速度勾配λ（＝（ε／ν）１／２）と衝突半径Ｒ（＝２ｒ）で無次元化された衝突頻度因子、横軸はコルモゴロフスケールで無次元化された粒子半径である。衝突によって粒子数が５％減少するまでを一つのランとして、それぞれのランにおける平均衝突頻度因子を算出した。比較のために、Onishi et al. (2009)（Onishi, R., Takahashi,K. and Komori, S., “Influence of Gravity on Collisions of MonodispersedDroplets in Homogeneous Isotropic Turbulence”, Physics of Fluids, Vol. 21(2009), 125108）で得られたＲｅλ＝５４．３の場合の結果も示す。Onishiet al. (2009)は、流体計算に擬スペクトル法を用いた。その点は、有限差分法を用いる本実施形態と異なるが、計算領域の大きさ、物性値、計算格子数、粒子数等の設定、そして標準偏差の算出法は本実施形態と同様である。本実施形態の結果とOnishi et al. (2009)の結果はよく一致することから、本実施形態の信頼性が確かめられた。
【００７６】
続いて、衝突頻度因子について述べる。図８にストークス数Ｓｔ（＝τｐλ）が１の粒子の衝突頻度因子を示す。縦軸は無次元化された衝突頻度因子、横軸は乱流レイノルズ数Ｒｅλである。黒塗りのプロットはＤＮＳの結果を表している。丸のプロットは本実施形態のよるものであり、三角のプロットはOnishi et al. (2009)によるものであり、四角のプロットはWoittiezet al.(2009)（非特許文献１）によるものである。白抜きのプロットは大西ら(2007)（大西領，小森悟，“乱流中における同一径粒子間の衝突因子のモデル化”, 日本機械学会論文集B編，Vol. 73, No. 730 (2007), pp. 1307-1314.）のモデル予測結果を表している。そして、モデルから伸びる矢印は、モデルが予測しようとした値を指している。
【００７７】
既往研究で得られてきたＲｅλ＜８５という低レイノルズ数の場合には、モデルは数割の誤差はあるものの、レイノルズ数の増加と共に衝突頻度因子が小さくなる傾向はよく予測できる。しかし、本実施形態により初めて得られたＲｅλ＞８５高レイノルズ数の場合は、モデルの予測誤差が大きい。モデルはレイノルズ数の増加と共に衝突頻度因子は急激に小さくなると予測するのに対して、ＤＮＳの結果はそのような明確な減少傾向を示さない。これは、大西ら(2007)のモデルが低レイノルズ数の場合のデータを基にして得られた経験パラメータを使用しているためである。
【００７８】
例えば、Falkovich et al. (2002)（Falkovich, G., Fouxon, A. and Stepanov, M. G., “Acceleration of RainInitiation by Cloud Turbulence”, Natrure, Vol.419 (2002), pp.151-154）は、乱流の間欠性が引き起こす“パチンコ効果（sling effect）”が高レイノルズ数の時の衝突頻度を増大すると予測している。そのような効果が考慮されていないため、大西ら(2007)のモデルは高レイノルズ数の場合に衝突頻度を過小評価すると考えられる。
【００７９】
続いて、並列計算性能について述べる。大規模並列計算には、高い線形拡張性（scalability）が求められる。線形拡張性の評価には、強線形拡張性（strong scaling）と弱線形拡張性（weak scaling）の２種類が用いられる。前者を評価するためには、全体の計算規模を固定した上で、演算プロセス数（本実施形態におけるノード１０の数）を増やした時に演算時間がどれだけ減少するかを調べる。例えば、演算プロセス数を２倍にしたときに演算時間が半分になれば理想的な強線形拡張性を持つといえる。一方、後者を評価するためには、演算プロセスに対する計算規模を固定した上で、演算プロセス数を増やした時に演算時間がどれだけ変化するかを調べる。例えば、全体の計算規模と演算プロセス数とを共に２倍にしても、演算時間が変化しなければ高い弱線形拡張性を持つといえる。
【００８０】
図９に、本実施形態におけるＮ２５６計算に対して演算プロセス数を変化させた場合の（ａ）１時間ステップに要した演算時間［ｓ］の結果、及び（ｂ）演算速度［ＭＦＬＯＰＳ］の結果を示す。両図の横軸は共に演算プロセス数である。全体の計算規模（２５６３格子、２６２，１４４粒子）を固定した上で演算プロセス数を変化させたので、上述の強線形拡張性を調べた結果である。この性能測定にはＳＧＩ ＡＬＴＩＸ４７００システム（６．４ＧＦＬＯＰＳ／ｃｏｒｅ）が用いられた。図９（ａ）から、演算時間が演算プロセス数に反比例して減少していくことがわかる。それに対応して、図９（ｂ）では、演算速度がプロセス数にほぼ比例して増加していくことがわかる。
【００８１】
例えば、ＩＢＭ ＢｌｕｅＧｅｎｅシステムを使った等方性乱流場に対する大規模並列ＤＮＳでは、理論ピーク性能の約１０％程度の性能であったとの報告がある。本実施形態においては、粒子運動の計算及び衝突判定までを入れた上でも、流体計算のみの場合と同程度の高い性能を達成できた。
【００８２】
上述したように、本実施形態は、乱流場を有限差分法によって計算し、粒子運動をラグラジアン追跡法によって計算し、粒子間の衝突データを取得する高効率な大規模並列計算法である。本実施形態は、以下のような特徴を有している。
（ｉ）定常等方性乱流場を得るための強制法として、有限差分法の高い並列化率を損なわないＲＣＦ法を用いた。
（ｉｉ）粒子間衝突を効率的に検出するためにセルインデックス法を用いた。
（ｉｉｉ）共有メモリ型並列（自動並列化ライブラリ）と分散メモリ型並列（ＭＰＩライブラリ）の両方に対応し、両方を同時に用いるハイブリッド計算に対応した。
（ｉｖ）ＭＰＩ並列計算領域からはみ出した粒子情報の通信やセルインデックス情報の更新を所定の周期の時間ステップ毎にだけ行うようにコードを設計し、分散メモリ型並列計算でのプロセス間通信コストを削減した。
【００８３】
実際に、本実施形態に基づいた処理を実行することによって、Ｒｅλ＝２０９という高レイノルズ数における慣性粒子の衝突頻度データを得ることができた。これは従来の粒子衝突計算における最大レイノルズ数Ｒｅλ＝８５を大幅に更新する貴重なデータである。また、本実施形態で得られた高レイノルズ数における衝突頻度データを用いて、既往の衝突頻度モデルの検証も行った。その結果、既往の衝突頻度モデルは高レイノルズ数の時に衝突頻度を過小評価することが明らかになった。
【符号の説明】
【００８４】
１…粒子状態計算システム、１０…ノード、１１…粒子状態計算部、１２…周期判断部、１３…自仮想粒子情報出力部、１４…他仮想粒子情報入力部、１５…自仮想粒子状態出力部、１６…他仮想粒子状態入力部、１７…出力部。
【技術分野】
【０００１】
本発明は、分散メモリ型並列計算により粒子の状態を計算する粒子状態計算装置及び粒子状態計算方法に関する。
【背景技術】
【０００２】
雲の発達初期において雲粒子（水粒子）の粒子系分布が急激に広がり（粒子系分布の分散が急激に大きくなり）、速やかな降雨開始をもたらす現象が長い間議論されている。その原因として、雲乱流による雲粒子の衝突促進効果、乱流エントレインメント、巨大凝結核、雲粒子の乱流拡散等が挙げられている。中でも、初めに挙げた雲乱流による衝突促進効果が最も長い間議論されてきた。それに伴い、乱流中での粒子の衝突頻度に関する研究が盛んに行われてきた。
【０００３】
これまでに粒子の乱流衝突頻度を予測するモデルがいくつか提案されてきた。例えば、非特許文献１に記載されたような、直接数値計算（Direct Numerical Simulation，ＤＮＳ）によって求められた衝突頻度因子データを基にしたものが提案されている。
【先行技術文献】
【非特許文献】
【０００４】
【非特許文献１】Woittiez, E.J.P., Jonker, H.J.J. and Portela, L.M., “On the Combined Effects of Turbulence and Gravity onDroplet Collisions in Clouds: a Numerical Study”, Journal of the AtmosphericSciences, Vol. 66 (2009), pp.1926-1943.
【発明の概要】
【発明が解決しようとする課題】
【０００５】
しかしながら、上記の従来の方法で得られる予測結果は、以下のように高レイノルズ数での信頼性が確かめられていない。従来の研究で得られたデータは、非特許文献１に示されるように最大でも８５というテイラーマイクロスケール基準レイノルズ数Ｒｅλ（＝ｕ´ｌλ／ν、ここでｕ´は速度変動強度、ｌλはテイラーマイクロスケール、νは動粘度係数である）で得られたデータである。そのレイノルズ数は、実際の雲乱流のＲｅλ＝１０３〜４という高レイノルズ数とはほど遠い。そのため、従来の乱流衝突モデルをそのまま雲乱流のシミュレーションに適用することの妥当性には疑問が残る。
【０００６】
これを解決するためには、できるだけ高いレイノルズ数における衝突頻度因子データが必要である。そして、そのためには大規模な並列計算が必須となる。並列計算には、各プロセスが同じ物理メモリを共有する共有メモリ型並列計算（ｏｐｅｎＭＰや自動並列化ライブラリなど）と、各プロセスがそれぞれ別個の物理メモリを使用する分散メモリ型並列計算（Message Passing Interface（ＭＰＩ）ライブラリが一般的）の２種類に分けられる。
【０００７】
前者は、せいぜい数十プロセスを使った並列計算しか実行できず、また使用できるメモリも限られる。従って、大規模な並列計算を実行できないという短所がある。これまでに、乱流中での粒子衝突に対する共有メモリ型（ｏｐｅｎＭＰ等）の並列コードは開発されてきた。しかし、上記の理由により、共有メモリ型の並列計算では、高レイノルズ数における粒子衝突データを取得することは難しい。
【０００８】
一方、後者の分散メモリ型並列計算では、ハードウェアを寄せ集めれば、多くのプロセスを使用した並列計算が可能である。分散メモリ型並列計算では、それぞれのプロセスが（分散した）メモリを使用して、必要な時にデータ通信を行いながら計算を進める。データ通信速度はメモリアクセス速度に比べて極端に遅いため、効率的なデータ通信が求められる。
【０００９】
しかし、粒子衝突に対する分散メモリ型並列計算コードは開発されてこなかった。その一つの原因として、気相乱流計算に対して主にスペクトル法を用いたＤＮＳが行われてきたことが挙げられる。分散メモリ型並列環境でスペクトル計算する場合には、全領域データの（全対全の）通信が発生する。膨大なデータ通信を効率的に処理する技術が必要となるため、現状では３次元領域分割に対応したスペクトル法は開発されていない。２次元領域分割による大規模並列計算では、一つ一つの領域は細長い、いわゆるペンシル型の形状になり、（表面積）／（体積）が大きくなる。粒子計算のような近傍プロセス間で情報の受け渡しが必要となる並列計算の場合、通信量は領域の表面積に比例するため、（表面積）／（体積）の増大は（通信量）／（計算量）の増大につながり、結果として並列効率の低下を引き起こす。そのため２次元分割では、粒子運動に対する大規模並列計算に向かない。これらのことから、２次元分割にしか対応していないスペクトル法を用いて乱流を計算し、高レイノルズ数における粒子衝突データを取得するには大きな困難があると考えられる。
【００１０】
本発明は、上記を鑑みてなされたものであり、効率的に分散メモリ型並列計算により粒子状態を計算し、その結果より大規模なモデルに対する計算を可能とする粒子状態計算装置及び粒子状態計算方法を提供することを目的とする。
【課題を解決するための手段】
【００１１】
上記目的を達成するために、本発明に係る粒子状態計算装置は、複数の分割された計算領域の一つである自領域に配置される粒子の状態を時間ステップ毎に繰り返し計算する粒子状態計算装置であって、各時間ステップにおける自領域の粒子の状態を計算する粒子状態計算手段と、粒子状態計算手段における対象の時間ステップが、予め設定された周期の時間ステップか否かを判断する周期判断手段と、周期判断手段によって予め設定された周期の時間ステップであると判断された時間ステップに、粒子状態計算手段によって計算された粒子の位置に基づいて、当該粒子から、自領域に隣接する分割された計算領域である隣接領域において計算に用いられる自領域仮想粒子を特定して、特定した自領域仮想粒子を示す情報を当該隣接領域の粒子の状態を計算する別の粒子状態計算装置に出力する自仮想粒子情報出力手段と、別の粒子状態計算装置から、自領域において計算に用いられる他領域仮想粒子を特定する情報を入力する他仮想粒子情報入力手段と、各時間ステップにおいて粒子状態計算手段によって算出された自領域仮想粒子の状態を示す情報を別の粒子状態計算装置に出力する自仮想粒子状態出力手段と、各時間ステップにおいて別の粒子状態計算装置から他領域仮想粒子の状態を示す情報を入力する他仮想粒子状態入力手段と、粒子状態計算手段による計算に応じた情報を出力する出力手段と、を備え、粒子状態計算手段は、他仮想粒子情報入力手段によって入力される他領域仮想粒子を特定する情報、及び他仮想粒子状態入力手段によって入力された他領域仮想粒子の状態を示す情報を用いて粒子の状態を計算することを特徴とする。
【００１２】
本発明に係る粒子状態計算装置では、予め設定された周期の時間ステップにおいて自領域における自領域仮想粒子が特定されて、特定された自領域仮想粒子を示す情報が隣接領域の状態を計算する別の粒子状態計算装置に出力される。その一方で、別の粒子状態計算装置から、自領域において計算に用いられる他領域仮想粒子を特定する情報が入力される。そして、各時間ステップにおいて、自領域仮想粒子の状態を示す情報が出力されると共に、他領域仮想粒子の状態を示す情報が入力される。入力された他領域仮想粒子を特定する情報、及び当該粒子の状態を示す情報を用いて自領域の粒子の状態が計算される。即ち、粒子状態計算装置間でやり取りする情報に係る粒子を限定された仮想粒子とすることで、各粒子状態計算装置で必要な情報を確実に入力させると共に粒子状態計算装置間の通信を効率的にする。また、仮想粒子を特定する情報のやりとりを、予め設定された周期の時間ステップとすることで粒子状態計算装置間の通信が更に効率的になる。これによって、本発明に係る粒子状態計算装置によれば、効率的に分散メモリ型並列計算により粒子状態を計算し、その結果より大規模なモデルに対する計算が可能となる。
【００１３】
粒子は、流体中の粒子であり、粒子状態計算手段は、有限差分法によって流体の状態を算出して、当該流体の影響を考慮して粒子の状態を計算する、ことが望ましい。この構成によれば、流体中の粒子に係る大規模なモデル計算が可能となる。
【００１４】
粒子状態計算手段は、他仮想粒子情報入力手段によって入力される他領域仮想粒子を特定する情報に基づいて、粒子間の近傍関係を特定して、当該近傍関係に基づいて粒子の状態を計算することが望ましい。この構成によれば、粒子間の衝突等の粒子間の相互作用が生じるモデルの計算を確実に行うことができる。
【００１５】
ところで、本発明は、上記のように粒子状態計算装置の発明として記述できる他に、以下のように粒子状態計算方法の発明としても記述することができる。これはカテゴリが異なるだけで、実質的に同一の発明であり、同様の作用及び効果を奏する。
【００１６】
即ち、本発明に係る粒子状態計算方法は、複数の分割された計算領域の一つである自領域に配置される粒子の状態を時間ステップ毎に繰り返し計算する粒子状態計算装置による粒子状態計算方法であって、各時間ステップにおける自領域の粒子の状態を計算する粒子状態計算ステップと、粒子状態計算ステップにおける対象の時間ステップが、予め設定された周期の時間ステップか否かを判断する周期判断ステップと、周期判断ステップにおいて予め設定された周期の時間ステップであると判断された時間ステップに、粒子状態計算ステップにおいて計算された粒子の位置に基づいて、当該粒子から、自領域に隣接する分割された計算領域である隣接領域において計算に用いられる自領域仮想粒子を特定して、特定した自領域仮想粒子を示す情報を当該隣接領域の粒子の状態を計算する別の粒子状態計算装置に出力する自仮想粒子情報出力ステップと、別の粒子状態計算装置から、自領域において計算に用いられる他領域仮想粒子を特定する情報を入力する他仮想粒子情報入力ステップと、各時間ステップにおいて粒子状態計算ステップにおいて算出された自領域仮想粒子の状態を示す情報を別の粒子状態計算装置に出力する自仮想粒子状態出力ステップと、各時間ステップにおいて別の粒子状態計算装置から他領域仮想粒子の状態を示す情報を入力する他仮想粒子状態入力ステップと、粒子状態計算ステップにおける計算に応じた情報を出力する出力ステップと、を含み、粒子状態計算ステップにおいて、他仮想粒子情報入力ステップにおいて入力される他領域仮想粒子を特定する情報、及び他仮想粒子状態入力ステップにおいて入力された他領域仮想粒子の状態を示す情報を用いて粒子の状態を計算することを特徴とする。
【発明の効果】
【００１７】
本発明によれば、粒子状態計算装置間でやり取りする情報に係る粒子を限定された仮想粒子とすることで、各粒子状態計算装置で必要な情報を確実に入力させると共に粒子状態計算装置間の通信を効率的にする。また、仮想粒子を特定する情報のやりとりを、予め設定された周期の時間ステップとすることで粒子状態計算装置間の通信が更に効率的になる。これによって、本発明によれば、効率的に分散メモリ型並列計算により粒子状態を計算し、その結果より大規模なモデルに対する計算が可能となる。
【図面の簡単な説明】
【００１８】
【図１】本発明の実施形態に係る粒子状態計算装置であるノードの機能構成を示す図である。
【図２】流体及び粒子の計算領域及び（流体計算における）計算領域の分割を示す図である。
【図３】計算領域のセルへの分割を示す図である。
【図４】複数のノードにより構成される粒子状態計算システム全体で実行される処理を示すフローチャートである。
【図５】本発明の実施形態に係る粒子状態計算装置である各ノードで実行される処理（粒子状態計算方法）を示すフローチャートである。
【図６】図５のフローチャートの各処理に応じた計算領域における粒子を示す図であり、説明のために平面図としている。
【図７】本実施形態に係る粒子状態計算装置及び粒子状態計算方法により得られる粒子の衝突頻度因子の結果を示すグラフである。
【図８】本実施形態に係る粒子状態計算装置及び粒子状態計算方法により得られるストークス数が１の粒子の衝突頻度因子の結果を示すグラフである。
【図９】Ｎ２５６計算に対して演算プロセス数を変化させた場合の（ａ）１時間ステップに要した演算時間、（ｂ）演算速度を示すグラフである。
【発明を実施するための形態】
【００１９】
以下、図面と共に本発明に係る粒子状態計算装置及び粒子状態計算方法の好適な実施形態について詳細に説明する。なお、図面の説明においては同一要素には同一符号を付し、重複する説明を省略する。
【００２０】
図１に本実施形態に係る粒子状態計算装置であるノード１０を示す。ノード１０は、別のノード１０と共に粒子状態計算システム１を構成する。粒子状態計算システム１は、例えば、雲を構成する水粒子等の粒子の状態を計算（シミュレーション）するシステムである。各ノード１０は、ＣＰＵ（Central Processing Unit）、メモリ、ハードディスク等のハードウェアを備えるコンピュータである。各ノード１０はイントラネットやインターネット等のネットワークで接続されており、各ノード１０間で情報の送受信が可能になっている。粒子の状態の計算は、各ノード１０が別個の物理メモリを用いて分散して計算を行う分散メモリ型並列計算によって行われる。ここで各ノード１０における、計算の処理を（演算）プロセスと呼ぶ。
【００２１】
計算対象の粒子は、予め大きさが設定される計算領域に配置される。計算領域は複数の領域に分割される。各ノード１０は、分割された領域の一つが計算対象として割り当てられて、割り当てられた領域である自領域に配置される粒子の状態を時間ステップ毎に繰り返し計算する。計算される粒子の状態は、具体的には粒子の位置及び速度である。また、計算領域には流体（より具体的には、非圧縮流体）が流れており、粒子は流体中にあり流体の影響を受けることとしてもよい（即ち、流体についての計算も行う）。流体は、例えば、水粒子の状態に影響を与える気流である。計算される流体の状態は、具体的には流体の速度及び圧力である。
【００２２】
計算対象の領域が隣接しているノード１０同士は、粒子の移動等に応じて粒子の情報を送受信して、粒子の計算に利用する。
【００２３】
なお、粒子状態計算システム１には、分散メモリ型並列計算を統括的に制御する装置も含まれている。その装置が計算のためのパラメータの入力、計算領域の分割、及び各ノード１０への計算の指示等の分散メモリ型並列計算を統括的に制御する処理を行う。当該装置が複数のノード１０のうちの何れかであってもよい。
【００２４】
ここで、粒子及び流体の計算について説明する。まず、流体の支配方程式と数値計算法について説明する。非圧縮流体の方程式は、次の連続の式と、ナビエ−ストークス方程式である。
【数１】
ここで、ｕｉは流体の速度ベクトル、ｐは流体の圧力を示す変数である。ｔは時間ステップを示す変数である。添え字ｉ，ｊは、ｘ，ｙ，ｚ何れかの方向を表しており、式（１）及び（２）は総和規約に従って表記されている。全ての変数は、代表長さＬ０と代表速度Ｕ０とによって無次元化されている。また、Ｒｅは動粘性係数νを用いてＲｅ＝Ｕ０Ｌ０／νと定義されるパラメータとして与えられる値である。例えば、νを１気圧２９８Ｋでの値１．５×１０−５ｍ２／ｓに設定される。式（２）の右辺最終項ｆｉは外力項である。
【００２５】
対流項の差分の計算には、Morinishi, Y., Lund, T.,Vasilyev, O. and Moin, P., “Fully Conservative Higher Order Finite DifferenceSchemes for Incompressible Flow”, Journal of Computational Physics, Vol. 143(1998), pp.90-124. (4th-order scheme for Non-linear term)に示される保存型４次中心差分スキーム、粘性項の差分の計算には４次中心差分法を用いることができる。時間進行（時間発展）の計算には、２次のルンゲ−クッタ法を用い、時間発展後の速度が連続の式を満たすように速度と圧力とを同時緩和させながら反復計算するＨＳＭＡＣ法（Hirt, C.W. and Cook, J.L., “Calculating Three-Dimensional Flowsaround Structures and over Rough Terrain”, Journal of Computational Physics,Vol. 10 (1972), pp.324-340）を用いることができる。その反復計算の中での空間差分法には２次中心差分法を用いることができる。反復計算の収束判定条件は、例えば、∂ｉｕｉ／（Ｕ０／Δ）＜１０−３（Δは計算格子幅）とすることができる。
【００２６】
外力項ｆｉには、Reduced-CommunicationForcing（ＲＣＦ）法を用いることができる（Onishi, R, Baba, Y. andTakahashi, K., “Efficient Large-Scale Forcing in Finite-Difference Simulationsof Steady Isotropic Turbulence”, Topics on Chaotic Systems -selected papersfrom CHAOS 2010 international conference, Skiadas, C.H. et al. eds.,World Scientific Publishing Co. (in press)）。これは、大スケール運動にのみエネルギーを注入する大スケール強制法の一つであり、有限差分法において高並列効率を達成可能な強制法である。例えば、波数ｋの大きさが｜ｋ｜＜２．５の大スケール運動のエネルギーが保たれるように設定される。
【００２７】
続いて、流体の計算領域と並列計算法について説明する。図２に示すように、一辺２πＬ０の立方体計算領域に対して、Ｎ３個のスタッガード格子を配する。つまり、計算格子幅Δ＝２πＬ０／Ｎとなる。ここで、Ｎは予め設定あるいは入力されるパラメータである。具体的な計算コードは、ＭＰＩライブラリを用いたＦｏｒｔｒａｎ９０コードで記述することができる。これは、３次元領域分割に対応しており、立方体計算領域の座標軸であるｘ、ｙ及びｚの各方向を等間隔に分割することができる。各方向の分割数Ｍｘ、Ｍｙ及びＭｚは、ＭＰＩプロセス総数ＭｐｒｏｃとＭｐｒｏｃ＝ＭｘＭｙＭｚの関係にある。このとき、各プロセスはｎｘ×ｎｙ×ｎｚ（ｎｘ＝Ｎ／Ｍｘ、ｎｙ＝Ｎ／Ｍｙ及びｎｚ＝Ｎ／Ｍｚ）の大きさを持った“ＭＰＩ分割領域”内の流体及び粒子運動の計算を担当する。Ｍｘ、Ｍｙ及びＭｚは、任意のＮの約数の値に予め設定あるいは入力されるパラメータである。なお、上記の例では、Ｎ３に等間隔に分割することとしているが、Ｎｘ，Ｎｙ，Ｎｚと各方向に分割数を変えることもできる。
【００２８】
例えば、以下の表に示す計算条件で気相乱流場の計算が行われる。全てのケースで、クーラン数（Ｕ０Δｔ／Δ）が０．２以下になるように時間刻みΔｔを設定した。ここで、時間刻みΔｔは繰り返し計算の時間ステップの間隔である。但し、粒子の大きさによっては、口述するようにΔｔを更に小さく設定しなければならない場合もある。
【表１】
【００２９】
引き続いて、粒子の状態の計算について説明する。粒子の運動は、例えば、ラグラジアン粒子追跡法を用いて計算することができる。線形ストークス抗力を仮定すると支配方程式は次式となる。
【数２】
ここで、ｖｐは粒子速度、ｕｉ（アンダーバーの記載は省略する）は粒子位置における流体速度を示す変数である。ここで、添え字ｉはｘ，ｙ，ｚ方向を表す。τｐは粒子緩和時間である。式（１）と同様に全ての変数は、代表長さＬ０と代表速度Ｕ０とによって無次元化されている。また、粒子に働く重力は無視する。粒子緩和時間τｐは、次式で計算される。
【数３】
ここで、ｒは粒子半径、ρｐ及びρｆはそれぞれ粒子及び流体の密度であり、予め設定あるいは入力されるパラメータである。例えば、雲粒を想定して１気圧２９８Ｋでの水及び空気の密度の値から、ρｐ／ρｆ＝８４０とする。線形ストークス抗力モデルは粒子レイノルズ数Ｒｅｐ（＝２ｒ｜ｖｐ−ｕ｜／ν）が１に比べて十分小さい場合に適用可能である。雲粒を想定した場合、およそｒ＝４０μｍのときＲｅｐ＝１である。それ以上大きな粒子の場合には非線形抗力モデルを用いるべきであるが、簡単のために全ての場合で線形ストークス抗力モデルを用いてもよい。また、現実の雲では水滴の空気に対する質量割合最大でも１０−３程度であるので、粒子は希薄であり流体運動に影響を与えないとしてもよい。粒子位置における流体速度ｕは、粒子近傍の８点の速度データから線形補間により算出される。時間進行の計算には、流体計算と同じく、２次のルンゲ−クッタ法を用いる。
【００３０】
続いて、粒子の衝突の計算について説明する。衝突は２粒子間衝突のみを考慮する。２粒子がΔｔ間に移動する軌跡を求めて、それらの軌跡の最小距離が粒子直径以下であるか否かを判断して、粒子直径以下であれば衝突とみなす。衝突後には一方の粒子を消滅させる。従って、一回衝突が起こると粒子の総数が１個減少する。これを利用して、粒子総数の減少率から衝突頻度Ｎｃ及び衝突頻度因子Ｋｃ（＝２Ｎｃ／ｎｐ２、ｎｐは粒子数密度）を算出する。ある時間ステップｎにおける粒子総数Ｎｐｎ、１ステップ前の粒子総数Ｎｐｎ−１及び計算領域の体積Ｖｄ（＝（２πＬ０）３）から、ｎ−１／２ステップにおける衝突頻度Ｎｃｎ−１／２は以下のように算出される。
【数４】
この衝突頻度とｎ−１／２ステップにおける粒子数密度ｎｐｎ−１／２（＝０．５（Ｎｐｎ−１＋Ｎｐｎ）／Ｖｄ）から、ｎ−１／２ステップにおける衝突頻度因子を次の式により算出する。
【数５】
この衝突頻度因子を時間平均することによって平均衝突因子を算出する。
【００３１】
計算領域内にＮｐ個の粒子が存在する場合、特殊な方法を使わなければ、Ｎｐ（Ｎｐ−１）／２〜０．５Ｎｐ２回の衝突判定計算をしなければならない。衝突判定回数を削減する方法として分子動力学法で使用されるセルインデックス法（Allen, M.P., and Tildesley, D.J., “Computer Simulation of Liquids”,Oxford University Press, (1987), p.408）がある。これは、図３に示すように衝突頻度領域をセルに小分けし、ある粒子を含むセルとそのセルとを囲むｘ，ｙ，ｚ各方向の前後、合計２７個のセルの中の粒子との衝突を判定する。但し、逐次的に衝突を判定していく場合には、２７個ではなく、８個のセル（例えば、対象粒子を含むセルとｘ，ｙ，ｚ各方向前方側のセル）を判定の対象とすれば、重複判定がなく効率がよい。このセルインデックス法では、粒子一つ一つに番号を付けた上で、何番と何番の粒子がどのセル内にあるかを管理するリスト（セルインデックス）を作成する。
【００３２】
図３に示すように計算領域をＭｃｅｌｌ３個のセルに分割した場合、各セル内の平均粒子数Ｎｐ，ｃｅｌｌはＮｐ，ｃｅｌｌ＝Ｎｐ／Ｍｃｅｌｌ３となる。ある粒子一つに対する衝突判定回数は、粒子を含むセル内の他粒子との衝突判定回数（Ｎｐ，ｃｅｌｌ−１）／２と、ｘ，ｙ，ｚ各方向後方側のセル内粒子との衝突判定回数７Ｎｐ，ｃｅｌｌの和、つまり（Ｎｐ，ｃｅｌｌ−１）／２＋７Ｎｐ，ｃｅｌｌと見積もられる。結局、衝突判定回数の総数は、Ｎｐ×（（Ｎｐ，ｃｅｌｌ−１）／２＋７Ｎｐ，ｃｅｌｌ）〜７．５Ｎｐ２／Ｍｃｅｌｌ３である。つまり、Ｍｃｅｌｌ＞２．５に設定すれば衝突判定回数の回数を削減できる。この見積もりでは、Ｍｃｅｌｌを大きくすればするほど衝突判定回数の回数は小さくなる。しかし、セルの中に粒子がほとんどないような状況では、セルの中に粒子があるかないかを調べる処理が無駄になる。そのため、Ｎｐ，ｃｅｌｌが１以上になるようにＭｃｅｌｌを設定するのが望ましい。
【００３３】
本実施形態では、有限差分法を用いて気相乱流を計算し、乱流中で運動する粒子同士の衝突を検出する。ＭＰＩライブラリを用いた分散メモリ型並列計算に対応させるためには、下記の３つの対応するＭＰＩ通信の機能を実現しなければならない。即ち、当該機能に対応するＭＰＩ通信プログラムを記述しなければならない。
（Ａ）差分計算に必要な袖領域の気相データの取得
（Ｂ）ＭＰＩ分割領域（個々のノード１０に割り当てられた計算領域）からはみ出した粒子を隣のプロセスに移動させる処理
（Ｃ）ＭＰＩ分割領域の境界近傍にある粒子情報の共有
ある粒子ペアがＭＰＩ分割領域の境界をまたいで近接する場合、その情報を両側のプロセスが共有していなければ衝突を検出できない可能性がある。この問題に対処するため上記の（Ｃ）の操作が必要になる。本発明は、セルインデックスを活用して、（Ｂ）（Ｃ）を効率的に処理するためものである。
【００３４】
引き続いて、ノード１０の本発明の係る機能を詳細に説明する。図１に示すようにノード１０は、粒子状態計算部１１と、周期判断部１２と、自仮想粒子情報出力部１３と、他仮想粒子情報入力部１４と、自仮想粒子状態出力部１５と、他仮想粒子状態入力部１６と、出力部１７とを備えて構成される。なお、各構成要素の更に具体的な機能、各構成要素によって行われる更に具体的な処理については、後述する処理フローの説明で説明する。これらの機能は、例えばノード１０においてプログラムが実行されることによって実現される。
【００３５】
粒子状態計算部１１は、各時間ステップにおける自領域の粒子の状態を計算する粒子状態計算手段である。自領域とは、自ノード１０に割り当てられた計算領域である。粒子状態計算部１１は、具体的には上述した方法によって粒子の状態を計算する。なお、計算に必要なパラメータについては、上述した分散メモリ型並列計算を統括的に制御する装置から入力される。また、計算に必要な計算式等は予め粒子状態計算部１１に記憶されており、それらが読み出されて計算が行われる（他の構成要素についても同様）。
【００３６】
粒子状態計算部１１は、粒子間の相互作用（本実施形態では、粒子同士の衝突）の判断、計算のために上述したようにセルインデックス法に応じた処理を行う。即ち、各粒子に粒子を特定する情報である番号を付与する。粒子状態計算部１１は、自領域にある粒子の番号を、当該粒子が含まれるセルを特定する情報に対応付けて記憶、管理する。セルは、上述したように自領域を更に分割した領域である。
【００３７】
粒子状態計算部１１は、後述するように、他仮想粒子情報入力部１４によって入力される他領域仮想粒子を特定する情報、及び他仮想粒子状態入力部１６によって入力された他領域仮想粒子の状態を示す情報を用いて粒子の状態を計算する。仮想粒子とは、自領域（自ノード１０）における状態の計算対象の粒子ではあるが、自領域に隣接する分割された計算領域である隣接領域（別のノード１０）における粒子の状態の計算に必要な粒子である。各ノード１０間で、仮想粒子の情報を送受信して、その情報を粒子の計算に利用する。自領域にある粒子のうち、隣接領域において計算に用いられる仮想粒子を自領域仮想粒子と呼ぶ。隣接領域にある粒子のうち、自領域において計算に用いられる仮想粒子を他領域仮想粒子と呼ぶ。
【００３８】
粒子状態計算部１１は、隣接領域（別のノード１０）との間で、自領域からはみ出して隣接領域に入った粒子の情報及び隣接領域からはみ出して自領域に入った粒子の情報の送受信も行う。
【００３９】
また、粒子状態計算部１１は、上述したように、有限差分法によって流体の状態を算出して、当該流体の影響を考慮して粒子の状態を計算する。
【００４０】
粒子状態計算部１１は、他仮想粒子情報入力部１４によって入力される他領域仮想粒子を特定する情報に基づいて、粒子間の近傍関係を特定して、当該近傍関係に基づいて粒子の状態を計算する。本実施形態では近傍関係とは衝突が判断される粒子のペアである。具体的には、後述するようにセルインデックスを用いる。
【００４１】
周期判断部１２は、粒子状態計算部１１における計算の時間ステップが、予め設定された周期の時間ステップか否かを判断する周期判断手段である。この周期を示す値は２以上の自然数であり、予め周期判断部１２に記憶されている。周期判断部１２は、例えば、粒子状態計算部１１によって管理される現時点での時間ステップの情報を参照することによって上記の判断を行う。周期判断部１２は、計算の時間ステップが予め設定された周期の時間ステップであると判断するとその旨を粒子状態計算部１１及び自仮想粒子情報出力部１３に通知する。
【００４２】
自仮想粒子情報出力部１３は、周期判断部１２によって予め設定された周期の時間ステップであると判断された時間ステップに、自領域仮想粒子を特定して、特定した自領域仮想粒子を示す情報を隣接領域の粒子の状態を計算する別のノード１０に出力する自仮想粒子情報出力手段である。自領域仮想粒子の特定は、粒子状態計算部１１によって計算された粒子の位置に基づいて行われる。自領域中のセルのうち、隣接領域に隣接するセルに位置する粒子が自領域仮想粒子とされる。自仮想粒子情報出力部１３は、隣接領域毎に自領域仮想粒子の特定し、自領域仮想粒子の番号を、隣接領域の粒子の状態を計算する別のノード１０に出力する。また、自仮想粒子情報出力部１３は、自領域仮想粒子の番号を自仮想粒子状態出力部１５に通知する。
【００４３】
他仮想粒子情報入力部１４は、隣接領域の粒子の状態を計算する別のノード１０から、他領域仮想粒子を特定する情報である他領域仮想粒子の番号を入力する他仮想粒子情報入力手段である。他仮想粒子情報入力部１４は、入力した他領域仮想粒子の番号を粒子状態計算部１１に通知する。
【００４４】
自仮想粒子状態出力部１５は、各時間ステップにおいて粒子状態計算部１１によって算出された自領域仮想粒子の状態を示す情報を別のノード１０に出力する自仮想粒子状態出力手段である。自仮想粒子状態出力部１５は、自仮想粒子情報出力部１３から通知された自領域仮想粒子の番号によって示される粒子の状態の情報を粒子状態計算部１１から取得して別のノード１０に出力する。出力対象のノード１０は、自仮想粒子情報出力部１３が番号の情報を出力したノード１０である。出力する粒子の状態の情報は、粒子の位置（座標）を示す情報と、粒子の速度（方向、大きさ）を示す情報である。
【００４５】
他仮想粒子状態入力部１６は、各時間ステップにおいて別のノード１０から他領域仮想粒子の状態を示す情報を入力する他仮想粒子状態入力手段である。自仮想粒子状態出力部１５は、入力した他領域仮想粒子の状態を示す情報を粒子状態計算部１１に出力する。粒子の情報を入力する他領域仮想粒子は、他仮想粒子情報入力部１４によって特定する情報が入力された他領域仮想粒子である。他仮想粒子状態入力部１６は、入力した他領域仮想粒子の状態を示す情報を粒子状態計算部１１に出力する。
【００４６】
出力部１７は、粒子状態計算部１１による計算に応じた情報を出力する出力手段である。具体的には例えば、出力部１７は、粒子状態計算部１１による計算により得られる、上述した平均衝突頻度因子を計算するための情報を、分散メモリ型並列計算を統括的に制御する装置に出力する。
【００４７】
引き続いて、図４及び図５のフローチャートを用いて、本実施形態に係るノード１０で実行される処理（粒子状態計算方法）を説明する。まず、図４のフローチャートを用いて、粒子状態計算システム１全体の処理の概要について説明し、続いて、図５のフローチャートを用いて各ノード１０による本発明に係る処理について説明する。
【００４８】
粒子状態計算システム１では、まず、分散メモリ型並列計算を統括的に制御する装置に計算用のパラメータが入力される（Ｓ０１）。このパラメータは、計算領域の大きさ、流体計算のための格子数、粒子計算のためのセル数、粒子数、各粒子の大きさ、密度、速度、分散メモリ型並列計算の分割数及び時間刻み（時間ステップの間隔）等の情報である。この入力は、例えば、ユーザによる当該装置の操作により行われる。
【００４９】
分散メモリ型並列計算を統括的に制御する装置では、入力されたパラメータに基づいて、各ノード１０が担当する計算領域を決定して、各ノードに対して計算に必要な情報を出力する（Ｓ０２）。各ノード１０では計算が行われる（Ｓ０３）。続いて、各ノード１０の計算結果に応じた情報（例えば、平均衝突頻度因子）が出力される（Ｓ０４）。結果の出力は、例えば、各ノード１０による計算結果が、分散メモリ型並列計算を統括的に制御する装置に集められて、当該装置によって最終結果が生成されて出力される。
【００５０】
引き続いて、図５のフローチャートを用いて、各ノード１０における処理（図４のフローチャートにおけるＳ０３に相当）について説明する。また、図６に各処理に対応する計算領域を示す。図６に示す図では、計算領域は２次元領域としており、あるプロセスが受け持つＭｃｅｌｌ，ｘ×Ｍｃｅｌｌ，ｙ個のセル（Ｍｃｅｌｌ，ｘ＝Ｍｃｅｌｌ／Ｍｘ及びＭｃｅｌｌ，ｙ＝Ｍｃｅｌｌ／Ｍｙ）が実線で描かれている。実線で描かれた領域は対象プロセスが担当するＭＰＩ分割領域（自領域）でもあり、点線部分は別のプロセス（近傍プロセス）が計算を担当している領域を示している。以下では、セルインデックス作成時にＭＰＩ分割領域内にある粒子を実粒子と呼ぶ。衝突判定のためには、Ｍｃｅｌｌ，ｘ＋１やＭｃｅｌｌ，ｙ＋１の領域にある粒子も考慮する必要がある。それらの粒子が、仮想粒子である。仮想粒子は衝突検出のためだけに用いられて、その運動は計算されない。図６において黒塗りの点は実粒子であり、白抜きの点は（他）仮想粒子である。
【００５１】
本処理では、まず、粒子状態計算部１１によって所定の時間ステップにおける流体の状態の計算が行われる（Ｓ１１、粒子状態計算ステップ）。計算される流体の状態は流体の速度及び圧力である。また、粒子状態計算部１１によって所定の時間ステップにおける実粒子の状態の計算が行われる（Ｓ１２、粒子状態計算ステップ）。計算される実粒子の状態は実粒子の位置及び速度である。これらの計算は、前の時間ステップ又は初期状態における流体及び実粒子の状態、並びに計算後の流体の状態に基づいて、上述したようにルンゲ−クッタ法により行われる。なお、流体及び実粒子の状態の計算は、必ずしも順番に行われるわけではなく、流体及び実粒子それぞれについてルンゲ−クッタ法の１次ステップ及び２次ステップ等のように交互あるいは同時に計算される。図６（ａ）に示す状態が前の時間ステップ又は初期状態（即ち、Ｓ１１の直前）であり、図６（ｂ）に示す状態が計算後の状態（即ち、Ｓ１２の直後）である。
【００５２】
Ｓ１２の後、図６（ｂ）に示すようにＭＰＩ分割計算領域からはみ出す（実）粒子が現れる。はみ出した粒子の処理（上述した（Ｂ）の処理）を毎時間ステップ実行すると、多大なノード１０間の通信が発生して、並列性能を低下させるおそれがある。そこで、以下のように予め設定された周期のみにはみ出した粒子の処理を行わせる。はみ出した実粒子の計算を続けるためにＭＰＩ分割計算領域外の流体速度データが必要になる。各プロセス（粒子状態計算部１１）は、上述（Ａ）の処理に伴って、４次中心差分計算のために４格子分の袖データを持つ。そのため、粒子が３．５Δ分（スタッガード格子を採用しているため４Δではなく３．５Δ）以上はみ出さない限りは、余分なＭＰＩ通信（ノード１０間の通信）をすることなくｕを算出することができる。つまり、流体計算を有限差分法によって計算するからこそ可能な効率化手法である。また、以下に示すようにセルインデックス作成も予め設定された周期毎に行う。
【００５３】
Ｓ１２の後、周期判断部１２によって、現時点の時間ステップが予め設定された周期の時間ステップか否かが判断される（Ｓ１３、周期判断ステップ）。現時点の時間ステップが予め設定された周期の時間ステップであると判断されると（Ｓ１３のＹｅｓ）、周期判断部１２から粒子状態計算部１１及び自仮想粒子情報出力部１３にその旨が通知され、その場合の処理（Ｓ１４〜Ｓ１９）が行われる。
【００５４】
まず、図６（ｃ）に示すように、粒子状態計算部１１において記憶されている他仮想粒子に関する情報が除去される（Ｓ１４）。続いて、粒子状態計算部１１によって自領域から出た実粒子の情報が入出力される（Ｓ１５、粒子状態計算ステップ）。具体的には、自領域から出た実粒子の情報が、当該実粒子が位置する領域の計算を担当するノードに出力され、自領域における実粒子から削除される。この実粒子は、他領域における実粒子となる。一方、他領域から出て自領域に入った粒子の情報が、他領域の計算を担当するノード１０から受け取られ、自領域の実粒子とされる（図６（ｄ）の矢印で示される粒子である）。即ち、実粒子のノード（領域）間での引き渡しが行われる。ここでノード１０間の通信（ＭＰＩ通信）が発生する。
【００５５】
続いて、粒子状態計算部１１によってセルインデックスが生成（更新）される。セルインデックスは、セルに含まれる（セル上に位置する）実粒子の番号をそのセルを特定する情報に対応付けた情報である（Ｓ１６、粒子状態計算ステップ）。図６（ｅ）に示す（２，２）のセルの実粒子のリストは、［４，１０２，６０５，６５１］である。セルインデックスは、粒子状態計算部１１によって記憶され管理される。なお、（２，２）のセルとは、ｘ軸において２番目、ｙ軸において２番目のセルであることを示す（以下同様）。
【００５６】
続いて、自仮想粒子情報出力部１３によって、自領域仮想粒子が特定されて、自領域仮想粒子を示す情報が隣接領域の粒子の状態を計算する別のノード１０に出力される（Ｓ１７、自仮想粒子情報出力ステップ）。具体的には、粒子状態計算部１１によってセルインデックスにおいて、そのセルに含むセルが自領域仮想粒子であるとされるセルに対応付けられた実粒子が自領域仮想粒子とされる。そのようなセルは、領域間の位置関係によって決定される。当該セルを特定する情報と自領域仮想粒子の番号とが自仮想粒子情報出力部１３から別のノード１０に出力される。また、自領域仮想粒子の番号は、自仮想粒子情報出力部１３から自仮想粒子状態出力部１５に通知される。
【００５７】
その一方で、別のノード１０から送信される他領域仮想粒子の番号と当該粒子を含むセルとが対応付けられた情報が、他仮想粒子情報入力部１４によって入力される（Ｓ１８、他仮想粒子情報入力ステップ）。この情報は、別のノード１０におけるＳ１７に相当する処理によって送信されたものである。上記の情報は、他仮想粒子情報入力部１４から粒子状態計算部１１に入力される。これにより、図６（ｆ）に示すように、粒子状態計算部１１によって自領域に隣接するセルに位置する他仮想粒子の番号が把握される。Ｓ１７及びＳ１８ではノード１０間の通信（ＭＰＩ通信）が発生する。
【００５８】
続いて、粒子状態計算部１１によって、自領域の各粒子について、衝突の有無をチェックする粒子のペアが生成される（Ｓ１９、粒子状態計算ステップ）。この粒子のペアは、対象粒子と当該対象粒子を含むセル及び当該セルに隣接しているセルに位置する粒子とのペアである。このペアの生成の際には、他仮想粒子も考慮される（実粒子と他仮想粒子のペア）。
【００５９】
（ｉ，ｊ）番目のセル（（ｉ，ｊ）セル）は、自分自身を含めて、（ｉ，ｊ）、（ｉ，ｊ±１）、（ｉ±１，ｊ）、（ｉ±１，ｊ±１）の計９個のセルに囲まれている（３次元の場合は２７個）。しかし、セル番号を順に増やして衝突をチェックする場合には、（ｉ，ｊ）、（ｉ，ｊ＋１）、（ｉ＋１，ｊ）、（ｉ＋１，ｊ＋１）の計４個（３次元の場合は８個）のセルを調べていけば、重複することなく全ての組み合わせをチェックできる。そのため粒子のペアのリストを作る際には、ｉ＝０やｊ＝０のセルを含める必要はない。図６（ｆ）に示す番号が１３０のペアとなるリストは、上記の４つのセルに含まれる粒子である［３３２，６９，８７１，７７７，１０９，４，６０５，１０２，６５１］である。
【００６０】
Ｓ１９の後、あるいはＳ１２の後、周期判断部１２によって現時点の時間ステップが予め設定された周期の時間ステップでないと判断された場合（Ｓ１３のＮｏ）、以下の処理が行われる。
【００６１】
まず、自仮想粒子状態出力部１５によって、自仮想粒子情報出力部１３から通知された自領域仮想粒子のその時点の時間ステップの状態を示す情報が取得される。取得される情報は、Ｓ１２において粒子状態計算部１１よって計算されたものである。取得された情報は、自仮想粒子状態出力部１５から上記の別のノード１０に出力される（Ｓ２０、自仮想粒子状態出力ステップ）。
【００６２】
その一方で、別のノード１０から送信される他領域仮想粒子の状態を示す情報が、他仮想粒子状態入力部１６によって入力される（Ｓ２１、他仮想粒子状態入力ステップ）。この情報は、別のノード１０におけるＳ２０に相当する処理によって送信されたものである。上記の情報は、他仮想粒子状態入力部１６から粒子状態計算部１１に入力される。これにより、図６（ｇ）に示すように、粒子状態計算部１１によって自領域に隣接するセルに位置する他仮想粒子の状態である位置及び速度が把握される。Ｓ２１及びＳ２２ではノード１０間の通信（ＭＰＩ通信）が発生する。
【００６３】
続いて、粒子状態計算部１１では、図６（ｈ）に示すように、上述したように粒子間の衝突の判断が行われて、その判断に基づく計算が行われる（Ｓ２２、粒子状態計算ステップ）。この判断は、Ｓ１９で生成された全てのペアの間で行われる。実粒子間の衝突の判断は、それぞれＳ１２で計算された粒子の状態に基づいて行われる。実粒子と他領域仮想粒子との衝突の判断は、実粒子についてはＳ１２で計算された粒子の状態に基づいて、他領域仮想粒子についてはＳ２１で入力された粒子の状態に基づいて行われる。
【００６４】
続いて、粒子状態計算部１１では、計算の終了の判断が行われる（Ｓ２３、粒子状態計算ステップ）。この判断は、粒子状態計算部１１に予め設定あるいは入力された判断基準に基づいて行われる。計算が終了しないと判断された場合（Ｓ２３のＮｏ）には、次の時間ステップに進められてＳ１１からの処理が繰り返される（Ｓ２４、粒子状態計算ステップ）。
【００６５】
計算が終了する判断された場合（Ｓ２３のＹｅｓ）には、出力部１７によって、粒子状態計算部１１による計算に応じた情報が出力される（Ｓ２５、出力ステップ）。以上で、ノード１０における処理が終了する。なお、上記の処理における時間ステップの更新（Ｓ２４）の後、毎ステップあるいは数ステップに一度、例えば、粒子の衝突、粒子の速度及び粒子のエネルギー等の統計量を算出して出力することとしてもよい。
【００６６】
上述した処理では、一度セルインデックスが生成された後、予め設定された周期の時間ステップを経過するまで更新されないので、その間に粒子は当初のセルから移動する可能性がある。しかし、上記の周期の時間ステップの間の粒子の移動距離がセル幅Δｃｅｌｌ（＝２πＬ０／Ｍｃｅｌｌ）よりも小さければ、衝突する可能性のある粒子ペアを漏らす可能性を無視できる。本実施形態では、例えば、Δｃｅｌｌ＝４Δ（つまり、Ｍｃｅｌｌ＝Ｎ／４）のとき、上記の周期を８に設定し、クーラン数を０．２以下に設定すると、粒子は流体速度程度で移動すると考えると上記の周期の時間ステップの間に粒子はせいぜい１．６Δ程度しか移動しない。つまり、この設定では、衝突検出の漏れは無視できるほど小さい。実際、いくつかのケースで上記の周期を１に設定した場合と、８に設定した場合とで比較して、衝突ペアに変化が無いことを確認した。
【００６７】
上述したように本実施形態では、ノード１０間でやり取りする情報に係る粒子を限定された仮想粒子とすることで、各ノード１０で必要な情報を確実に入力させると共にノード１０間の通信を効率的にする。また、仮想粒子を特定する情報のやりとりを、予め設定された周期の時間ステップとすることでノード１０間の通信が更に効率的になる。これによって、本発明に係る本実施形態によれば、効率的に分散メモリ型並列計算により粒子状態を計算し、その結果より大規模なモデルに対する計算が可能となる。
【００６８】
また、本実施形態のように流体中の粒子を計算することとすれば、流体中の粒子に係る大規模なモデル計算が可能となる。但し、必ずしも流体中の粒子の計算でなくてもよく、分散メモリ型並列計算により粒子状態を計算するものに対して適用可能である。
【００６９】
また、本実施形態のように、セルインデックスを用いて粒子間のペアの相互作用を判断することとすれば、粒子間の衝突等の粒子間の相互作用が生じるモデルの計算を確実に行うことができる。
【００７０】
なお、上述した実施形態における粒子の計算については、定常な気相乱流場が形成されてから、粒子を混入させて行うこととしてもよい。また、衝突の判定については、粒子の混入から一定時間経過後（例えば、Ｔ０（＝Ｌ０／Ｕ０）以上）としてもよい。
【００７１】
なお、本実施形態においては雲をシミュレーションする例を示したが、必ずしも雲を対象とした計算に限られず、粒子の状態を計算するものであれば適用が可能である。例えば、小惑星の生成シミュレーション等にも適用可能である。
【００７２】
以下に、上述した実施形態による具体的な計算結果を示す。なお、定常な気相乱流相場が形成されたことを確認後、粒子を混入させた。粒子運動を計算する場合、時間刻みΔｔを０．５τｐ以下に設定した。この条件と流体計算のクーラン数による条件を同時に満たすΔｔを流体計算と粒子計算との両方に用いた。また、粒子混入からＴ０（＝Ｌ０／Ｕ０）以上経過後から衝突判定を開始した。
【００７３】
以下の表に粒子の計算条件を示す。最左列に示されているケース名は、上記の表と対応している。全てのケースで粒子数密度ｎｐが雲の代表値である１０８オーダーになるように設定した。セルに含まれる平均粒子数はＮ６４及びＮ２５６の場合は１、Ｎ５１２の場合は４．８であった。
【表２】
【００７４】
以下の表に、統計的に定常状態に達した気相等方性乱流場の乱流特性量を示す。表中のｕ´は速度変動強度、Ｒｅλ（＝ｕ´ｌλ／ν、ｌλはテイラーマイクロスケール）はテイラーマイクロスケール基準乱流レイノルズ数、ｋｍａｘ（＝Ｎ／２）は最大有効波数、ｌη（＝（ν３／ε）１／４、εはエネルギー離散率）はコルモゴロフスケールである。ｋｍａｘｌηは計算解像度の指標であり、運動エネルギー等の２次元量を議論する場合には１程度、スキューネスやフラットネス等の高次元を議論する場合には２程度必要だと言われている。従来の方法においては粒子衝突の計算では１．３〜１．４程度に設定されることが多いが、本例では通常よりも高い解像度である２程度になるように設定した。
【表３】
【００７５】
図７にＮ６４の場合（Ｒｅλ＝５２．３）の衝突頻度因子の結果を示す。縦軸は局所速度勾配λ（＝（ε／ν）１／２）と衝突半径Ｒ（＝２ｒ）で無次元化された衝突頻度因子、横軸はコルモゴロフスケールで無次元化された粒子半径である。衝突によって粒子数が５％減少するまでを一つのランとして、それぞれのランにおける平均衝突頻度因子を算出した。比較のために、Onishi et al. (2009)（Onishi, R., Takahashi,K. and Komori, S., “Influence of Gravity on Collisions of MonodispersedDroplets in Homogeneous Isotropic Turbulence”, Physics of Fluids, Vol. 21(2009), 125108）で得られたＲｅλ＝５４．３の場合の結果も示す。Onishiet al. (2009)は、流体計算に擬スペクトル法を用いた。その点は、有限差分法を用いる本実施形態と異なるが、計算領域の大きさ、物性値、計算格子数、粒子数等の設定、そして標準偏差の算出法は本実施形態と同様である。本実施形態の結果とOnishi et al. (2009)の結果はよく一致することから、本実施形態の信頼性が確かめられた。
【００７６】
続いて、衝突頻度因子について述べる。図８にストークス数Ｓｔ（＝τｐλ）が１の粒子の衝突頻度因子を示す。縦軸は無次元化された衝突頻度因子、横軸は乱流レイノルズ数Ｒｅλである。黒塗りのプロットはＤＮＳの結果を表している。丸のプロットは本実施形態のよるものであり、三角のプロットはOnishi et al. (2009)によるものであり、四角のプロットはWoittiezet al.(2009)（非特許文献１）によるものである。白抜きのプロットは大西ら(2007)（大西領，小森悟，“乱流中における同一径粒子間の衝突因子のモデル化”, 日本機械学会論文集B編，Vol. 73, No. 730 (2007), pp. 1307-1314.）のモデル予測結果を表している。そして、モデルから伸びる矢印は、モデルが予測しようとした値を指している。
【００７７】
既往研究で得られてきたＲｅλ＜８５という低レイノルズ数の場合には、モデルは数割の誤差はあるものの、レイノルズ数の増加と共に衝突頻度因子が小さくなる傾向はよく予測できる。しかし、本実施形態により初めて得られたＲｅλ＞８５高レイノルズ数の場合は、モデルの予測誤差が大きい。モデルはレイノルズ数の増加と共に衝突頻度因子は急激に小さくなると予測するのに対して、ＤＮＳの結果はそのような明確な減少傾向を示さない。これは、大西ら(2007)のモデルが低レイノルズ数の場合のデータを基にして得られた経験パラメータを使用しているためである。
【００７８】
例えば、Falkovich et al. (2002)（Falkovich, G., Fouxon, A. and Stepanov, M. G., “Acceleration of RainInitiation by Cloud Turbulence”, Natrure, Vol.419 (2002), pp.151-154）は、乱流の間欠性が引き起こす“パチンコ効果（sling effect）”が高レイノルズ数の時の衝突頻度を増大すると予測している。そのような効果が考慮されていないため、大西ら(2007)のモデルは高レイノルズ数の場合に衝突頻度を過小評価すると考えられる。
【００７９】
続いて、並列計算性能について述べる。大規模並列計算には、高い線形拡張性（scalability）が求められる。線形拡張性の評価には、強線形拡張性（strong scaling）と弱線形拡張性（weak scaling）の２種類が用いられる。前者を評価するためには、全体の計算規模を固定した上で、演算プロセス数（本実施形態におけるノード１０の数）を増やした時に演算時間がどれだけ減少するかを調べる。例えば、演算プロセス数を２倍にしたときに演算時間が半分になれば理想的な強線形拡張性を持つといえる。一方、後者を評価するためには、演算プロセスに対する計算規模を固定した上で、演算プロセス数を増やした時に演算時間がどれだけ変化するかを調べる。例えば、全体の計算規模と演算プロセス数とを共に２倍にしても、演算時間が変化しなければ高い弱線形拡張性を持つといえる。
【００８０】
図９に、本実施形態におけるＮ２５６計算に対して演算プロセス数を変化させた場合の（ａ）１時間ステップに要した演算時間［ｓ］の結果、及び（ｂ）演算速度［ＭＦＬＯＰＳ］の結果を示す。両図の横軸は共に演算プロセス数である。全体の計算規模（２５６３格子、２６２，１４４粒子）を固定した上で演算プロセス数を変化させたので、上述の強線形拡張性を調べた結果である。この性能測定にはＳＧＩ ＡＬＴＩＸ４７００システム（６．４ＧＦＬＯＰＳ／ｃｏｒｅ）が用いられた。図９（ａ）から、演算時間が演算プロセス数に反比例して減少していくことがわかる。それに対応して、図９（ｂ）では、演算速度がプロセス数にほぼ比例して増加していくことがわかる。
【００８１】
例えば、ＩＢＭ ＢｌｕｅＧｅｎｅシステムを使った等方性乱流場に対する大規模並列ＤＮＳでは、理論ピーク性能の約１０％程度の性能であったとの報告がある。本実施形態においては、粒子運動の計算及び衝突判定までを入れた上でも、流体計算のみの場合と同程度の高い性能を達成できた。
【００８２】
上述したように、本実施形態は、乱流場を有限差分法によって計算し、粒子運動をラグラジアン追跡法によって計算し、粒子間の衝突データを取得する高効率な大規模並列計算法である。本実施形態は、以下のような特徴を有している。
（ｉ）定常等方性乱流場を得るための強制法として、有限差分法の高い並列化率を損なわないＲＣＦ法を用いた。
（ｉｉ）粒子間衝突を効率的に検出するためにセルインデックス法を用いた。
（ｉｉｉ）共有メモリ型並列（自動並列化ライブラリ）と分散メモリ型並列（ＭＰＩライブラリ）の両方に対応し、両方を同時に用いるハイブリッド計算に対応した。
（ｉｖ）ＭＰＩ並列計算領域からはみ出した粒子情報の通信やセルインデックス情報の更新を所定の周期の時間ステップ毎にだけ行うようにコードを設計し、分散メモリ型並列計算でのプロセス間通信コストを削減した。
【００８３】
実際に、本実施形態に基づいた処理を実行することによって、Ｒｅλ＝２０９という高レイノルズ数における慣性粒子の衝突頻度データを得ることができた。これは従来の粒子衝突計算における最大レイノルズ数Ｒｅλ＝８５を大幅に更新する貴重なデータである。また、本実施形態で得られた高レイノルズ数における衝突頻度データを用いて、既往の衝突頻度モデルの検証も行った。その結果、既往の衝突頻度モデルは高レイノルズ数の時に衝突頻度を過小評価することが明らかになった。
【符号の説明】
【００８４】
１…粒子状態計算システム、１０…ノード、１１…粒子状態計算部、１２…周期判断部、１３…自仮想粒子情報出力部、１４…他仮想粒子情報入力部、１５…自仮想粒子状態出力部、１６…他仮想粒子状態入力部、１７…出力部。
【特許請求の範囲】
【請求項１】
複数の分割された計算領域の一つである自領域に配置される粒子の状態を時間ステップ毎に繰り返し計算する粒子状態計算装置であって、
各時間ステップにおける前記自領域の粒子の状態を計算する粒子状態計算手段と、
前記粒子状態計算手段における対象の時間ステップが、予め設定された周期の時間ステップか否かを判断する周期判断手段と、
前記周期判断手段によって予め設定された周期の時間ステップであると判断された時間ステップに、前記粒子状態計算手段によって計算された粒子の位置に基づいて、当該粒子から、前記自領域に隣接する前記分割された計算領域である隣接領域において計算に用いられる自領域仮想粒子を特定して、特定した自領域仮想粒子を示す情報を当該隣接領域の粒子の状態を計算する別の粒子状態計算装置に出力する自仮想粒子情報出力手段と、
前記別の粒子状態計算装置から、前記自領域において計算に用いられる他領域仮想粒子を特定する情報を入力する他仮想粒子情報入力手段と、
各時間ステップにおいて粒子状態計算手段によって算出された前記自領域仮想粒子の状態を示す情報を前記別の粒子状態計算装置に出力する自仮想粒子状態出力手段と、
各時間ステップにおいて前記別の粒子状態計算装置から前記他領域仮想粒子の状態を示す情報を入力する他仮想粒子状態入力手段と、
前記粒子状態計算手段による計算に応じた情報を出力する出力手段と、を備え、
前記粒子状態計算手段は、前記他仮想粒子情報入力手段によって入力される前記他領域仮想粒子を特定する情報、及び前記他仮想粒子状態入力手段によって入力された前記他領域仮想粒子の状態を示す情報を用いて粒子の状態を計算する粒子状態計算装置。
【請求項２】
前記粒子は、流体中の粒子であり、
前記粒子状態計算手段は、有限差分法によって流体の状態を算出して、当該流体の影響を考慮して粒子の状態を計算する、
ことを特徴とする請求項１に記載の粒子状態計算装置。
【請求項３】
前記粒子状態計算手段は、前記他仮想粒子情報入力手段によって入力される前記他領域仮想粒子を特定する情報に基づいて、粒子間の近傍関係を特定して、当該近傍関係に基づいて粒子の状態を計算することを特徴とする請求項１又は２に記載の粒子状態計算装置。
【請求項４】
複数の分割された計算領域の一つである自領域に配置される粒子の状態を時間ステップ毎に繰り返し計算する粒子状態計算装置による粒子状態計算方法であって、
各時間ステップにおける前記自領域の粒子の状態を計算する粒子状態計算ステップと、
前記粒子状態計算ステップにおける対象の時間ステップが、予め設定された周期の時間ステップか否かを判断する周期判断ステップと、
前記周期判断ステップにおいて予め設定された周期の時間ステップであると判断された時間ステップに、前記粒子状態計算ステップにおいて計算された粒子の位置に基づいて、当該粒子から、前記自領域に隣接する前記分割された計算領域である隣接領域において計算に用いられる自領域仮想粒子を特定して、特定した自領域仮想粒子を示す情報を当該隣接領域の粒子の状態を計算する別の粒子状態計算装置に出力する自仮想粒子情報出力ステップと、
前記別の粒子状態計算装置から、前記自領域において計算に用いられる他領域仮想粒子を特定する情報を入力する他仮想粒子情報入力ステップと、
各時間ステップにおいて粒子状態計算ステップにおいて算出された前記自領域仮想粒子の状態を示す情報を前記別の粒子状態計算装置に出力する自仮想粒子状態出力ステップと、
各時間ステップにおいて前記別の粒子状態計算装置から前記他領域仮想粒子の状態を示す情報を入力する他仮想粒子状態入力ステップと、
前記粒子状態計算ステップにおける計算に応じた情報を出力する出力ステップと、を含み、
前記粒子状態計算ステップにおいて、前記他仮想粒子情報入力ステップにおいて入力される前記他領域仮想粒子を特定する情報、及び前記他仮想粒子状態入力ステップにおいて入力された前記他領域仮想粒子の状態を示す情報を用いて粒子の状態を計算する粒子状態計算方法。
【請求項１】
複数の分割された計算領域の一つである自領域に配置される粒子の状態を時間ステップ毎に繰り返し計算する粒子状態計算装置であって、
各時間ステップにおける前記自領域の粒子の状態を計算する粒子状態計算手段と、
前記粒子状態計算手段における対象の時間ステップが、予め設定された周期の時間ステップか否かを判断する周期判断手段と、
前記周期判断手段によって予め設定された周期の時間ステップであると判断された時間ステップに、前記粒子状態計算手段によって計算された粒子の位置に基づいて、当該粒子から、前記自領域に隣接する前記分割された計算領域である隣接領域において計算に用いられる自領域仮想粒子を特定して、特定した自領域仮想粒子を示す情報を当該隣接領域の粒子の状態を計算する別の粒子状態計算装置に出力する自仮想粒子情報出力手段と、
前記別の粒子状態計算装置から、前記自領域において計算に用いられる他領域仮想粒子を特定する情報を入力する他仮想粒子情報入力手段と、
各時間ステップにおいて粒子状態計算手段によって算出された前記自領域仮想粒子の状態を示す情報を前記別の粒子状態計算装置に出力する自仮想粒子状態出力手段と、
各時間ステップにおいて前記別の粒子状態計算装置から前記他領域仮想粒子の状態を示す情報を入力する他仮想粒子状態入力手段と、
前記粒子状態計算手段による計算に応じた情報を出力する出力手段と、を備え、
前記粒子状態計算手段は、前記他仮想粒子情報入力手段によって入力される前記他領域仮想粒子を特定する情報、及び前記他仮想粒子状態入力手段によって入力された前記他領域仮想粒子の状態を示す情報を用いて粒子の状態を計算する粒子状態計算装置。
【請求項２】
前記粒子は、流体中の粒子であり、
前記粒子状態計算手段は、有限差分法によって流体の状態を算出して、当該流体の影響を考慮して粒子の状態を計算する、
ことを特徴とする請求項１に記載の粒子状態計算装置。
【請求項３】
前記粒子状態計算手段は、前記他仮想粒子情報入力手段によって入力される前記他領域仮想粒子を特定する情報に基づいて、粒子間の近傍関係を特定して、当該近傍関係に基づいて粒子の状態を計算することを特徴とする請求項１又は２に記載の粒子状態計算装置。
【請求項４】
複数の分割された計算領域の一つである自領域に配置される粒子の状態を時間ステップ毎に繰り返し計算する粒子状態計算装置による粒子状態計算方法であって、
各時間ステップにおける前記自領域の粒子の状態を計算する粒子状態計算ステップと、
前記粒子状態計算ステップにおける対象の時間ステップが、予め設定された周期の時間ステップか否かを判断する周期判断ステップと、
前記周期判断ステップにおいて予め設定された周期の時間ステップであると判断された時間ステップに、前記粒子状態計算ステップにおいて計算された粒子の位置に基づいて、当該粒子から、前記自領域に隣接する前記分割された計算領域である隣接領域において計算に用いられる自領域仮想粒子を特定して、特定した自領域仮想粒子を示す情報を当該隣接領域の粒子の状態を計算する別の粒子状態計算装置に出力する自仮想粒子情報出力ステップと、
前記別の粒子状態計算装置から、前記自領域において計算に用いられる他領域仮想粒子を特定する情報を入力する他仮想粒子情報入力ステップと、
各時間ステップにおいて粒子状態計算ステップにおいて算出された前記自領域仮想粒子の状態を示す情報を前記別の粒子状態計算装置に出力する自仮想粒子状態出力ステップと、
各時間ステップにおいて前記別の粒子状態計算装置から前記他領域仮想粒子の状態を示す情報を入力する他仮想粒子状態入力ステップと、
前記粒子状態計算ステップにおける計算に応じた情報を出力する出力ステップと、を含み、
前記粒子状態計算ステップにおいて、前記他仮想粒子情報入力ステップにおいて入力される前記他領域仮想粒子を特定する情報、及び前記他仮想粒子状態入力ステップにおいて入力された前記他領域仮想粒子の状態を示す情報を用いて粒子の状態を計算する粒子状態計算方法。
【図１】
【図２】
【図３】
【図４】
【図５】
【図７】
【図８】
【図９】
【図６】
【図２】
【図３】
【図４】
【図５】
【図７】
【図８】
【図９】
【図６】
【公開番号】特開２０１２−１２８７９３（Ｐ２０１２−１２８７９３Ａ）
【公開日】平成２４年７月５日（２０１２．７．５）
【国際特許分類】
【出願番号】特願２０１０−２８１８０８（Ｐ２０１０−２８１８０８）
【出願日】平成２２年１２月１７日（２０１０．１２．１７）
【出願人】（５０４１９４８７８）独立行政法人海洋研究開発機構 (110)
【公開日】平成２４年７月５日（２０１２．７．５）
【国際特許分類】
【出願日】平成２２年１２月１７日（２０１０．１２．１７）
【出願人】（５０４１９４８７８）独立行政法人海洋研究開発機構 (110)
[ Back to top ]