特開 2023-181746 シミュレーション方法、シミュレーション装置、及びプログラム

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】公開特許公報(A)

(11)【公開番号】P2023181746

(43)【公開日】2023-12-25

(54)【発明の名称】シミュレーション方法、シミュレーション装置、及びプログラム

(51)【国際特許分類】

   G16C 10/00 20190101AFI20231218BHJP

【ＦＩ】

G16C10/00

【審査請求】未請求

【請求項の数】5

【出願形態】ＯＬ

(21)【出願番号】P 2022095057

(22)【出願日】2022-06-13

(71)【出願人】

【識別番号】000002107

【氏名又は名称】住友重機械工業株式会社

(74)【代理人】

【識別番号】100105887

【弁理士】

【氏名又は名称】来山 幹雄

(72)【発明者】

【氏名】松村 裕也

(57)【要約】

【課題】複数のスレッドを並列実行する機能を有効に利用することが可能なシミュレーション方法を提供する。
【解決手段】計算領域に複数の粒子を配置し、各粒子に作用する力に基づいて運動方程式を数値的に解き、粒子の運動を解析する。この解析は、計算領域に複数の粒子を挿入する処理を含む。計算領域に複数の粒子を挿入する処理は、複数のスレッドを並列実行する機能を持つ処理部の複数のスレッドのそれぞれが、挿入すべき複数の粒子のうち少なくとも一部の粒子を担当して、計算領域に複数の粒子のそれぞれの挿入位置を決定する処理を実行する。
【選択図】図４

【特許請求の範囲】

【請求項1】

計算領域に複数の粒子を配置し、各粒子に作用する力に基づいて運動方程式を数値的に解き、粒子の運動を解析するシミュレーション方法において、
前記計算領域に複数の粒子を挿入する処理を含み、
前記計算領域に複数の粒子を挿入する処理は、複数のスレッドを並列実行する機能を持つ処理部の前記複数のスレッドのそれぞれが、挿入すべき複数の粒子のうち少なくとも一部の粒子を担当して、前記計算領域に複数の粒子のそれぞれの挿入位置を決定する処理を実行するシミュレーション方法。

【請求項2】

前記複数の粒子を挿入する処理において、
前記複数のスレッドのそれぞれが、
前記計算領域に新たに挿入する少なくとも１つの粒子の挿入位置の候補となる候補位置を決定し、
自スレッドが決定した前記候補位置と、他のスレッドが決定した前記候補位置との間で、粒子に課される位置条件が満たされるまで前記候補位置を修正する、
その後、自スレッドが決定した前記候補位置と、前記計算領域に存在する既存の粒子の位置との間で、前記位置条件が満たされるように、前記候補位置を修正する請求項１に記載のシミュレーション方法。

【請求項3】

シミュレーション条件が入力される入力部と、
複数のスレッドを並列実行する機能を持つ処理部と、
前記処理部が実行するプログラムが記憶された記憶部と
を備え、
前記処理部は、
前記入力部に入力されたシミュレーション条件に基づいて、計算領域に複数の粒子を配置し、各粒子に作用する力に基づいて運動方程式を数値的に解くことによって粒子の運動を解析する処理と、
前記計算領域に複数の粒子を挿入する処理と
を実行し、
前記計算領域に複数の粒子を挿入する処理において、前記複数のスレッドのそれぞれが、挿入すべき複数の粒子のうちの少なくとも１つを担当して、前記複数の粒子のそれぞれの挿入位置を決定するシミュレーション装置。

【請求項4】

前記計算領域に複数の粒子を挿入する処理において、
前記処理部で並列実行される前記複数のスレッドのそれぞれは、
前記複数の粒子のそれぞれの挿入位置の候補となる候補位置を決定し、
自スレッドが決定した前記候補位置と、他のスレッドが決定した前記候補位置との間で、粒子に課される位置条件が満たされるまで、前記候補位置を修正し、
その後、前記候補位置と、前記計算領域に存在する既存の粒子の位置との間で、前記位置条件が満たされるまで、前記候補位置を修正する請求項３に記載のシミュレーション装置。

