再表 2021-24295 データ構造、数値計算方法、および数値計算プログラム

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

【公報種別】再公表特許(A1)

(11)【国際公開番号】WO/0

(43)【国際公開日】2021年2月11日

【発行日】2021年9月13日

(54)【発明の名称】データ構造、数値計算方法、および数値計算プログラム

(51)【国際特許分類】

   G06F 30/23 20200101AFI20210816BHJP

   G06F 17/13 20060101ALI20210816BHJP

   G16Z 99/00 20190101ALI20210816BHJP

【ＦＩ】

   G06F17/50 612H

   G06F17/13

   G16Z99/00

【審査請求】有

【予備審査請求】未請求

【全頁数】22

【出願番号】特願2019-569977(P2019-569977)

(21)【国際出願番号】PCT/0/0

(22)【国際出願日】2019年8月2日

(11)【特許番号】特許第6713626号(P6713626)

(45)【特許公報発行日】2020年6月24日

(81)【指定国】 AP(BW,GH,GM,KE,LR,LS,MW,MZ,NA,RW,SD,SL,ST,SZ,TZ,UG,ZM,ZW),EA(AM,AZ,BY,KG,KZ,RU,TJ,TM),EP(AL,AT,BE,BG,CH,CY,CZ,DE,DK,EE,ES,FI,FR,GB,GR,HR,HU,IE,IS,IT,LT,LU,LV,MC,MK,MT,NL,NO,PL,PT,RO,RS,SE,SI,SK,SM,TR),OA(BF,BJ,CF,CG,CI,CM,GA,GN,GQ,GW,KM,ML,MR,NE,SN,TD,TG),AE,AG,AL,AM,AO,AT,AU,AZ,BA,BB,BG,BH,BN,BR,BW,BY,BZ,CA,CH,CL,CN,CO,CR,CU,CZ,DE,DJ,DK,DM,DO,DZ,EC,EE,EG,ES,FI,GB,GD,GE,GH,GM,GT,HN,HR,HU,ID,IL,IN,IR,IS,JO,JP,KE,KG,KH,KN,KP,KR,KW,KZ,LA,LC,LK,LR,LS,LU,LY,MA,MD,ME,MG,MK,MN,MW,MX,MY,MZ,NA,NG,NI,NO,NZ,OM,PA,PE,PG,PH,PL,PT,QA,RO,RS,RU,RW,SA,SC,SD,SE,SG,SK,SL,SM,ST,SV,SY,TH,TJ,TM,TN,TR,TT

(71)【出願人】

【識別番号】518446075

【氏名又は名称】ＤｅｅｐＦｌｏｗ株式会社

(74)【代理人】

【識別番号】100193378

【弁理士】

【氏名又は名称】野口 明生

(72)【発明者】

【氏名】深川 宏樹

【テーマコード（参考）】

5B046

5B056

5B146

5L049

【Ｆターム（参考）】

5B046CA01

5B046FA06

5B046JA09

5B046JA10

5B056BB03

5B146BA00

5B146DJ03

5B146DJ04

5B146DJ07

5B146EA08

5L049DD02

(57)【要約】

数値計算すべき領域をセル複体に分割し、前記セル複体の双対複体を定め、前記数値計算に用いる物理量を前記セル複体および前記双対複体に割り当てたデータ構造を用いて非構造格子を用いても安定性が高い数値計算をする。また、このデータ構造を用いた数値計算方法および数値計算プログラムを提供する。

【特許請求の範囲】

【請求項1】

数値計算すべき領域をセル複体に分割し、前記セル複体の双対複体を定め、前記数値計算に用いる物理量を前記セル複体および前記双対複体に割り当てたデータ構造。

【請求項2】

前記セル複体および前記双対複体における、ｋ次のセル複体と（ｎ−ｋ）次の双対複体を対にして管理した請求項１に記載のデータ構造。

【請求項3】

前記セル複体および前記双対複体のいずれか一方に割り当てた前記物理量のデータのみを保持し、他方に割り当てた前記物理量のデータを参照する際にホッジスター演算子で変換して前記数値計算に用いる請求項２に記載のデータ構造。

【請求項4】

数値計算すべき領域をセル複体に分割し、前記セル複体の双対複体を定め、前記数値計算に用いる物理量を前記セル複体および前記双対複体に割り当てたデータ構造と、
前記数値計算に用いる支配方程式を前記セル複体および前記双対複体上の微分形式の方程式に変換したものと、
を用いて行う数値計算方法。

【請求項5】

前記方程式は、場の時間微分を、前記場および前記場の一階空間微分ならびに二階以上の高階空間微分で表現したものであることを特徴とする請求項４に記載の数値計算方法。

【請求項6】

前記微分形式の方程式を、前記セル複体と前記双対複体の間の双対性を用いて計算される近傍における前記方程式の値を用いて近似して行う請求項４または請求項５に記載の数値計算方法。

【請求項7】

前記方程式は、適切な関数Ｆを用いて下記式の形で表現されることを特徴とする請求項４から請求項６のいずれか１項に記載の数値計算方法。ただし、φは場であり、ｍは時刻の離散化のインデックスであり、ｌは前記セル複体または前記双対複体のインデックスであり、（ｌ，ｊ）はインデックスｌの前記セル複体または前記双対複体の近傍のセル複体または双対複体を表すインデックスである。

【数1】

【請求項8】

