(57)【特許請求の範囲】
【請求項１】
３次元の解析対象領域において接触判定格子を用いて接触判定計算装置が行う粒子の接触判定計算方法であって、
前記接触判定計算装置は、粒子の接触判定計算を行う接触判定計算部を備え、
前記接触判定計算部は、
格子番号が（２ｉ−１，２ｉ−１，２ｉ−１）番の接触判定格子における粒子の接触判定計算を並列に行うステップと、
格子番号が（２ｉ，２ｉ−１，２ｉ−１）番の接触判定格子における粒子の接触判定計算を並列に行うステップと、
格子番号が（２ｉ−１，２ｉ，２ｉ−１）番の接触判定格子における粒子の接触判定計算を並列に行うステップと、
格子番号が（２ｉ，２ｉ，２ｉ−１）番の接触判定格子における粒子の接触判定計算を並列に行うステップと、
格子番号が（２ｉ−１，２ｉ−１，２ｉ）番の接触判定格子における粒子の接触判定計算を並列に行うステップと、
格子番号が（２ｉ，２ｉ−１，２ｉ）番の接触判定格子における粒子の接触判定計算を並列に行うステップと、
格子番号が（２ｉ−１，２ｉ，２ｉ）番の接触判定格子における粒子の接触判定計算を並列に行うステップと、
格子番号が（２ｉ，２ｉ，２ｉ）番の接触判定格子における粒子の接触判定計算を並列に行うステップと、
を所定の順序で繰り返し行う（ｉは１以上の整数）ことを特徴とする、粒子の接触判定計算方法。
【請求項２】
３次元の解析対象領域において接触判定格子を用いて行う粒子の接触判定計算装置であって、
粒子の接触判定計算を行う接触判定計算部を備え、
前記接触判定計算部は、
格子番号が（２ｉ−１，２ｉ−１，２ｉ−１）番の接触判定格子における粒子の接触判定計算を並列に行うステップと、
格子番号が（２ｉ，２ｉ−１，２ｉ−１）番の接触判定格子における粒子の接触判定計算を並列に行うステップと、
格子番号が（２ｉ−１，２ｉ，２ｉ−１）番の接触判定格子における粒子の接触判定計算を並列に行うステップと、
格子番号が（２ｉ，２ｉ，２ｉ−１）の接触判定格子における粒子の接触判定計算を並列に行うステップと、
格子番号が（２ｉ−１，２ｉ−１，２ｉ）番の接触判定格子における粒子の接触判定計算を並列に行うステップと、
格子番号が（２ｉ，２ｉ−１，２ｉ）番の接触判定格子における粒子の接触判定計算を並列に行うステップと、
格子番号が（２ｉ−１，２ｉ，２ｉ）番の接触判定格子における粒子の接触判定計算を並列に行うステップと、
格子番号が（２ｉ，２ｉ，２ｉ）番の接触判定格子における粒子の接触判定計算を並列に行うステップと、
を所定の順序で繰り返し行う（ｉは１以上の整数）ことを特徴とする、粒子の接触判定計算装置。
【請求項３】
コンピュータを、３次元の解析対象領域において接触判定格子を用いて行う粒子の接触判定計算装置として機能させるためのコンピュータプログラムであって、
前記接触判定計算装置は、粒子の接触判定計算を行う接触判定計算部を備え、
前記接触判定計算部は、
格子番号が（２ｉ−１，２ｉ−１，２ｉ−１）番の接触判定格子における粒子の接触判定計算を並列に行うステップと、
格子番号が（２ｉ，２ｉ−１，２ｉ−１）番の接触判定格子における粒子の接触判定計算を並列に行うステップと、
格子番号が（２ｉ−１，２ｉ，２ｉ−１）番の接触判定格子における粒子の接触判定計算を並列に行うステップと、
格子番号が（２ｉ，２ｉ，２ｉ−１）番の接触判定格子における粒子の接触判定計算を並列に行うステップと、
格子番号が（２ｉ−１，２ｉ−１，２ｉ）番の接触判定格子における粒子の接触判定計算を並列に行うステップと、
格子番号が（２ｉ，２ｉ−１，２ｉ）番の接触判定格子における粒子の接触判定計算を並列に行うステップと、
格子番号が（２ｉ−１，２ｉ，２ｉ）番の接触判定格子における粒子の接触判定計算を並列に行うステップと、
格子番号が（２ｉ，２ｉ，２ｉ）番の接触判定格子における粒子の接触判定計算を並列に行うステップと、
を所定の順序で繰り返し行う（ｉは１以上の整数）ことを特徴とする、コンピュータプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、粒子の接触判定計算方法、粒子の接触判定計算装置、及びコンピュータプログラムにかかり、特に、所定の解析対象領域において接触判定格子を用いて行う、粒子の接触判定計算方法、粒子の接触判定計算装置、及びコンピュータプログラムに関する。
【背景技術】
【０００２】
粉体計算やＭＰＳ（Moving Particle Semi-implicit）法による接触を考慮した粒子解析において、粒子接触判定計算は、全体の時間の多くを占めると言われており、高速化が必須な処理である。特に近年の計算機においては、ノード内に８スレッド以上の並列実行可能なＣＰＵが普及している。したがって、並列化による高速化が課題となっている。
【０００３】
従来、粒子接触判定計算においては、接触判定計算の計算量を減らす方法が使われてきた。接触判定計算を全粒子について行うのではなく、接触する可能性のある近傍の粒子についてのみ物理量の計算を行う。この近傍の粒子を絞り込むために、接触判定格子を用いる方法がある（例えば、非特許文献１参照）。また、キャッシュを有効利用することで並列化時の台数効果を上げるため、粒子に番号を振ったうえで、粒子を接触判定時の演算が局所化するように並び替える方法が知られている（例えば、非特許文献２参照）。
【先行技術文献】
【非特許文献】
【０００４】
【非特許文献１】酒井幹夫編著：粉体の数値シミュレーション，丸善出版 (2012)
【非特許文献２】西浦泰介，阪口秀:ＧＰＵを用いたＤＥＭの高速化アルゴリズム, Transactions ofJSCES, paper No.20100007 (2010)
【発明の概要】
【発明が解決しようとする課題】
【０００５】
ところで、粉体計算で用いられる離散要素法（Discrete Element Method：ＤＥＭ)で必要となる、接触する各粒子に働く接触力の足しこみ計算において、従来は、作用−反作用の法則を用いて演算量の削減を行っている。この従来法をスレッド並列化する際、接触する相手の情報も同時に更新することになるが、この操作は他スレッドが更新する情報の書き換えにほかならない。したがって、複数のスレッド間でデータの読み書きが競合し、排他制御が必要となる。この排他制御はＯｐｅｎＭＰを利用する場合、critical指示文やatomic指示子を入れるのが普通である（例えば、非特許文献１参照）。
【０００６】
一方、critical指示文等を入れない並列化では、接触時に相手の情報を更新しない方法により排他制御を回避してきた。しかしこの方法は、critical指示文を入れる方法に対して、約２倍の演算量の増大になる。つまり、排他制御する方法では演算量が１／２になるメリットがある反面、排他制御によるスレッド並列化のオーバーヘッドが生じる。一方で、排他制御がない方法では演算量が増える反面、スレッド並列化時のオーバーヘッドが少ない。つまり、両者は欠点と利点があり、状況に応じて最適な方法が異なる。
【０００７】
そこで本発明は、このような問題に鑑みてなされたもので、その目的とするところは、接触判定計算において更なる高速化を図ることの可能な、粒子の接触判定計算方法、粒子の接触判定計算装置、及びコンピュータプログラムを提供することにある。
【課題を解決するための手段】
【０００８】
上記課題を解決するために、本発明によれば、以下のような粒子の接触判定計算方法、粒子の接触判定計算装置、及びコンピュータプログラムが提供される。
【０００９】
本発明の粒子の接触判定計算方法は、３次元の解析対象領域において接触判定格子を用いて接触判定計算装置が行う粒子の接触判定計算方法であって、前記接触判定計算装置は、粒子の接触判定計算を行う接触判定計算部を備え、前記接触判定計算部は、格子番号が（２ｉ−１，２ｉ−１，２ｉ−１）番の接触判定格子における粒子の接触判定計算を並列に行うステップと、格子番号が（２ｉ，２ｉ−１，２ｉ−１）番の接触判定格子における粒子の接触判定計算を並列に行うステップと、格子番号が（２ｉ−１，２ｉ，２ｉ−１）番の接触判定格子における粒子の接触判定計算を並列に行うステップと、格子番号が（２ｉ，２ｉ，２ｉ−１）番の接触判定格子における粒子の接触判定計算を並列に行うステップと、格子番号が（２ｉ−１，２ｉ−１，２ｉ）番の接触判定格子における粒子の接触判定計算を並列に行うステップと、格子番号が（２ｉ，２ｉ−１，２ｉ）番の接触判定格子における粒子の接触判定計算を並列に行うステップと、格子番号が（２ｉ−１，２ｉ，２ｉ）番の接触判定格子における粒子の接触判定計算を並列に行うステップと、格子番号が（２ｉ，２ｉ，２ｉ）番の接触判定格子における粒子の接触判定計算を並列に行うステップと、を所定の順序で繰り返し行う（ｉは１以上の整数）ことを特徴とする。
【００１０】
本発明の他の粒子の接触判定計算方法は、１次元の解析対象領域において接触判定格子を用いて行う粒子の接触判定計算方法であって、格子番号が奇数番の接触判定格子における粒子の接触判定計算を並列に行うステップと、格子番号が偶数番の接触判定格子における粒子の接触判定計算を並列に行うステップと、を交互に行うことを特徴とする。
【００１１】
また、本発明の他の粒子の接触判定計算方法は、２次元の解析対象領域において接触判定格子を用いて行う粒子の接触判定計算方法であって、格子番号が（２ｉ−１，２ｉ−１）番の接触判定格子における粒子の接触判定計算を並列に行うステップと、格子番号が（２ｉ，２ｉ−１）番の接触判定格子における粒子の接触判定計算を並列に行うステップと、格子番号が（２ｉ−１，２ｉ）番の接触判定格子における粒子の接触判定計算を並列に行うステップと、格子番号が（２ｉ，２ｉ）番の接触判定格子における粒子の接触判定計算を並列に行うステップと、を所定の順序で繰り返し行う（ｉは１以上の整数）ことを特徴とする。
【００１２】
また、本発明の他の粒子の接触判定計算方法は、ｎ次元（ｎは１以上の整数）の解析対象領域において接触判定格子を用いて行う粒子の接触判定計算方法であって、前記接触判定格子を、２ｎ個の格子からなる単位格子にグループ分けし、前記各単位格子において、前記単位格子内の同じ位置の格子に同じ値を付与することで、前記単位格子内の２ｎ個の格子に１〜２ｎの値を付与し、同じ値の格子における粒子の接触判定計算を並列に行うステップを、所定の順序で繰り返し行うことを特徴とする。
【００１３】
また、本発明にかかる粒子の接触判定計算装置は、粒子の接触判定計算を行う接触判定計算部を備え、接触判定計算部は、上述した粒子の接触判定計算方法と同様の技術的特徴を有する。
【００１４】
また、本発明にかかるコンピュータプログラムは、コンピュータを、上記粒子の接触判定計算装置として機能させるためのコンピュータプログラムである。
【００１５】
本発明によれば、接触判定格子を、粒子が接触しないようにグループ化した上で並列計算を行うことができる。このため、接触する相手の情報も同時に更新する、いわゆる作用−反作用の法則を用いて演算量の削減を行うことができる。このような並列計算を行うことで、演算量はcritical指示文による方法と同じく１／２にしつつも、排他制御を必要としない。このようにして、接触判定計算において更なる高速化を図ることが可能である。
【００１６】
本発明の基本的な技術的特徴は、上述したように、接触判定格子を、粒子が接触しないようにグループ化した上で並列計算を行うことである。このことは、接触判定格子を複数のカラーに色分けすることとしてイメージすることができる。そこで、本明細書において、本発明の粒子の接触判定計算方法を「マルチカラー接触判定法」とも称する。また、このマルチカラー接触判定法で用いる接触判定格子のことを、「マルチカラー接触判定格子」とも称する。
【発明の効果】
【００１７】
本発明によれば、作用−反作用の法則を用いて演算量の削減を行いつつも、排他制御を必要としない。このようにして、接触判定計算において更なる高速化を図ることが可能である。本発明のその他の効果については、以下の発明を実施するための形態の項でも説明する。
【図面の簡単な説明】
【００１８】
【図１】接触判定格子のイメージを示す図である。
【図２】１次元のマルチカラー接触判定格子を概念的に示す図である。
【図３】マルチカラー接触判定法のプログラムを示す図である。
【図４】１次元のマルチカラー接触判定法の流れを示す図である。
【図５】２次元のマルチカラー接触判定格子を概念的に示す図である。
【図６】２次元のマルチカラー接触判定法の流れを示す図である。
【図７】３次元のマルチカラー接触判定格子を概念的に示す図である。
【図８】３次元のマルチカラー接触判定法の流れを示す図である。
【図９】本実施形態の方法（Method３）と従来法（Method１、Method２）との効果の比較を示す図である。
【発明を実施するための形態】
【００１９】
以下に添付図面を参照しながら、本発明の好適な実施の形態について詳細に説明する。さらに、従来の方法と本実施形態にかかるマルチカラー接触判定法とのスレッド並列実行性能を比較することで、マルチカラー接触判定法の有効性を検証する。なお、本明細書及び図面において、実質的に同一の機能構成を有する構成要素については、同一の符号を付することにより重複説明を省略する。
【００２０】
まず、本実施形態の前提となる接触判定計算について説明する。
【００２１】
はじめに、粒子を登録する。ここでは一例として、疑似乱数（以下、単に乱数ともいう）で発生させて粒子を登録する。次に粒子番号の並び替えを行う。本実施形態では、解析対象となる領域（例えば、３次元領域）に任意の数だけ粒子を発生させる。粒子の位置は、乱数で決められる。そして、各粒子について接触判定計算を行う。この接触判定計算には、接触する可能性のある粒子を判定するために、以下に述べる接触判定格子を用いる。
【００２２】
図１は、接触判定格子のイメージを示す図である。図１に示したように、解析対象領域を格子状に分割する。その後、格子１００内に入っているすべての粒子１１０を、対応する判定格子上に登録する。格子幅を粒子直径と同じかそれ以上にすることにより、接触する可能性のある粒子を同一格子内と隣接格子内に存在する粒子に限定することができるため、計算量を減らすことができる。なお、格子幅は、粒子直径以上であれば、必ずしも粒子直径と同じでなくてもよい。
【００２３】
接触する粒子に対する接触力を更新する際、自分自身の情報だけでなく、相手の情報も同時に更新する。このとき、右隣り（又は上隣りなど）の格子のみの粒子と接触判定を行う。これは、作用−反作用の法則を用いて演算量を１／２に削減するための処理である。
【００２４】
この処理について、ＯｐｅｎＭＰ並列化することを考える。このとき、例えば粒子番号０番と粒子番号１番の粒子が、同一格子内にあるとする。また、粒子番号０番の粒子はスレッド０、粒子番号１番の粒子はスレッド１で処理されるとする。このとき、スレッド０において更新を行っている最中に、スレッド１において更新を行う状況がありうる。すなわち、スレッド０とスレッド１の間で、同時に値を更新しようとするタイミングがありうる。この状態では、スレッド０もしくはスレッド１の結果のみが反映され、逐次の結果と一致しなくなる。したがって、この処理は、排他制御が必要である。以下、衝突判定計算において、スレッド間で値を同時に書き込むことを「競合」と呼ぶことにする。
【００２５】
従来、競合を回避するには、ＯｐｅｎＭＰで提供されるcritical指示子を用いたディレクティブであるcritical指示文やatomic文により、排他制御を行うことで競合の問題を制御してきた（例えば、非特許文献１を参照）。また、critical指示文を使わない方法として、接触時に相手の情報を更新しないようにすることで、競合を回避してきた。ただしこの方法では、接触力計算に関して冗長計算を行うことから、演算量が２倍に増加する。なお、本明細書において、作用−反作用の原理を利用せずに接触時に相手の情報を更新しない接触力計算のことを、「冗長計算」と称する。
【００２６】
以上のような問題を解消する方法として、以下、本実施形態にかかるマルチカラー接触判定法について、詳細に説明する。
【００２７】
マルチカラー接触判定法は、マルチカラー接触判定格子を用いた接触判定計算の方法である。マルチカラー接触判定格子とは、従来の接触判定計算に用いる接触判定格子を用いるが、接触判定計算のためのアクセスの仕方が異なる接触判定格子のことである。以下ではまず、１次元の場合における例を示すことで、マルチカラー接触判定法のアルゴリズムを説明する。
【００２８】
（１）１次元の解析対象領域
図２は、１次元のマルチカラー接触判定格子を概念的に示す図である。図２（ａ）は、接触判定格子の格子番号を示し、図２（ｂ）は、マルチカラー接触判定格子を示し、図２（ｃ）は、マルチカラーを付すための単位格子を示す。なお、本実施形態のマルチカラー接触判定法の適用の前提は、接触格子の格子幅を粒子直径、あるいはそれ以上に取ることである。また、接触判定格子を用いた接触判定計算には、格子番号ごとに並列計算する方法（格子並列）と、粒子番号ごとに計算する方法（粒子並列）が知られているが、本実施形態では、格子並列を前提とする。すなわち、一つの接触格子内に存在する粒子の接触判定計算は同一のスレッドで行う。
【００２９】
ＯｐｅｎＭＰ並列化を考慮すると、接触判定格子の性質上、対象格子番号ｘの両隣格子である、格子番号（ｘ−１）と格子番号（ｘ＋１）には依存関係があるが、それ以外の格子には依存関係がない。１次元の場合にこれを一般化すると、格子番号が奇数番の接触判定格子どうしには依存関係がなく、また、格子番号が偶数番の接触判定格子どうしには依存関係がない。ここで、依存関係がないとは、データの競合が起き得ないために並列実行が可能であることをいい、依存関係があるとは、データ競合が起き得るために並列実行ができないことをいう。
【００３０】
そこで、図２（ｂ）に示したマルチカラー接触判定格子を採用する。図２（ｂ）に示したように、並列実行が可能である奇数番の格子番号ｘ(ｘ＝１、３、５、７、９、・・・)に対して、値“１”を付与する。また、並列実行が可能である偶数番の格子番号ｘ（ｘ＝２、４、６、８、・・・）に対して、値“２”を付与する。この値“１”と値“２”は、接触判定格子を２色（２１色）のカラーで色分けすることを意味する。このようにしてマルチカラー接触判定格子を定義することができる。なお、値“１”と値“２”は、単に両者を区別するものに過ぎず、番号自体に特別の意味はない。奇数番の格子に対して値“２”を付与し、偶数番の格子に対して値“１”を付与してもよいことは言うまでもない。
【００３１】
接触判定格子への値の付与の仕方として、次のように考えることもできる。まず、図２（ａ）の接触判定格子を２個（２１個）ずつの単位格子（１次元単位格子）にグループ分けする。図２（ｃ）は、１次元単位格子を示す（ｉは１以上の整数）。そして、各１次元単位格子に対して、単位格子内の同じ位置の格子に同じ値を付与する。図２（ｃ）の一例では、１次元単位格子の左位置（格子番号２ｉ−１）の格子に対して値“１”が付与されている。また、右位置（格子番号２ｉ）の格子に対して値“２”が付与されている。これをすべての１次元単位格子に対して行う。このようにして、図２（ｂ）に示したマルチカラー接触判定格子を定義することができる。
【００３２】
以上のように番号付けされたマルチカラー接触判定格子を用いた、粒子の接触判定計算方法について、図３及び図４を参照しながら説明する。図３は、マルチカラー接触判定法のプログラムを示す図である。ここでは１次元であるため、ｎ＝１であり、ＭＡＸ＿ＣＯＬＯＲ＝２１＝２である。図４に示したように、まず、値“１”が付与された格子の接触判定計算を並列に行う（ステップＳ１０２）。その後、値“２”が付与された格子の接触判定を行う（ステップＳ１０４）。この接触判定を必要に応じて繰り返し行う。例えば、個々の粒子の位置と速度を小さな時間刻み巾（ΔＴ）でもって、少しずつ更新するような粉体解析の場合には、各時間ステップにおいて、粒子の接触判定を行い、接触している粒子に対しては接触力を求めるようにしてもよい。
【００３３】
上述したように、本実施形態のマルチカラー接触判定法の適用の前提は、接触格子の格子幅を粒子直径、あるいはそれ以上に取ることである。そうすることで、同じ値（カラー）の格子内にある粒子どうしは接触することがない。また、本実施形態では、格子並列を前提とし、一つの接触格子内に存在する粒子の接触判定計算は同一のスレッドで行う。そのためマルチカラー接触判定法では、自分と相手の接触力を更新する際に、原理上、競合が起きない。その結果、並列化をしても複数のスレッド間でデータの競合が起きない。なお、マルチカラー接触判定法では、逐次の接触判定計算とは加減算の順序が異なる。ただし、数学上における演算結果の一致は保証される。以上のようにして、マルチカラー接触判定法では、並列実行（計算）が可能である。
【００３４】
以上は１次元の解析対象に対するマルチカラー接触判定法の説明である。２次元、３次元に解析対象を広げる場合においても、マルチカラー接触判定法は、マルチカラー接触判定格子を２次元、３次元に拡張することで拡張が可能となる。以下、順に説明する。
【００３５】
（２）２次元の解析対象領域
まず、２次元の解析対象領域について、図５、図６を参照しながら説明する。図５は、２次元のマルチカラー接触判定格子を概念的に示す図である。図５（ａ）は、２次元のマルチカラー接触判定格子を示し、図５（ｂ）は、マルチカラーを付すための単位格子を示す。図６は、２次元のマルチカラー接触判定法の流れを示す図である。
【００３６】
マルチカラー接触判定法を２次元に拡張するには、図５（ａ）に示したように、２次元にマルチカラー接触判定格子に拡張する。図５（ａ）に示したように、接触判定格子上のデータ依存を回避するため、４色（２２色）に色分けする。このことで、１次元と同様の議論で、マルチカラー接触判定法が実現できる。すなわち、図５（ａ）に示したように、マルチカラー接触判定格子においては、値“１”が付与された格子どうしには依然関係がない。同様に、値“２”が付与された格子どうしには依存関係がなく、値“３”が付与された格子どうしには依存関係がなく、値“４”が付与された格子どうしには依存関係がない。このため、上記１次元と同様の議論で並列実行が可能である。
【００３７】
接触判定格子への値の付与の仕方について説明する。まず、図５（ａ）を２×２個ずつの単位格子（２次元単位格子）に分割する。図５（ｂ）は、２次元単位格子の一つを示す（ｉは１以上の整数）。そして、各２次元単位格子に対して、単位格子内の同じ位置の格子に同じ値を付与する。図５（ｂ）の一例では、格子番号が（２ｉ−１，２ｉ−１）番の格子に値“１”が付与されている。また、格子番号が（２ｉ，２ｉ−１）番の格子に値“２”が付与されている。また、格子番号が（２ｉ−１，２ｉ）番の格子に値“３”が付与されている。また、格子番号が（２ｉ，２ｉ）番の格子に値“４”が付与されている。これをすべての２次元単位格子に対して行う。このようにして、図５（ａ）に示したマルチカラー接触判定格子を定義することができる。
【００３８】
以上のように番号付けされたマルチカラー接触判定格子を用いた、粒子の接触判定計算方法について、図６を参照しながら説明する。まず、格子番号が（２ｉ−１，２ｉ−１）番の格子、すなわち、値“１”が付与された格子の接触判定計算を並列に行う（ステップＳ２０２）。次に、格子番号が（２ｉ，２ｉ−１）番の格子、すなわち、値“２”が付与された格子の接触判定を行う（ステップＳ２０４）。次に、格子番号が（２ｉ−１，２ｉ）番の格子、すなわち、値“３”が付与された格子の接触判定を行う（ステップＳ２０６）。次に、格子番号が（２ｉ，２ｉ）番の格子、すなわち、値“４”が付与された格子の接触判定を行う（ステップＳ２０８）。当該処理は、図３に示したマルチカラー接触判定法のプログラムにおいて、ｎ＝２、ＭＡＸ＿ＣＯＬＯＲ＝２２＝４として実行することができる。
【００３９】
（３）３次元の解析対象領域
次に、３次元の解析対象領域について、図７、図８を参照しながら説明する。図７は、３次元のマルチカラー接触判定格子を概念的に示す図である。図７（ａ）は、３次元のマルチカラー接触判定格子を示し、図７（ｂ）は、マルチカラーを付すための単位格子を示す。図８は、３次元のマルチカラー接触判定法の流れを示す図である。
【００４０】
マルチカラー接触判定法を３次元に拡張するには、図７（ａ）に示したように、３次元にマルチカラー接触判定格子を拡張する。図７（ａ）に示したように、接触判定格子上のデータ依存を回避するため、８色（２３色）色に色分けする。このことで、１次元、２次元と同様の議論で、マルチカラー接触判定法が実現できる。すなわち、図７（ａ）に示したように、マルチカラー接触判定格子においては、値“１”が付与された格子どうしには依然関係がない。同様に、値“２”が付与された格子どうしには依存関係がなく、値“３”が付与された格子どうしには依存関係がなく、値“４”が付与された格子どうしには依存関係がなく、値“５”が付与された格子どうしには依存関係がなく、値“６”が付与された格子どうしには依存関係がなく、値“７”が付与された格子どうしには依存関係がなく、値“８”が付与された格子どうしには依存関係がない。このため、上記１次元、２次元と同様の議論で並列実行が可能である。
【００４１】
接触判定格子への値の付与の仕方について説明する。まず、図７（ａ）を２×２×２個ずつの単位格子（３次元単位格子）に分割する。図７（ｂ）は、３次元単位格子の一つを示す（ｉは１以上の整数）。そして、各３次元単位格子に対して、単位格子内の同じ位置の格子に同じ値を付与する。
【００４２】
図７（ｂ）の一例では、格子番号が（２ｉ−１，２ｉ−１，２ｉ−１）番の格子に値“１”が付与されている。また、格子番号が（２ｉ，２ｉ−１，２ｉ−１）番の格子に値“２”が付与されている。また、格子番号が（２ｉ−１，２ｉ，２ｉ−１）番の格子に値“３”が付与されている。また、格子番号が（２ｉ，２ｉ，２ｉ−１）番の格子に値“４”が付与されている。また、格子番号が（２ｉ−１，２ｉ−１，２ｉ）番の格子に値“５”が付与されている。また、格子番号が（２ｉ，２ｉ−１，２ｉ）番の格子に値“６”が付与されている。また、格子番号が（２ｉ−１，２ｉ，２ｉ）番の格子に値“７”が付与されている。また、格子番号が（２ｉ，２ｉ，２ｉ）番の格子に値“８”が付与されている。これをすべての３次元単位格子に対して行う。このようにして、図７（ａ）に示したマルチカラー接触判定格子を定義することができる。
【００４３】
以上のように番号付けされたマルチカラー接触判定格子を用いた、粒子の接触判定計算方法について、図８を参照しながら説明する。まず、格子番号が（２ｉ−１，２ｉ−１，２ｉ−１）番の格子、すなわち、値“１”が付与された格子の接触判定計算を並列に行う（ステップＳ３０２）。次に、格子番号が（２ｉ，２ｉ−１，２ｉ−１）番の格子、すなわち、値“２”が付与された格子の接触判定を行う（ステップＳ３０４）。次に、格子番号が（２ｉ−１，２ｉ，２ｉ−１）番の格子、すなわち、値“３”が付与された格子の接触判定を行う（ステップＳ３０６）。次に、格子番号が（２ｉ，２ｉ，２ｉ−１）番の格子、すなわち、値“４”が付与された格子の接触判定を行う（ステップＳ３０８）。
【００４４】
以下同様に、格子番号が（２ｉ−１，２ｉ−１，２ｉ）番の格子、すなわち、値“５”が付与された格子の接触判定計算を並列に行う（ステップＳ３１０）。次に、格子番号が（２ｉ，２ｉ−１，２ｉ）番の格子、すなわち、値“６”が付与された格子の接触判定を行う（ステップＳ３１２）。次に、格子番号が（２ｉ−１，２ｉ，２ｉ）番の格子、すなわち、値“７”が付与された格子の接触判定を行う（ステップＳ３１４）。次に、格子番号が（２ｉ，２ｉ，２ｉ）番の格子、すなわち、値“８”が付与された格子の接触判定を行う（ステップＳ３１６）。当該処理は、図３に示したマルチカラー接触判定法のプログラムにおいて、ｎ＝３、ＭＡＸ＿ＣＯＬＯＲ＝２３＝８として実行することができる。
【００４５】
以上、本実施形態にかかる粒子の接触判定計算方法について説明した。本実施形態の方法は、専用の装置（粒子の接触判定計算装置）によって実行してもよい。かかる粒子の接触判定計算装置は、粒子の接触判定計算部を備えることができる。そして、かかる粒子の接触判定計算部は、本実施形態で説明した粒子の接触判定計算方法を実行することができる。
【００４６】
また、専用の装置ではなく、汎用のコンピュータに所定のコンピュータプログラムを組み込むことによって、コンピュータを、上述したような粒子の接触判定計算装置として機能させてもよい。かかるコンピュータプログラムは、ＣＲ−ＲＯＭ、ＤＶＤ−ＲＯＭ、半導体メディアなどの各種記憶媒体に格納された形で市場を流通させることができる。また、様々なネットワークを介した伝達、例えばインターネットを介したダウンロードによって市場を流通させるようにしてもよい。
【００４７】
以上、本実施形態に特徴的なマルチカラー接触判定法について説明した。次いで、性能評価について説明する。
【００４８】
本性能評価では、以下の仕様により評価を行った。
・インテル(R)  Core(TM)  i7-2600プロセッサー(クアッドコア/3.40GHz/TB時最大3.80GHz/8MBスマートキャッシュ/HT対応)
・8GBメモリ[4GB×2 (DDR3  SDRAM  PC3-10600)/デュアルチャネル]
・Windows(R)7  Professional正規版(OEM)［64bit］
【００４９】
本性能評価では、粒子径は０．００２、粒子数は６２５０万、ステップ数は２０に設定した。また判定時に各粒子は接触数の情報を更新するが、３次元の接触判定格子を利用するため、１粒子あたりの接触数情報の上限は１２とする。
【００５０】
本実施形態にかかるマルチカラー接触判定法の有効性を見るために、以下の表１のように３つのプログラムを用意した。表１は、Methodによる実装方式の違いを示す。
【００５１】
【表１】
0005874155.tif 000002
【００５２】
表１では、接触判定格子並列化時のデータの競合を防ぐためにcritical指示文を入れたプログラムをmethod１とする。接触判定格子並列において、相手の情報を書き込まず、計算量が２倍になるプログラムをmethod２とする。以上のmethod１、method２は従来法である。最後に、本実施形態の方法であるマルチカラー接触判定格子を用いて並列化するプログラムをmethod３とする。
【００５３】
１ノードを占有し、ＯｐｅｎＭＰにより最大８スレッドで並列実行したときの計算時間と台数効果について性能を評価する。ここで台数効果とは、並列処理の効果を示す一般的な指標である。式（１）のように「逐次実行における計算時間」に対する「並列実行における実行時間」の比で表される。
【００５４】
Ｓ＝逐次実行に要する計算時間／並列計算に要する計算時間・・・式（１）
【００５５】
なお、式（１）における計算時間とは、接触判定計算に要する計算時間と粒子の接触力計算時間を合計した総計算時間のことを指す。
【００５６】
各methodにおいて、逐次計算（スレッド数１）と並列計算（スレッド数２、４、６、及び８）での計算時間の比較を行う。ここで計算時間とは、５ステップ目以降の計算時間を測定し、その値を平均したものである。計算時間の結果を、以下の表２及び図９に示す。表２は、各methodの実行時間（秒）を示す。
【００５７】
【表２】
0005874155.tif 000003
【００５８】
表２から、８スレッド時において、本実施形態にかかるマルチカラー接触判定法を利用することで最大で１．６２倍（method１の５．６５６[秒]に対する、method３の３．４９６[秒]）の速度向上が得られることがわかる。
【００５９】
以上、本実施形態では、粉体解析等で用いられる接触判定計算において、マルチカラー接触判定格子を用いた新しいスレッド並列化手法である、マルチカラー接触判定法について説明した。また、いくつかの従来法に対する性能評価を行った。従来法では、格子並列を採用しているが、計算の効率化のために、粒子番号のリナンバリングを行っている。本実施形態のマルチカラー接触判定法では、このリナンバリングに加え、格子のアクセスの仕方を工夫することで、スレッド並列化時にデータ競合が起こらないようにしている。そのため、従来法に対しより高い並列化効率が期待できる。そして、上述したように、本実施形態のマルチカラー接触判定法は、計算時間と台数効果の両方において、従来手法より優れた結果が得られた。
【００６０】
以上、添付図面を参照しながら本発明の好適な実施形態について説明したが、本発明はかかる例に限定されないことは言うまでもない。当業者であれば、特許請求の範囲に記載された範疇内において、各種の変更例または修正例に想到し得ることは明らかであり、それらについても当然に本発明の技術的範囲に属するものと了解される。
【００６１】
例えば、上記実施形態では、解析対象空間として、１次元〜３次元空間について説明したが、本発明はこれに限定されず、任意の次元数の解析対象空間に適用可能である。
【産業上の利用可能性】
【００６２】
本発明は、粒子の接触判定計算方法、粒子の接触判定計算装置、及びコンピュータプログラムに利用可能であり、特に、所定の解析対象領域において接触判定格子を用いて行う、粒子の接触判定計算方法、粒子の接触判定計算装置、及びコンピュータプログラムに利用可能である。
【符号の説明】
【００６３】
１０  接触判定格子
２０  粒子
１００  １次元接触判定格子
１１０  マルチカラー接触判定格子（１次元）
１２０  １次元単位格子
２１０  マルチカラー接触判定格子（２次元）
２２０  ２次元単位格子
３１０  マルチカラー接触判定格子（３次元）
３２０  ３次元単位格子
【図１】
0005874155.tif 000004
【図２】
0005874155.tif 000005
【図３】
0005874155.tif 000006
【図４】
0005874155.tif 000007
【図５】
0005874155.tif 000008
【図６】
0005874155.tif 000009
【図７】
0005874155.tif 000010
【図８】
0005874155.tif 000011
【図９】
0005874155.tif 000012