【請求項5】

複数のスレッドを並列実行する機能を持つ処理部に、
シミュレーション条件に基づいて、計算領域に複数の粒子を配置し、各粒子に作用する力に基づいて運動方程式を数値的に解くことによって粒子の運動を解析する処理と、
前記計算領域に複数の粒子を挿入する処理と
を実行させるプログラムであって、
前記計算領域に複数の粒子を挿入する処理において、前記複数のスレッドのそれぞれが、挿入すべき複数の粒子のうちの少なくとも１つを担当して、前記複数の粒子のそれぞれの挿入位置を決定するプログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、シミュレーション方法、シミュレーション装置、及びプログラムに関する。

【背景技術】

【0002】

分子動力学法を用いたシミュレーションにおいて、計算領域に粒子を挿入することによって圧力や流量の制御が行われる。この処理を並列で実行する場合、一般的には、計算領域を空間的に複数の領域に分割し、並列実行される複数のスレッドのそれぞれが、分割された複数の領域のそれぞれにおける粒子の挿入処理を実行する。例えば、下記の特許文献１には、計算対象の中に複数の小さな部分領域を設け、部分領域で分子動力学法による詳細な解析を行うシミュレーション方法が開示されている。

【先行技術文献】

【特許文献】

【0003】

【特許文献1】特開平１０－９７５１８号公報

【発明の概要】

【発明が解決しようとする課題】

【0004】

分子動力学法を用いて計算する場合には、計算対象となる領域にある程度の大きさが必要であるため、計算領域を複数の領域に分割するときの分割数には限界がある。グラフィックス・プロセシング・ユニット（ＧＰＵ）のように、多数のスレッドを並列実行することができるプロセッサを用いる場合、多数のスレッドのそれぞれに１つの分割領域を割り振って計算を行うと、分割領域が割り振られない多くのスレッドが発生する。このため、プロセッサの計算能力を十分生かすことができない。

【0005】

本発明の目的は、複数のスレッドを並列実行する機能を有効に利用することが可能なシミュレーション方法、シミュレーション装置、及びプログラムを提供することである。

【課題を解決するための手段】

【0006】

本発明の一観点によると、
計算領域に複数の粒子を配置し、各粒子に作用する力に基づいて運動方程式を数値的に解き、粒子の運動を解析するシミュレーション方法において、
前記計算領域に複数の粒子を挿入する処理を含み、
前記計算領域に複数の粒子を挿入する処理は、複数のスレッドを並列実行する機能を持つ処理部の前記複数のスレッドのそれぞれが、挿入すべき複数の粒子のうち少なくとも一部の粒子を担当して、前記計算領域に複数の粒子のそれぞれの挿入位置を決定する処理を実行するシミュレーション方法が提供される。

【0007】

本発明の他の観点によると、
シミュレーション条件が入力される入力部と、
複数のスレッドを並列実行する機能を持つ処理部と、
前記処理部が実行するプログラムが記憶された記憶部と
を備え、
前記処理部は、
前記入力部に入力されたシミュレーション条件に基づいて、計算領域に複数の粒子を配置し、各粒子に作用する力に基づいて運動方程式を数値的に解くことによって粒子の運動を解析する処理と、
前記計算領域に複数の粒子を挿入する処理と
を実行し、
前記計算領域に複数の粒子を挿入する処理において、前記複数のスレッドのそれぞれが、挿入すべき複数の粒子のうちの少なくとも１つを担当して、前記複数の粒子のそれぞれの挿入位置を決定するシミュレーション装置が提供される。

【0008】

本発明のさらに他の観点によると、
複数のスレッドを並列実行する機能を持つ処理部に、
シミュレーション条件に基づいて、計算領域に複数の粒子を配置し、各粒子に作用する力に基づいて運動方程式を数値的に解くことによって粒子の運動を解析する処理と、
前記計算領域に複数の粒子を挿入する処理と
を実行させるプログラムであって、
前記計算領域に複数の粒子を挿入する処理において、前記複数のスレッドのそれぞれが、挿入すべき複数の粒子のうちの少なくとも１つを担当して、前記複数の粒子のそれぞれの挿入位置を決定するプログラムが提供される。