複数の演算器と前記演算器が演算に用いるデータを記憶するための階層型の記憶装置とを有するコンピュータに数値計算を実行させる数値計算プログラムであって、
前記数値計算すべき領域をセル複体に分割し、前記セル複体の双対複体を定め、前記セル複体および前記双対複体に割り当てた前記数値計算に用いる物理量を前記階層型の記憶装置に格納するステップと、
前記数値計算に用いる支配方程式を前記セル複体および前記双対複体上の微分形式の方程式に変換したものを用いて、前記セル複体および前記双対複体に割り当てた物理量に関する数値計算を行うステップと、
を含む数値計算プログラム。

【請求項9】

前記階層型の記憶装置に格納されている物理量をホッジスター演算子で変換するステップを含む請求項８に記載の数値計算プログラム。

【請求項10】

前記セル複体および前記双対複体上の微分形式の方程式をそれぞれ前記複数の演算器に割り当てて並列処理することを特徴とする請求項８または請求項９に記載の数値計算プログラム。

【請求項11】

前記方程式は、場の時間微分を、前記場および前記場の一階空間微分ならびに二階以上の高階空間微分で表現したものであることを特徴とする請求項８から請求項１０のいずれか１項に記載の数値計算プログラム。

【請求項12】

前記微分形式の方程式を、前記セル複体と前記双対複体の間の双対性を用いて計算される近傍における前記方程式の値を用いて近似して行う請求項８から請求項１１のいずれか１項に記載の数値計算プログラム。

【請求項13】

前記方程式は、適切な関数Ｆを用いて下記式の形で表現されることを特徴とする請求項８から請求項１２のいずれか１項に記載の数値計算プログラム。ただし、φは場であり、ｍは時刻の離散化のインデックスであり、ｌは前記セル複体または前記双対複体のインデックスであり、（ｌ，ｊ）はインデックスｌの前記セル複体または前記双対複体の近傍のセル複体または双対複体を表すインデックスである。

【数2】

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、データ構造、数値計算方法、および数値計算プログラムに関する。

【背景技術】

【0002】

現在の設計工程ではＣＡＤ（Ｃｏｍｐｕｔｅｒ−Ａｉｄｅｄ Ｄｅｓｉｇｎ）が一般的に用いられている。コンピュータを用いた作図は、手作業で作図することに比べて作業負担が少ないのはもちろんのこと、修正作業や編集作業も大幅に効率化できるからである。また、デジタル化された図面は、データの共有にも優れており、共同作業を支援することができる。

【0003】

一方、研究開発の工程でもＣＡＥ（Ｃｏｍｐｕｔｅｒ−Ａｉｄｅｄ ｅｎｇｉｎｅｅｒｉｎｇ）も活用が進められている。研究開発の工程では、多くの試作品の強度などの性能が検証されるが、ＣＡＥを活用すると、試作品の性能をコンピュータ上のシミュレーションで検証することが可能だからである。

【先行技術文献】

【非特許文献】

【0004】

【非特許文献1】深川宏樹，陽解法を用いた軸受の流体構造連成解析，スーパーコンピューティングニュース，２０１８，０８，Ｖｏｌ．２０，Ｎｏ．Ｓｐｅｃｉａｌ Ｉｓｓｕｅ ２，ｐｐ． ４−９

【発明の概要】

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

【0005】

上記のように、製品の研究開発や設計で活用されているＣＡＤとＣＡＥであるが、両者のデータを流用するには幾つかの問題がある。ＣＡＤに好ましいデータ構造とＣＡＥに好ましいデータ構造があり、一般にこれらは一致しないからである。

【0006】

例えば、ＣＡＥに用いるためには、構造格子とよばれる６面体で構成されたデータ構造が好ましい。構造格子で構成されたデータ構造では、数値計算との相性がよいだけではなく、スタッガード格子のような構造格子に特有の計算手法を用いることができる。

【0007】

スタッガード格子とは、構造格子を半周期ずらした位置に計算点を設けて、構造格子の中心にスカラー量を割り当て、構造格子の面の中心に当該面に垂直なベクトル量を割り当てることで、計算の安定性を高める手法である。例えば、スタッガード格子を用いた流体解析では、構造格子の中心に圧力（スカラー量）という物理量を割り当て、構造格子の面の中心に流速（ベクトル量）という物理量を割り当てることで数値計算が行われる。

【0008】

一方、ＣＡＤでは、非構造格子が求められることもある。構造格子では、形状がある程度単純なものにしか対応することができないからである。マップドメッシュと呼ばれる手法を用いて立方体を変換した構造格子を用いれば、ある程度柔軟な形状に対応でき、ＣＡＥでも利用し得る（非特許文献１参照）。しかしながら、この手法を用いることにも限界がある。

【0009】

そこで、非構造格子でも用いることができる数値計算の手法が求められており、その一つの例として有限体積法が知られている。有限体積法は、コントロールボリュームと呼ばれる領域で解析する領域を分割し、各コントロールボリュームにおいて保存方程式を適用する。すると、コントロールボリューム同士が重ならない限り、解析する領域全体でも保存性が満足される。

【0010】

しかしながら、有限体積法は、非構造格子でも用いることができるものの、計算安定性の観点では、スタッガード格子を用いた計算手法に劣る。

【0011】

本発明は、上記に鑑みてなされたものであり、その目的は、非構造格子を用いても安定性が高い数値計算をすることができるデータ構造、数値計算方法、および数値計算プログラムを提供することである。

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

【0012】

上記課題を解決するために、本発明の一態様にかかるデータ構造は、数値計算すべき領域をセル複体に分割し、前記セル複体の双対複体を定め、前記数値計算に用いる物理量を前記セル複体および前記双対複体に割り当てる。

