JP2015530636A - 粒子流動のシミュレーションシステム及びその方法 - Google Patents

粒子流動のシミュレーションシステム及びその方法 Download PDF

Info

Publication number
JP2015530636A
JP2015530636A JP2015521951A JP2015521951A JP2015530636A JP 2015530636 A JP2015530636 A JP 2015530636A JP 2015521951 A JP2015521951 A JP 2015521951A JP 2015521951 A JP2015521951 A JP 2015521951A JP 2015530636 A JP2015530636 A JP 2015530636A
Authority
JP
Japan
Prior art keywords
gpu
particle
particles
calculation
node
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
JP2015521951A
Other languages
English (en)
Other versions
JP6009075B2 (ja
Inventor
磊 楊
磊 楊
記 斎
記 斎
園 田
園 田
笑菲 高
笑菲 高
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Institute of Modern Physics of CAS
Original Assignee
Institute of Modern Physics of CAS
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Institute of Modern Physics of CAS filed Critical Institute of Modern Physics of CAS
Publication of JP2015530636A publication Critical patent/JP2015530636A/ja
Application granted granted Critical
Publication of JP6009075B2 publication Critical patent/JP6009075B2/ja
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2210/00—Indexing scheme for image generation or computer graphics
    • G06T2210/52—Parallel processing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2210/00—Indexing scheme for image generation or computer graphics
    • G06T2210/56—Particle system, point based geometry or rendering

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Physics (AREA)
  • Algebra (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本発明が開示したＧＰＵに基づく粒子流動のシミュレーションシステム及び方法は、クライアントから入力される粒子のモデリング情報に応じて粒子の情報を生成し、幾何体の情報を生成し、粒子の情報及び幾何体の情報を受信し、粒子数及び各計算ノードにフリーなＧＰＵの数に応じて、どれらの計算ノードにおけるどれらのＧＰＵを使用するかを確定し、そして、確定したＧＰＵの数及び粒子の空間における分布状況に応じてどれらの粒子がどの計算ノードのどのＧＰＵによって算出されるかを確定し、確定した結果に応じて割り当てを行い、粒子の衝突による各粒子のストレスを複数のＧＰＵにおいて並列に算出して加速度を更に算出し、粒子の流動のシミュレーションを行い、シミュレーションの結果を表示することを含む。本発明の実施例によって、高密度の粒子の擬似実験シミュレーションを実現でき、エネルギー消耗を低下すると共に演算の効率を向上することができる。【選択図】図４

Description

本発明は、粒子流動のシミュレーションの技術分野に関する。具体的には、粒子物質又は固体構造の研究に適用できる、ＧＰＵに基づく粒子シミュレーションシステム及びその方法に関する。
粒子システムは、ずっと注目されている研究内容である。例えば、食品制御、化学、土木工事、オイルガス、鉱物採掘、製薬、粉末冶金、エネルギー等の産業分野に多く応用されている。理論的な研究において、如何にして積み上げて最も密着な堆積に達するか、砂の堆がどのような状況で崩れるかを研究して雪崩等の課題への研究が行われている。人々が、関連的な粒子システムを研究するために、大型の実験用粒子システムを設立する必要があり、手間がかかる。そして、粒子システムの一部は、コストが高く、極端な条件下で運行する必要があるため、実験で完成する可能性がない。しかしながら、虚構の実験に基づくシミュレーションシステムには、類似の問題が存在していない。
現在、粒子システム模擬の算出方法は、主にＤＥＭ（離散要素法：Ｄｉｓｃｒｅｔｅ Ｅｌｅｍｅｎｔ Ｍｅｔｈｏｄ）方法である。ＤＥＭ方法は、有限要素法、数値流体力学（ＣＦＤ）に継いで、物質システム問題を分析するためのもう１種の数値算出方法である。ＤＥＭ方法は、微小的な体系のパラメータ化モデルを構築したことによって、粒子行為の模擬及び分析を行い、粒子、構造、流体、電磁及びその結合等に関するたくさんの綜合的な問題を解決するために、プラットフォームを提供し、科学過程の分析、製品設計の最適化及び研究開発への力強いツールともなっている。現在、ＤＥＭ方法は、科学研究における適用に加え、科学技術及び工業分野においても熟しつつ、粒子物質の研究、岩土工事及び地質工事などの科学及び応用分野から、工業過程及び工業製品の設計、研究開発の分野まで広げ、たくさんの工業分野において重要な成果を収めた。
ＤＥＭ方法は、シミュレーション精度が高いが、計算量が大きいとの特徴を有する。現在、ＤＥＭ方法は、主にＣＰＵを用いて実現される。これらの方法は、ＣＰＵの算出能力が不足であることによる算出規模が不足となり、納得できる機器時間内において、非常に小さい空間サイズ及び時間範囲のみしか算出できない。或いは、大規模ひいては超大規模なＣＰＵコンピュータのクラスタを建設する必要とし、建設のコストが高く、そして、電力消費量が大きすぎて、使用及びメンテナンスのコストも非常に高くなっている。また、現在、ＣＰＵで実現したＤＥＭ方法は、粒子数が少ない或いは低密度粒子の衝突を実現できたとしても、高密度の大量の粒子衝突の模擬を完全に実現することができない。
ＧＰＵ（Ｇｒａｐｈｉｃｓ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ、グラフィックス・プロセッシング・ユニット）で汎用計算を行う技術は、ますます熟に成りつつである。例えば、ｎＶＩＤＩＡ及びＡＭＤという現在の２つの表示カードメーカーは、いずれもＧＰＵ汎用計算をサポートできる。本願の発明人は、上記の課題に鑑みて、ＧＰＵに基づく粒子流動のシミュレーションシステム及び方法を提供した。
本発明によれば、高密度の粒子擬似実験シミュレーションを実現し、エネルギー消耗量を低下させると共に算出効率を向上することができるＧＰＵに基づく粒子流動のシミュレーションシステム及びその方法を提供した。
本発明の一局面によれば、ＧＰＵに基づく粒子流動のシミュレーション方法を提供した。該ＧＰＵに基づく粒子流動のシミュレーション方法は、並列の複数のＧＰＵ上で離散要素法（ＤＥＭ）方法を実行して粒子流動のシミュレーションを行うＧＰＵに基づく粒子流動のシミュレーション方法であって、ＤＥＭ方法で粒子をモデリングし、作られたＤＥＭモデルを複数の粒子として割り当て、該複数の粒子を複数の計算ノードに割り当て処理を行い、各計算ノードのＣＰＵ及びＧＰＵに記憶空間がそれぞれ割り当てされ、ＣＰＵにおいてデータの初期化を行い、初期化されたデータをＣＰＵの記憶空間から前記ＧＰＵの記憶空間へコピーするステップａと、
上記各計算ノードのＧＰＵは、各粒子の処理を行い、各計算ノードのＧＰＵの各ストリーミングプロセッサは、１つの粒子を処理し、ＧＰＵの記憶空間に記憶される粒子の座標及び粒子の速度を更新するステップｂと、
ステップｂの処理過程において、各計算ノードが制御する粒子を確定し、各計算ノードが制御する粒子の個数をＣＰＵの記憶空間へコピーして、各計算ノードがどれらの粒子を算出するかを均衡負荷の原則に従って動的に確定できるように、ＧＰＵの記憶空間における粒子の数に応じて動的分割を行うステップｃと、
ＭＰＩインターフェースプロトコルによって、上記データが動的分割された粒子を各計算ノード間に遷移させるステップｄと、
ステップｃで取得した各計算ノードが制御する粒子に応じて、ＧＰＵにおいて重畳領域を算出し、データをＣＰＵメモリへコピーしてから、ＭＰＩインターフェースプロトコルによってデータのやり取りを行うステップeと、
各計算ノードのＧＰＵにおける各ストリーミングプロセッサは、各粒子の座標に応じて、各粒子が位置するＧＰＵの記憶空間におけるグリッドの番号を算出するステップｆと、
各計算ノードのＧＰＵにおける各ストリーミングプロセッサは、各粒子の運動中のストレス及び加速度を処理して算出するステップｇと、
各計算ノードのＧＰＵにおける各ストリーミングプロセッサは、各粒子の速度を処理するステップｈと、
指定された歩数に達するまでに、ステップｂに戻すステップｉと、
マスターノード及び計算ノードの記憶空間を釈放するステップｊと、を含む。
一実施例においては、ステップｂ、ステップｆ、ステップｇ及びステップｈは、ＧＰＵにおいて各粒子に対して並列なデータ処理を行う。すなわち、各ＧＰＵが粒子に対する処理は、同期的に行われる。
一実施例においては、ステップｄにおいて、前記粒子が各ノード間で遷移し、粒子がノード間で伝送遷移する方法を用いて、すなわち、ＭＰＩインターフェースで関数を送受信し、粒子の各物理量の送受信を実現し、そして、粒子のノード間における伝送遷移を実現した。
一実施例においては、ステップｅにおいて、前記ＧＰＵにおいて重畳領域を算出することは、ＧＰＵにおいて重畳領域（Ｏｖｅｒｌａｐ領域）を算出することを利用し、ＧＰＵの１つのストリーミングプロセッサが１つのグリッドを処理することを含む。三次元の場合において、それぞれのグリッドは、２６個のグリッドに隣接し、隣接するグリッドが現在の計算ノード内に位置するか否かを判断し、位置しなければ、ｏｖｅｒｌａｐ領域とし、他のノードから遷移して取得する。
本発明の別の局面によれば、ＧＰＵに基づく粒子流動のシミュレーション方法を提供した。該方法は、
粒子の材料、粒子のパラメーター、境界の条件、幾何体の形状、及び粒子の初期分布の領域を確定し、予め定められた粒子の分布領域及び数量に応じて粒子を生成するモデリングステップと、
粒子の総数及び複数の計算ノードにフリーなＧＰＵ数に応じて、最適なＧＰＵ数を確定し、最適なＧＰＵ数及び現在フリーなＧＰＵ数に応じて、算出関与のＧＰＵを確定し、算出関与のＧＰＵの状態をビジーに設置するタスク管理ステップと、
算出ステップと、を含み、
該算出ステップは、
各計算ノードの算出関与のＧＰＵを初期化し、算出に必要な粒子の情報を各ＧＰＵに発信するステップと、
各ＧＰＵが、予め定められた速度を並列に更新し、受信した粒子の情報をソートして各自のソートセルリストを生成するステップと、
各ＧＰＵが、現在各自のコースにおける非ゼローのグリッドの番号及びグリッドにおける粒子数を並列に算出し、マスターノードに発信して、マスターノードによって各ＧＰＵの最適な粒子数に応じて、グリッドを動的分割し、各ＧＰＵが並列に算出するグリッドの数及びグリッドの番号を確定するステップと、
マスターノードが確定した結果に応じて、各ＧＰＵが粒子情報を並列に送受信し、各ＧＰＵにおいて各自のソートセルリストを生成し直すステップと、
各ＧＰＵにおいて現在時刻の衝突リストを生成し、現在時刻の衝突リストと一つ前の時刻の衝突リストと接線相対変位に応じて、各ＧＰＵにおいて接線相対変位の位置を並列に調整するステップと、
接触力学モデルによって、各ＧＰＵにおいて各粒子のストレス及び加速度を並列に算出するステップと、
現在の算出結果を記憶するステップと、
算出が完成していなければ、各ＧＰＵが予め定められた速度及び座標を並列に更新するステップに戻し、そうでなければ、算出ステップを終了させるステップと、を含む。
一実施例において、前記方法は、更に表示ステップを含み、上記表示ステップは、境界の条件を確定し、幾何体の境界を透明な曲面で作るステップと、粒子の位置及び粒子の直径に応じて、粒子を同じ色又は異なる色のペレットで描くステップと、階調画像でスカラ場を表示し、粒子情報を重み付けしてグリッドにマッピングすることによって、ベクトル場を流線描き方法で描くステップとを含む。
一実施例において、全ての粒子の物理情報を外部の記憶装置に格納する。
一実施例において、各ＧＰＵは、関連の物理統計量を並列に算出する。
一実施例において、予め定められた粒子の分布領域及び数量によって粒子を生成することは、粒子の数量条件を満たすまでに、比較的に小さい空間内でいくつかの粒子を生成して、これらの粒子を平行移動してコピーを行い、他の空間を充填する。
一実施例において、ソートセルリストは、粒子が位置するグリッドに従って、全ての粒子をソートする。
一実施例において、動的分割方法を採用して、非ゼローのグリッドの番号及びグリッドにおける数粒子数をＧＰＵにおいて並列に算出する。
一実施例において、各ＧＰＵにおいて、１つのスレッド（thread）が１つの粒子に対応するとの方式を用いて算出を行う。
一実施例において、接線相対変位を算出することは、一つ前の時刻の接線相対変位を記録し、現在時刻の衝突リストに応じてそれを更新することを含む。
一実施例において、コピー又はポインターやり取り技術を用いて、現在の算出結果をアレイに記憶する。
本発明のもう一つの局面によれば、ＧＰＵに基づく粒子流動のシミュレーションシステムを提供した。該ＧＰＵに基づく粒子流動のシミュレーションシステムは、
粒子の材料、粒子のパラメーター、境界の条件、幾何体の形状、及び粒子の初期分布の領域を確定し、予め定められた粒子の分布領域及び数量に応じて粒子を生成するように構成されるモデリングモジュールと、
粒子の総数及び複数の計算ノードにフリーなＧＰＵ数に応じて、最適なＧＰＵ数を確定し、最適なＧＰＵ数及び現在フリーなＧＰＵ数に応じて、算出関与のＧＰＵを確定し、算出関与のＧＰＵの状態をビジーに設置するように構成されるタスク管理モジュールと、
算出モジュールと、を含み、
該算出モジュールは、
各計算ノードの算出関与のＧＰＵを初期化し、算出に必要な粒子の情報を各ＧＰＵに発信し、
各ＧＰＵが、予め定められた速度及び座標を並列に更新し、受信した粒子の情報をソートして各自のソートセルリストを生成し、
各ＧＰＵが、現在各自のコースにおける非ゼローのグリッドの番号及びグリッドにおける粒子数を並列に算出し、マスターノードに発信して、マスターノードによって各ＧＰＵの最適な粒子数に応じて、グリッドを動的分割し、各ＧＰＵが並列に算出するグリッドの数及びグリッドの番号を確定し、
マスターノードが確定した結果に応じて、各ＧＰＵが粒子情報を並列に送受信し、各ＧＰＵにおいて各自のソートセルリストを生成し直し、
各ＧＰＵにおいて現在時刻の衝突リストを生成し、現在時刻の衝突リストと一つ前の時刻の衝突リストと接線相対変位に応じて、各ＧＰＵにおいて接線相対変位の位置を並列に調整し、
接触力学モデルによって、各ＧＰＵにおいて各粒子のストレス及び加速度を並列に算出し、
現在の算出結果を記憶し、
算出が完成していなければ、各ＧＰＵが予め定められた速度及び座標を並列に更新するステップに戻し、そうでなければ、算出ステップを終了させるように構成される。
一実施例において、前記システムは、更に表示モジュールを含み、上記表示モジュールは、境界の条件を確定し、幾何体の境界を透明な曲面で作り、粒子の位置及び粒子の直径に応じて、粒子を同じ色又は異なる色のペレットで描き、階調画像でスカラ場を表示し、粒子情報を重み付けしてグリッドにマッピングすることによって、ベクトル場を流線描き方法で描くように構成される。
本発明のまた一局面によれば、ＧＰＵに基づく粒子流動のシミュレーションシステムを提供した。該ＧＰＵに基づく粒子流動のシミュレーションシステムは、
クライアントから入力される粒子のモデリング情報に応じて、粒子の情報を生成すると共に、幾何体の情報を生成するように構成される前端サーバと、
前端サーバから粒子の情報及び幾何体の情報を受信し、粒子の数及び各計算ノードにフリーなＧＰＵの数に応じて、どれらの計算ノードにおけるどれらのＧＰＵを使用するかを確定し、そして、確定したＧＰＵの数及び粒子の空間における分布状況に応じてどれらの粒子がどの計算ノードのどのＧＰＵによって算出されるかを確定し、確定した結果によって割り当てるように構成される管理ノードと、
それぞれが複数のＧＰＵを含み、複数のＧＰＵにおいて粒子の衝突による各粒子のストレスを並列に算出し、加速度を更に算出して、粒子の流動をシミュレーションするように構成される複数の計算ノードと、
シミュレーションの結果を表示するように構成される後端サーバと、を備える。
一実施例において、前端サーバは、幾何体を有限の曲面に分解し、これらの曲面に番号をつけることによって、幾何体の情報を生成する。
一実施例において、後端サーバは、表示されるシミュレーション結果において、幾何体の境界を透明な曲面で作り、粒子の位置及び粒子の直径によって、粒子を同じ色又は異なる色のペレットで描き、且つ、階調画像でスカラ場を表示し、粒子情報を重み付けしてグリッドにマッピングすることによって、ベクトル場を流線描き方法で描く。
一実施例において、前端サーバ、管理ノード、計算ノード及び後端サーバは、ＩＢ（ＩｎｆｉｎｉＢａｎｄ）ネットワークによって通信する。
本発明によれば、複数のＧＰＵに基づく、モデリングから結果表示までのシミュレーションシステムを実現し、複数のＧＰＵのハードウェア特徴を利用して、複数のＧＰＵの粒子流動のシミュレーション方法を実現した。本発明の実施例によれば、ＧＰＵの強い浮動小数点演算能力、広い帯域幅及び複数の軽量計算コアという特徴によって、ＧＰＵ内の大量のストリーミングプロセッサを十分に利用し、分子動力学の加速アルゴリズムをＤＥＭアルゴリズムに合理的に引き入れ、ＤＥＭアルゴリズムをＧＰＵのハードウェア構造により適応できる。複数のＧＰＵで実現する場合、該アルゴリズムが、データを動的分割して負荷均衡を実現する方法を採用し、Ｏｖｅｒｌａｐ領域及び通信量を低減し、ＧＰＵ及びＣＰＵの利用率及び演算効率を大きく向上できた。納得できるエネルギー消耗及び時間の条件下で、非常によい算出効果を取得し、エネルギー消耗が小さく、メンテナンスコストが低く、且つ演算の効率を向上する効果を奏した。
以下、図面及び実施例によって、本発明の技術案をより詳細に説明する。
図１は、本発明の実施例に係るＧＰＵに基づく粒子流動のシミュレーションシステムの構造模式図である。 図２は、本発明の一実施例に係るＧＰＵに基づく粒子流動のシミュレーション方法のフローチャートである。 図３は、本発明の別の実施例に係るＧＰＵに基づく粒子流動のシミュレーションシステムのモジュール構造模式図である。 図４は、本発明の実施例に係る計算モジュールの操作フローチャートである。
以下、本発明の好ましい実施例について、図面を参照して説明する。ここで記述された好ましい実施例が本発明を説明して解釈することのみに用いられ、本発明を限定するものでないことを理解すべきである。
図１は、本発明の実施例に係る、ＧＰＵに基づく粒子流動のシミュレーションシステムの構造模式図である。図１に示すように、該システムは、前端サーバ１０、後端サーバ２０、管理ノード３０、複数の計算ノード４０−１，…，４０−Ｎ（Ｎは１よりも大きい整数である）、ＩＢスイッチ装置５０及びイーサネット（登録商標）・スイッチ装置６０を含む。また、図１は、該システムがクライアント及び記憶装置を含むことを更に示している。クライアントは、インターネットを介して前端サーバ１０と通信可能であり、現場の実験員に粒子流動のシミュレーション実験を遠隔距離で行わせることを可能にした。例えば、ユーザは、クライアントで、例えば粒子数、大きさや材料などの情報（ヤング率、ポアソン（Ｐｏｉｓｓｏｎ）比、密度、回復係数等）及び粒子の分布範囲、摩擦係数、境界条件などのパラメーターというモデリングに必要な情報またはパラメーターを入力すると共に、粒子ボールと接触する幾何体の材料情報を提供し、これらの情報またはパラメーターを、前端サーバに伝送することができる。該外部の記憶装置は、フリーズ、停電などの意外状況の発生によるデータの紛失を防止するよう、例えば各計算ノードの算出結果を記憶することができる。ここで、クライアント及び外部の記憶装置は選択可能であり、例えば、ユーザは、直接に前端サーバで入力を行っても良い。あるいは、計算ノードの算出結果は前端または後端サーバなどに記憶されることができる。
図１において、前端サーバ１０、後端サーバ２０及び計算ノード４０は、ＩＢスイッチ装置５０を介して互いに接続されている。そして、前端サーバ１０、管理ノード３０及び計算ノード４０は、イーサネット（登録商標）スイッチ装置６０を介して互いに接続されている。しかしながら、本発明の実施例は、ほかの任意の適宜な接続方式を採用することも可能である。一実施例において、計算ノード４０は、ＧＰＵ加速カードを有する高性能クラスタであっても良い。また、一実施例において、それぞれの計算ノードは、いずれもＧＦ１１０コア以上のＮＶＩＤＩＡ汎用計算カードを有する。一実施例において、計算ノードは、４０ＧｂのＩＢ（ＩｎｆｉｎｉＢａｎｄ）ネットワーク接続を使用する。一実施例において、前後端サーバは、それぞれＱｕａｄｒｉｏ６０００表示カードを有する一台のグラフィックスワークステーションである。例えば、ワークステーションのメモリは、３２Ｇよりも大きく、ＩＢネットワークカードを有する。
本発明の実施例のＧＰＵに基づく粒子流動のシミュレーションシステムにおいて、前端サーバ１０は、クライアントから入力された粒子モデリング情報に応じて粒子情報を生成すると共に、幾何体情報を生成する。例えば、前端サーバ１０は、粒子のサイズ、材料及び幾何構造に関する入力を受信でき、粒子を交互的に増加・削除し、粒子の位置を移動させることもできる。前端サーバ１０は、幾何体を有限的な曲面に分解し、これらの曲面に番号をつけることによって、幾何体情報を生成することができる。管理ノード３０は、現在の各計算ノードの運行状態、ＧＰＵの作業状態、記憶状況などを任意に観察できると共に、各タスク間に衝突が発生しないことを保証するように、提出されたタスクを中止することもできる。例えば、管理ノード３０は、粒子情報及び幾何体情報を、前端サーバ１０から受信し、粒子数及び各計算ノードにおいてフリーなＧＰＵの数に応じて、どれらの計算ノードのどれらのＧＰＵを使用するかを確定する。その後、確定されたＧＰＵの数及び粒子が空間における分布状況によって、どれらの粒子が、どの計算ノードのどのＧＰＵによって算出されるかを確定して、確定した結果に応じて割り当てを行う。計算モジュール全体が各々の計算ノード４０から構成され、複雑な境界問題を処理することができ、複数のＧＰＵを並列に運行し、中断（例えば、停電）機能を有し、中断前の状態に引き続いて演算することができる。該計算モジュールは、データの動的分割方法及びポインターやり取り技術を用いて、データの動的平衡を保証する。例えば、各計算ノード４０は、それぞれのＧＰＵにおいて粒子衝突による各粒子のストレスを並列に算出し、加速度を更に算出して、粒子流動のシミュレーションを行う。後端サーバ２０は、シミュレーションの結果を表示する。例えば、現在の粒子の構造、温度場、ストリーム場、圧力場などのパラメーターを動的に表示する。また、交互の方式で観察の角度を調整し、粒子グループを任意に拡大縮小させても良い。例えば、後端サーバ２０は、ディスプレイなどの出力装置を含んでも良い。後端サーバ２０は、幾何体の境界を透明的な曲面で作り出し、粒子の位置及び粒子の直径に応じて、粒子を、同じ色又は異なる色のペレットで描き、そして、階調画像で温度場などのスカラ場を表示し、粒子情報を重み付けてグリッドにマッピングすることによって、ストリーム場、圧力場などのベクトル場を、流線描き方法で描くことができる。
以上のシステムは、本発明の基本的な構想の一種の表現のみである。当業者は、上記各部品の機能を更に割り当てて組み合わせることによって他のシステムを構築して形成できることを理解すべきである。また、機能が十分に強ければ、上記各部品の機能は、１つのコンピュータ又はワークステーションに集積しても良い。
図２は、本発明の実施例のシミュレーションシステムに実行される、ＧＰＵに基づく粒子流動のシミュレーション方法のフローチャートである。図２に示すように、該シミュレーション方法は、以下のステップを含む。
ステップ２０１：ＤＥＭ方法を用いて粒子のモデリングを行い、作成したＤＥＭモデルを複数の粒子として割り当て、該複数の粒子を複数の計算ノードに割り当て処理を行う。それぞれの計算ノードのＣＰＵ及びＧＰＵに記憶空間が割り当てされており、ＣＰＵにおいてデータの初期化を行い、初期化されたデータを、ＣＰＵの記憶空間から前記ＧＰＵの記憶空間へコピーする。
ステップ２０２：上記それぞれの計算ノードのＧＰＵは、各粒子を処理する。その中、それぞれの計算ノードのＧＰＵの各ストリーミングプロセッサは、１つの粒子を処理し、ＧＰＵの記憶空間に記憶された粒子の座標及び粒子の速度を更新する。
ステップ２０３：ＧＰＵの記憶空間に記憶された粒子の座標が変化することによって、負荷の均衡を保証するために、毎回の算出において各ノードが算出する粒子は異なっている。まず、それぞれの計算ノードのＧＰＵは、該ノードが制御する粒子の数を算出し、各ＧＰＵが制御する粒子の数をＣＰＵの記憶空間へコピーして、ＧＰＵの記憶空間におけるグリッドの粒子数に応じてデータの動的分割を行う。すなわち、負荷均衡の原則に従い、各ノードがどれらの粒子を算出するかを算出する。
ステップ２０４：ＭＰＩインターフェースプロトコルによって、データが動的分割された上記粒子を、それぞれの計算ノード間に遷移させる。
ステップ２０５：ステップ２０３で取得した各計算ノードが制御する粒子によって、ＧＰＵにおいて重畳領域を算出し、データをＣＰＵのメモリ内へコピーし、その後、ＭＰＩインターフェースプロトコルによってデータのやり取りを行う。
ステップ２０６：それぞれの計算ノードのＧＰＵにおける各ストリーミングプロセッサは、それぞれの粒子の座標に応じて、各粒子がＧＰＵの記憶空間に位置するグリッドの番号を算出する。
ステップ２０７：それぞれの計算ノードのＧＰＵにおける各ストリーミングプロセッサは、それぞれの粒子運動中のストレス及び加速度の算出処理を行う。
ステップ２０８：それぞれの計算ノードのＧＰＵにおける各ストリーミングプロセッサは、それぞれの粒子速度の算出処理を行う。
ステップ２０９：指定の歩数に達するまでにステップ２０２に戻り、ＤＥＭ方法を完成させる。
ステップ２１０：マスターノード及び計算ノードの記憶空間を釈放する。
その中、前記ステップ２０２、ステップ２０６、ステップ２０７及びステップ２０８においては、ＧＰＵを用いてそれぞれの粒子に対して並列なデータ処理を行う。すなわち、それぞれのＧＰＵが粒子に対する処理は、同期的に行われるものである。
ステップ２０４において、前記粒子が各ノード間における遷移は、粒子がノード間に伝送して遷移する方法を用いる。すなわち、ＭＰＩインターフェースを用いて関数を送受信し、粒子の各物理量の発信及び受信を実現し、そして粒子がノード間における伝送及び遷移を実現する。受信関数は、ＭＰＩ_Ｓｅｎｄ()及びＭＰＩ_Ｒｅｃｖ()関数である。
ステップ２０５において、前記ＧＰＵにおいて重畳領域（Ｏｖｅｒｌａｐ区）を算出することは、ＧＰＵにおいてＯｖｅｒｌａｐ領域を算出する方法を利用している。すなわち、ＧＰＵの１つのストリーミングプロセッサは、１つのグリッドの処理を行う。三次元の場合には、それぞれのグリッドは、２６個のグリッドに隣接し、そして、隣接のグリッドが現在の計算ノード中に位置するか否かを判断し、位置しなければ、ｏｖｅｒｌａｐ領域として算出し、他のノードから遷移して取得する。
具体的には、以下の通りである。
ステップ１：それぞれの計算ノードは、ＣＰＵ及びＧＰＵにおいて記憶空間を設け、ＣＰＵにおいてデータを初期化して、ＧＰＵへコピーする。
ステップ２：
計算ノードのＧＰＵの各ストリーミングプロセッサは、１つの粒子の処理を行い、１歩の粒子座標及び１／２歩の粒子速度を並列に更新する。ＣＵＤＡのＫｅｒｎｅｌ関数：
__global__void UpdateP(double＊x1, double ＊x2, double ＊x3,
double ＊vx, double ＊vy, double ＊vz,
double ＊ax, double ＊ay, double ＊az,
unsigned int NumParticles);
__global__void UpdateV (double ＊vx, double ＊vy, double ＊vz,
double ＊wx, double ＊wy, double ＊wz,
double ＊ax, double ＊ay, double ＊az,
double ＊bx, double ＊by, double ＊bz,
unsigned int NumParticles);
が含まれている。呼び出す際に、ＣＵＤＡのシンタックス（syntax）条件に従って、以下の方式：
UpdateV <<<gridsize, blocksize>>>(vx, vy, vz,
wx, wy, wz,
ax, ay, az,
bx, by, bz,
NumParticles);
を用いて呼び出す。この２つの関数の「ｂｌｏｃｋ」及び「ｇｒｉｄ」は、いずれも一次元の方式を採用し、異なる粒子数に対して、ｂｌｏｃｋ及びｇｒｉｄの値を調整でき、算出時間に対して一定の影響を与えている。
ステップ３：それぞれの計算ノードのＧＰＵにおいて、該ノードが制御する粒子を算出し、ＣＰＵへコピーして、グリッド内の粒子数に従ってデータの動的分割を行う。
算出過程において、粒子が、異なるノード間に遷移することによって、負荷が不均衡の場合を避けるのに、本発明は、データを動的に分割する方式を用いて、それぞれのノードの計算量を平衡させる。
初期状態では、仮に、Ｍ個のグリッドを有し、各グリッドにおける粒子数Ｘが同じであり、Ｍ個のグリッド（Ｇ０〜ＧＭ−１）は、それぞれ、均等にＮ段に分割され、それぞれ、Ｎ個のノード（Ｐ０〜ＰＮ−１）で処理される。これにより、それぞれのノードが算出する粒子数は、(Ｍ／Ｎ)＊Ｘとなる。反復算出された後で、各ノードＰｉが算出するグリッド範囲内の粒子総数が変化する。このため、それぞれノードが算出するグリッドの範囲を調整することによって、算出粒子の総数を変更させることができる。データの動的分割は、以下のように実現されることになる。即ち：
（１）それぞれのノードは、グリッド全体の数量Ｍと同じであるｉｎｔ型のアレイｉＣｅｌｌＣｏｕｎｔを維持する。ＣＵＤＡのコア関数ｃａｌｃＰａｒｔｉｃｌｅＮｕｍＰｅｒＣｅｌｌ（）を呼び出してそれぞれのグリッドにおける粒子の個数を算出し、それをｉＣｅｌｌＣｏｕｎｔに格納する。この時、ｉＣｅｌｌＣｏｕｎｔ中の粒子個数は、局所的なことであり、現在のノードが算出する粒子が、グリッドにおける個数のみを記録した。
（２）ＰＩＤ=０のノードをＲＯＯＴノードとし、ＭＰＩ減少関数ＭＰＩ_Ｒｅｄｕｃｅ（）を呼び出して、全てのノードｉＣｅｌｌＣｏｕｎｔの情報を、加算の操作によってＲＯＯＴノードのｉＧｌｏｂａｌＣｅｌｌＣｏｕｎｔアレイに集める。このとき、ｉＧｌｏｂａｌＣｅｌｌＣｏｕｎｔアレイに記録される各グリッドの粒子個数は、全体的なことであり、全ての粒子が各グリッドにおける個数である。
（３）ｉＧｌｏｂａｌＣｅｌｌＣｏｕｎｔアレイを用いて、各ノードの算出グリッド範囲の分割を行う。この分割は、ＣＰＵ+ＧＰＵの方式を採用している。分割のステップは、以下の通りである。
ノード個数Ｎに応じて、アレイｉＧｌｏｂａｌＣｅｌｌＣｏｕｎｔをＮ段に均等に分割し、それぞれのノードの算出グリッド範囲が同じであるとする。各ノードの算出グリッド範囲が、アレイｉＤｉｖｉｄｅｄＲｅｓｕｌｔに格納される。初期状態の場合、ｉＤｉｖｉｄｅｄＲｅｓｕｌｔにおける各元素の値は、{0,M/N-1,M/N,2M/N-1,...,(N-1)M/N,M-1}であり、ノードｉの範囲は、ｉＤｉｖｉｄｅｄＲｅｓｕｌｔ[ｉ＊２]及びｉＤｉｖｉｄｅｄＲｅｓｕｌｔ[ｉ＊２+１]によって取得することができる。
ＣＵＤＡコア関数ｄＲｅｄｕｃｅＰｅｒＳｅｇ（）を呼び出して、各段の粒子個数をそれぞれ求めて、アレイｉＰａｒｔｉｃｌｅｓＣｏｕｎｔＰｅｒＳｅｇ={Ｘ０,Ｘ１,...,ＸＮ−１}に格納する。
ＣＰＵにより、ｉＤｉｖｉｄｅｄＲｅｓｕｌｔ、ｉＰａｒｔｉｃｌｅｓＣｏｕｎｔＰｅｒＳｅｇ及びｉＧｌｏｂａｌＣｅｌｌＣｏｕｎｔに基づいて最終的な分割結果を確定する。まず、理想的な状況の下での各ノード算出粒子の個数ｉＰａｒｔｉｃｌｅｓＰｅｒＮｏｄｅＩｄｅａｌを確定して、ｉＰａｒｔｉｃｌｅｓＣｏｕｎｔＰｅｒＳｅｇ[０]の値を読み出す。若し、ｉＰａｒｔｉｃｌｅｓＣｏｕｎｔＰｅｒＳｅｇ[０] > ｉＰａｒｔｉｃｌｅｓＰｅｒＮｏｄｅＩｄｅａｌであれば、ノード０が処理する範囲は、大きすぎると分かる。このため、
iParticlesCountPerSeg[0] - iGlobalCellCount[iDividedResult[0*2+1]],
iParticlesCountPerSeg[1]+iGlobalCellCount[iDividedResult[0*2+1]],
iDividedResult[1*2] = iDividedResult[0*2+1],
iDividedResult[0*2+1]-1，
ｉＰａｒｔｉｃｌｅｓＣｏｕｎｔＰｅｒＳｅｇ[０]は、ｉＰａｒｔｉｃｌｅｓＰｅｒＮｏｄｅＩｄｅａｌと同じあるいは近接になるまでに、上記過程を繰り返して行う。若し、ｉＰａｒｔｉｃｌｅｓＣｏｕｎｔＰｅｒＳｅｇ[０] < ｉＰａｒｔｉｃｌｅｓＰｅｒＮｏｄｅＩｄｅａｌであれば、ノード０が処理する範囲は小さすぎると分かる。このため、上記過程を反対の方向への処理を行う。ｉＰａｒｔｉｃｌｅｓＣｏｕｎｔＰｅｒＳｅｇ[０]は、ｉＰａｒｔｉｃｌｅｓＰｅｒＮｏｄｅＩｄｅａｌと同じあるいは近接になった際に、ｉＤｉｖｉｄｅｄＲｅｓｕｌｔ[０], ｉＤｉｖｉｄｅｄＲｅｓｕｌｔ[０＊２+１]は、ノード０の算出範囲となる。
（３）の過程を繰り返して行い、全ての分段に対して処理を行った後、各ノードの処理するグリッドの範囲を取得することができる。
（４）ＲＯＯＴノードは、ＭＰＩ_ＢＣａｓｔ（）関数を呼び出して、分割結果を全てのノードにブロードキャストする。
ステップ４：
ＭＰＩインターフェースプロトコルを用いて、データが分割された粒子を各ノード間に遷移させる。
各ノードは、グリッドの分割結果ｉＤｉｖｉｄｅｄＲｅｓｕｌｔに応じて、ｉＳｅｎｄＧｒｉｄＩｎｆｏアレイ及びｉＳｅｎｄＰａｒｔｉｃｌｅｓＯｆｆｓｅｔアレイを確定する。アレイｉＳｅｎｄＧｒｉｄＩｎｆｏ及びｉＳｅｎｄＰａｒｔｉｃｌｅｓＯｆｆｓｅｔの大きさは、グリッド全体の数と同じであり、その中、ｉＳｅｎｄＧｒｉｄＩｎｆｏは、各グリッドがどのノードに位置するかを記録するものである。ｉＳｅｎｄＰａｒｔｉｃｌｅｓＯｆｆｓｅｔは、各グリッドにおいて１番目の粒子が粒子アレイに位置する位置を記録するものである。
連結リストｇｒｉｄＩｎｆｏの長さに応じて、現在のノードが粒子をｉＳｅｎｄＮｏｄｅＣｏｕｎｔ個のノードに発信することを確定し、発信情報をアレイｉＳｅｎｄＩｎｆｏに書き込む。アレイｉＳｅｎｄＩｎｆｏの長さはｉＳｅｎｄＮｏｄｅＣｏｕｎｔ＊３である。ここで、ｉＳｅｎｄＩｎｆｏ[ｉ＊３]は、受信粒子のノードの番号ＰＩＤＲであり、ｉＳｅｎｄＩｎｆｏ[ｉ＊３+１]は発信粒子の個数であり、ｉＳｅｎｄＩｎｆｏ[ｉ＊３+２]は、発信ノードの番号ＰＩＤＳである。
ＲＯＯＴノードは、ＭＰＩ_Ｇａｔｈｅｒｖ（）関数を呼び出して、全てのノードのｉＳｅｎｄＩｎｆｏアレイをｉＧｌｏｂａｌＳｅｎｄＩｎｆｏアレイに集めさせる。ｉＧｌｏｂａｌＳｅｎｄＩｎｆｏ[ｉ＊３]の値に応じて、小さい順にソートし、更にＭＰＩ_Ｓｃａｔｔｅｒｖ（）関数を呼び出し、ｉＧｌｏｂａｌＳｅｎｄＩｎｆｏ[ｉ＊３]の値に応じて、トリプルを、対応のノードに発信する。
各ノードは、ＲＯＯＴから発信されたトリプルを受信し、アレイｉＲｅｃｖＩｎｆｏに格納してから、粒子の発信及び受信を開始する。
ステップ５：ステップ３で取得した各ノードが制御する粒子に応じて、ＧＰＵにおいてＯｖｅｒｌａｐ領域を算出し、データをＣＰＵメモリへコピーする。そして、ＭＰＩインターフェースプロトコルに基づいてデータのやり取りを行う。
三次元のＤＥＭは、算出過程において、各グリッドが、隣接の２６個のグリッド（ｏｖｅｒｌａｐグリッド）における粒子データを必要となっているため、各ノードのグリッド算出範囲及び伝送粒子を分割し直した後で、各ノードは、算出が正確に行うことを確保するように、ｏｖｅｒｌａｐグリッドを必ず取得する。Ｏｖｅｒｌａｐのやり取り過程は、以下のように実現される。
受信した粒子を粒子アレイに格納すると共に、発信した粒子を粒子アレイから削除する。位置するグリッドの番号に従って、新たな粒子アレイを小さい順にソートして、ｉＣｅｌｌＣｏｕｎｔ及びｉＳｅｎｄＰａｒｔｉｃｌｅｓＯｆｆｓｅｔアレイを算出し直す。
ｉＤｉｖｉｄｅｄＲｅｓｕｌｔアレイが記録された現在のノード処理グリッド範囲に応じて、それぞれの範囲内のグリッドに隣接する隣接グリッドを算出し、現在のノードに位置しない隣接グリッドの番号及びそれが位置するノードの番号を確定する。
ＲＯＯＴノードは、ＭＰＩ_Ｇａｔｈｅｒｖを呼び出して、各ノードのｉＳｅｎｄＩｎｆｏアレイをＲＯＯＴノードのｉＧｌｏｂａｌＳｅｎｄＩｎｆｏアレイに集めさせる。ｉＧｌｏｂａｌＳｅｎｄＩｎｆｏ[ｉ＊３]に従って小さい順にソートした後で、ＭＰＩ_Ｓｃａｔｔｅｒｖ（）関数を呼び出し、ｉＧｌｏｂａｌＳｅｎｄＩｎｆｏ[ｉ＊３]の値に応じて、トリプルを、対応のノードに発信する。
各ノードは、ＲＯＯＴから発信されたトリプルを、アレイｉＲｅｃｖＩｎｆｏに格納させる。ｉＣｅｌｌＣｏｕｎｔ[ｉＲｅｃｖＩｎｆｏ[ｉ＊３+１]]に応じて、何個の粒子を、番号がｉＲｅｃｖＩｎｆｏ[ｉ＊３+２]であるノードに発信するかを確定すると共に、ｉＳｅｎｄＧｒｉｄＩｎｆｏ[ｉＲｅｃｖＩｎｆｏ[ｉ＊３＋１]]=ｉＲｅｃｖＩｎｆｏ[ｉ＊３＋２]とさせる。
ステップ２中の方法を用いて、ｏｖｅｒｌａｐグリッドにおける粒子を、指定されたノードに発信する。
ステップ６：各計算ノードのＧＰＵにおける１つのストリーミングプロセッサにおいて１つの粒子の処理を行う。粒子の座標に応じて、各粒子が位置するグリッドの番号を算出する。
グリッドの番号は、記憶空間を節約するために、行毎に１次元に記憶される。ＣＵＤＡコア関数
calcHash<<<gridsize, blocksize>>> (ParticleHash, ParticleIndex,
x1, x2, x3,
NumParticles);
を呼び出して、粒子が位置するグリッドの番号ＰａｒｔｉｃｌｅＨａｓｈが取得される。算出領域外の粒子に対して、その粒子が位置するグリッドを算出する時、それを算出領域内のあるグリッドに人為的に入れて、算出に影響しない。
そして、Ｃｅｌｌ−ｌｉｓｔの条件に従って、以下のｋｅｒｎｅｌを用いて、ＰａｒｔｉｃｌｅＨａｓｈによりｃｅｌｌ−ｌｉｓｔを生成する。
CalcCellStartEnd<<<gridsize, blocksize>>> (cellStart, cellEtart,
ParticleHash, ParticleIndex,
NumParticles)
上記の結果に応じて、以下のｋｅｒｎｅｌ関数、即ち、
nbrlstgen<<<gridsize, blocksize>>>(NbrLst, NbrLstcnt,
x1, x2, x3,
ParticleIndex, ParticleHash,
CellStart, CellEnd, NumParticles);
を呼び出して、各粒子の隣接リストＮｂｒＬｓｔを生成する。新たに生成したＮｂｒＬｓｔによって、新たな接線相対変位Ｕを算出する。
ステップ７：各計算ノードのＧＰＵにおける１つのストリーミングプロセッサは、１つの粒子の処理を行い、その粒子のストレス及び加速度を算出する。
ステップ６で取得したＮｂｒＬｓｔ及びＵに応じて、粒子の座標、速度、角速度と合わせ、ＤＥＭ方法の条件に従って、各粒子のストレス及びトルクを算出する。ニュートンの第二法則に従って、各粒子の加速度及び角加速度を算出する。
ステップ８：ステップ７で算出した加速度及び角加速度に応じて、１／２歩の粒子の速度を更新する。具体的な方式はステップ２と同じである。
ステップ９：条件を満たすまでに、ステップ２に戻って循環し、引き続きの算出を行う。
ステップ１０：ＧＰＵ装置のメモリに必要なデータをＣＰＵメモリへコピーして、マスターノード及び計算ノードの記憶空間を釈放する。
以下の表１には、上記シミュレーション方法で実行した結果が示されている。プログラムは、ｎＶＩＤＩＡのＧＰＵにおいて異なる歩数で運行される。なお、異なるｂｌｏｃｋ及びＴｈｒｅａｄの数をそれぞれ採用して実行される。
Figure 2015530636
図３は、本発明の別の実施例に係る、ＧＰＵに基づく粒子流動シミュレーションシステムのモジュール構造を示す模式図である。図３に示すように、該モジュール化されたシミュレーションシステムは、モデリングモジュール３０２と、タスク管理モジュール３０４と、計算モジュール３０６と表示モジュール３０８とを含む。図１を参照して、例えば、モデリングモジュール３０２は、前端サーバ１０で実現されることができ、タスク管理モジュール３０４は、管理ノード３０で実現されることができ、計算モジュール３０６は、計算ノード４０のクラスタで実現されることができ、表示モジュール３０８は、後端サーバ２０で実現されることができる。しかしながら、これらのモジュールは、適宜な方式で、例えば１つ又は複数のコンピューターで実現されることもできる。
モデリングモジュール３０２は、粒子を生成するために必要な情報、例えば、粒子の数、大きさ、材料などの情報（ヤング率、ポアソン(Ｐｏｉｓｓｏｎ)比、密度、回復係数など）及び粒子の分布範囲、摩擦係数、境界条件などのパラメーターを受信し、粒子と接触する幾何体の材料の情報を提供する。
モデリングモジュール３０２は、受信した情報に応じて、必要な粒子モデル（単に「粒子」ということもできる）を生成する。生成した粒子間に重畳しない、あるいは、重畳が小さいことを確保するために、以下のいくつかの種類の方法を用いて粒子モデルを生成することが可能である。（１）規則生成法、すなわち、所定の範囲内で規則的な粒子を生成する。ただし、粒子の半径の０.１％〜１％に相当する変動を加える必要がある。（２）１つの粒子を生成する毎に、その粒子を、以前に生成された全ての粒子と比較し、重畳するか否かを検出する。若し、重畳すれば、その粒子を生成し直すことになる。そうでなければ、生成が成功することと見なす。（３）まず、小さい空間内で方法（２）を用いていくつかの粒子を生成し、そして、粒子の数の条件を満たすまでに、これらの粒子をコピーして平行移動させ、他の空間を充填する。これは、粒子分布のランダム性を向上できると共に、算出の時間を節約することもできる。上記３つの方法以外に、粒子の数が比較的に少ない場合について、空間の範囲が確定された後で、交互的な方法によって、マウスでクリックして生成してもよい。
粒子が生成された後で、モデリングモジュール３０２は、幾何体の情報に対する処理を行う。幾何体を有限的な曲面に分解し、これらの曲面に番号をつける。次に、生成された粒子、幾何体及び他の材料の情報をタスク管理モジュール３０４に供給する。
タスク管理モジュール３０４は、まず、伝送される粒子の数及びフリーなＧＰＵの数に応じて、現在のタスクに対してノード及びＧＰＵを割り当てる。若し、リソースが不足であれば、それをユーザに通知し、あるいは、待ちや放棄をユーザに選択させる。ＧＰＵを確定した後で、初期の粒子の位置情報を管理ノード３０のＧＰＵに記憶し、ＧＰＵの数及び粒子が空間中の分布状況に応じてどれらの粒子がどの計算ノード４０のどのＧＰＵカードによって算出されるかを確定する。タスク管理モジュール３０４は、確定した結果を、計算モジュール３０６へ伝送し、各計算ノード４０に割り当てる。
各計算ノード４０が自身に必要な粒子を取得した後で、まず、現在の加速度に応じて１／２歩を積分し、１／２歩後の速度を取得し、そして、この速度及び現在の粒子の座標値に応じて全ての粒子の位置の更新を行う。
位置を更新した後で衝突の検出を行う。このとき、空間をいくつかのグリッドに分割する必要がある。いずれか一つの粒子のストレス状況を算出する時、該粒子と隣接するグリッド内の粒子がその粒子に衝突するか否かの算出のみを行えば良い。若し、衝突が生じれば、衝突粒子を衝突リストに入れ、衝突粒子の個数に１を加算する。
粒子ボールのストレスを算出する時、まず、該衝突粒子の座標、速度、角速度の情報を抽出して、ストレスの算出を行う。その後、全ての衝突粒子に対して合力を求めると共に、粒子の加速度を算出する。粒子周辺の幾何体のストレスについて、まず、粒子と幾何体との間の距離を算出し、該距離が粒子の半径よりも小さい場合、該粒子が幾何体に衝突していると見なす。幾何体を、質量が無限大であり且つ速度及び角速度場が０である粒子とし、粒子が幾何体から受ける力を同様に算出することができる。
中断後で引き続きの算出を保証するために、実際な需要に応じて、一歩の算出データを一時間毎に格納することができる。該計算モジュール３０６は、需要に応じて、堆積係数、平均堆積密度、温度粘性係数などの物理量を算出して記憶してもよい。算出完成後、若し、ユーザが結果を可視化にしたければ、データを表示モジュール３０８に発信することができる。
以下、図４を参照して、計算モジュール３０６の操作フローを記述する。該実施例において、計算モジュール３０６の算出過程は、「ソートされるセルリスト」を採用することができる。該方法は、全ての粒子に対して、粒子が位置するグリッドに従ってソートし、ｃｅｌｌＳｔａｒｔ及びＣｅｌｌＥｎｄとの２つのアレイの優勢を十分に利用する。該方法は、構造が簡単で、実現しやすく、効率が高いという特徴を有する。そのため、該方法は、各種類の高密度の粒子衝突に適用し、粒子の高い速度によるクロス・ノードの伝送という問題を解決することができる。
粒子を記述する物理量は、座標ｐｏｓ、速度ｖｅｌ、角速度ｗ、加速度ａ、角加速度ｂｅｔａ、粒子の接線相対変位Ｕを有する。これらの変数は、いずれも三次元の変数である。また、粒子が位置するグリッドの番号ｈａｓｈ、粒子の永久全体番号ｐｉｄ及び一時局所番号ｉｎｄｅｘ、粒子の衝突リストＣｏｌｌｉｄｅＬｉｓｔ、及び衝突の粒子数ＣｏｌｌｉｄｅＬｉｓｔＣｎｔを更に有する。
セルとは、上記分割によって取得したグリッドである。本明細書では、「セル」と「グリッド」との意味は同じであり、両者を互換して使用することができる。セルｉを記述する変数は、ｃｅｌｌＳｔａｒｔ[ｉ]、ｃｅｌｌＥｎｄ[ｉ]、ｃｅｌｌＣｏｕｎｔ[ｉ]を有し、ただし、ｉはセルの番号を示し、ｃｅｌｌＳｔａｒｔ[ｉ]はセルｉの開始粒子の番号を示し、ｃｅｌｌＥｎｄ[ｉ]は、セルｉの終了粒子の番号を示し、ｃｅｌｌＣｏｕｎｔ[ｉ]は、セルｉの粒子総数を示す。
コース通信を記述するための二次元のアレイは、ＰａｒｔｉｃｌｅｓＳｅｎｄＴｏＥａｃｈＮｏｄｅと称しても良い。ｉ行目ｊ列目のエレメント[ｉ][ｊ]は、ｉ番目のノードからｊ番目のノードへ発信する粒子の総数を示す。
本発明の採用した時間積分アルゴリズムは、速度ｖｅｒｌｅｔアルゴリズムである（従来の積分アルゴリズムであり、例えば、http://en.wikipedia.org/wiki/Verlet_integrationを参照する）。
図４に示すように、ステップ４０１において、初期化を行う。ＧＰＵ及びＣＰＵの記憶空間を設け、算出した粒子の情報を各計算ノードのＧＰＵに発信することを含む。
ステップ４０２において、予め定められた速度及び座標を更新する。例えば、加速度（又は角加速度）に応じて１／２歩の速度（又は角速度）を更新した直後、速度に応じて粒子の座標の更新を行う。以下の式に示すようになる。
Figure 2015530636
以上の２つのステップは、いずれもそれぞれの計算ノードのＧＰＵにおいて並列に完成されるものである。ＧＰＵ中の一つのスレッド（ｔｈｒｅａｄ）は、一つの粒子に対応し、ＧＰＵの最も高い効率に達した。
このように、新たな座標を取得した。新たな座標及び新たな速度（角速度）における加速度（角加速度）を算出する必要がある。
粒子の座標が変わったため、もともとＡコース（又はＧＰＵ）で算出すべきの粒子は、このときにＢコースで算出すべきとなる可能性がある。このように、Ａコースの該粒子の全ての情報を、Ｂコースに発信する必要がある。
まず、各計算ノードのＧＰＵにおいて、各粒子が位置するグリッドの番号Ｈａｓｈを算出する。各粒子が位置するグリッドの番号Ｈａｓｈ及び粒子の局所の自然番号ｉｎｄｅｘでｋｅｙ−ｖａｌｕｅのソートを行う。このステップは、ｔｈｒｕｓｔライブラリ（従来熟したライブラリであり、ｃｕｄａに集積され、例えばhttp://code.google.com/p/thrust/を参照する）で完成する。ソートされたｈａｓｈに応じて、ＧＰＵにおいて、並列に算出を行い、各グリッドｉのｃｅｌｌＳｔａｒｔ[ｉ]、ｃｅｌｌＥｎｄ[ｉ]及びｃｅｌｌＣｏｕｎｔ[ｉ]を取得する。すなわち、ステップ４０３を実行する。
ソートのｉｎｄｅｘに従って、粒子の全ての物理量のソートを行う。
ここまで、粒子が位置するグリッドの番号に従って、粒子の全ての物理量をソートし直し、各グリッドｉのｃｅｌｌＳｔａｒｔ[ｉ]、ｃｅｌｌＥｎｄ[ｉ]、ｃｅｌｌＣｏｕｎｔ[ｉ]は「ソートされるセルリスト」と総称する。
そして、ステップ４０４において、動的分割を行う。具体的には、各計算ノードは、自身が有する粒子のグリッド及び粒子の数を、複数の計算ノードにおけるマスターノードに発信する。すなわち、各計算ノードにおいて、ｃｅｌｌＣｏｕｎｔ[ｉ]!=０のとき、ｉ及びｃｅｌｌＣｏｕｎｔ[ｉ]をマスターノードに発信する。それぞれの計算ノードが発信したｃｅｌｌＣｏｕｎｔ[ｉ]を、マスターノードによって累積して、空間全体のｃｅｌｌＣｏｕｎｔ[ｉ]を取得する。マスターノードは、空間全体のｃｅｌｌＣｏｕｎｔ[ｉ]に応じて、それぞれのＧＰＵの算出する粒子を分割し直す。分割の原則は、グリッドを単位とし、それぞれのＧＰＵがいずれも連続的なグリッドを算出し、且つ、グリッドにおける粒子の総数が、各ＧＰＵの平均の粒子の数に近接する、ということである。このように、それぞれのＧＰＵはいずれも、粒子の座標変化による粒子の算出範囲を取得した。
新たな算出範囲及び現在の各ＧＰＵの算出範囲に応じて、関連の粒子情報を送受信する。ＧＰＵが送受信する必要とする粒子の総数を確定するために、二次元のアレイ：ＰａｒｔｉｃｌｅｓＳｅｎｄＴｏＥａｃｈＮｏｄｅを作成する。該アレイのそれぞれの一次元の大きさは、いずれもコースの数（又はＧＰＵの数）である。ＰａｒｔｉｃｌｅｓＳｅｎｄＴｏＥａｃｈＮｏｄｅ[ｉ][ｊ]の意味は、ｉ番目のＧＰＵがｊ番目のＧＰＵへ発信する必要とする粒子の総数であり、すなわち、ｊ番目のＧＰＵがｉ番目のＧＰＵから受信する粒子の総数である。該アレイの対角線のエレメントは、全てゼローである。該アレイのｉ行目に対して求めた和は、ｉ番目のＧＰＵが発信する粒子の総数である。ｊ列目に対して求めた和は、ｊ番目のＧＰＵが受信する粒子の総数である。ｃｅｌｌＳｔａｒｔ及びｃｅｌｌＣｏｕｎｔを入力として、アレイＰａｒｔｉｃｌｅｓＳｅｎｄＴｏＥａｃｈＮｏｄｅを算出する。その同時に、ＳｅｎｄＳｔａｒｔを算出する。ＳｅｎｄＳｔａｒｔも二次元のアレイであり、ＳｅｎｄＳｔａｒｔ[ｉ][ｊ]は、ｉ番目のＧＰＵがｊ番目のＧＰＵへ発信する一番目の粒子のアレイにおける位置である。このように、発信のために、発信しようとする粒子の情報をＧＰＵから取得して、発信粒子のバッファーに伝送することができる。次に、受信のために、アレイの列に対して和を求めることによって、それぞれのＧＰＵが受信する粒子の総数を確定でき、対応のバッファーを設けることができる。全ての送受信が完成するまでに、ＭＰＩの標準関数における、例えば、非同期の送受信方式ＭＰＩ_Ｉｒｅｃｖ関数及びＭＰＩ_Ｉｓｅｎｄ関数などによって、対応の粒子の物理情報を送受信する。
ｃｕｄａＭｅｍｃｐｙＨｏｓｔＴｏＤｅｖｉｃｅ関数（既知の関数であり、ＧＰＵにおいてホストメモリとのやり取りデータを記憶する）によって、受信したアレイを、直接にＧＰＵの各アレイの末端へコピーし、送受信バッファーを釈放する。
このとき、それぞれのＧＰＵに対して算出すべきの新たな粒子の情報は、全て取得されたが、新たに加入した粒子及び発信した粒子を考慮する必要があり、「ソートされるセルリスト」を算出し直すことによって、ソートされた物理量のアレイを取得することができる。
各ＧＰＵの算出する粒子が独立せず、すなわち、ＧＰＵ間に重畳（Ｏｖｅｒｌａｐ）領域があるため、ステップ４０５において、各ＧＰＵの算出するグリッドの番号に応じて、該ＧＰＵが必要なＯｖｅｒｌａｐ領域を算出することができる。動的分割と類似する方法を用いて、それぞれのＧＰＵは、必要なＯｖｅｒｌａｐ領域の粒子の物理情報を取得して、それぞれのアレイの末端に記憶する。このように、Ｏｖｅｒｌａｐ領域を加えた物理情報のアレイは、完全的にソートされていないが、同一のグリッドにおける粒子は、連続的に記憶されている。その同時、それぞれのグリッドのｃｅｌｌＳｔａｒｔ及びｃｅｌｌＥｎｄを算出する。
ステップ４０６において、粒子の情報及びｃｅｌｌＳｔａｒｔ、ｃｅｌｌＥｎｄに応じて、現在の全ての粒子の衝突リストを算出する。その方法は、以下のようである。即ち、いずれか一つの粒子ｉに対して、まず、ｔｅｘｔｕｒｅ ｍｅｍｏｒｙ（テクスチャーメモリ）によってその座標を取得し、それが位置するグリッドの番号を算出し、その自身を含む周辺の２７つのグリッドにおける他の全ての粒子をスキャンする。若し、他の粒子と該粒子とのセントロイド距離が両者の半径の合計よりも小さければ、この粒子を、該粒子の衝突リストＣｏｌｌｉｄｅＬｉｓｔ[ｉ][ＣｏｌｌｉｄｅＬｉｓｔＣｎｔ[ｉ]]にマークして、衝突リストの数ＣｏｌｌｉｄｅＬｉｓｔＣｎｔ[ｉ]に１を加算する。
接線相対変位は、２つの粒子が接触するときのみに存在する。現在時刻のいずれか一つの粒子ｉのストレスを算出するために、一つ前の時刻の接線相対変位を必要とする。該接線相対変位を記憶するアレイＵの次元の大きさは、ＣｏｌｌｉｄｅＬｉｓｔの次元の大きさと同じである。Ｕ[ｉ][ｊ]は、粒子ｉと粒子ＣｏｌｌｉｄｅＬｉｓｔ[ｉ][ｊ]との接線相対変位を記憶する。したがって、算出結果の正確性を確保するために、粒子の現在時刻のストレスを算出する前に、現在時刻のＣｏｌｌｉｄｅＬｉｓｔ、一つ前の時刻のＣｏｌｌｉｄｅＬｉｓｔＯｌｄ及びＵＯｌｄに応じて、アレイＵをソートし直しなければならない。このソート過程は、ＧＰＵにおいて実現される。具体的には、入力された一つ前の時刻の衝突リストＣｏｌｌｉｄｅＬｉｓｔＯｌｄと、ＣｏｌｌｉｄｅＬｉｓｔＣｎｔと、ＣｏｌｌｉｄｅＬｉｓｔＯｌｄに対応するＵＯｌｄとを用い、現在時刻の衝突リストＣｏｏｌｉｄｅＬｉｓｔ[ＣｏｌｌｉｄｅＬｉｓｔＣｎｔ]を入力とし、ＵＯｌｄの順序を調整して、現在時刻のアレイＵを取得する。
このように、力を算出するための全ての正しいアレイを取得した。ステップ４０７において、ＨＭ接触力学モデルによって、それぞれの粒子のストレスを算出する。具体的には、座標ｐｏｓ、速度ｖｅｌ、角速度ｗ、粒子の接線相対変位Ｕ、衝突リストＣｏｌｌｉｄｅＬｉｓｔ[ＣｏｌｌｉｄｅＬｉｓｔＣｎｔ]を用いて、ＨＭ接触力学の式に従って、それぞれの粒子の加速度ａ及び角加速度ｂｅｔａを算出することができる。
新な加速度ａ（角加速度ｂ）を取得した後で、ステップ４０８において、以上の速度に従って１／２歩の速度を再びに更新する。
ここまで、計算モジュールにおける完全な一歩の演算を完成した。
現在の全ての粒子の物理情報のアレイを格納し、次回のアレイのために準備する。ここで、ステップ４０９においては、コピー又はポインターやり取り技術を採用することができる。ポインターやり取り技術は、現在のアレイと次回に算出するアレイとの最初のアドレスをやり取りして、データのコピーが必要な比較的に長い時間を低減することができる。
ステップ４１０において、外部への記憶を行うか否かを判断する。必要であれば、ステップ４１１において、全ての粒子の全ての物理情報を外部の記憶装置に格納して、停電した後で算出し直すというリスクを防止することができる。ステップ４１２において、統計するか否かを判断する。必要であれば、ステップ４１３において、例えば、平均値、分散などの関連的な統計物理量を算出する。ステップ４１４において、算出の終了条件を満たすか否かを判断する。例えば、予め定められた回数の算出を実行したか否かを判断する。算出が完成していなければ、ステップ４０２に戻す。そうでなければ、算出を終了させ、結果を格納して、記憶の空間を釈放する。
国際的に著名なソフトウェアｌａｍｍｐｓ（広く適用されるオープンソース・ソフトウェアであり、http://lammps.sandia.gov/を参照することができる）の８コアのＣＰＵに基づく実施と比べると、本発明のＧＰＵ （例えば、Ｔｅｌｓａ Ｍ２０９０）に基づくシミュレーション方法の演算速度は、１０倍ほど向上できた。
当業者は、本発明の主旨及び範囲を逸脱しない限り、本発明に対する変更や変形をすることができる。このように、本発明のこれらの補正及び変形が本発明の特許請求の範囲及びそれと同様な技術範囲に属すれば、本発明もこれらの変更及び変形を含む。

Claims (20)

  1. 並列の複数のＧＰＵ上で離散要素法（ＤＥＭ）方法を実行して粒子流動のシミュレーションを行うＧＰＵに基づく粒子流動のシミュレーション方法であって、
    ＤＥＭ方法で粒子をモデリングし、作られたＤＥＭモデルを複数の粒子として割り当てし、該複数の粒子を複数の計算ノードに割り当て処理を行い、各計算ノードのＣＰＵ及びＧＰＵに記憶空間がそれぞれ割り当てされており、ＣＰＵにおいてデータの初期化を行い、初期化されたデータをＣＰＵの記憶空間から前記ＧＰＵの記憶空間へコピーするステップａと、
    上記各計算ノードのＧＰＵは、各粒子の処理を行い、各計算ノードのＧＰＵの各ストリーミングプロセッサは、１つの粒子を処理し、ＧＰＵの記憶空間に記憶される粒子の座標及び粒子の速度を更新するステップｂと、
    ステップｂの処理過程において、各計算ノードが制御する粒子を確定し、各計算ノードが制御する粒子の個数をＣＰＵの記憶空間へコピーして、各計算ノードがどれらの粒子を算出するかを均衡負荷の原則に従って動的に確定できるように、ＧＰＵの記憶空間における粒子の数に応じて動的分割を行うステップｃと、
    ＭＰＩインターフェースプロトコルによって、上記データが動的分割された粒子を各計算ノード間に遷移させるステップｄと、
    ステップｃで取得した各計算ノードが制御する粒子に応じて、ＧＰＵにおいて重畳領域を算出し、データをＣＰＵメモリへコピーしてから、ＭＰＩインターフェースプロトコルによってデータのやり取りを行うステップeと、
    各計算ノードのＧＰＵにおける各ストリーミングプロセッサは、各粒子の座標に応じて、各粒子が位置するＧＰＵの記憶空間におけるグリッドの番号を算出するステップｆと、
    各計算ノードのＧＰＵにおける各ストリーミングプロセッサは、各粒子の運動中のストレス及び加速度を処理して算出するステップｇと、
    各計算ノードのＧＰＵにおける各ストリーミングプロセッサは、各粒子の速度を処理するステップｈと、
    指定された歩数に達するまでに、ステップｂに戻すステップｉと、
    マスターノード及び計算ノードの記憶空間を釈放するステップｊと、を含むことを特徴とするＧＰＵに基づく粒子流動のシミュレーション方法。
  2. ステップｂ、ステップｆ、ステップｇ及びステップｈは、
    ＧＰＵにおいて各粒子に対して並列なデータ処理を行う請求項１に記載のＧＰＵに基づく粒子流動のシミュレーション方法。
  3. ステップｄにおいて、前記粒子が各ノード間で遷移し、ＭＰＩインターフェースで関数を送受信し、粒子の各物理量の送受信を実現する請求項１に記載のＧＰＵに基づく粒子流動のシミュレーション方法。
  4. ステップｅにおいて、前記ＧＰＵにおいて重畳領域を算出することは、
    ＧＰＵにおいて重畳領域を算出することを利用し、ＧＰＵの１つのストリーミングプロセッサが１つのグリッドを処理することを含み、
    三次元の場合において、各グリッドは２６個のグリッドに隣接し、隣接のグリッドが現在の計算ノードに位置するか否かを判断し、位置しなければ、重畳領域とし、他のノードから遷移して取得する請求項１に記載のＧＰＵに基づく粒子流動のシミュレーション方法。
  5. ＧＰＵに基づく粒子流動のシミュレーション方法であって、
    粒子の材料、粒子のパラメーター、境界の条件、幾何体の形状、及び粒子の初期分布の領域を確定し、予め定められた粒子の分布領域及び数量に応じて粒子を生成するモデリングステップと、
    粒子の総数及び複数の計算ノード上にフリーなＧＰＵ数に応じて、最適なＧＰＵ数を確定し、最適なＧＰＵ数及び現在フリーなＧＰＵ数に応じて、算出関与のＧＰＵを確定し、算出関与のＧＰＵの状態をビジーに設置するタスク管理ステップと、
    算出ステップと、を含み、
    該算出ステップは、
    各計算ノードの算出関与のＧＰＵを初期化し、算出に必要な粒子の情報を各ＧＰＵに発信するステップと、
    各ＧＰＵが、予め定められた速度及び座標を並列に更新し、受信した粒子の情報をソートして各自のソートセルリストを生成するステップと、
    各ＧＰＵが、現在各自のコースにおける非ゼローのグリッドの番号及びグリッドにおける粒子数を並列に算出し、マスターノードに発信して、マスターノードによって各ＧＰＵの最適な粒子数に応じて、グリッドを動的分割し、各ＧＰＵが並列に算出するグリッドの数及びグリッドの番号を確定するステップと、
    マスターノードが確定した結果に応じて、各ＧＰＵが粒子情報を並列に送受信し、各ＧＰＵにおいて各自のソートセルリストを生成し直すステップと、
    各ＧＰＵにおいて現在時刻の衝突リストを生成し、現在時刻の衝突リストと一つ前の時刻の衝突リストと接線相対変位に応じて、各ＧＰＵにおいて接線相対変位の位置を並列に調整するステップと、
    接触力学モデルによって、各ＧＰＵにおいて各粒子のストレス及び加速度を並列に算出するステップと、
    現在の算出結果を記憶するステップと、
    算出が完成していなければ、各ＧＰＵが予め定められた速度及び座標を並列に更新するステップに戻し、そうでなければ、算出ステップを終了させるステップと、を含むことを特徴とするＧＰＵに基づく粒子流動のシミュレーション方法。
  6. 更に、表示ステップを含み、
    上記表示ステップは、
    境界の条件を確定し、幾何体の境界を透明な曲面で作るステップと、
    粒子の位置及び粒子の直径に応じて、粒子を同じ色又は異なる色のペレットで描くステップと、
    階調画像でスカラ場を表示し、粒子情報を重み付けしてグリッドにマッピングすることによって、ベクトル場を流線描き方法で描くステップとを含む請求項５に記載のＧＰＵに基づく粒子流動のシミュレーション方法。
  7. 算出結果とする全ての粒子情報を外部の記憶装置に格納するステップとを更に含む請求項５に記載のＧＰＵに基づく粒子流動のシミュレーション方法。
  8. 各ＧＰＵが、粒子に関する物理統計量を並列に算出することを更に含む請求項５に記載のＧＰＵに基づく粒子流動のシミュレーション方法。
  9. 予め定められた粒子の分布領域及び数量に応じて粒子を生成することは、
    粒子の数量条件を満たすまでに、比較的に小さい空間内でいくつかの粒子を生成して、これらの粒子を平行移動してコピーを行い、他の空間を充填することを含む請求項５に記載のＧＰＵに基づく粒子流動のシミュレーション方法。
  10. ソートセルリストは、粒子が位置するグリッドに応じて、全ての粒子をソートするリストである請求項５に記載のＧＰＵに基づく粒子流動のシミュレーション方法。
  11. 動的分割方法を採用して、非ゼローのグリッドの番号及びグリッドにおける粒子数をＧＰＵにおいて並列に算出する請求項５に記載のＧＰＵに基づく粒子流動のシミュレーション方法。
  12. 各ＧＰＵにおいては、
    １つのスレッドが１つの粒子に対応するとの方式を用いて算出を行う請求項５に記載のＧＰＵに基づく粒子流動のシミュレーション方法。
  13. 接線相対変位を算出することは、
    一つ前の時刻の接線相対変位を記録し、現在時刻の衝突リストに応じてそれを更新することを含む請求項５に記載のＧＰＵに基づく粒子流動のシミュレーション方法。
  14. コピー又はポインターやり取り技術を用いて、現在の算出結果をアレイに記憶する請求項５に記載のＧＰＵに基づく粒子流動のシミュレーション方法。
  15. ＧＰＵに基づく粒子流動のシミュレーションシステムであって、
    粒子の材料、粒子のパラメーター、境界の条件、幾何体の形状、及び粒子の初期分布の領域を確定し、予め定められた粒子の分布領域及び数量に応じて粒子を生成するように構成されるモデリングモジュールと、
    粒子の総数及び複数の計算ノードにフリーなＧＰＵ数に応じて、最適なＧＰＵ数を確定し、最適なＧＰＵ数及び現在フリーなＧＰＵ数に応じて、算出関与のＧＰＵを確定し、算出関与のＧＰＵの状態をビジーに設置するように構成されるタスク管理モジュールと、
    算出モジュールと、を含み、
    該算出モジュールは、
    各計算ノードの算出関与のＧＰＵを初期化し、算出に必要な粒子の情報を各ＧＰＵに発信し、
    各ＧＰＵが、予め定められた速度及び座標を並列に更新し、受信した粒子の情報をソートして各自のソートセルリストを生成し、
    各ＧＰＵが、現在各自のコースにおける非ゼローのグリッドの番号及びグリッドにおける粒子数を並列に算出し、マスターノードに発信して、マスターノードによって各ＧＰＵの最適な粒子数に応じて、グリッドを動的分割し、各ＧＰＵが並列に算出するグリッドの数及びグリッドの番号を確定し、
    マスターノードが確定した結果に応じて、各ＧＰＵが粒子情報を並列に送受信し、各ＧＰＵにおいて各自のソートセルリストを生成し直し、
    各ＧＰＵにおいて現在時刻の衝突リストを生成し、現在時刻の衝突リストと一つ前の時刻の衝突リストと接線相対変位に応じて、各ＧＰＵにおいて接線相対変位の位置を並列に調整し、
    接触力学モデルによって、各ＧＰＵにおいて各粒子のストレス及び加速度を並列に算出し、
    現在の算出結果を記憶し、
    算出が完成していなければ、各ＧＰＵが予め定められた速度及び座標を並列に更新するステップに戻し、そうでなければ、算出ステップを終了させるように構成されることを特徴とするＧＰＵに基づく粒子流動のシミュレーションシステム。
  16. 更に、表示モジュールを含み、
    上記表示モジュールは、
    境界の条件を確定し、幾何体の境界を透明な曲面で作り、
    粒子の位置及び粒子の直径に応じて、粒子を同じ色又は異なる色のペレットで描き、
    階調画像でスカラ場を表示し、粒子情報を重み付けしてグリッドにマッピングすることによって、ベクトル場を流線描き方法で描くように構成される請求項１５に記載のＧＰＵに基づく粒子流動のシミュレーションシステム。
  17. ＧＰＵに基づく粒子流動のシミュレーションシステムであって、
    クライアントから入力される粒子のモデリング情報に応じて、粒子の情報を生成すると共に、幾何体の情報を生成するように構成される前端サーバと、
    前端サーバから粒子の情報及び幾何体の情報を受信し、粒子の数及び各計算ノード上にフリーなＧＰＵの数に応じて、どれらの計算ノードにおけるどれらのＧＰＵを使用するかを確定し、そして、確定したＧＰＵの数及び粒子の空間における分布状況に応じてどれらの粒子がどの計算ノードのどのＧＰＵによって算出されるかを確定し、確定した結果によって割り当てるように構成される管理ノードと、
    それぞれが複数のＧＰＵを含み、複数のＧＰＵにおいて粒子の衝突による各粒子のストレスを並列に算出し、加速度を更に算出して、粒子の流動をシミュレーションするように構成される複数の計算ノードと、
    シミュレーションの結果を表示するように構成される後端サーバと、を備えることを特徴とするＧＰＵに基づく粒子流動のシミュレーションシステム。
  18. 前端サーバは、
    幾何体を有限の曲面に分解し、これらの曲面に番号をつけることによって、幾何体の情報を生成する請求項１７に記載のＧＰＵに基づく粒子流動のシミュレーションシステム。
  19. 後端サーバは、
    表示されるシミュレーション結果において、幾何体の境界を透明な曲面で作り、
    粒子の位置及び粒子の直径に応じて、粒子を同じ色又は異なる色のペレットで描き、
    且つ、階調画像でスカラ場を表示し、粒子情報を重み付けしてグリッドにマッピングすることによって、ベクトル場を流線描き方法で描く請求項１７に記載のＧＰＵに基づく粒子流動のシミュレーションシステム。
  20. 前端サーバ、管理ノード、計算ノード及び後端サーバは、
    ＩＢ（ＩｎｆｉｎｉＢａｎｄ）ネットワークによって通信する請求項１７に記載のＧＰＵに基づく粒子流動のシミュレーションシステム。
JP2015521951A 2012-12-20 2013-05-22 粒子流動のシミュレーションシステム及びその方法 Active JP6009075B2 (ja)

Applications Claiming Priority (3)

Application Number Priority Date Filing Date Title
CN201210558134.6 2012-12-20
CN201210558134 2012-12-20
PCT/CN2013/076045 WO2014094410A1 (zh) 2012-12-20 2013-05-22 颗粒流动仿真系统和方法

Publications (2)

Publication Number Publication Date
JP2015530636A true JP2015530636A (ja) 2015-10-15
JP6009075B2 JP6009075B2 (ja) 2016-10-19

Family

ID=49193522

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2015521951A Active JP6009075B2 (ja) 2012-12-20 2013-05-22 粒子流動のシミュレーションシステム及びその方法

Country Status (5)

Country Link
US (1) US10007742B2 (ja)
JP (1) JP6009075B2 (ja)
CN (1) CN103324780B (ja)
GB (1) GB2523640B (ja)
WO (1) WO2014094410A1 (ja)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2017091262A (ja) * 2015-11-11 2017-05-25 富士通株式会社 粒子シミュレーションプログラム、粒子シミュレーション装置、及び計算機資源配分方法
KR20190048790A (ko) * 2017-10-31 2019-05-09 한국과학기술원 공간 데이터를 분산 처리하는 시스템 및 방법
KR20200023954A (ko) * 2018-08-27 2020-03-06 이에이트 주식회사 복수의 프로세서를 이용한 입자 기반의 유체 시뮬레이션 방법 및 유체 시뮬레이션 장치
JP2021009568A (ja) * 2019-07-01 2021-01-28 富士通株式会社 情報処理装置、粒子シミュレーション方法および粒子シミュレーションシステム
KR20210136327A (ko) * 2020-05-07 2021-11-17 국방과학연구소 도심지 유동해석 방법 및 장치

Families Citing this family (50)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103617088B (zh) * 2013-11-29 2018-07-24 深圳中微电科技有限公司 在不同类型线程中分配内核资源的方法、装置及其处理器
US10210287B2 (en) 2013-12-31 2019-02-19 Disney Enterprises, Inc. Augmented material point method for simulating phase changes and varied materials
US10019827B2 (en) * 2013-12-31 2018-07-10 Disney Enterprises, Inc. Material point method for simulation of granular materials
WO2015108980A1 (en) * 2014-01-17 2015-07-23 Conocophillips Company Advanced parallel "many-core" framework for reservoir simulation
CN105389855B (zh) * 2014-08-26 2019-11-01 三星电子株式会社 对对象进行建模的方法和设备
JP6547547B2 (ja) * 2015-09-25 2019-07-24 富士通株式会社 粒子シミュレーションプログラム、計算機資源配分方法、および粒子シミュレーション装置
CN107016180A (zh) * 2017-03-30 2017-08-04 中国石油大学(北京) 一种颗粒流动仿真方法
US10262390B1 (en) * 2017-04-14 2019-04-16 EMC IP Holding Company LLC Managing access to a resource pool of graphics processing units under fine grain control
US10275851B1 (en) 2017-04-25 2019-04-30 EMC IP Holding Company LLC Checkpointing for GPU-as-a-service in cloud computing environment
CN107230242B (zh) * 2017-06-07 2020-09-25 广州酷狗计算机科技有限公司 粒子映射方法和装置
US10325343B1 (en) 2017-08-04 2019-06-18 EMC IP Holding Company LLC Topology aware grouping and provisioning of GPU resources in GPU-as-a-Service platform
CN107766148B (zh) * 2017-08-31 2021-02-19 北京百度网讯科技有限公司 一种异构集群及任务处理方法和装置
US10698766B2 (en) 2018-04-18 2020-06-30 EMC IP Holding Company LLC Optimization of checkpoint operations for deep learning computing
CN110795228B (zh) 2018-08-03 2023-08-25 伊姆西Ip控股有限责任公司 用于训练深度学习模型的方法和制品、以及计算系统
US10776164B2 (en) 2018-11-30 2020-09-15 EMC IP Holding Company LLC Dynamic composition of data pipeline in accelerator-as-a-service computing environment
CN109753726A (zh) * 2019-01-04 2019-05-14 东南大学 一种基于边界盒搜索方法和gpu的球磨机介质运动仿真方法
CN109977505B (zh) * 2019-03-13 2024-02-09 南京大学 基于gpu矩阵计算的离散元孔隙系统快速搜索方法
CN110275732B (zh) * 2019-05-28 2023-02-21 上海交通大学 在ARMv8处理器上质点网格法的并行实现方法
US11875101B2 (en) * 2019-06-20 2024-01-16 Asml Netherlands B.V. Method for patterning process modelling
CN110348156B (zh) * 2019-07-18 2022-10-14 河南理工大学 一种炉料颗粒在高炉旋转溜槽内部运动的仿真方法
CN112035995B (zh) * 2019-07-19 2024-07-09 交通运输部天津水运工程科学研究所 基于gpu计算技术的非结构网格潮汐潮流数值模拟方法
CN112528456B (zh) * 2019-09-18 2024-05-07 曙光信息产业(北京)有限公司 一种异构节点计算系统及方法
CN110766709B (zh) * 2019-10-30 2021-03-30 成都理工大学 基于图像识别的滑坡颗粒堆积特征识别方法
CN111222262B (zh) * 2019-10-30 2023-09-29 中国中元国际工程有限公司 基于质量比影响的气浮隔振平台性能优化设计方法
CN110992456B (zh) * 2019-11-19 2021-09-07 浙江大学 一种基于位置动力学的雪崩模拟方法
CN111027244B (zh) * 2019-12-03 2024-03-12 天津大学 一种百亿级颗粒模型的构建方法
CN112001108B (zh) * 2020-07-08 2024-02-02 中国人民解放军战略支援部队信息工程大学 锥束ct蒙特卡洛仿真集群并行加速方法及系统
CN112100939B (zh) * 2020-09-14 2023-06-16 福建天晴在线互动科技有限公司 一种基于Compute Shader的实时流体仿真方法及其系统
CN112380788B (zh) * 2020-11-06 2022-03-01 天津大学 一种超椭球颗粒与流场双向耦合的半解析计算方法
CN113177346B (zh) * 2020-11-17 2022-06-10 西北工业大学 一种锅炉煤粉运输安全性判断方法及系统
CN112380793B (zh) * 2020-11-18 2024-02-13 上海交通大学 基于gpu的湍流燃烧数值模拟并行加速实现方法
CN112464543B (zh) * 2021-01-28 2021-04-06 北京科技大学 一种计算vim冶炼过程中的夹杂物运动的方法
CN113283048B (zh) * 2021-03-12 2024-09-24 南京航空航天大学 一种多枝晶运动相场法并行模拟的碰撞检测和合并方法
CN113221200B (zh) * 2021-04-15 2022-10-25 哈尔滨工程大学 一种适用于堆芯颗粒分布不确定性分析的三维高效随机排布方法
CN113239522B (zh) * 2021-04-20 2022-06-28 四川大学 一种基于计算机集群的大气污染物扩散模拟方法
CN113282976B (zh) * 2021-04-30 2023-04-11 重庆大学 基于comsol的粉床构建方法
CN112948643B (zh) * 2021-05-13 2021-08-06 中国空气动力研究与发展中心计算空气动力研究所 一种基于线程并行的结构化网格流线积分方法
CN113343605B (zh) * 2021-06-29 2022-03-25 华南理工大学 一种汽车涂层抗石击性标准实验的仿真对标方法
CN114186468B (zh) * 2021-12-16 2024-03-29 西安交通大学 一种基于gpu加速的二次电子发射模拟方法及系统及设备
CN114254572B (zh) * 2021-12-16 2024-01-02 西北工业大学太仓长三角研究院 考虑污染物沉积的航发压气机流场性能预测方法及系统
CN114491864B (zh) * 2022-01-26 2022-12-13 哈尔滨工程大学 一种具有参数化、可重构特征的核动力管网模型预处理方法
CN114444362B (zh) * 2022-02-18 2024-09-17 中国石油大学(北京) 一种球形颗粒表面液相信息识别方法、系统和介质
CN115600449B (zh) * 2022-07-05 2024-06-18 浙江大学 基于细长柔性物件组成的大规模颗粒系统的模拟预测方法
CN115719047B (zh) * 2022-11-14 2023-06-06 沐曦集成电路(上海)有限公司 基于波形gpu联合仿真系统
CN116738806A (zh) * 2023-04-18 2023-09-12 西北农林科技大学 含丁坝水流交汇区微塑料输运规律的模拟方法
CN116562065B (zh) * 2023-07-12 2023-09-12 北京凌云智擎软件有限公司 一种网格拓扑的转换方法、设备及装置
CN117272769B (zh) * 2023-09-27 2024-06-07 长江大学 一种颗粒混合过程仿真数据的处理方法
CN117217132B (zh) * 2023-11-08 2024-02-09 四川新能源汽车创新中心有限公司 一种基于cfd-dem耦合的颗粒运动模拟方法
CN118069374B (zh) * 2024-04-18 2024-06-18 清华大学 数据中心智能训练仿真事务加速方法、装置、设备及介质
CN118095144B (zh) * 2024-04-23 2024-06-25 临沂大学 一种基于计算流体力学的高碱煤灰沉积模拟方法及系统

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2009069930A (ja) * 2007-09-11 2009-04-02 Prometech Software Inc 粒子法シミュレーションのためのスライスデータ構造、およびスライスデータ構造を利用した粒子法シミュレーションのｇｐｕへの実装方法
JP2010238030A (ja) * 2009-03-31 2010-10-21 Japan Agengy For Marine-Earth Science & Technology 粒子シミュレーション装置及び粒子シミュレーション方法
US20110032255A1 (en) * 2009-08-10 2011-02-10 Dem Solutions Ltd Method and apparatus for discrete element modeling with a virtual geometry object
WO2012034176A1 (en) * 2010-09-15 2012-03-22 Commonwealth Scientific And Industrial Research Organisation Discrete element method
US20120259602A1 (en) * 2011-04-08 2012-10-11 Japan Agency For Marine-Earth Science And Technology Particle simulator and method of simulating particles

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7952583B2 (en) * 2000-06-19 2011-05-31 Mental Images Gmbh Quasi-monte carlo light transport simulation by efficient ray tracing
JP5124845B2 (ja) 2007-09-14 2013-01-23 本田技研工業株式会社 金属リング周長補正方法
CN101727653B (zh) * 2008-10-31 2012-03-07 中国科学院过程工程研究所 一种基于图形处理器的多组分系统离散模拟计算方法
US9038088B2 (en) * 2011-03-10 2015-05-19 Nec Laboratories America, Inc. Load balancing on hetrogenous processing cluster based on exceeded load imbalance factor threshold determined by total completion time of multiple processing phases
US9760376B1 (en) * 2016-02-01 2017-09-12 Sas Institute Inc. Compilation for node device GPU-based parallel processing

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2009069930A (ja) * 2007-09-11 2009-04-02 Prometech Software Inc 粒子法シミュレーションのためのスライスデータ構造、およびスライスデータ構造を利用した粒子法シミュレーションのｇｐｕへの実装方法
JP2010238030A (ja) * 2009-03-31 2010-10-21 Japan Agengy For Marine-Earth Science & Technology 粒子シミュレーション装置及び粒子シミュレーション方法
US20110032255A1 (en) * 2009-08-10 2011-02-10 Dem Solutions Ltd Method and apparatus for discrete element modeling with a virtual geometry object
WO2012034176A1 (en) * 2010-09-15 2012-03-22 Commonwealth Scientific And Industrial Research Organisation Discrete element method
US20120259602A1 (en) * 2011-04-08 2012-10-11 Japan Agency For Marine-Earth Science And Technology Particle simulator and method of simulating particles

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
JPN6016005219; 都築 怜理 外３名: '複数ノードのＧＰＵによる大規模パッシブスカラー粒子計算の強スケーリングと動的負荷分散' 情報処理学会研究報告 ２０１１（平成２３）年度▲６▼ ［ＤＶＤ-ＲＯＭ］ , 20120415, ｐｐ．１-１０, 一般社団法人情報処理学会 *
JPN6016005220; 白幡 晃一 外２名: 'ＧＰＵを考慮したＭａｐＲｅｄｕｃｅのタスクスケジューリング' 情報処理学会研究報告 平成２２年度▲３▼ ［ＣＤ-ＲＯＭ］ , 20101015, ｐｐ．１-８, 一般社団法人情報処理学会 *
JPN6016005221; 原田 隆宏 外３名: '複数のＧＰＵを用いた粒子法シミュレーションの並列化' 情報処理学会論文誌 論文誌ジャーナル Ｖｏｌ．４９ Ｎｏ．１２ 【修正版】 ［ＣＤ-ＲＯＭ］ 第４９巻 第１２号, 20081215, ｐｐ．４０６７-４０７９, 社団法人情報処理学会 *

Cited By (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2017091262A (ja) * 2015-11-11 2017-05-25 富士通株式会社 粒子シミュレーションプログラム、粒子シミュレーション装置、及び計算機資源配分方法
KR20190048790A (ko) * 2017-10-31 2019-05-09 한국과학기술원 공간 데이터를 분산 처리하는 시스템 및 방법
KR101980699B1 (ko) * 2017-10-31 2019-05-22 한국과학기술원 공간 데이터를 분산 처리하는 시스템 및 방법
KR20200023954A (ko) * 2018-08-27 2020-03-06 이에이트 주식회사 복수의 프로세서를 이용한 입자 기반의 유체 시뮬레이션 방법 및 유체 시뮬레이션 장치
KR102181989B1 (ko) * 2018-08-27 2020-11-23 이에이트 주식회사 복수의 프로세서를 이용한 입자 기반의 유체 시뮬레이션 방법 및 유체 시뮬레이션 장치
JP2021009568A (ja) * 2019-07-01 2021-01-28 富士通株式会社 情報処理装置、粒子シミュレーション方法および粒子シミュレーションシステム
JP7244757B2 (ja) 2019-07-01 2023-03-23 富士通株式会社 情報処理装置、粒子シミュレーション方法および粒子シミュレーションシステム
US11892387B2 (en) 2019-07-01 2024-02-06 Fujitsu Limited Information processing apparatus and method and system for particle simulation
KR20210136327A (ko) * 2020-05-07 2021-11-17 국방과학연구소 도심지 유동해석 방법 및 장치
KR102371345B1 (ko) * 2020-05-07 2022-03-07 국방과학연구소 도심지 유동해석 방법 및 장치

Also Published As

Publication number Publication date
JP6009075B2 (ja) 2016-10-19
US10007742B2 (en) 2018-06-26
CN103324780A (zh) 2013-09-25
WO2014094410A1 (zh) 2014-06-26
GB201500658D0 (en) 2015-03-04
CN103324780B (zh) 2016-03-16
GB2523640A (en) 2015-09-02
GB2523640B (en) 2020-05-27
US20150213163A1 (en) 2015-07-30

Similar Documents

Publication Publication Date Title
JP6009075B2 (ja) 粒子流動のシミュレーションシステム及びその方法
CN110706341B (zh) 一种城市信息模型的高性能渲染方法、装置及存储介质
US10262392B2 (en) Distributed and parallelized visualization framework
JP5227103B2 (ja) 近傍粒子探索に用いるデータ構造の構築方法、そのプログラム、およびそのプログラムを格納した記憶媒体
Valero-Lara et al. Accelerating fluid–solid simulations (Lattice-Boltzmann & Immersed-Boundary) on heterogeneous architectures
Hernandez et al. Simulating and visualizing real-time crowds on GPU clusters
Bruckschen et al. Real-time out-of-core visualization of particle traces
CN103345580A (zh) 基于格子Boltzmann方法的并行CFD方法
CN107016180A (zh) 一种颗粒流动仿真方法
Guan et al. pRPL 2.0: Improving the parallel raster processing library
CN108701275A (zh) 粒子模拟装置、粒子模拟方法以及粒子模拟程序
Akkurt et al. An efficient edge based data structure for the compressible Reynolds‐averaged Navier–Stokes equations on hybrid unstructured meshes
Kuester et al. Visualization of particle traces in virtual environments
CN100492371C (zh) 分布型cad装置
Hung et al. Automatic clustering method for real-time construction simulation
Chen et al. Flexible and rapid animation of brittle fracture using the smoothed particle hydrodynamics formulation
Plimpton et al. Rendezvous algorithms for large-scale modeling and simulation
CN113628318B (zh) 一种基于光线追踪的分布式实时神经元渲染方法和系统
CN112149287B (zh) 一种面向负载均衡的交通仿真路网图形化分割方法及系统
Yong et al. GVM based intuitive simulation web application for collision detection
Pacevič et al. Visualization of cracks by using the local Voronoi decompositions and distributed software
Fehling Algorithms for massively parallel generic hp-adaptive finite element methods
Oh et al. Practical simulation of hierarchical brittle fracture
Niebling et al. Integrated simulation workflows in computer aided engineering on hpc resources
Fofanov et al. Optimization of load balancing algorithms in parallel modeling of objects using a large number of grids

Legal Events

Date Code Title Description
A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20160216

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20160516

TRDD Decision of grant or rejection written
A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20160830

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20160913

R150 Certificate of patent or registration of utility model

Ref document number: 6009075

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250