【発明の効果】

【0009】

複数のスレッドのそれぞれが、挿入すべき複数の粒子のうち少なくとも一部の粒子を担当して粒子の挿入位置を決定することにより、複数のスレッドを有効に利用することが可能である。

【図面の簡単な説明】

【0010】

【図1】図１は、実施例によるシミュレーション装置のブロック図である。

【図2】図２は、実施例によるシミュレーション方法の手順を示すフローチャートである。

【図3】図３は、計算領域、既に存在する既存の複数の粒子、及び新たに挿入する複数の粒子を示す模式図である。

【図4】図４は、複数の粒子を挿入する手順を示すフローチャートである。

【図5】図５Ａは、複数のスレッドで決定された候補位置の一例を示す模式図であり、図５Ｂは、候補位置を修正した後の候補位置の一例を示す模式図である。

【図6】図６Ａは、修正前の１つのスレッドの候補位置及び既存の粒子の模式図であり、図６Ｂは、修正後の候補位置及び既存の粒子の模式図である。

【図7】図７は、計算領域内の粒子数の時間変化を示すグラフである。

【図8】図８は、計算開始時点における粒子の分布を示す模式図である。

【図9】図９は、１４００タイムステップ目における粒子１１の分布を示す模式図である。

【発明を実施するための形態】

【0011】

図１～図９を参照して、一実施例によるシミュレーション方法、及びシミュレーション装置について説明する。

【0012】

図１は、本実施例によるシミュレーション装置のブロック図である。実施例によるシミュレーション装置は、入力部５０、処理部５１、出力部５２、及び記憶部５３を含む。入力部５０から処理部５１にシミュレーション条件等が入力される。さらに、オペレータから入力部５０に各種指令（コマンド）等が入力される。入力部５０は、例えば通信装置、リムーバブルメディア読取装置、キーボード等で構成される。

【0013】

処理部５１は、入力されたシミュレーション条件及び指令に基づいて分子動力学法またはくりこみ群分子動力学法（以下、単に分子動力学法という。）を用いたシミュレーションを行う。処理部５１は、複数のスレッドを並列実行する機能を持つ。処理部５１として、例えば、グラフィックス・プロセシング・ユニット（ＧＰＵ）が用いられる。

【0014】

処理部５１によるシミュレーションの結果が出力部５２に出力される。シミュレーション結果には、例えば、シミュレーション対象物を再現した解析モデルを構成する複数の粒子の状態、対象物の物理量の空間的変化、時間的変化等を表す情報が含まれる。分子動力学法によるシミュレーションを処理部５１に実行させるためのプログラムが、記憶部５３に記憶されている。出力部５２は、通信装置、リムーバブルメディア書込み装置、ディスプレイ等を含む。

【0015】

図２は、本実施例によるシミュレーション方法の手順を示すフローチャートである。まず、処理部５１（図１）は、入力部５０（図１）に入力されたシミュレーション条件を取得する（ステップＳＡ１）。シミュレーション条件には、シミュレーション対象物の形状及び大きさ、物性値、粒子間相互作用を定義する情報、外力を定義する情報、初期条件、境界条件、時間刻み幅、解析終了条件等が含まれる。粒子間相互作用は、例えばレナードジョーンズポテンシャルを用いて定義される。

【0016】

次に、処理部５１は、シミュレーション条件に基づいて計算の初期化を行う（ステップＳＡ２）。具体的には、シミュレーション対象物の形状及び大きさに応じて計算領域を設定し、初期条件に応じて計算領域内に複数の粒子を配置する。さらに、各粒子に初期速度を付与する。

【0017】

次に、処理部５１は、シミュレーション条件に基づいて、各粒子に作用する力を計算する（ステップＳＡ３）。次に、処理部５１は、各粒子について運動方程式を解き、各粒子の位置及び速度を更新する（ステップＳＡ４）。処理を終了するか否かを判定し（ステップＳＡ５）、処理を終了する場合には、解析結果を出力部５２（図１）に出力する（ステップＳＡ７）。

【0018】