【0013】

また、本発明の一態様にかかる数値計算方法は、数値計算すべき領域をセル複体に分割し、前記セル複体の双対複体を定め、前記数値計算に用いる物理量を前記セル複体および前記双対複体に割り当てたデータ構造と、前記数値計算に用いる支配方程式を前記セル複体および前記双対複体上の微分形式の方程式に変換したものと、を用いて行う。

【0014】

また、本発明の一態様にかかる数値計算プログラムは、複数の演算器と前記演算器が演算に用いるデータを記憶するための階層型の記憶装置とを有するコンピュータに数値計算を実行させる数値計算プログラムであって、前記数値計算すべき領域をセル複体に分割し、前記セル複体の双対複体を定め、前記セル複体および前記双対複体に割り当てた前記数値計算に用いる物理量を前記階層型の記憶装置に格納するステップと、前記数値計算に用いる支配方程式を前記セル複体および前記双対複体上の微分形式の方程式に変換したものを用いて、前記セル複体および前記双対複体に割り当てた物理量に関する数値計算を行うステップと、を含む。

【発明の効果】

【0015】

本発明によれば、非構造格子を用いても安定性が高い数値計算をすることができるデータ構造、数値計算方法、および数値計算プログラムを提供することができる。

【図面の簡単な説明】

【0016】

【図1】図１は、構造格子と非構造格子の例を示す図である。

【図2】図２は、４面体の複体構造を示す図である。

【図3】図３は、バウンダリ準同型を例示する図である。

【図4】図４は、セル複体と双対複体の関係を例示する図である。

【図5】図５は、セル複体と双対複体のペア管理を例示する図である。

【図6】図６は、共有メモリ型の階層構造の概念を示す図である。

【図7】図７は、分散メモリ型の階層構造の概念を示す図である。

【図8】図８は、近傍と双対性の関係を示す図である。

【図9】図８は、参照の局所性が高い陽解法の概念を示す図である。

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

【0017】

以下、図面を参照しながら、本発明の実施形態にかかるデータ構造、数値計算方法、および数値計算プログラムを詳細に説明する。ただし、以下の説明で参照される図面は模式的なものであり、寸法またはその比率が実際のものとは異なる場合がある。

【0018】

本実施形態にかかるデータ構造は、構造格子に限定せず、非構造格子を含む一般構造格子を用いる。図１は、構造格子と非構造格子の例を示す図である。構造格子とは、図１（ａ）に示すように、６面体を一単位として計算領域を分割する計算格子である。一方、非構造格子とは、構造格子以外のものをいう。非構造格子とは、したがって、図１（ｂ）に示すように、４面体や５面体を含む任意の多面体を組み合わせて計算領域を分割する計算格子のことをいう。

【0019】

ただし、本実施形態にかかるデータ構造における多面体とは、数学の定義でいうところのセル複体である。つまり、単なる形状としての多面体を考えるのではなく、複体構造を有する多面体である。セル複体とは、ｎ次元セルと呼ばれるｋ次元ディスクと同相であるユニットが所定の階層構造を有しているものをいう。一般的な多面体は、面や辺といったｋ次元ディスクと同相なユニットで構成されており、セル複体は多面体の一般化に相当している。

【0020】

例えば、図２のように、４面体は、４個の０次複体（頂点）と６個の１次複体（辺）と４個の２次複体（面）と１個の３次複体（体；バルク）とに分解することができる。本実施形態にかかるデータ構造では、計算領域を任意の多面体で分割するが、このときの多面体には、上記のような複体としての構成を有しているもの考える。言い換えると、本実施形態にかかるデータ構造では、数値計算すべき領域をセル複体に分割することになる。なお、６面体も、複体構造を有する多面体である。したがって、本実施形態にかかるデータ構造は、形状としての６面体の使用を排除するものではなく、複体構造を有するセル複体に数値計算すべき領域を分割することで、一般の多面体の使用を許容するものである。また、面は曲面であっても計算可能である。

【0021】

ここで、各複体には向き付けを行う。向き付けられたｋ次複体が生成する自由アーベル群Ｃ_ｋ（Ｘ）の列（チェイン複体）には、バウンダリ準同型と呼ばれるチェイン準同型∂_ｋ：Ｃ_ｋ（Ｘ）→Ｃ_ｋ−１（Ｘ）が定義される。図３は、向き付けられた２次複体（面）から向き付けられた１次複体（辺）へのバウンダリ準同型∂_２：Ｃ_２（Ｘ）→Ｃ_１（Ｘ）の例を示す図である。

【0022】

このバウンダリ準同型は、チェイン複体Ｃ_ｋ（Ｘ）に対して、以下の様に作用し、∂_ｋ∂_ｋ−１＝０を満たす。

【0023】

【数1】

【0024】

一方、複体の構造を有した多面体の分割には、双対複体を定義することができる。通常の数値計算では、数値計算を行うべき領域を分割したものをメッシュと呼ぶので、本実施形態にかかるデータ構造で用いる双対複体は、双対メッシュと呼び得るものである。ただし、本実施形態にかかるデータ構造では、領域を分割するメッシュは単に形状だけではなく複体としての構造を有しており、同様に、双対メッシュも単に形状だけではなく複体としての構造を有しているものを考える。

【0025】

