JP6261130B2 - 粒子シミュレーション装置、粒子シミュレーション方法及び粒子シミュレーションプログラム - Google Patents

粒子シミュレーション装置、粒子シミュレーション方法及び粒子シミュレーションプログラム Download PDF

Info

Publication number
JP6261130B2
JP6261130B2 JP2014115567A JP2014115567A JP6261130B2 JP 6261130 B2 JP6261130 B2 JP 6261130B2 JP 2014115567 A JP2014115567 A JP 2014115567A JP 2014115567 A JP2014115567 A JP 2014115567A JP 6261130 B2 JP6261130 B2 JP 6261130B2
Authority
JP
Japan
Prior art keywords
particle
pair
particles
setting
matrix
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.)
Expired - Fee Related
Application number
JP2014115567A
Other languages
English (en)
Other versions
JP2015230535A (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.)
Japan Agency for Marine Earth Science and Technology
Original Assignee
Japan Agency for Marine Earth Science and Technology
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 Japan Agency for Marine Earth Science and Technology filed Critical Japan Agency for Marine Earth Science and Technology
Priority to JP2014115567A priority Critical patent/JP6261130B2/ja
Priority to PCT/JP2015/065609 priority patent/WO2015186633A1/ja
Priority to US15/315,254 priority patent/US10354099B2/en
Priority to EP15803426.4A priority patent/EP3153982A4/en
Publication of JP2015230535A publication Critical patent/JP2015230535A/ja
Application granted granted Critical
Publication of JP6261130B2 publication Critical patent/JP6261130B2/ja
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B5/00—ICT specially adapted for modelling or simulations in systems biology, e.g. gene-regulatory networks, protein interaction networks or metabolic networks
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06G—ANALOGUE COMPUTERS
    • G06G7/00—Devices in which the computing operation is performed by varying electric or magnetic quantities
    • G06G7/48—Analogue computers for specific processes, systems or devices, e.g. simulators
    • G06G7/485—Analogue computers for specific processes, systems or devices, e.g. simulators for determining the trajectory of particles, e.g. of electrons
    • 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
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
    • G16C10/00—Computational theoretical chemistry, i.e. ICT specially adapted for theoretical aspects of quantum chemistry, molecular mechanics, molecular dynamics or the like

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • General Physics & Mathematics (AREA)
  • General Health & Medical Sciences (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Computing Systems (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Mathematical Physics (AREA)
  • Computer Hardware Design (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Molecular Biology (AREA)
  • Biophysics (AREA)
  • Medical Informatics (AREA)
  • Evolutionary Biology (AREA)
  • Biotechnology (AREA)
  • Physiology (AREA)
  • Mathematical Optimization (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Algebra (AREA)
  • Geometry (AREA)
  • Evolutionary Computation (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Chemical & Material Sciences (AREA)
  • High Energy & Nuclear Physics (AREA)
  • Dispersion Chemistry (AREA)
  • Analytical Chemistry (AREA)
  • Biochemistry (AREA)
  • Immunology (AREA)
  • Pathology (AREA)

Description

本発明は、粒子シミュレーション装置、粒子シミュレーション方法及び粒子シミュレーションプログラムに関する。
粒子シミュレーションは、基本的シミュレーション技術の一つであり、流体あるいは固体といった基礎的物質のほか、土砂あるいは紛体のような工学的重要な物質、またまたタンパク質のような生物及び医療において重要な物質のシミュレーションに用いられる。例えば、土砂崩れ、液状化、津波等の災害シミュレーションや、タンパク質設計等の創薬シミュレーションに用いられる。これら様々な物質シミュレーションは、きわめて応用上価値の高いものであり、シミュレーション技術の進歩は、これらの応用を迅速に発展させるために重要なものである。
粒子シミュレーションは、粒子一個一個の位置と速度とを変数として保持し、その変化をモデルに基づいて追跡することにより粒子群全体の運動を実現する。特に粒子シミュレーションとして広く使用されているＤＥＭ（Discrete Element Method）（土、砂）、ＳＰＨ（SmoothedParticle Hydrodynamics）（流体）、ＭＤ（Molecular Dynamics）（分子、タンパク質）といった手法では、粒子間の相互作用を想定し、特にその相互作用が対称であることを仮定する（作用反作用の法則）。以下の非特許文献１では、この対称性を用いて粒子シミュレーションを高速化するという手法が提案されている。この手法は、並列計算が可能なようにＧＰＵ（Graphics Processing Unit）計算向けに作成されたものである。この手法において、対称性の利用はペアリストと呼ばれる粒子間での相互作用ペアのリストが作成されて使用されるという特徴を持つ。このペアリストに基づいた相互作用の計算によって全体のシミュレーション時間が短縮される。
D. Nishiura, H. Sakaguchi,Parallel-vector algorithms for particle simulations on shared-memorymultiprocessors, J. Comp. Phys. 230 (2011) 1923−1938.
非特許文献１に記載された方法では、各粒子に対して粒子の位置に基づいた粒子番号が付与され、粒子毎に別の粒子とのペア（接触候補となる粒子とのペア）が構成される。このペアに対しては、ペアを特定する番号が付与され、ペア毎に接触力が計算される。上記の方法では、粒子毎の接触力の総和を演算するため、粒子毎に自身が構成要素となっているペアの番号のリストが作成される。
このリストは、自身よりも粒子番号が大きな粒子とのペアについてのリストと、自身よりも粒子番号が小さな粒子とのペアについてのリストとの二つのリストが生成される。このうち、自身よりも粒子番号が小さな粒子とのペアについてのリストを生成するため、非特許文献１に記載された方法では、ｉｆ文（条件の判断）を含む、二重のｆｏｒループ（ダブルｆｏｒループ）によって探索を行うアルゴリズムが用いられる。このアルゴリズムは、多くのメモリアクセスが必要な多数のトライアルアンドエラー処理を行うものである。
また、このような多数のトライアルアンドエラー処理は、ＧＰＵの計算性能を大幅に低下させるワープダイバージェンスを引き起こすため、ＧＰＵで実行する場合には問題が大きいものとなる。
本発明は、上記の問題点に鑑みてなされたものであり、演算効率を向上させることができる粒子シミュレーション装置、粒子シミュレーション方法及び粒子シミュレーションプログラムを提供することを目的とする。
上記の目的を達成するために、本発明に係る粒子シミュレーション装置は、作業空間内の複数の粒子について他の粒子との相互作用力に基づき位置及び速度を算出し、粒子の挙動をシミュレーションする粒子シミュレーション装置であって、複数の粒子それぞれについて、粒子の位置を示す位置情報を取得する位置情報取得手段と、複数の粒子それぞれに対してソート可能な粒子番号を設定する粒子番号設定手段と、位置情報取得手段によって取得された位置情報に基づいて、互いに近傍に位置する粒子のペアを選択すると共に、粒子番号設定手段によって設定された、当該ペアを構成する粒子の一方の粒子番号に基づき当該ペアに対してペア番号を設定するペア設定手段と、ペア設定手段によって設定されたペア番号、及び当該ペアを構成する粒子の粒子番号を行の構成要素とする行列を生成し、当該ペアを構成する粒子のもう一方の粒子番号に基づいて当該行列の列の順序をソートし、ソートした行列に基づき、粒子の粒子番号から当該粒子が構成するペアのペア番号を参照するための参照用情報を生成する参照用情報生成手段と、ペア設定手段によって選択されたペアそれぞれに係る粒子同士の相互作用判定を行い、相互作用していると判定された粒子同士の相互作用力を計算する相互作用力演算手段と、参照用情報生成手段によって生成された参照用情報に基づいて、相互作用力演算手段によって計算された相互作用力から粒子毎の相互作用力の総和を計算する総和演算手段と、総和演算手段によって計算された粒子毎の相互作用力の総和に基づいて、粒子の位置及び速度を算出する粒子情報算出手段と、を備える。
本発明に係る粒子シミュレーション装置では、ソートした行列を用いることで、二重のｆｏｒループを用いずに上記のリストに相当する参照用情報が生成される。これにより、多数のトライアルアンドエラー処理が不要となる。即ち、本発明に係る粒子シミュレーション装置によれば演算効率を向上させることできる。
参照用情報生成手段は、行列の一つの列の構成要素を、ペアを構成する粒子の粒子番号のうち小さい方の粒子番号とし、別の列の構成要素を、ペアを構成する粒子の粒子番号のうち大きい方の粒子番号とすることとしてもよい。この構成によれば、適切かつ確実にソートした行列に基づき、参照用情報を生成することができる。これにより、適切かつ確実に本発明を実施することができる。
粒子番号設定手段は、位置情報取得手段によって取得された位置情報に基づいて、複数の粒子それぞれに対してソート可能な粒子番号を設定することとしてもよい。この構成によれば、適切かつ確実に粒子に粒子番号を付与することができる。これにより、適切かつ確実に本発明を実施することができる。
作業空間は、複数のセルに分割されており、位置情報取得手段は、位置情報として、複数の粒子それぞれについて、粒子が位置するセルを示す情報を取得する、こととしてもよい。この構成によれば、効率的かつ容易に近傍に位置する粒子のペアを選択することができ、演算効率を更に向上させることできる。
ところで、本発明は、上記のように粒子シミュレーション装置の発明として記述できる他に、以下のように粒子シミュレーション方法及び粒子シミュレーションプログラムの発明としても記述することができる。これはカテゴリが異なるだけで、実質的に同一の発明であり、同様の作用及び効果を奏する。
即ち、本発明に係る粒子シミュレーション方法は、作業空間内の複数の粒子について他の粒子との相互作用力に基づき位置及び速度を算出し、粒子の挙動をシミュレーションする粒子シミュレーション装置の動作方法である粒子シミュレーション方法であって、複数の粒子それぞれについて、粒子の位置を示す位置情報を取得する位置情報取得ステップと、複数の粒子それぞれに対してソート可能な粒子番号を設定する粒子番号設定ステップと、位置情報取得ステップにおいて取得された位置情報に基づいて、互いに近傍に位置する粒子のペアを選択すると共に、粒子番号設定ステップにおいて設定された、当該ペアを構成する粒子の一方の粒子番号に基づき当該ペアに対してペア番号を設定するペア設定ステップと、ペア設定ステップにおいて設定されたペア番号、及び当該ペアを構成する粒子の粒子番号を行の構成要素とする行列を生成し、当該ペアを構成する粒子のもう一方の粒子番号に基づいて当該行列の列の順序をソートし、ソートした行列に基づき、粒子の粒子番号から当該粒子が構成するペアのペア番号を参照するための参照用情報を生成する参照用情報生成ステップと、ペア設定ステップにおいて選択されたペアそれぞれに係る粒子同士の相互作用判定を行い、相互作用していると判定された粒子同士の相互作用力を計算する相互作用力演算ステップと、参照用情報生成ステップにおいて生成された参照用情報に基づいて、相互作用力演算ステップにおいて計算された相互作用力から粒子毎の相互作用力の総和を計算する総和演算ステップと、総和演算ステップにおいて計算された粒子毎の相互作用力の総和に基づいて、粒子の位置及び速度を算出する粒子情報算出ステップと、を含む。
また、本発明に係る粒子シミュレーションプログラムは、コンピュータを、作業空間内の複数の粒子について他の粒子との相互作用力に基づき位置及び速度を算出し、粒子の挙動をシミュレーションする粒子シミュレーション装置として機能させる粒子シミュレーションプログラムであって、コンピュータを、複数の粒子それぞれについて、粒子の位置を示す位置情報を取得する位置情報取得手段と、複数の粒子それぞれに対してソート可能な粒子番号を設定する粒子番号設定手段と、位置情報取得手段によって取得された位置情報に基づいて、互いに近傍に位置する粒子のペアを選択すると共に、粒子番号設定手段によって設定された、当該ペアを構成する粒子の一方の粒子番号に基づき当該ペアに対してペア番号を設定するペア設定手段と、ペア設定手段によって設定されたペア番号、及び当該ペアを構成する粒子の粒子番号を行の構成要素とする行列を生成し、当該ペアを構成する粒子のもう一方の粒子番号に基づいて当該行列の列の順序をソートし、ソートした行列に基づき、粒子の粒子番号から当該粒子が構成するペアのペア番号を参照するための参照用情報を生成する参照用情報生成手段と、ペア設定手段によって選択されたペアそれぞれに係る粒子同士の相互作用判定を行い、相互作用していると判定された粒子同士の相互作用力を計算する相互作用力演算手段と、参照用情報生成手段によって生成された参照用情報に基づいて、相互作用力演算手段によって計算された相互作用力から粒子毎の相互作用力の総和を計算する総和演算手段と、総和演算手段によって計算された粒子毎の相互作用力の総和に基づいて、粒子の位置及び速度を算出する粒子情報算出手段と、として機能させる。
本発明では、ソートした行列を用いることで、二重のｆｏｒループを用いずに上記のリストに相当する参照用情報が生成されるため、多数のトライアルアンドエラー処理が不要となる。即ち、本発明によれば演算効率を向上させることできる。
本発明の実施形態に係る粒子シミュレーション装置の機能構成を示す図である。 本実施形態における作業領域、当該作業領域に含まれる粒子、及び当該粒子によって構成されるペアを示す図である。 本実施形態において生成等される行列、配列を示す図である。 本発明の実施形態に係る粒子シミュレーション装置で実行される処理全体（粒子シミュレーション方法）を示すフローチャートである。 参照用情報であるＲ（２）を生成する処理を示すフローチャートである。 非特許文献１に記載された方法において生成等される行列、配列を示す図である。 非特許文献１に記載された方法の（ａ）性能試験に用いたモデル（ｂ）性能試験の結果を示す図である。 本実施形態、及び非特許文献１に記載された方法における、性能試験での、粒子密度とペアインデックスマトリクスを生成するためにかかった経過時間との関係を示すグラフである。 粒子シミュレーションにおける粒子密度の例を示す図である。 本実施形態、及び非特許文献１に記載された方法における、性能試験での、ステップ毎の経過時間及びメモリアクセス数を示すグラフである。 性能試験に用いたモデルのシミュレーションの実行中の状態（スナップショット）を示す図である。 本実施形態、及び非特許文献１に記載された方法における、性能試験での、時刻（ステップ）と有効粒子密度及びペアインデックスマトリクスを生成するためにかかった経過時間との関係を示すグラフである。 本発明の実施形態に係る粒子シミュレーションプログラムの構成を、記録媒体と共に示す図である。
以下、図面と共に本発明に係る粒子シミュレーション装置、粒子シミュレーション方法及び粒子シミュレーションプログラムの実施形態について詳細に説明する。なお、図面の説明においては同一要素には同一符号を付し、重複する説明を省略する。
図１に本実施形態に係る粒子シミュレーション装置１０を示す。粒子シミュレーション装置１０は、作業空間内の複数の球形の粒子の挙動をシミュレーション（解析）する装置である。具体的には、粒子シミュレーション装置１０は、シミュレーション上の時刻（ステップ）毎の各粒子の位置及び速度に基づき各粒子に働く力を算出する。各粒子に働く力には、粒子間での相互作用である接触（衝突）による相互作用力である接触力が含まれる。粒子シミュレーション装置１０は、算出した力に基づき次の時刻における各粒子の位置及び速度を算出する。
本実施形態に係る粒子シミュレーション装置１０によるシミュレーションの対象となる粒子は、従来の粒子シミュレーションの対象となっていた任意の粒子を含む。例えば、上述したような土砂や粉体を対象とすることができる。あるいは、流体や個体を複数の粒子からなるものと仮定して対象とすることとしてもよい。本実施形態に係る粒子シミュレーション装置１０によるシミュレーションにより、物理的な問題をシミュレーションすることができる。
なお、本実施形態に係る粒子シミュレーション装置１０の機能は、上述した非特許文献１及び特開２０１０−２３８０３０号公報（特許文献１）に記載されたシミュレーションの機能を改良したものである。従って、本実施形態に係る粒子シミュレーション装置１０の本発明に係る機能以外の部分は、非特許文献１及び特許文献１に記載の内容で実現されていてもよい。なお、特段の記載をしていない部分については、粒子シミュレーション装置１０の機能は、特許文献１に記載されているものと同様である。
粒子シミュレーション装置１０は、例えば、ＣＰＵ（CentralProcessing Unit）、ＧＰＵ（Graphics Processing Unit）、メモリ、ハードディスク、ディスプレイ等のハードウェアを備えるコンピュータとして構成される。これらの構成要素がプログラム等により動作することによって、後述する粒子シミュレーション装置１０としての機能が発揮される。粒子シミュレーション装置１０は、並列演算が可能な装置において特に効果的に動作する。なお、粒子シミュレーション装置１０は、演算装置としては必ずしもＧＰＵを備えている必要はなく、ＣＰＵのみを備えた構成（スカラー機）であってもよい。
図１に示すように粒子シミュレーション装置１０は、機能的な構成要素として、粒子情報保持部１１と、位置情報取得部１２と、粒子番号設定部１３と、ペア設定部１４と、参照用情報生成部１５と、接触力演算部１６と、総和演算部１７と、粒子情報算出部１８とを備えて構成される。
本実施形態において粒子が運動する領域である作業領域は、三次元の空間であり、図２に示すように一辺の大きさが予め設定された立方体のセルに分割（区分）されている。粒子シミュレーション装置１０は、シミュレーションの処理を行う前に作業空間を予めセルに分割しており、作業空間がどのようにセルに分割されているか予め把握している。上記の一辺の大きさは、例えば、Ｄｍａｘ×（１．０＋α）とする。Ｄｍａｘは、シミュレーションの対象となる複数の粒子の粒子径のうち、最大の値である。αは、後述する接触候補リストの更新頻度を調整するパラメータであり、例えば、α＝０．２である。また、一辺の大きさは、後述するカットオフ長より大きい値としてもよい。また、作業領域内の各セルにはセルを特定するセル番号が付されている。セル番号は、例えば、作業空間内のセルの位置に応じて順番に付されている。
粒子情報保持部１１は、作業領域内の複数の粒子それぞれについての粒子情報を保持する手段である。粒子情報は、粒子の座標、粒子の速度及び粒子半径を示す情報を含む。粒子の座標は、作業空間における粒子の位置を示す三次元座標である。粒子の速度は、並進速度及び回転速度を含む。粒子の座標及び粒子の速度については、シミュレーションの開始時の情報（初期情報）は、予め粒子情報保持部１１に粒子シミュレーション装置１０のユーザ等により入力されており、また、シミュレーション中の情報は、後述する粒子情報算出部１８によって更新される。粒子半径は、予め粒子情報保持部１１に粒子シミュレーション装置１０のユーザ等により入力されている。また、粒子情報保持部１１は、粒子情報以外のシミュレーションに利用される情報を、予め入力して保持していてもよい。このような情報としては、摩擦係数、弾性係数、粘性減衰係数、反発係数等である。
位置情報取得部１２は、複数の粒子それぞれについて、粒子の位置を示す位置情報を取得する位置情報取得手段である。位置情報取得部１２は、位置情報として、複数の粒子それぞれについて、粒子が位置するセルを示す情報であるセル番号を取得する。具体的には、位置情報取得部１２は、粒子情報保持部１１に保持された各粒子の現時刻（現ステップ）の粒子情報を取得する。位置情報取得部１２は、取得した粒子情報によって示される粒子の座標がどのセルに含まれるか否かを判断し、当該座標を含むセルのセル番号を当該粒子の位置情報とする。位置情報取得部１２は、取得した粒子の位置情報を粒子番号設定部１３及びペア設定部１４に出力する。また、位置情報取得部１２は、取得した粒子情報を接触力演算部１６に入力する。
なお、位置情報取得部１２による位置情報の取得及び位置情報の取得に基づく処理は、各時刻で行われる必要はなく、後述するペア設定部１４による粒子のペアの設定（更新）が必要な場合に行われることとしてもよい。そこで、位置情報取得部１２は、ペアの設定（更新）が行われるべきか否かを判定することとしてもよい。この判定は、例えば、以下のように行う。各粒子について現時刻までの各時刻における粒子の座標に基づき積算移動距離を算出して、当該粒子の半径をｒとしたときに当該積算移動距離がｒαよりも大きいか否かを判断する。何れかの粒子において、積算移動距離がｒαよりも大きいと判断された場合、ペアの設定（更新）が行われるべきと判定する。なお、積算移動距離は、ペアの設定（更新）を行った場合には初期化される。また、最初の時刻では、上記の判定は行わずにペアの設定を行うものとする。
位置情報取得部１２は、ペアの設定（更新）が必要であると判定した場合には、取得した粒子の位置情報を粒子番号設定部１３及びペア設定部１４に出力し、ペアの設定（更新）の処理を行われる。位置情報取得部１２は、ペアの設定（更新）が必要ではないと判定した場合には、接触力演算部１６に対してその旨を通知し、当該時刻ではペアの設定（更新）は行わず、接触力演算部１６による処理を行うこととしてもよい。
粒子番号設定部１３は、複数の粒子それぞれに対してソート可能な粒子番号（粒子インデックス）を設定する粒子番号設定手段である。粒子番号設定部１３は、位置情報取得部１２によって取得された位置情報に基づいて、複数の粒子それぞれに対してソート可能な粒子番号を設定する。粒子番号設定部１３は、粒子を位置情報（粒子が所属するセルのセル番号順）に並べ、その順に粒子番号を設定する（付け直す）。粒子番号としては、例えば、１から昇順の整数が設定される。この粒子番号の設定によって、例えば、図２（ａ）に示すように各粒子に粒子の位置に応じて順序付けられた粒子番号（図２（ａ）の例では、１〜１０）が設定される。なお、粒子情報の設定の際、特許文献１に示されているように粒子番号から、粒子情報保持部１１に保持された粒子情報を取得できるようにされる。粒子番号設定部１３は、設定した粒子番号をペア設定部１４に出力する。
ペア設定部１４は、位置情報取得部１２によって取得された位置情報に基づいて、互いに近傍に位置する粒子のペアを選択すると共に、粒子番号設定部１３によって設定された、当該ペアを構成する粒子の一方の粒子番号に基づき当該ペアに対してペア番号（ペアインデックス）を設定するペア設定手段である。互いに近傍に位置する粒子とは、接触する可能性がある２つの粒子である。
本実施形態では、図３（ａ）に示す相互作用行列Ｕの上三角部分を構成する（ことを想定する）。ここで、ｉ，ｊはそれぞれ粒子番号である。相互作用行列Ｕは、粒子ｉと粒子ｊとが近傍に位置する粒子のペアを構成するものである場合、Ｕｉ，ｊ＝１であり、粒子ｉと粒子ｊとが近傍に位置する粒子のペアを構成するものでない場合、Ｕｉ，ｊ＝０である行列である。
ペア設定部１４は、位置情報取得部１２及び粒子番号設定部１３から入力された情報に基づき、セル毎にセルに属する粒子の粒子番号の最大値と最小値とを粒子シミュレーション装置１０のメモリ上に記憶する。ペア設定部１４は、同一のセル及び隣接するセル（例えば、本実施形態のように三次元のセルであれば合計で２７個のセル）に属する粒子同士を近傍に位置する粒子のペアとする。
セルの隣接関係は、予めペア設定部１４に記憶されている。ペア設定部１４は、上記のように記憶されたセルに属する粒子の粒子番号の最大値と最小値とを参照して、各粒子に対して、当該粒子の粒子番号よりも大きい粒子番号が設定された粒子からペアとなる粒子を決定する。ペア設定部１４は、各粒子に対して、当該粒子が所属するセル及び当該セルに隣接するセルに属している粒子のうち、自身よりも大きな粒子番号が設定された粒子をペアとなる粒子として決定する。なお、粒子番号は粒子が所属するセルのセル番号順に付されているので、自身の粒子番号より小さい粒子番号の粒子しか含まれないセルは参照する必要がない。従って、粒子ｉに対して、ｉ＞ｊとなる粒子ｊが含まれる１４個の隣接するセルのみを探索すればよい。
このとき、単に隣接するセルに属する粒子同士をペアにするのではなく、粒子中心間距離に基づいて粒子同士がペアになるか否かを判断してもよい。例えば、粒子中心間距離が（ｒｉ＋ｒｊ）（１．０＋α）以下であるか否かを判断し、当該条件を満たす場合のみにペアとすることとしてもよい。ここで、ｒｉ，ｒｊは、それぞれ粒子ｉ，ｊの粒子半径である。
図３（ａ）に示す相互作用行列Ｕは、メモリ上のスペースが多く必要となるため、ペア設定部１４は、各粒子の粒子番号ｉ毎に自身の粒子番号よりも大きい粒子番号を持つと共にペアとなる粒子の粒子番号を格納した図３（ｂ）に示す配列Ｕ（１）を生成する。また、ペア設定部１４は、図３（ｂ）に示すように、Ｕ（１）から各粒子の粒子番号ｉ毎に自身の粒子番号よりも大きい粒子番号を持つと共にペアとなる粒子の数（配列Ｕ（１）のｉ毎の要素の数）ｎｉ＜ｊを算出する。
ペア設定部１４は、上記のように設定された粒子のペアに対してペア番号を設定する。ペア設定部１４は、ペア番号を、ペアとなる２つの粒子番号のうち小さい粒子番号が小さい順に設定する。ペア番号としては、例えば、１から昇順の整数が設定される。ペアとなる２つの粒子番号のうち小さい粒子番号が同一である場合、大きい方の粒子番号が小さい順に設定する。このように、ペア番号は、当該ペアを構成する粒子の一方の粒子番号（２つの粒子の粒子番号のうち小さい方の粒子番号）に基づき設定される。ペア設定部１４は、具体的には、特許文献１に示されるように、粒子番号ｉ毎の自身の粒子番号よりも大きい粒子番号を持つと共にペアとなる粒子の数のプレフィックス和ｓｉ＜ｊに基づき粒子番号を設定する。このように設定されたペア番号は、図２（ｂ）に示されるようになる。ペア設定部１４は、粒子番号ｉ毎の自身の粒子番号ｉよりも大きい粒子番号ｊを持つ粒子とのペアのペア番号を格納した図３（ｂ）に示す行列Ｒ（１）を生成する。ペア設定部１４は、生成したＵ（１）、ｓｉ＜ｊ、行列Ｒ（１）を参照用情報生成部１５に出力する。
また、ペア設定部１４は、特許文献１に示されているように現在のペアが１ステップ前にもペアであったかを調べる。ペア設定部１４は、特許文献１に示されているように現在のペアが１ステップ前にもペアであった場合には、１ステップ前の当該ペアの接触力を現在のペアにも引き継ぐ処理を行う。
参照用情報生成部１５は、ペア設定部１４によって設定されたペア番号、及び当該ペアを構成する粒子の粒子番号を行の構成要素とする行列を生成し、当該ペアを構成する粒子のもう一方の粒子番号に基づいて当該行列の列の順序をソートし、ソートした行列に基づき、粒子の粒子番号ｉから当該粒子が構成するペアのペア番号を参照するための参照用情報を生成する参照用情報生成手段である。参照用情報生成部１５は、行列の一つの列の構成要素を、ペアを構成する粒子の粒子番号のうち小さい方の粒子番号ｉとし、別の列の構成要素を、ペアを構成する粒子の粒子番号のうち大きい方の粒子番号ｊとする。
参照用情報（ペアリスト、ペアインデックスマトリクス）は、各ペアについて粒子間に発生する接触力が算出された後、粒子毎の接触力の総和を算出する際に参照されるものである。自身の粒子番号ｉよりも大きい粒子番号（ｉ＜ｊ）を持つ粒子ｊとのペアのペア番号を参照するための参照用情報は、ペア設定部１４によって生成された行列Ｒ（１）である。参照用情報生成部１５は、自身の粒子番号ｉよりも小さい粒子番号（ｉ＞ｊ）を持つ粒子ｊとのペアのペア番号を参照するための参照用情報である行列Ｒ（２）を生成する。Ｒ（２）は、粒子番号ｉ毎の自身の粒子番号よりも小さい粒子番号を持つ粒子とのペアのペア番号を格納した行列である。
参照用情報生成部１５は、ペア設定部１４から入力されたＵ（１）及びＲ（１）から、図３（ｃ）に示す３つの列を有する行列Ｔ（１）（トリプレット）を生成する。第１列はペア番号ｐである。第２列及び第３列は、第１列のペア番号に示されるペアを構成する粒子の粒子番号である。第２列は当該粒子番号のうち小さい方の粒子番号ｉであり、第３列は大きい方の粒子番号ｊである。参照用情報生成部１５は、要素をペア番号ｐ及び粒子番号ｉ（一方の粒子番号）が昇順となるように行を並べて、行列Ｔ（１）を生成する。
続いて、参照用情報生成部１５は、第３列の粒子番号ｊ（もう一方の粒子番号）が昇順となるように行列Ｔ（１）の行をソートして（行を並べ替えて）、図３（ｃ）に示す行列Ｔ（２）を生成する。参照用情報生成部１５は、行列Ｔ（２）の第３列の粒子番号ｊ毎の要素数をカウントして、図３（ｄ）に示すｎｉ＞ｊを算出する（ここで、ｎｉ＞ｊのインデックス「ｉ＞ｊ」におけるｉが、行列Ｔ（２）の第３列の粒子番号ｊに相当する）。例えば、行列Ｔ（２）の第３列の粒子番号ｊにおいてｊ＝３の行は１つであるのでｉ＝３のときのｎｉ＞ｊは１となり、ｊ＝９の行は３つであるのでｉ＝９のときのｎｉ＞ｊは３となる。ｎｉ＞ｊは、各粒子の粒子番号ｉ毎に自身の粒子番号よりも小さい粒子番号を持つと共にペアとなる粒子の数である。なお、ｎｉ＞ｊは、非特許文献１の手法でも算出されているが、非特許文献１の手法では２７の隣接セルの検索によって得られるものである。
続いて、参照用情報生成部１５は、ｎｉ＞ｊからプレフィックス和ｓｉ＞ｊ（ｉ）＝Σｉｎｉ＞ｊ（ｋ）を算出する。ここで、ｋは粒子番号を示す（総和計算で変動される）変数である。参照用情報生成部１５は、行列Ｔ（２）、ｎｉ＞ｊ及びｓｉ＞ｊ（ｉ）から、図３（ｄ）に示す行列Ｒ（２）を生成する。具体的には、参照用情報生成部１５は、Ｔ（２）の第１列のうち、ｓｉ＞ｊ（ｉ−１）＋１，…，ｓｉ＞ｊ（ｉ−１）＋ｎｉ＞ｊ番目の行のペア番号ｐを行列Ｒ（２）のｉ番目（粒子番号ｉ）の行の要素とする。図３（ｄ）に示すように、行列Ｒ（２）のｉ番目（粒子番号ｉ）の要素数は、ｎｉ＞ｊである。
参照用情報生成部１５による行列Ｒ（２）の生成は、以下に示すアルゴリズムで実現することができる。
参照用情報生成部１５は、生成した参照用情報であるＲ（１）及びＲ（２）を総和演算部１７に出力する。また、参照用情報生成部１５は、ペア番号ｐと当該ペア番号ｐによって示されるペアを構成する粒子の粒子番号ｉ，ｊとの対応が特定可能な情報（例えば、行列Ｔ（１））を接触力演算部１６に出力する。
接触力演算部１６は、ペア設定部１４によって選択されたペアそれぞれに係る粒子同士の接触判定（相互作用判定）を行い、接触（相互作用）していると判定された粒子同士の接触力を計算する相互作用力演算手段である。接触力演算部１６は、参照用情報生成部１５から入力される情報によって、ペア設定部１４によって選択されたペアのペア番号ｐと当該ペア番号ｐによって示されるペアを構成する粒子の粒子番号ｉ，ｊとの対応を特定する。
接触力演算部１６は、位置情報取得部１２から、粒子情報保持部１１に保持されているペアを構成する粒子の粒子情報を入力する。また、接触力演算部１６は、粒子情報保持部１１に保持されている、接触力を算出に必要なパラメータを取得し、接触力の計算に用いてもよい。接触力演算部１６は、粒子情報に基づいて粒子間距離（あるいは、粒子中心間距離）を算出して、当該粒子間距離が予め接触力演算部１６に記憶された閾値（カットオフ長）よりも小さいか否かを判断することで粒子同士の接触判定を行う。粒子間距離が閾値よりも小さい場合には、粒子同士が接触していると判定する。
接触力演算部１６は、粒子同士が接触していると判定した場合には、当該ペアについて接触力の計算を行う。この接触力の計算では、接触力の並進成分と回転成分とがそれぞれ計算される。接触力の計算は、特許文献１に記載されているようにＶｏｉｇｔモデルに基づいて行われる。接触力演算部１６は、粒子同士が接触していないと判定した場合には、当該ペアについて接触力をゼロとする。
粒子間に働く接触力は、粒子間で対称的なものである。即ち、粒子番号ｉの粒子が粒子番号ｊの粒子から受ける接触力をＦｉ，ｊとし、粒子番号ｊの粒子が粒子番号ｉの粒子から受ける接触力をＦｊ，ｉとすると、Ｆｉ，ｊ＝−Ｆｊ，ｉの関係が成り立つ。そこで、各ペアｐに対して算出する接触力Ｆｐは、Ｆｐ＝σｐ ｉ，ｊＦｉ，ｊとすることができる。σｐ ｉ，ｊは、ペア＜ｉ，ｊ＞への射影演算子であり、インバースを有する。具体的には例えば、接触力演算部１６は、ペアｐに対して算出する接触力Ｆｐを、小さい粒子番号ｉを有する粒子が大きい粒子番号ｊを有する粒子から受ける接触力として算出する。接触力演算部１６は、特許文献１に記載されているように、計算した接触力Ｆｐをペア番号ｐに対応付けられた配列で記憶する。この情報は、総和演算部１７によって参照される。
シミュレーションの高速化のため、特許文献１に記載されているように接触力演算部１６による接触判定及び接触力の計算は、ペア番号（特許文献１では接触候補ペアリスト番号）をスレッド化して行われる。
総和演算部１７は、参照用情報生成部１５によって生成された参照用情報に基づいて、接触力演算部１６によって計算された接触力から粒子毎の接触力の総和を計算する総和演算手段である。粒子番号ｉを有する各粒子に生じる粒子の接触力の総和Ｆｉは、以下の式によって算出できる。

上記の式においてチルダ付きのσｐ ｉ，ｊは、σｐ ｉ，ｊのインバースである。σｐ ｉ，ｊのインバースは、上記のようにＦｐを小さい粒子番号ｉを有する粒子が大きい粒子番号ｊを有する粒子から受ける接触力とした場合には、ｉ＜ｊであれば１（Ｆｐの順符号）であり、ｊ＜ｉであれば−１（Ｆｐの逆符号）である。
具体的には、総和演算部１７は、参照用情報生成部１５によって生成された参照用情報Ｒ（１），Ｒ（２）を参照して、総和の計算対象となる粒子のペアを特定する。続いて、総和演算部１７は、当該ペアについて接触力演算部１６によって計算された接触力Ｆｐを取得する（読み出す）。総和演算部１７は、取得した接触力から、上記の式に基づいて当該粒子に生じる接触力の総和Ｆｉを算出する。なお、σｐ ｉ，ｊのインバースの値は、参照用情報Ｒ（１），Ｒ（２）の何れかからペアが特定されたかによって特定する。参照用情報Ｒ（１）からペアが特定された場合、ｉ＜ｊであるのでσｐ ｉ，ｊのインバースは１となる。参照用情報Ｒ（２）からペアが特定された場合、ｉ＞ｊであるのでσｐ ｉ，ｊのインバースは−１となる。総和演算部１７は、算出した各粒子の接触力の総和Ｆｉを示す情報を粒子情報算出部１８に出力する。
粒子情報算出部１８は、総和演算部１７によって計算された粒子毎の接触力の総和に基づいて、次の時刻（ステップ）における粒子の位置及び速度を算出する粒子情報算出手段である。具体的には、粒子毎に粒子情報保持部１１に保持された粒子情報を取得し、粒子情報に示される現時刻（現ステップ）の座標及び速度と接触力の総和とから、次の時刻（ステップ）の座標及び速度を算出する。この算出は、例えば、特許文献１に示されるようにｌｅａｐ−ｆｌｏｇ法を用いて行うことができる。粒子情報算出部１８は、算出した次の時刻（ステップ）における粒子の位置及び速度で、粒子毎に粒子情報保持部１１に保持された粒子情報を更新する。粒子情報算出部１８によって、全ての粒子についての更新が行われると、次の時刻（ステップ）の処理が行われる。
粒子シミュレーション装置１０では、１つの時刻（ステップ）での演算が完了する毎に、シミュレーションの終了条件を満たしているか否かが判断される。例えば、予め設定した回数（時刻（ステップ））の演算が終了した場合、終了条件を満たしていると判断される。終了条件を満たしていると判断された場合には、粒子シミュレーション装置１０では、シミュレーションが終了される。この場合、例えば、表示装置や他の装置への演算結果の出力等が行われる。終了条件を満たしていないと判断された場合には、次の時刻（ステップ）の演算が繰り返し行われる。以上が、粒子シミュレーション装置１０の構成である。
引き続いて、図４及び図５のフローチャートを用いて、本実施形態に係る粒子シミュレーション装置１０の動作方法である、粒子シミュレーション装置１０で実行される処理（粒子シミュレーション方法）を説明する。図４を用いて本処理全体の処理を説明し、図５を用いて参照用情報Ｒ（２）の生成を説明する。本処理は、例えば、粒子シミュレーション装置１０のユーザが、粒子シミュレーション装置１０に対して解析を開始する操作を行うことで開始される。
粒子シミュレーション装置１０では、位置情報取得部１２によって、粒子情報保持部１１に保持された各粒子の現時刻の粒子情報が取得される（Ｓ０１、位置情報取得ステップ）。粒子情報によって示される粒子の座標を含むセルのセル番号が、当該粒子の位置情報とされる。続いて、位置情報取得部１２によって、ペアの設定（更新）が必要であるか否かが判定される（Ｓ０２）。なお、最初の時刻の処理では、この判定は行われず、必ずペアの設定が行われる。
ペアの設定（更新）が必要であると判定された場合（Ｓ０２のＹＥＳ）には、位置情報取得部１２から、粒子の位置情報が粒子番号設定部１３及びペア設定部１４に出力される。また、位置情報取得部１２から、粒子情報が接触力演算部１６に出力される。続いて、粒子番号設定部１３によって、各粒子に対して、位置情報に基づいて粒子番号が設定される（２度目以降の設定の場合、付け直される）（Ｓ０３、粒子番号設定ステップ）。設定された粒子番号は、粒子番号設定部１３からペア設定部１４に出力される。
引き続いて、ペア設定部１４によって、各粒子に対してペアとなる粒子が決定され、各粒子の粒子番号ｉ毎に自身の粒子番号よりも大きい粒子番号を持つと共にペアとなる粒子の粒子番号を格納した配列Ｕ（１）が生成される（Ｓ０４、ペア設定ステップ）。続いて、ペア設定部１４によって、粒子のペアに対してペア番号が設定され、粒子番号ｉ毎の自身の粒子番号よりも大きい粒子番号を持つ粒子とのペアのペア番号を格納した行列Ｒ（１）が生成される（Ｓ０５、ペア設定ステップ）。ペア設定部１４によって生成された情報は、参照用情報生成部１５に出力される。
続いて、参照用情報生成部１５によって、自身の粒子番号ｉよりも小さい粒子番号（ｉ＞ｊ）を持つ粒子ｊとのペアのペア番号を参照するための参照用情報である行列Ｒ（２）が生成される（Ｓ０６、参照用情報生成ステップ）。
参照用情報生成部１５による行列Ｒ（２）の生成処理について、図６のフローチャートを用いて、より詳細に説明する。本処理では、ペア設定部１４から入力されたＵ（１）及びＲ（１）から行列Ｔ（１）が生成される（Ｓ６１）。続いて、行列Ｔ（１）の第３列の粒子番号ｊに基づいて、行列Ｔ（１）の行がソートされて、行列Ｔ（２）が生成される（Ｓ６２）。続いて、行列Ｔ（２）の第３列の粒子番号ｊ毎の要素数ｎｉ＞ｊ、及び行列Ｔ（２）の第３列（Ｔ（２） ３）のプレフィックス和ｓｉ＞ｊが算出される（Ｓ６３）。続いて、行列Ｔ（２）、ｎｉ＞ｊ及びｓｉ＞ｊ（ｉ）から行列Ｒ（２）が生成される（Ｓ６４）。以上が、行列Ｒ（２）の生成処理である。
図４に戻り、続いて、生成された参照用情報であるＲ（１）及びＲ（２）は、参照用情報生成部１５から総和演算部１７に出力される。また、ペア番号ｐと当該ペア番号ｐによって示されるペアを構成する粒子の粒子番号ｉ，ｊとの対応が特定可能な情報（例えば、行列Ｔ（１））が、参照用情報生成部１５から接触力演算部１６に出力される。
Ｓ０２において、ペアの設定（更新）が必要でないと判定された場合（Ｓ０２のＮＯ）、及びＳ０６の後、続いて、接触力演算部１６によって、各ペアについて、ペアを構成する粒子同士の接触判定が行われる（Ｓ０７、相互作用力演算ステップ）。なお、ペアの更新がされなかった場合には、本処理では一つ前の時刻（ステップ）のペアの情報が用いられる。続いて、接触力演算部１６によって、粒子同士が接触していると判定されたペアについて接触力の計算が行われる（Ｓ０８、相互作用力演算ステップ）。なお、粒子同士が接触していないと判定されたペアについては、接触力はゼロとされる。計算された接触力はペア番号ｐに対応付けられた配列で記憶され、総和演算部１７によって参照される。
続いて、総和演算部１７によって、参照用情報生成部１５によって生成された参照用情報Ｒ（１）、Ｒ（２）が参照されて、接触力演算部１６によって計算された接触力から粒子毎の接触力の総和が計算される（Ｓ０９、総和演算ステップ）。算出された各粒子の接触力の総和Ｆｉを示す情報は、総和演算部１７から粒子情報算出部１８に出力される。
続いて、粒子情報算出部１８によって、総和演算部１７によって計算された粒子毎の接触力の総和に基づいて、次の時刻（ステップ）における粒子の位置及び速度が算出される（Ｓ１０、粒子情報算出ステップ）。続いて、粒子情報算出部１８によって算出された粒子の位置及び速度で粒子毎に粒子情報保持部１１に保持された粒子情報が更新される。
続いて、粒子シミュレーション装置１０では、シミュレーションの終了条件を満たしているか否かが判断される（Ｓ１１）。終了条件を満たしていると判断された場合（Ｓ１１のＹＥＳ）には、処理（シミュレーション）が終了される。終了条件を満たしていないと判断された場合（Ｓ１１のＮＯ）には、時刻（ステップ）が一つ進められて、次の時刻（ステップ）での上述した処理（Ｓ０１〜Ｓ１１）が行われる。以上が、本実施形態に係る粒子シミュレーション装置１０で実行される処理である。
ここで、本発明の比較対象となる非特許文献１に記載されたシミュレーション方法を説明する。非特許文献１に記載された方法においても、本実施形態と同様に作業領域のセルへの分割が行われており、また、各粒子には各粒子の位置に応じた粒子番号が設定されている。
非特許文献１に記載された方法では、図６（ａ）に示す相互作用行列Ｕを構成する（ことを想定する）。この行列Ｕは、本実施形態の図３（ａ）に示す相互作用行列Ｕに相当するものである。但し、非特許文献１に記載された方法では、相互作用行列Ｕには、本実施形態とは異なり下三角部分も含まれる。
非特許文献１に記載された方法では、格納するメモリスペースを考慮し、相互作用行列Ｕを図６（ｂ）に示す配列Ｕ（１）及びＵ（２）を生成する。配列Ｕ（１）は、本実施形態におけるＵ（１）と同じものである。配列Ｕ（２）は、各粒子の粒子番号ｉ毎に自身の粒子番号よりも小さい粒子番号を持つと共にペアとなる粒子の粒子番号を格納した配列である。これらの行列は、隣接する２７セルを探索することで構成される。また、同時に各粒子の粒子番号ｉ毎に自身の粒子番号よりも大きい粒子番号を持つと共にペアとなる粒子の数（配列Ｕ（１）のｉ毎の要素の数）ｎｉ＜ｊ、及び自身の粒子番号よりも小さい粒子番号を持つと共にペアとなる粒子の数（配列Ｕ（２）のｉ毎の要素の数）ｎｉ＞ｊを算出する。
非特許文献１に記載された方法では、続いて、例えば、図６（ａ）に示す相互作用行列Ｕの１の要素に順番の番号を振ることで図６（ｃ）に示すように粒子間のペアにペア番号を設定する（ことを想定する）。しかしながら上記のＵは実際には生成せず、本実施形態と同様に配列Ｕ（１）とプレフィックス和から、図６（ｄ）に示す行列Ｒ（１）を生成する。行列Ｒ（１）は、本実施形態におけるＲ（１）と同じものである。
続いて、Ｒ（２）が、ｉ番目のスレッドがＲ（２）に対応するＲ（１） ｉ，ｊの要素を探索することで生成される。このとき、ｉ番目のスレッドは、Ｕ（１） ｉ，ｊの値を読み出し、Ｕ（２）のＵ（１） ｉ，ｊ番目の行を探索する。スレッドは、値がｉとなる要素を見つけた場合にＲ（１） ｉ，ｊの要素をＲ（２）の要素とする。このプロセスがｊについて繰り返されて、Ｒ（２）が生成される。
即ち、非特許文献１での、行列Ｒ（２）の生成は、以下に示すアルゴリズムで実現することができる。行列Ｒ（２）は、本実施形態におけるＲ（２）と同じものである。
非特許文献１に記載された方法では、本実施形態と同様にＲ（１）、Ｒ（２）が用いられて粒子シミュレーションが行われる。
引き続いて、非特許文献１に記載された方法に対して行った性能試験の結果を示す。この性能試験では、作業領域であるシミュレーションボックス内に均一に粒子を配置し、Ｒ（１）、Ｒ（２）であるペアインデックスマトリクスを求めた。ここで、シミュレーションボックスはＬ３のサイズの三次元空間とし、図７（ａ）に示すように均一に粒子を配置した。シミュレーションボックスは、一辺ｌ＝Ｌ／４である立方体のセルに分割される。カットオフ長ｒｃは、ｌよりもわずかに小さくした。
この性能試験では、処理を３つのステップに分けている。ステップ１はセル分割及びソート、ステップ２はペアインデックスマトリクスの上三角部分の構築、ステップ３はペアインデックスマトリクスの下三角部分の構築である。図７（ｂ）に性能試験の結果を示す。図７（ｂ）のグラフはステップ毎の演算処理の実行にかかった経過時間である。この性能試験に示されるように他のステップに比べてステップ３に多くの演算処理が必要になっている。
上述したようにステップ３では、ｉ番目のスレッドがＲ（２）に対応するＲ（１） ｉ，ｊの要素を探索する。この処理では、ｉ番目のスレッドはＵ（２）の一つの行を探索し、この探索がＲ（１）の１つの行に対して繰り返される。即ち、非特許文献１に示す方法のアルゴリズムでは、ｉｆ文（条件の判断）を含む、二重のｆｏｒループ（ダブルｆｏｒループ）によって探索を行っている。これには、上述したように、ＧＰＵの計算性能を大幅に低下させるワープダイバージェンスを引き起こす等の問題があった。以上が、非特許文献１に記載されたシミュレーション方法である。
上記の非特許文献１に記載されたシミュレーション方法を踏まえて、本実施形態の効果について説明する。上述したように非特許文献１に記載された方法では、多数のトライアルアンドエラー処理が性能の低下の原因となっている。一方で、本実施形態では、ソートした行列を用いることで、二重のｆｏｒループを用いずに参照用情報が生成される。これにより、多数のトライアルアンドエラー処理が不要となる。即ち、本実施形態によれば、演算効率を向上させることできる。なお、本実施形態のようにペアリストを用いる方法は、アトミック演算が不要な方法であり、接触力の並列計算が可能なものである。
また、本実施形態のように、参照用情報Ｒ（２）を生成するための、３つの列を有する行列Ｔ（１）の列を、ペアを構成する粒子の粒子番号のうち小さい方の粒子番号ｉとし、もう一方の列を、ペアを構成する粒子の粒子番号のうち大きい方の粒子番号ｊとすることとしてもよい。この構成によれば、適切かつ確実にソートした行列Ｔ（２）に基づき、参照用情報Ｒ（２）を生成することができる。これにより、適切かつ確実に本発明を実施することができる。
また、本実施形態のように、位置情報に基づいて、複数の粒子それぞれに対してソート可能な粒子番号を設定することとしてもよい。この構成によれば、適切かつ確実に粒子に粒子番号を付与することができる。これにより、適切かつ確実に本発明を実施することができる。
また、本実施形態のように作業領域をセルに分割してシミュレーションを行うこととしてもよい。この構成によれば、効率的かつ容易に近傍に位置する粒子のペアを選択することができ、演算効率を更に向上させることできる。
なお、本実施形態では、粒子に働く相互作用力として、粒子間の接触による接触力を用いることとしたが、本発明における相互作用は接触には限られない。粒子間に力（相互作用力）が働くものであれば、任意の相互作用を対象にすることとしてもよい。
引き続いて、本実施形態に係る性能試験の結果を示す。当該試験を行うために、粒子位置の静的配列を用意し、それからペアインデックスマトリクスを構築するための計算コストを測定した。システムは、Ｌ３のサイズの三次元のシミュレーションボックスとした。このシミュレーションボックスでは、Ｎ＝１０３個の粒子が不規則に配置される。シミュレーションボックスは、一辺の長さがｌの立方体のセルに分割される。相互作用のカットオフ長ｒｃは、ｌよりわずかに小さく設定される。
本性能試験では、いくつかのｌの値が設定される。即ち、性能における、相互作用マトリクスのサイズの効果が調べられる。本性能試験では、非特許文献１に準じたコード（アルゴリズム）と、本実施形態に準じたコード（アルゴリズム）を実装し、それらを比較した。２つのコードは、それぞれに本質的な部分を除いては同一である。本コードは、ＣＵＤＡにより実装され、ＴＥＳＬＡ Ｃ２０７５上で試験された。本実施形態における重要な一連の動作は、ソーティング、カウンティング及びプレフィックス和である。それらの一連の動作を実現するため、ソーティング及びプレフィックス和については既存のＴｈｒｕｓｔライブラリを用いて実装した。性能は、ペアインデックスマトリクスを構築するためにかかった経過時間として測定された。
図８に、２つのコードの測定結果を示す。図５は、粒子密度と、ペアインデックスマトリクスを生成するためにかかった経過時間との関係を示す図である。経過時間は、セル分割、ソーティング、上三角部分及び下三角部分のペアインデックスマトリクスの構築にかかった時間の合計である。図５に示されるように、２つのコードの性能は粒子密度が小さい場合には同程度であるが、粒子密度が大きくなると本実施形態によるものの性能が極めて高くなる。この結果は、行列のサイズが大きくなると、本実施形態の性能を改善する効果が大きくなることを示している。本発明の手法は、オブジェクトが大きくなるとより効果が大きくなる分割統治技術を用いる。粒子密度が１である場合、ソーティングのオーバーヘッドが生じるため、性能の優位性はわずかに逆転する。
シミュレーションのタイプを理解することで、本発明の手法による改善をおこなうことができる。図９にＤＥＭ及びＳＰＨシミュレーションで仮定される特有なシミュレーションを示す。図９（ａ）に示すように、ＤＥＭでは、粒子が殻周辺のみで作用する強い接触力で別の粒子との間で相互作用するため、各セルはわずかの粒子しか含まない。一方で、図９（ｂ）に示すように、ＳＰＨの粒子は、非常に柔軟に作用する殻を有しており、多数の他の粒子との間で相互作用する。安定計算のためには、水力学の粒子のカットオフ半径は、平均粒子間距離の２〜４倍に設定されることが望ましい。これは、ＳＰＨシミュレーションで仮定される粒子密度は約十〜数十であり、本実施形態のアルゴリズムはこの範囲の密度で非常に効率的である。従って、本発明は、ＤＥＭよりもＳＰＨやＭＤのような多くの相互作用ペアを含む粒子システムを改善する。
更に、ＤＥＭにおいても、粒子の半径が大きく分散している場合には、本発明が効果的である。図９（ｃ）に示すように、例えば、ブラジルナッツ問題の研究では、最も大きな粒子のサイズが、最も小さな粒子のサイズの数倍とされる。最大の半径をｒｍａｘ、最小の半径をｒｍｉｎとすると、セル中の粒子密度はたかだか（ｒｍａｘ／ｒｍｉｎ）３となる。均一のセルを有する粒子システムをシミュレーションする際に、セルサイズは、最大の粒子を収容できるようにしなければならないためである。本実施形態のアルゴリズムは、このタイプのシミュレーションにおいて明らかに効率的である。
本実施形態に係る性能をより詳細に示す。通常、ＳＰＨシミュレーションで仮定されるように、ここではＬ／ｌ＝４とした。図１０（ａ）に、非特許文献１に準じたコードと、本実施形態に準じたコードとにおける、ペアインデックスマトリクスの構築に係る経過時間を示す。ステップ２及びステップ３において、性能が改善されている。ステップ２における性能改善は、探索範囲の減少によるものである。検索される近傍セルの数は、非特許文献１に準じたコードでは２７、本実施形態に準じたコードでは１４である。
極めて大きい改善がステップ３でなされている。これは、多数のトライアルアンドエラー処理によるものであると考えられ、これは、図１０（ｂ）に示すプロファイリング結果による裏付けられる。図１０（ｂ）は、ロード（メモリ読出）及びストア（格納）指示の数、並びにプロファイラによってカウントされた分岐（ダイバージェントブランチ）の数である。処理全体における、メモリ読出の指示の数は明らかに減少している。また、分岐の数の劇的な減少は、本実施形態に準じたコードがワープダイバージェンスを適切に回避していることを示している。上記の通り、本実施形態では、ペアインデックスマトリクスの構築の高速化に成功している。
また、以下のように、粒子シミュレーション全体に対して本発明の有効性を検証した。具体的には、ペアリストの構築、力の演算、力の合算、及び運動方程式の積算を含むＳＰＨシミュレーションを実行した。粒子は運動方程式に従い、粒子には粒子間の接触力ｆｉ，ｊに加えて、粒子に働く外力ｆｅｘｔも考慮する。本シミュレーションでは、粒子は重力による自由落下を行うものとした。即ち、ｆｅｘｔ＝ρｇである。ここで、ρはＳＰＨと合わせて計算される質量密度であり、ｇは重力加速度ベクトルである。応力と密度のラグラジアン形式での表現に、M. Muller, D. Charypar, M. Gross, Particle-based fluid simulationfor interactive applications, in: Proceedings of 2003 ACM SIGGRAPHSymposium,2003, pp. 154−159.に示される枠組みを用いた。上記の枠組みは、流体のシミュレーションの実行には不向きであるが、当該枠組みは極めてシンプルであるためアルゴリズムの改良の検証に便利である。
Ｌ＝１０ｍ、ｌ＝０．１ｍ、カーネル関数のカットオフ長ｒ＝０．０８ｍとした。図１１に示すように初期状態において、Ｎ＝１０６の粒子を柱状に配置し、粒子の間隔をｒｃ以上とした。運動方程式は、シンプレクティックスキームにより離散化され、時間間隔Δｔ＝５．０×１０−４ｓとした。ｔ＝５．０ｓとなるまで繰り返した。
システムの特徴を示す重要な指標の一つである有効粒子密度ρ＊は以下の式によって定義される。

本シミュレーションでは、粒子の質量等のパラメータは、安定有効密度が〜３．５／セルとなるように選択された。図１２（ａ）に示すように、シミュレーション中の有効密度の動きを測定した。図１２（ａ）に示すように、有効密度は、ｔが１．３程度で最大値をとり、流体は、ほぼ重力落下を介して圧縮された。ペアリストは、最大積算移動距離Δｘが臨界値を超えた場合、ｍａｘ（Δｘ）＞（ｌ−ｒｃ）／２で更新した。
また、図１２（ｂ）に示すように、各時刻（ステップ）における、力の計算、総和及び時間積算の実行にかかった経過時間を測定した。また、図１２（ｂ）のグラフに、ペアリストの構築にかかった経過時間を測定し、１００時刻（ステップ）の移動平均を示す。図１２（ｂ）のグラフは、有効密度が最大値を取るｔが１．３程度では、ペアリストの構築に最も時間がかかっており、本実施形態では、ペアリスト構築にかかる時間が、非特許文献１に記載された方法よりも短い時間となっている。このように、本実施形態は、特に重力の衝撃を介した流体の圧縮において、特に性能を改善している。
また、本実施形態は、大きな性能の低下をもたらさずに、粒子のペアリストを頻繁に更新可能とするものである。上記の結果は、本実施形態が、衝撃を介した大きな圧縮を含むシミュレーションを高速化できることを示すものである。これは、シミュレーションが、乱流といった粒子の激しい動きを含むものであっても、本実施形態が効率的であることを示している。更に、本例では、安定有効密度をいくらか小さい値（ρ＊〜３．５／セル）を選択していた。安定したＳＰＨシミュレーションを実行するために、より大きな安定密度をとるべきである。大きな密度は、大きな圧力振動等の不安定さを避けることができるためである。もし、大きな安定密度を用いれば、本実施形態の改善はより明らかになる。このように、本実施形態は、粒子シミュレーションの高速化に有効である。
引き続いて、上述した一連の粒子シミュレーション装置１０による処理をコンピュータに実行させるための粒子シミュレーションプログラムを説明する。図１３に示すように、粒子シミュレーションプログラム３０は、コンピュータに挿入されてアクセスされる、あるいはコンピュータが備える記録媒体２０に形成されたプログラム格納領域２１内に格納される。
粒子シミュレーションプログラム３０は、粒子情報保持モジュール３１と、位置情報取得モジュール３２と、粒子番号設定モジュール３３と、ペア設定モジュール３４と、参照用情報生成モジュール３５と、接触力演算モジュール３６と、総和演算モジュール３７と、粒子情報算出モジュール３８とを備えて構成される。粒子情報保持モジュール３１と、位置情報取得モジュール３２と、粒子番号設定モジュール３３と、ペア設定モジュール３４と、参照用情報生成モジュール３５と、接触力演算モジュール３６と、総和演算モジュール３７と、粒子情報算出モジュール３８とを実行させることにより実現される機能は、上述した粒子シミュレーション装置１０の粒子情報保持部１１と、位置情報取得部１２と、粒子番号設定部１３と、ペア設定部１４と、参照用情報生成部１５と、接触力演算部１６と、総和演算部１７と、粒子情報算出部１８とそれぞれ同様である。
なお、粒子シミュレーションプログラム３０は、その一部又は全部が、通信回線等の伝送媒体を介して伝送され、他の機器により受信されて記録（インストールを含む）される構成としてもよい。また、粒子シミュレーションプログラム３０の各モジュールは、１つのコンピュータでなく、複数のコンピュータのいずれかにインストールされてもよい。その場合、当該複数のコンピュータによるコンピュータシステムよって上述した一連の粒子シミュレーションプログラム３０の処理が行われる。
１０…粒子シミュレーション装置、１１…粒子情報保持部、１２…位置情報取得部、１３…粒子番号設定部、１４…ペア設定部、１５…参照用情報生成部、１６…接触力演算部、１７…総和演算部、１８…粒子情報算出部、２０…記録媒体、２１…プログラム格納領域、３０…粒子シミュレーションプログラム、３１…粒子情報保持モジュール、３２…位置情報取得モジュール、３３…粒子番号設定モジュール、３４…ペア設定モジュール、３５…参照用情報生成モジュール、３６…接触力演算モジュール、３７…総和演算モジュール、３８…粒子情報算出モジュール。

Claims (6)

  1. 作業空間内の複数の粒子について他の粒子との相互作用力に基づき位置及び速度を算出し、粒子の挙動をシミュレーションする粒子シミュレーション装置であって、
    前記複数の粒子それぞれについて、粒子の位置を示す位置情報を取得する位置情報取得手段と、
    前記複数の粒子それぞれに対してソート可能な粒子番号を設定する粒子番号設定手段と、
    前記位置情報取得手段によって取得された位置情報に基づいて、互いに近傍に位置する粒子のペアを選択すると共に、前記粒子番号設定手段によって設定された、当該ペアを構成する粒子の一方の粒子番号に基づき当該ペアに対してペア番号を設定するペア設定手段と、
    前記ペア設定手段によって設定されたペア番号、及び当該ペアを構成する粒子の粒子番号を行の構成要素とする行列を生成し、当該ペアを構成する粒子のもう一方の粒子番号に基づいて当該行列の列の順序をソートし、ソートした行列に基づき、粒子の粒子番号から当該粒子が構成するペアのペア番号を参照するための参照用情報を生成する参照用情報生成手段と、
    前記ペア設定手段によって選択されたペアそれぞれに係る粒子同士の相互作用判定を行い、相互作用していると判定された粒子同士の相互作用力を計算する相互作用力演算手段と、
    前記参照用情報生成手段によって生成された参照用情報に基づいて、前記相互作用力演算手段によって計算された相互作用力から粒子毎の相互作用力の総和を計算する総和演算手段と、
    前記総和演算手段によって計算された粒子毎の相互作用力の総和に基づいて、粒子の位置及び速度を算出する粒子情報算出手段と、
    を備える粒子シミュレーション装置。
  2. 前記参照用情報生成手段は、前記行列の一つの列の構成要素を、前記ペアを構成する粒子の粒子番号のうち小さい方の粒子番号とし、別の列の構成要素を、前記ペアを構成する粒子の粒子番号のうち大きい方の粒子番号とする請求項１に記載の粒子シミュレーション装置。
  3. 前記粒子番号設定手段は、前記位置情報取得手段によって取得された位置情報に基づいて、前記複数の粒子それぞれに対してソート可能な粒子番号を設定する請求項１又は２に記載の粒子シミュレーション装置。
  4. 前記作業空間は、複数のセルに分割されており、
    前記位置情報取得手段は、前記位置情報として、前記複数の粒子それぞれについて、粒子が位置するセルを示す情報を取得する、請求項１〜３の何れか一項に記載の粒子シミュレーション装置。
  5. 作業空間内の複数の粒子について他の粒子との相互作用力に基づき位置及び速度を算出し、粒子の挙動をシミュレーションする粒子シミュレーション装置の動作方法である粒子シミュレーション方法であって、
    前記複数の粒子それぞれについて、粒子の位置を示す位置情報を取得する位置情報取得ステップと、
    前記複数の粒子それぞれに対してソート可能な粒子番号を設定する粒子番号設定ステップと、
    前記位置情報取得ステップにおいて取得された位置情報に基づいて、互いに近傍に位置する粒子のペアを選択すると共に、前記粒子番号設定ステップにおいて設定された、当該ペアを構成する粒子の一方の粒子番号に基づき当該ペアに対してペア番号を設定するペア設定ステップと、
    前記ペア設定ステップにおいて設定されたペア番号、及び当該ペアを構成する粒子の粒子番号を行の構成要素とする行列を生成し、当該ペアを構成する粒子のもう一方の粒子番号に基づいて当該行列の列の順序をソートし、ソートした行列に基づき、粒子の粒子番号から当該粒子が構成するペアのペア番号を参照するための参照用情報を生成する参照用情報生成ステップと、
    前記ペア設定ステップにおいて選択されたペアそれぞれに係る粒子同士の相互作用判定を行い、相互作用していると判定された粒子同士の相互作用力を計算する相互作用力演算ステップと、
    前記参照用情報生成ステップにおいて生成された参照用情報に基づいて、前記相互作用力演算ステップにおいて計算された相互作用力から粒子毎の相互作用力の総和を計算する総和演算ステップと、
    前記総和演算ステップにおいて計算された粒子毎の相互作用力の総和に基づいて、粒子の位置及び速度を算出する粒子情報算出ステップと、
    を含む粒子シミュレーション方法。
  6. コンピュータを、作業空間内の複数の粒子について他の粒子との相互作用力に基づき位置及び速度を算出し、粒子の挙動をシミュレーションする粒子シミュレーション装置として機能させる粒子シミュレーションプログラムであって、
    前記コンピュータを、
    前記複数の粒子それぞれについて、粒子の位置を示す位置情報を取得する位置情報取得手段と、
    前記複数の粒子それぞれに対してソート可能な粒子番号を設定する粒子番号設定手段と、
    前記位置情報取得手段によって取得された位置情報に基づいて、互いに近傍に位置する粒子のペアを選択すると共に、前記粒子番号設定手段によって設定された、当該ペアを構成する粒子の一方の粒子番号に基づき当該ペアに対してペア番号を設定するペア設定手段と、
    前記ペア設定手段によって設定されたペア番号、及び当該ペアを構成する粒子の粒子番号を行の構成要素とする行列を生成し、当該ペアを構成する粒子のもう一方の粒子番号に基づいて当該行列の列の順序をソートし、ソートした行列に基づき、粒子の粒子番号から当該粒子が構成するペアのペア番号を参照するための参照用情報を生成する参照用情報生成手段と、
    前記ペア設定手段によって選択されたペアそれぞれに係る粒子同士の相互作用判定を行い、相互作用していると判定された粒子同士の相互作用力を計算する相互作用力演算手段と、
    前記参照用情報生成手段によって生成された参照用情報に基づいて、前記相互作用力演算手段によって計算された相互作用力から粒子毎の相互作用力の総和を計算する総和演算手段と、
    前記総和演算手段によって計算された粒子毎の相互作用力の総和に基づいて、粒子の位置及び速度を算出する粒子情報算出手段と、
    として機能させる粒子シミュレーションプログラム。
JP2014115567A 2014-06-04 2014-06-04 粒子シミュレーション装置、粒子シミュレーション方法及び粒子シミュレーションプログラム Expired - Fee Related JP6261130B2 (ja)

Priority Applications (4)

Application Number Priority Date Filing Date Title
JP2014115567A JP6261130B2 (ja) 2014-06-04 2014-06-04 粒子シミュレーション装置、粒子シミュレーション方法及び粒子シミュレーションプログラム
PCT/JP2015/065609 WO2015186633A1 (ja) 2014-06-04 2015-05-29 粒子シミュレーション装置、粒子シミュレーション方法及び粒子シミュレーションプログラム
US15/315,254 US10354099B2 (en) 2014-06-04 2015-05-29 Particle simulation device, particle simulation method, and particle simulation program
EP15803426.4A EP3153982A4 (en) 2014-06-04 2015-05-29 Particle simulation device, particle simulation method, and particle simulation program

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
JP2014115567A JP6261130B2 (ja) 2014-06-04 2014-06-04 粒子シミュレーション装置、粒子シミュレーション方法及び粒子シミュレーションプログラム

Publications (2)

Publication Number Publication Date
JP2015230535A JP2015230535A (ja) 2015-12-21
JP6261130B2 true JP6261130B2 (ja) 2018-01-17

Family

ID=54766706

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2014115567A Expired - Fee Related JP6261130B2 (ja) 2014-06-04 2014-06-04 粒子シミュレーション装置、粒子シミュレーション方法及び粒子シミュレーションプログラム

Country Status (4)

Country Link
US (1) US10354099B2 (ja)
EP (1) EP3153982A4 (ja)
JP (1) JP6261130B2 (ja)
WO (1) WO2015186633A1 (ja)

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP6671064B2 (ja) * 2016-03-03 2020-03-25 国立研究開発法人海洋研究開発機構 粒子シミュレーション装置、粒子シミュレーション方法及び粒子シミュレーションプログラム
CN109325250A (zh) * 2018-07-26 2019-02-12 四川大学 一种高速滑坡-碎屑流运动侵蚀效应的数值模拟方法及系统
CN109977505B (zh) * 2019-03-13 2024-02-09 南京大学 基于gpu矩阵计算的离散元孔隙系统快速搜索方法
JP7244757B2 (ja) * 2019-07-01 2023-03-23 富士通株式会社 情報処理装置、粒子シミュレーション方法および粒子シミュレーションシステム
CN113420410B (zh) * 2021-05-20 2024-08-27 华北水利水电大学 基于pfc的岩体非均质模型的构建与应用
CN113656656B (zh) * 2021-07-21 2024-01-26 南京南力科技有限公司 用于宽级配离散元颗粒系统的高效邻居检索方法及系统
JP2023127881A (ja) 2022-03-02 2023-09-14 富士通株式会社 シミュレーションプログラム、シミュレーション方法及びシミュレーション装置
WO2023209869A1 (ja) * 2022-04-27 2023-11-02 富士通株式会社 構造抽出プログラム、構造抽出方法および情報処理装置

Family Cites Families (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5596511A (en) 1993-01-27 1997-01-21 Fuji Xerox Co., Ltd. Computing method and apparatus for a many-body problem
EP2447852A3 (en) 2005-04-19 2015-01-07 D.E. Shaw Research, LLC Scalable method for the evaluation of distance-limited pairwise particle interactions
US7948485B1 (en) * 2005-12-12 2011-05-24 Sony Computer Entertainment Inc. Real-time computer simulation of water surfaces
JP5408611B2 (ja) * 2009-03-31 2014-02-05 独立行政法人海洋研究開発機構 粒子シミュレーション装置及び粒子シミュレーション方法
US8762117B2 (en) 2009-06-03 2014-06-24 Algoryx Simulation Ab Method, an apparatus and computer program product for simulating dynamic fluids
BR112013006284B1 (pt) * 2010-09-15 2021-06-15 Commonwealth Scientific And Industrial Research Organisation Método de elemento discreto e sistema de simulação de partícula
EP2509009A1 (en) 2011-04-05 2012-10-10 Japan Agency for Marine-Earth Science and Technology Particle simulator and method of simulating particles
US8554527B2 (en) * 2011-04-08 2013-10-08 Japan Agency For Marine-Earth Science And Technology Particle simulator and method of simulating particles
JP5844165B2 (ja) 2012-01-26 2016-01-13 住友重機械工業株式会社 解析装置
US8855982B2 (en) 2012-02-06 2014-10-07 Sumitomo Heavy Industries, Ltd. Analysis device and simulation method

Also Published As

Publication number Publication date
EP3153982A4 (en) 2018-10-10
EP3153982A1 (en) 2017-04-12
WO2015186633A1 (ja) 2015-12-10
US20170193251A1 (en) 2017-07-06
US10354099B2 (en) 2019-07-16
JP2015230535A (ja) 2015-12-21

Similar Documents

Publication Publication Date Title
JP6261130B2 (ja) 粒子シミュレーション装置、粒子シミュレーション方法及び粒子シミュレーションプログラム
Ashari et al. Fast sparse matrix-vector multiplication on GPUs for graph applications
Harmon et al. Subspace integration with local deformations
Wu et al. A comparative measurement study of deep learning as a service framework
Brown et al. Implementing molecular dynamics on hybrid high performance computers–Particle–particle particle-mesh
Schäfer et al. A smooth particle hydrodynamics code to model collisions between solid, self-gravitating objects
Rahimian et al. Petascale direct numerical simulation of blood flow on 200k cores and heterogeneous architectures
Esteves et al. Competitive k-means, a new accurate and distributed k-means algorithm for large datasets
AU2017227323B2 (en) Particle simulation device, particle simulation method, and particle simulation program
Lubbe et al. Analysis of parallel spatial partitioning algorithms for GPU based DEM
Rubio-Largo et al. Granular gas of ellipsoids: analytical collision detection implemented on GPUs
EP2808812A1 (en) Analysis device and simulation method
Shamseddine et al. A novel spatio-temporally adaptive parallel three-dimensional DSMC solver for unsteady rarefied micro/nano gas flows
Weatherley et al. Scaling benchmark of esys-particle for elastic wave propagation simulations
Hylton et al. Tuning the performance of a computational persistent homology package
US20230145348A1 (en) Force-directed graph layout
JP2011145943A (ja) 粒子シミュレーション装置及び粒子シミュレーション方法
Chappell et al. Wind-driven gas networks and star formation in galaxies: reaction-advection hydrodynamic simulations
Al-Madi et al. Scaling genetic programming for data classification using mapreduce methodology
Bülow et al. A scalable parallel Stokesian Dynamics method for the simulation of colloidal suspensions
Playne et al. Benchmarking multi-GPU communication using the shallow water equations
JP7285223B2 (ja) シミュレーション装置、シミュレーション方法、及びプログラム
Liu et al. A bubble‐inspired algorithm for finite element mesh partitioning
Lin et al. A parallel and scalable CAST-based clustering algorithm on GPU
EP4131084A1 (en) Program, data processing method, and data processing apparatus

Legal Events

Date Code Title Description
A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20170426

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: 20171205

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20171211

R150 Certificate of patent or registration of utility model

Ref document number: 6261130

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

LAPS Cancellation because of no payment of annual fees