処理を継続する場合には、所定の条件に基づいて、計算領域に複数の粒子を挿入する（ステップＳＡ６）。例えば、シミュレーション対象物の圧力を制御する場合には、所定の圧力に維持されるように粒子を挿入する。シミュレーション対象物の流量を制御する場合には、所定の流量が確保されるように粒子を挿入する。

【0019】

図３は、計算領域１０、既に存在する既存の複数の粒子１１、及び新たに挿入する複数の粒子１２を示す模式図である。図３では、既存の複数の粒子１１を、黒色で塗りつぶした円形記号で表し、新たに追加する粒子１２を中空の円形記号で表している。処理部５１（図１）の複数のスレッドのそれぞれに、１つの粒子１２の挿入位置を決定する処理を割り振る。なお、本実施例では、処理部５１のスレッドの数が、挿入すべき粒子１２の個数以上であると仮定する。一部のスレッドには、粒子１２を挿入する処理が割り振られない場合もある。粒子を挿入する手順については、後に図４～図６Ｂを参照して説明する。計算領域１０に粒子１２を挿入した後、処理部５１はステップＳＡ３からの手順を繰り返す。

【0020】

図４は、複数の粒子を挿入する手順（ステップＳＡ６）を示すフローチャートである。処理部５１（図１）は、複数のスレッドのそれぞれに、挿入する粒子一つを割り振る。複数のスレッドが並列実行する手順は、複数のスレッドで同一であるため、図４では、複数のスレッドのうち１つのスレッドに着目してその手順を示し、他のスレッドが実行する手順は省略して、一連の手順を１本の破線で示している。

【0021】

まず、複数のスレッドのそれぞれは、粒子を挿入すべき位置である候補位置を決定する（ステップＳＢ０１）。例えば、スレッドのそれぞれは、計算領域１０（図３）全体を対象として、計算領域１０の中からランダムに候補位置を選択する。次に、すべてのスレッドの間で、候補位置を共有する（ステップＳＢ０２）。スレッドのそれぞれは、自スレッドが決定した候補位置と、他スレッドが決定した候補位置とが、所定の位置条件を満たすか否かを判定する（ステップＳＢ０３）。

【0022】

例えば、自スレッドの候補位置から他のいずれかのスレッドの候補位置までの距離が、所定の最低許容距離未満である場合には、位置条件が満たされないと判定する。自スレッドの候補位置から他のいずれのスレッドの候補位置までの距離も、最低許容距離以上である場合には、位置条件が満たされると判定する。所定の位置条件が満たされる場合には、すべてのスレッドの候補位置が位置条件を満たすまで待機する（ステップＳＢ０４）。

【0023】

次に、最低許容距離について説明する。２つの粒子の間の距離が近すぎると、２つの粒子に作用する斥力が大きくなりすぎて、計算が破綻してしまう。最低許容距離は、計算が破綻してしまう距離よりもやや長く設定するとよい。例えば、最低許容距離を、粒子間の相互作用ポテンシャルから定義される粒子の直径よりもやや長い値にするとよい。

【0024】

現在の候補位置が位置条件を満たしていない場合は、処理をステップＳＢ０１に戻し、候補位置を修正する。この際にも、例えば候補位置を計算領域１０（図３）の範囲内でランダムに選択する。

【0025】

図５Ａは、複数のスレッドで決定された候補位置１３の一例を示す模式図である。ｉ番目のスレッドの候補位置１３ｉとｊ番目のスレッドの候補位置１３ｊとの間の距離ｒ_ｉｊは、最低許容距離以上であり、候補位置１３ｉと候補位置１３ｊとは、位置条件を満たす。ｉ番目のスレッドの候補位置１３ｉとｋ番目のスレッドの候補位置１３ｋとの間の距離ｒ_ｉｋは、最低許容距離未満であり、候補位置１３ｉと候補位置１３ｋとは、位置条件を満たさない。この場合、ｉ番目のスレッドとｋ番目のスレッドとが、それぞれ候補位置１３ｉ、１３ｋを修正する（ステップＳＢ０１）。

【0026】