具体的には、セル複体（メッシュ）の３次複体（体）の中心を定め、これを双対複体（双対メッシュ）の０次複体（頂点）にする。そして、隣り合う（２次複体（面）を共有する）セル複体（メッシュ）の中心として定められた双対複体（双対メッシュ）の０次複体（頂点）を結ぶラインを双対複体（双対メッシュ）の１次複体（辺）とする。このように定められた双対複体（双対メッシュ）の１次複体（辺）によって張られた面を双対複体（双対メッシュ）の２次複体（面）とし、さらに、双対複体（双対メッシュ）の２次複体（面）によって閉ざされた空間を双対複体（双対メッシュ）の３次複体（体）とする。図４は、セル複体と双対複体の関係を例示する図である。なお、図４では図示の都合上、２次元の多面体（つまり多角形）における双対複体を図示している。図中では、セル複体を実線で、双対複体を破線で示している。

【0026】

このように、セル複体（メッシュ）と双対複体（双対メッシュ）では、ｋ次の複体とｎ−ｋ次の複体の間には双対性が存在する。具体的には、セル複体（メッシュ）と双対複体（双対メッシュ）では、それぞれ、０次複体（頂点）と３次複体（体）が一対一に対応し、１次複体（辺）と２次複体（面）が一対一に対応し、２次複体（面）と１次複体（辺）が一対一に対応し、３次複体（体）と０次複体（頂点）が一対一に対応する。したがって、セル複体（メッシュ）と双対複体（双対メッシュ）の間に同型写像★（一対一対応関係）を定義することができる（例えば、後に詳述する図５も参照。）。

【0027】

また、双対複体はそれ自体も複体である。したがって、双対複体の列にもバウンダリ準同型を定義することが可能である。

【0028】

ところで、上記の話は純粋な数学的技術である。一方、上記の数学的技術を数値計算に用いるために、本実施形態にかかるデータ構造では双対複体上に物理量を割り当てる。すなわち、本実施形態にかかるデータ構造では、複体構造を有する多面体で数値計算を行うべき領域を分割し、当該多面体の複体構造に関する双対複体に数値計算に用いる物理量を割り当てる。

【0029】

また、ここで物理量とは、例えば流体解析の場合、圧力や流速であり、電磁場解析の場合、電磁ポテンシャルである。そして、一般に圧力や流速などの物理量は連続量である。これら連続量である物理量は、数値計算を行うべき領域の分割に対応したセル複体および双対複体に割り当てることで離散化される。

【0030】

この離散化の方法は簡単に言うとド・ラーム複体を離散化することに対応するので、ここで簡単に連続版のド・ラーム複体について説明する。Ｍを微分可能多様体としΩ^ｋ（Ｍ）をＭ上のｋ次微分形式の空間とする。なお、Ω^０（Ｍ）はＭ上の滑らかな関数の空間である。ｋ次微分形式には外微分ｄ^ｋという作用素が定義され、Ｍ上のｋ次微分形式の空間Ω^ｋ（Ｍ）の列は、チェイン複体を構成しｄ^ｋｄ^ｋ＋１＝０となる。

【0031】

【数2】

【0032】

一方、ド・ラーム複体の離散化の基本的な考え方は、上記のように離散化された複体上で微分形式を積分することで値を決める。本実施形態にかかるデータ構造では、数値計算に応用することを目的としているので、数値計算に用いる物理量（物理的な場）を双対複体上で積分することで、双対複体上に物理量を割り当てる。なお、このことは後に詳述するホッジスターの性質からも解るように、数値計算に用いる物理量（物理的な場）をセル複体上で積分することで、セル複体上に物理量を割り当てることと双対性を介して同値である。

【0033】

具体的には、σをｋ次の双対複体の元とし、ωをｋ次の微分形式とする場合、離散化された微分形式の値を下式で定める。

【数3】

【0034】

例えば、０次の複体とは点のことであり、０次の微分形式とはスカラー場とみなすことができるので、上記定義に従う離散化された微分形式は、当該点におけるスカラー場の値である。本実施形態にかかるデータ構造では、数値計算に用いる物理量（物理的な場）を双対複体上に割り当てるので、例えば流体解析における圧力分布の場合、数値計算すべき領域をセル複体に分割したものの双対複体の各頂点での圧力をサンプリングすることになる。このことは、数値計算すべき領域のセル複体での分割（メッシュ）と双対複体での分割（双対メッシュ）との双対性を考えると、数値計算すべき領域の分割を双対複体の各頂点で代表させて圧力をサンプリングすることになる。

【0035】

また、１次の複体とは辺のことであり、１次の微分形式は音楽準同型を介してベクトル場とみなすことができる。よって、辺上にベクトル場をサンプリングする場合、上記定義に従い、１次の複体（辺）上の１次の微分形式の積分としてサンプリングすることができる。本実施形態にかかるデータ構造では、例えば流体解析における流速分布の場合、数値計算すべき領域をセル複体に分割したものの双対複体の各辺に沿った流速をサンプリングすることになる。

【0036】

既に指摘したように、セル複体での分割（メッシュ）と双対複体での分割（双対メッシュ）との間には双対性があり、しかも、数値計算上も関連しているので、本実施形態にかかるデータ構造では、セル複体での分割（メッシュ）と双対複体での分割（双対メッシュ）とをペアにして管理する。ここで、ペアにして管理するとは、各複体に対するインデックスを共通化ないし関連化することを意味する。

【0037】

すなわち、図５に示すように、セル複体における３次の複体（体：バルク）と双対複体における０次の複体（点）をペアにし、セル複体おける２次の複体（面）と双対複体における１次の複体（辺）をペアにし、セル複体における１次の複体（辺）と双対複体における２次の複体（面）をペアにし、セル複体における０次の複体（点）と双対複体における３次の複体（体；バルク）をペアにして管理する。このペアは、上記したセル複体（メッシュ）と双対複体（双対メッシュ）の間の同型写像★によって定める。

【0038】

このように、セル複体での分割（メッシュ）と双対複体での分割（双対メッシュ）をペアにして管理し、インデックスを共通化ないし関連化することは、後述するように、コンピュータ上の数値計算プログラムの実行における参照の局所性を高める。また、双対複体での分割（双対メッシュ）は、構造格子におけるスタッガード格子と同様の機能を有しているので、本実施形態にかかるデータ構造安定性を用いれば、安定性が高い数値計算をすることができる。

【0039】

本実施形態にかかるデータ構造では、数値計算に用いる物理量（物理的な場）を双対複体上で積分することで、双対複体上に物理量を割り当てる。このように、双対複体上に物理量を割り当てる大きなメリットは、離散版のド・ラーム複体では連続版のド・ラーム複体で成り立つ微分形式の性質の多くが継承されることにある。したがって、離散版のド・ラーム複体において定義された微分形式を離散微分形式と呼ぶことがある。

【0040】

例えば、ストークスの定理の類似の関係式がなりたつ。

【0041】

【数4】

【0042】

また、ホッジスター演算子＊の離散化版が、セル複体（メッシュ）と双対複体（双対メッシュ）との間の同型写像★を用いて以下の様に定義される。

【0043】

【数5】

【0044】

上記定義から明らかなように、離散化されたホッジスター演算子は、セル複体（メッシュ）上の値と双対複体（双対メッシュ）上の値を相互に変換することを可能にする。

【0045】

また、微分形式では、一般的なベクトル解析における演算子が外微分を用いて記述することができる。そして、これら微分形式を用いた演算が、離散版のド・ラーム複体でも成立する。以下、勾配ｇｒａｄ、発散ｄｉｖ、回転ｒｏｔについて説明する。

【0046】

勾配ｇｒａｄは、音楽同型（Ｍｕｓｉｃａｌ ｉｓｏｍｏｒｐｈｉｓｍ）♯を用いて以下の様に記述できる。

【0047】

【数6】

【0048】

発散ｄｉｖは、音楽同型（Ｍｕｓｉｃａｌ ｉｓｏｍｏｒｐｈｉｓｍ）♭とホッジスター演算子＊を用いて以下の様に記述できる。

【0049】

【数7】

【0050】

回転ｒｏｔは、音楽同型（Ｍｕｓｉｃａｌ ｉｓｏｍｏｒｐｈｉｓｍ）＃，♭とホッジスター演算子＊を用いて以下の様に記述できる。

【0051】

【数7】

【0052】

その他、上記例示した以外にも、離散微分形式では、ラプラシアンやリー微分や内部積などの演算子も定義することが可能である。

【0053】

以上のように、本実施形態にかかるデータ構造は、構造格子に限定せず、非構造格子を含む一般構造格子を用いる。一方、本実施形態にかかるデータ構造は、数値計算すべき領域をセル複体に分割し、当該セル複体の双対複体を定める。そして、数値計算に用いる物理量（物理的な場）をセル複体および双対複体上に物理量を割り当てる。

【0054】

本実施形態にかかるデータ構造では、離散微分形式を用いた演算によって、３次元ユークリッド空間上のベクトル解析で用いられる演算子のすべてを実行可能である。したがって、本実施形態にかかるデータ構造は、多くの数値計算に応用可能である。以下、本実施形態にかかるデータ構造の応用例として、流体解析および電磁場解析に対する応用について説明する。

【0055】

〔流体解析〕
流体解析では、支配方程式として、連続の方程式とナビエ・ストークス方程式が用いられる。連続の方程式とナビエ・ストークス方程式は、以下の様に微分形式を用いて記述することが可能である。

【0056】

【数8】

【0057】

ここで、ρは流体内の質量密度であり、*ρは、その双対ホッジであり、要素内の総質量を表す。速度場u^#に関するＬｉｅ微分をL_u#と書く。計量テンソルをg( , )とし、u:=g(u^#, )と定める。圧力Pは質量密度ρの関数とする。なお、Δはホッジラプラシアンとし、すなわちΔ＝δｄ＋ｄδ，δ＝（−１）^ｋ＊^−１ｄ＊である。

【0058】

上記式の形からも解るように、連続の方程式とナビエ・ストークス方程式は、本実施形態にかかるデータ構造上の演算として実行可能である。

【0059】

〔電磁場解析〕
電磁場解析では、支配方程式として、電磁ポテンシャルが用いられる。アンペールの法則およびガウスの法則は、ローレンツゲージを用いてゲージ固定をすると以下の様に微分形式を用いて記述することが可能である。

【0060】

【数9】

【0061】

ここで、Ａはベクトルポテンシャル、φはスカラーポテンシャル、Jは電流密度、ρは電荷である。上記式の形からも解るように、電磁ポテンシャルは、本実施形態にかかるデータ構造上の演算として実行可能である。

【0062】

以上の様に、本実施形態にかかるデータ構造は、多くの数値計算に応用可能である。以下、本実施形態にかかる数値計算および数値計算プログラムをコンピュータ上に実現することに関し説明する。

【0063】