図５Ｂは、候補位置１３を修正した後の候補位置の一例を示す模式図である。破線で示した候補位置１３ｉ、１３ｋが除去され、実線で示された位置に新たに候補位置１３ｉ、１３ｋが決定される。

【0027】

図４に示すように、すべてのスレッドの候補位置が位置条件を満たしたら、スレッドのそれぞれは、自スレッドの候補位置と、既存の粒子１１（図３）の位置とが、位置条件を満たすか否かを判定する（ステップＳＢ０５）。位置条件が満たされる場合は、現在の候補位置を、粒子の挿入位置として確定する（ステップＳＢ０８）。

【0028】

自スレッドの候補位置と、既存の粒子１１（図３）の位置とが、位置条件を満たさない場合は、候補位置を修正する（ステップＳＢ０６）。以下、図６Ａ及び図６Ｂを参照して、ステップＳＢ０６で候補位置を修正する処理について説明する。図６Ａは、修正前のｉ番目のスレッドの候補位置１３ｉ及び既存の粒子１１の模式図であり、図６Ｂは、修正後の候補位置１３ｉ及び既存の粒子１１の模式図である。図６Ａ及び図６Ｂにおいて、既存の粒子１１を黒く塗りつぶした円形記号で表し、候補位置１３ｉを中空の円形記号で表す。

【0029】

図６Ａに示すように、候補位置１３ｉから既存の粒子１１ｊの位置までの距離ｒ_ｉｊが最低許容距離未満である。このため、候補位置１３ｉと、既存の粒子１１ｊの位置とは、位置条件を満たしていない。このとき、ｉ番目のスレッドは、既存の粒子１１ｊによる相互作用ポテンシャルから、候補位置１３ｉに位置する粒子に作用する力を計算する。候補位置１３ｉと既存の粒子１１ｊの位置とが、位置条件を満たしていない場合には、２つの粒子に斥力が作用する。

【0030】

図６Ｂに示すように、位置条件を満たしていない候補位置１３ｉを、位置条件を満たす方向に微小距離だけ移動させる。移動させる方向は、例えば、相手側の既存の粒子１１ｊからの斥力の方向とすればよい。移動させる距離は、予め設定しておくとよい。例えば、粒子直径の１倍～３倍程度に設定しておくとよい。

【0031】

図４に示すように、ステップＳＢ０６で候補位置を修正した後、修正後の候補位置１３ｉ（図６Ｂ）と、すべての既存の粒子１１の位置とが、位置条件を満たすか否かを判定する（ステップＳＢ０７）。修正後の候補位置１３ｉ（図６Ｂ）と、すべての既存の粒子１１の位置とが、位置条件を満たしている場合は、現在の候補位置１３ｉを、粒子の挿入位置として確定する（ステップＳＢ０８）。

【0032】

修正後の候補位置１３ｉ（図６Ｂ）と、いずれかの既存の粒子１１の位置とが、位置条件を満たしていない場合は、候補位置の修正を繰り返す。規定回数の修正を行っても位置条件が満たされない場合は、処理をステップＳＢ０１に戻し、粒子を挿入する候補位置を新たに決定する。規定回数は、例えば１０以上１００以下に設定するとよい。

【0033】

複数のスレッドのそれぞれは、ステップＳＢ０８で候補位置を確定した後、すべてのスレッドの候補位置が確定するまで待機する（ステップＳＢ０９）。すべてのスレッドの候補位置が確定すると、処理部５１（図１）は、確定した複数の候補位置に、それぞれ粒子を挿入する（ステップＳＢ１０）。

【0034】

次に、上記実施例の優れた効果について説明する。
本実施例では、複数のスレッドのそれぞれに、挿入する粒子の一つの挿入位置を決定する処理を割り振るため、挿入する粒子の個数が多い場合に、処理部５１の計算能力を効率よく利用することができる。

【0035】