〔装置構成〕
図６は、共有メモリ型の階層構造の概念を示す図であり、図７は、分散メモリ型の階層構造の概念を示す図である。本発明の実施形態にかかる数値計算方法および数値計算プログラムは、複数の演算器とこれらが演算に用いるデータを記憶するための階層型の記憶装置とを有するコンピュータにおいて好適に実施される。そこで、図６および図７を参照しながら、複数の演算器と階層構造の記憶装置を備えるコンピュータにおける計算処理について説明する。なお、複数の演算器を備えるコンピュータは並列計算機や並列コンピュータとも呼ばれ、この並列計算機上の計算処理は並列処理や並列コンピューティングとも呼ばれる。

【0064】

図６に示される概念図には、複数の演算器１０_１〜１０_ｎとこれらが演算に用いるデータを記憶するための階層型の記憶装置１１とを備える並列計算機１００が記載されている。階層型の記憶装置１１とは、各演算器１０_１〜１０_ｎと主記憶装置１２との間に階層化された記憶装置を配したものである。図１に示される例では、演算器１０_１〜１０_ｎから順に、レジスタ１３_１〜１３_ｎ、Ｌ１キャッシュメモリ１４_１〜１４_ｎ、Ｌ２キャッシュメモリ１５_１〜１５_ｍ、主記憶装置１２の構成を有している。なお、図６に示される構成は、一つの例であり、階層型の記憶装置１１の構成はこれに限定されるものではなく、階層の程度としても、Ｌ２キャッシュメモリ１５_１〜１５_ｍと主記憶装置１２との間にＬ３キャッシュメモリを備えるなどの変形例がある。

【0065】

これらの記憶装置の処理速度は、主記憶装置１２、Ｌ２キャッシュメモリ１５_１〜１５_ｍ、Ｌ１キャッシュメモリ１４_１〜１４_ｎ、レジスタ１３_１〜１３_ｎの順で高速であり、メモリのデータ転送能力であるメモリバンド幅は、演算器１０_１〜１０_ｎの演算性能に比べて低い（Ｂｙｔｅ／Ｆｌｏｐが低い）。このような特性から、参照の局所性を高めることが演算性能を発揮させるために重要となる。

【0066】

なお、参照の局所性が高い状態とは、一度参照されたデータが再度参照されることや当該データに近いアドレスのデータが参照される可能性が高いという状態である。言い換えると、演算器１０_１〜１０_ｎがレジスタ１３_１〜１３_ｎまたはＬ１キャッシュメモリ１４_１〜１４_ｎもしくはＬ２キャッシュメモリ１５_１〜１５_ｍなど上位の記憶装置に記憶されているデータを用いて計算を行うことができる状態である。

【0067】

逆に言えば、演算器１０_１〜１０_ｎが主記憶装置１２に記憶されているデータを用いなければいけない状況を少なくする戦略が必要である。例えば、そのような状況は、演算器１０_１が計算した結果を用いて演算器１０_ｎが計算をしなければいけないときに起きる。また、演算器１０_１が計算した結果を用いて演算器１０_ｎが計算をしなければいけない場合には、演算器１０_１の計算が終わらなければ演算器１０_ｎの計算を開始することができず、待ち時間が発生してしまうという問題も発生してしまう。

【0068】

図７は、分散メモリ型の階層構造の概念を示す図であり、大型計算機システムなどで採用されている。図７に示される並列計算機２００では、各ノード１６_１〜１６_ｎが複数の演算器とこれらが演算に用いるデータを記憶するための階層型の記憶装置とを備えており、さらに、各ノード１６_１〜１６_ｎにおける主記憶装置１２_１〜１２_ｎがネットワーク１７を介して互いに共有されている。ここで、主記憶装置１２_１〜１２_ｎの共有とは、主記憶装置１２_１〜１２_ｎに単一のアドレスを付与することで、各ノード１６_１〜１６_ｎの演算器が、主記憶装置１２_１〜１２_ｎに記憶されているデータを相互に参照することができることをいう。

【0069】

上記のような構成のコンピュータでも、参照の局所性を高めることが演算性能を発揮させるために重要となる。特に、あるノード１６_１〜１６_ｍの演算器が、他のノード１６_１〜１６_ｍの主記憶装置１２_１〜１２_ｎに記憶されているデータを参照するような状況を可能な限り少なくすることが重要である。

【0070】

上記観点から、現在の主流となっている計算機のアーキテクチャでは参照の局所性が高い数値計算方法および数値計算プログラムが求められている。そこで、本発明の実施形態にかかる数値計算方法および数値計算プログラムは、陽解法を用いた数値計算を行う。以下、陽解法に関する説明を行う。

【0071】

〔陽解法〕
陽解法とは、第１時刻tにおける状態量のみから微小時間後の第2時刻t+Δtにおける位置xの状態量を計算する方法である。一方、陰解法は、第２時刻t+Δtにおける位置xの状態量を計算するために、第２時刻t+Δtにおける他の位置x_iの状態量も用いる方法である。端的には、陰解法は、第２時刻t+Δtにおける位置xの状態量を計算するために、第２時刻t+Δtにおける他の位置x_iの状態量も含んだ連立方程式を解く必要があり、陽解法では、このような連立方程式を解く必要がないことを意味する。

【0072】

このように、陽解法では、連立方程式を解く必要がないので、時間ステップ当たりの計算量およびメモリ転送量が少ない。さらに、第２時刻t+Δtにおける位置xの状態量を計算するために、第１時刻tにおける位置xの近傍のみの状態量から計算するようにすれば、より一層、参照の局所性が高まる。そこで、以下で参照の局所性が高い陽解法を適用し得る条件について説明する。

【0073】

結論としては、ある領域の場φ(n,t)の時間発展が、以下の式（１）に示すように、場φ(n,t)と当該場の空間に関する一階微分φ′(n,t)と二階微分φ″(n,t)の関数fで表現することができればよい。なお、位置をnとし、時刻をtとしている。

【0074】

【数10】

【0075】

上記式のように表現される場φ(n,t)を時刻および各複体に関して離散化する。時刻tの離散化は、インデックスmを用いてt(m+1)=t(m)+Δtとする。一方、位置ｎの離散化は、各複体のインデックスｌを用いてn(l)とする。すると、各複体のインデックスlに対応した頂点n(l)の近傍の頂点n(l,j)は、例えば図８のように表すことができる。なお、図８は、図示の都合上、２次元版を記載しているが、３次元でも同様に近傍を表すことができる。また、ここでは、インデックスlに対応した頂点n(l)の近傍の頂点n(l,j)のインデックスを(l,j)としているが、表示の都合上であり双対複体のインデックスが２変数であることを意味するものではない。これは、セル複体での分割（メッシュ）と双対複体での分割（双対メッシュ）のインデックスを共通化ないし関連化して管理していることを意味している。

【0076】

図８に示すように、セル複体の頂点n(l)の近傍には、双対複体の頂点n(l,j)が存在している。そして、双対複体の頂点n(l,j)は、次のように計算できる。まず、セル複体の頂点n(l)の双対★n(l)を求める。この★は、先述した同型写像である。なお、★n(l)は、図８中の斜線で示される領域である。そして、セル複体の頂点n(l)の周辺の双対複体の頂点n(l,j)は、バウンダリ準同型∂を用いて、∂(★n(l))の要素となっている。

【0077】

ここで、頂点n(l)から近傍の頂点n(l,j)ためには、同型写像★とバウンダリ準同型∂を用いればよい。そして、先述したように、本実施形態のデータ構造では、セル複体での分割（メッシュ）と双対複体での分割（双対メッシュ）とを同型写像★を用いてペアにして管理し、頂点n(l)と頂点n(l,j)に割り当てたデータは、コンピュータ上のメモリにおいても近い位置に保持する。このように、本実施形態のデータ構造の利用した数値計算プログラムは、セル複体での分割（メッシュ）と双対複体での分割（双対メッシュ）とをペアにして管理することで、コンピュータ上の実行において参照の局所性が高いデータ配置を実現する。

【0078】

そして、この近傍を用いて、場の一階微分φ′(n(l),t(m))と二階微分φ″(n(l),t(m))は、近傍の双対複体の頂点n(l,j)における場φ(n(l,j),t(m))の関数で近似できる。したがって、f(φ(n,t),φ′(n,t),φ″(n,t))は、F(φ(n(l,j),t(m)))と近似できる。

【0079】

結果、上記式は、以下のように表現される。

【0080】

【数11】

【0081】

また、上記式をインデックス(l,m)で表示すれば、以下のように表現される。

【0082】

【数12】

【0083】

上記式から解るように、φ(l,m+1)がφ(l,m)とφ((l,j),m)で定まる。ここで重要なことは、上記式の右辺には、m+1が含まれないことである。このことは、第１時刻t(m)における状態量のみから微小時間後の第2時刻t(m+1)における頂点n(l)の状態量を計算することができることを意味している。

【0084】

また、上記式の右辺には、インデックス(l,j)が含まれているが、これは複体の頂点n(l)の近傍の双対複体の頂点をn(l,j)と表したのだから、φ((l,j),m)は、第１時刻t(m)における頂点n(l)の近傍の頂点の場φの値である。図９は、参照の局所性が高い陽解法の概念を示す図である。上記式に示される方程式は、図９に示されるように、時刻のインデックスm+1における位置インデックスlのφ(l,m+1)を求める際に、時刻のインデックスmにおける位置のインデックスlのφ(l,m)とその近傍のインデックス(l,j)におけるφ((l,j),m)を用いるだけでよいことを意味している。インデックスlは、φの値をメモリに記憶する際のアドレスの近さに対応するので、φ(l,m+1)の計算に際し、φ(l,m)とφ((l,j),m)という参照の局所性が高いデータを用いていることになる。

【0085】

以上の議論から解るように、上記式（１）のように表現できる場合、参照の局所性が高い陽解法を用いることができる。また、上記式（１）の右辺に当該場の空間に関する任意の高階微分が含まれていても、同様にして本方法を用いることができる。そして、上記したナビエ・ストークス方程式、電磁ポテンシャルの方程式は、上記式（１）のように表現できる。したがって、複数の演算器と演算器が演算に用いるデータを記憶するための階層型の記憶装置とを有するコンピュータを用いて、階層型の記憶装置に記憶されている第１時刻における物理量から微小時間後の第２時刻における物理量を数値計算する際に参照の局所性が高い陽解法を適用し得る。

【0086】

つまり、本発明の実施形態にかかる数値計算プログラムは、図６に示されるような共有メモリ型の並列計算機に対する指令として好適であり、本発明の実施形態にかかる数値計算方法は、図７に示されるような共有メモリ型の並列計算機における処理として好適である。また、分散メモリ型の並列計算機における各ノードが共有メモリ型の並列計算機として構成されている、共有分散メモリ型の並列計算機に対しても、本発明の実施形態にかかる数値計算方法および数値計算プログラムが有効であることは言うまでもない。

【0087】