また、上記実施例では、候補位置同士が位置条件を満たすか否かを判定し（ステップＳＢ０３）、候補位置同士の位置条件が満たされた状態が確定した後、候補位置と既存粒子の位置との位置条件が満たされるか否かを判定している（ステップＳＢ０５）。ステップＳＢ０３で候補位置同士の位置条件が満たされないと判定された場合には、位置条件を満たしていない少なくとも２つの候補条件を決定した複数のスレッドにおいて、候補位置を決定する処理（ステップＳＢ０１）を再度実行しなければならない。

【0036】

これに対して、ステップＳＢ０５で候補位置と既存粒子の位置との位置条件が満たされないと判定された場合には、位置条件を満たさない候補位置を決定したスレッドに閉じて、候補位置の修正を行うことができる。ステップＳＢ０５をステップＳＢ０３の後に実行するため、計算の手戻りが少なくなるという優れた効果が得られる。

【0037】

次に、上記実施例の変形例について説明する。
上記実施例では、ステップＳＢ０５において判断基準となる位置条件として、最低許容距離を定義している。その他に、候補位置におけるポテンシャルエネルギを判断基準として用いてもよい。例えば、周囲の既存粒子によって、候補位置に発生するポテンシャルエネルギを計算し、候補位置のポテンシャルエネルギが計算領域のポテンシャルエネルギの平均値より大きい場合に、位置条件が満たされないと判定するようにしてもよい。

【0038】

上記実施例では、ステップＳＢ０６で候補位置を修正する際に、候補位置を所定の微小距離移動させている。その他に、ニュートンラフソン法のような反復法を用いて、候補位置を修正してもよい。

【0039】

上記実施例では、計算領域１０（図３）の全体を対象として、計算領域１０の中から粒子を挿入する候補位置をランダムに選択している。その他に、計算領域１０の中に、粒子を挿入する特定の範囲を定義し、この特定の範囲内に複数の粒子を挿入するようにしてもよい。

【0040】

上記実施例では、処理部５１（図１）のスレッドの数が、挿入する粒子の個数以上である場合について説明した。挿入する粒子の個数がスレッドの数より多い場合には、スレッドの数以下の個数の粒子を挿入する処理を複数回繰り返すとよい。また、上記実施例では、複数のスレッドのそれぞれに、挿入する１つの粒子の候補位置を決定する処理を割り振っている。その他の方法として、複数のスレッドのそれぞれに、挿入する複数の粒子の候補位置を決定する処理を割り振り、１つのスレッドが同時に複数の候補位置を決定するようにしてもよい。

【0041】

次に、図７～図９を参照して、上記実施例による方法を用いて複数の粒子を挿入する計算を実際に行った結果について説明する。この計算では、立方体形状の計算領域に１００タイムステップごとに５０個の粒子を挿入する処理を実行した。粒子間相互作用には、無次元化したレナードジョーンズポテンシャルを用い、全方向に周期境界条件を適用した。ステップＳＢ０３及びステップＳＢ０５の位置条件に用いる最小許容距離を、粒子直径の１．１２倍とした。

【0042】

図７は、計算領域内の粒子数の時間変化を示すグラフである。横軸はタイムステップ数を表し、縦軸は計算領域内の粒子数を表す。図８は、計算開始時点における粒子１１の分布を示す模式図であり、図９は、１４００タイムステップを実行した時点における粒子１１の分布を示す模式図である。図７に示すように、時間の経過とともに粒子数が予定通りに増加していることがわかる。また、図８及び図９に示すように、計算領域１０内に粒子１１が想定通りに挿入されていることがわかる。

【0043】

上記実施例による方法を用いて、計算が破綻することなく、想定通りに粒子が挿入されていることが確認された。

【0044】

上記実施例は例示であり、実施例及び変形例に示した構成の部分的な置換または組み合わせが可能であることは言うまでもない。実施例及び変形例の同様の構成による同様の作用効果については実施例及び変形例ごとには逐次言及しない。さらに、本発明は上述の実施例及び変形例に制限されるものではない。例えば、種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。

【符号の説明】

【0045】

１０ 計算領域
１１、１１ｊ 既存の粒子
１２ 挿入する粒子
１３、１３ｉ、１３ｊ、１３ｋ 候補位置
５０ 入力部
５１ 処理部
５２ 出力部
５３ 記憶部

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】