以上、図面を参照しながら本発明を実施形態に基づいて説明してきたが、本発明は上記の実施形態よって限定されるものではない。上記説明した数値計算方法は、適切に数値計算プログラムとして実装することができ、逆に上記説明した数値計算プログラムの実行は、数値計算方法の実施とみなし得る。また、上記説明した数値計算プログラムは、コンピュータが読み取り可能な非トランジェント（non-transient）な媒体に記録された状態で生産等の実施を行うことができる。

【符号の説明】

【0088】

１００，２００ 並列計算機
１０_１〜１０_ｎ 演算器
１１ 階層型の記憶装置
１２ 主記憶装置
１３_１〜１３_ｎ レジスタ
１４_１〜１４_ｎ Ｌ１キャッシュメモリ
１５_１〜１５_ｍ Ｌ２キャッシュメモリ
１６_１〜１６_ｎ ノード
１７ ネットワーク

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【手続補正書】

【提出日】2020年3月30日

【手続補正1】

【補正対象書類名】特許請求の範囲

【補正対象項目名】全文

【補正方法】変更

【補正の内容】

【特許請求の範囲】

【請求項1】

複数の演算器と前記演算器が演算に用いるデータを記憶するための階層型の記憶装置とを有するコンピュータにおいて形状および物理量を計算するためのデータ構造であって、 前記計算すべき領域を、前記形状に対応させてセル複体と双対複体に分解し、前記物理量を前記セル複体および前記双対複体に割り当て、
前記演算器が前記セル複体と前記双対複体に定められたホッジスター演算子とバウンダリ準同型を用いて前記記憶装置のデータを参照するデータ構造。

【請求項2】

前記セル複体におけるｋ次のセル複体と、前記双対複体における（ｎ−ｋ）と対にして、前記計算に用いるインデックスを共通化ないし関連化した請求項１に記載のデータ構造。

【請求項3】

前記セル複体および前記双対複体のいずれか一方に割り当てた前記物理量のデータのみを保持し、他方に割り当てた前記物理量のデータを参照する際に前記ホッジスター演算子で変換して前記計算に用いる請求項２に記載のデータ構造。

【請求項4】

計算すべき領域を、セル複体と双対複体に分割し、前記計算に用いる物理量を前記セル複体および前記双対複体に割り当て、前記セル複体と前記双対複体に定められたホッジスター演算子とバウンダリ準同型を用いてデータを参照するデータ構造と、
前記計算に用いる支配方程式を前記セル複体および前記双対複体上の微分形式の方程式に変換したものと、
を用いてコンピュータ上で行う数値計算方法。

【請求項5】

前記方程式は、場の時間発展を、前記場および前記場の一階空間微分ならびに二階以上の高階空間微分で表現したものであることを特徴とする請求項４に記載の数値計算方法。

【請求項6】

前記微分形式の方程式を、前記セル複体と前記双対複体の間の双対性を用いて計算される近傍における前記方程式の値を用いて近似して行う請求項４または請求項５に記載の数値計算方法。

【請求項7】

前記方程式は、適切な関数Ｆを用いて下記式の形で表現されることを特徴とする請求項４から請求項６のいずれか１項に記載の数値計算方法。ただし、φは場であり、ｍは時刻の離散化のインデックスであり、ｌは前記セル複体または前記双対複体のインデックスであり、（ｌ，ｊ）はインデックスｌの前記セル複体または前記双対複体の近傍のセル複体または双対複体を表すインデックスである。

【数1】

【請求項8】

複数の演算器と前記演算器が演算に用いるデータを記憶するための階層型の記憶装置とを有するコンピュータにおいて形状および物理量の数値計算を実行させる数値計算プログラムであって、
前記計算すべき領域を、前記形状に対応させてセル複体と双対複体に分割し、前記セル複体および前記双対複体に割り当てた前記物理量を前記階層型の記憶装置に格納するステップと、
前記演算器が前記セル複体と前記双対複体に定められたホッジスター演算子とバウンダリ準同型を用いて前記記憶装置のデータを参照するステップと、
前記計算に用いる支配方程式を前記セル複体および前記双対複体上の微分形式の方程式に変換したものを用いて前記形状および前記物理量に関する数値計算を行うステップと、
を含む数値計算プログラム。

【請求項9】

前記階層型の記憶装置に格納されている物理量を前記ホッジスター演算子で変換するステップを含む請求項８に記載の数値計算プログラム。

【請求項10】

前記セル複体および前記双対複体上の微分形式の方程式をそれぞれ前記複数の演算器に割り当てて並列処理することを特徴とする請求項８または請求項９に記載の数値計算プログラム。

【請求項11】

前記方程式は、場の時間発展を、前記場および前記場の一階空間微分ならびに二階以上の高階空間微分で表現したものであることを特徴とする請求項８から請求項１０のいずれか１項に記載の数値計算プログラム。

【請求項12】

前記微分形式の方程式を、前記セル複体と前記双対複体の間の双対性を用いて計算される近傍における前記方程式の値を用いて近似して行う請求項８から請求項１１のいずれか１項に記載の数値計算プログラム。

【請求項13】

前記方程式は、適切な関数Ｆを用いて下記式の形で表現されることを特徴とする請求項８から請求項１２のいずれか１項に記載の数値計算プログラム。ただし、φは場であり、ｍは時刻の離散化のインデックスであり、ｌは前記セル複体または前記双対複体のインデックスであり、（ｌ，ｊ）はインデックスｌの前記セル複体または前記双対複体の近傍のセル複体または双対複体を表すインデックスである。

【数2】

【国際調査報告】