(57)【要約】超大規模構造問題を解く方法として、領域分割法（ＤＤＭ）、Ｎｅｕｍａｎｎ前処理つきＤＭＭ、ＢＤＤ法、並列ＣＧ法が知られているが、発散して解が求まらない、計算時間が極めて長くなるという問題がある。本発明は、自由度が百万以上の超大規模構造問題を解く際、解が発散することなく求まり、しかも短時間の計算時間で解が求まる、計算パーフォンマンスの高い超大規模構造問題を解くシステムを提供するものである。自由度１００万以上の超大規模構造問題を解く並列有限要素法計算システムであって、領域分割をする手段、部分領域を各プロセッサーの担当する部分へ分配する手段、剛性マトリックスを作成する手段、部分領域全体重ね合わせ運動を設定する手段、全自由度の前処理つき射影ＣＧ法の初期設定を行う手段、全自由度の前処理つき射影ＣＧ法の反復計算をする手段及び変位解を出力する手段を備えたことを特徴とする並列有限要素法計算システムである。また、該システムを機能させるためのプログラム及び該プログラムを記録したコンピュータ読み取り可能な記録媒体を含む。
2002093412.tif 000002
【特許請求の範囲】
【請求項１】
自由度１００万以上の超大規模構造問題を解く並列有限要素法計算システムであって、領域分割を行う手段、部分領域を各プロセッサーの担当する部分へ分配する手段、剛性マトリックスを作成する手段、部分領域全体重ね合わせ運動を設定する手段、全自由度の前処理つき射影ＣＧ法の初期設定を行う手段、全自由度の前処理つき射影ＣＧ法の反復計算をする手段及び変位解を出力する手段を備えたことを特徴とする並列有限要素法計算システム。
【請求項２】
前記部分領域全体重ね合わせ運動を設定する手段が、全自由度表示の射影子を作成する手段、部分領域全体重ね合わせ運動マトリックスを作成する手段及び部分領域全体重ね合わせ運動マトリックスをＬＵ分解する手段を備えたことを特徴とする請求の範囲１に記載の並列有限要素法計算システム。
【請求項３】
前記全自由度の前処理つき射影ＣＧ法の初期設定を行う手段が、全自由度の初期変位を設定する手段、全自由度の初期残差を計算する手段、対角スケーリング前処理の計算をする手段、全自由度のＣＯＡＲＳＥ  ＧＲＩＤ前処理の計算をする手段及び全自由度のＣＧ法探索方向のベクトル初期値を設定する手段を備えたことを特徴とする請求の範囲１に記載の並列有限要素法計算システム。
【請求項４】
前記全自由度の前処理つき射影ＣＧ法の反復計算をする手段が、全自由度の変位を更新する手段、全自由度の残差を更新する手段、対角スケーリング前処理の計算をする手段、全自由度のＣＯＡＲＳＥ  ＧＲＩＤ前処理の計算をする手段、全自由度のＣＧ法探索方向のベクトルを更新する手段及び収束を判定する手段を備えたことを特徴とする請求の範囲１に記載の並列有限要素法計算システム。
【請求項５】
自由度１００万以上の超大規模構造問題を解く並列有限要素法計算システムであって、領域分割を行う手段、部分領域を各プロセッサーの担当する部分へ分配する手段、剛性マトリックスを作成する手段、部分領域全体重ね合わせ運動を設定する手段、全自由度の前処理つき射影ＣＧ法の初期設定を行う手段、全自由度の前処理つき射影ＣＧ法の反復計算をする手段及び変位解を出力する手段として機能させるための並列有限要素法計算プログラム。
【請求項６】
前記部分領域全体重ね合わせ運動を設定する手段が、全自由度表示の射影子を作成する手段、部分領域全体重ね合わせ運動マトリックスを作成する手段及び部分領域全体重ね合わせ運動マトリックスをＬＵ分解する手段として機能させるための請求の範囲５に記載の並列有限要素法計算プログラム。
【請求項７】
前記全自由度の前処理つき射影ＣＧ法の初期設定を行う手段が、全自由度の初期変位を設定する手段、全自由度の初期残差を計算する手段、対角スケーリング前処理の計算をする手段、全自由度のＣＯＡＲＳＥ  ＧＲＩＤ前処理の計算をする手段及び全自由度のＣＧ法探索方向のベクトル初期値を設定する手段として機能させるための請求の範囲５に記載の並列有限要素法計算プログラム。
【請求項８】
前記全自由度の前処理つき射影ＣＧ法の反復計算をする手段が、全自由度の変位を更新する手段、全自由度の残差を更新する手段、対角スケーリング前処理の計算をする手段、全自由度のＣＯＡＲＳＥ  ＧＲＩＤ前処理の計算をする手段、全自由度のＣＧ法探索方向のベクトルを更新する手段及び収束を判定する手段として機能させるための請求の範囲５に記載の並列有限要素法計算プログラム。
【請求項９】
自由度１００万以上の超大規模構造問題を解く並列有限要素法計算システムであって、領域分割を行う手段、部分領域を各プロセッサーの担当する部分へ分配する手段、剛性マトリックスを作成する手段、部分領域全体重ね合わせ運動を設定する手段、全自由度の前処理つき射影ＣＧ法の初期設定を行う手段、全自由度の前処理つき射影ＣＧ法の反復計算をする手段及び変位解を出力する手段として機能させるための並列有限要素法計算プログラムを記録したコンピュータ読み取り可能な記録媒体。
【請求項１０】
前記部分領域全体重ね合わせ運動を設定する手段が、全自由度表示の射影子を作成する手段、部分領域全体重ね合わせ運動マトリックスを作成する手段及び部分領域全体重ね合わせ運動マトリックスをＬＵ分解する手段として機能させるための請求の範囲９に記載の並列有限要素法計算プログラムを記録したコンピュータ読み取り可能な記録媒体。
【請求項１１】
前記全自由度の前処理つき射影ＣＧ法の初期設定を行う手段が、全自由度の初期変位を設定する手段、全自由度の初期残差を計算する手段、対角スケーリング前処理の計算をする手段、全自由度のＣＯＡＲＳＥ  ＧＲＩＤ前処理の計算をする手段及び全自由度のＣＧ法探索方向のベクトル初期値を設定する手段として機能させるための請求の範囲９に記載の並列有限要素法計算プログラムを記録したコンピュータ読み取り可能な記録媒体。
【請求項１２】
前記全自由度の前処理つき射影ＣＧ法の反復計算をする手段が、全自由度の変位を更新する手段、全自由度の残差を更新する手段、対角スケーリング前処理の計算をする手段、全自由度のＣＯＡＲＳＥ  ＧＲＩＤ前処理の計算をする手段、全自由度のＣＧ法探索方向のベクトルを更新する手段及び収束を判定する手段として機能させるための請求の範囲９に記載の並列有限要素法計算プログラムを記録したコンピュータ読み取り可能な記録媒体。
【発明の詳細な説明】
技術分野
本発明は、自由度１００万以上の超大規模構造問題を効果的に解く、超大規模並列有限要素法ソルバアルゴリズムに関する。並列ＣＧアルゴリズムを基盤とし、領域分割（Ｄｏｍａｉｎ  Ｄｅｃｏｍｐｏｓｉｔｉｏｎ）に基づく共役射影勾配法アルゴリズムを取入れたものである。当該アルゴリズムを「ＣＧＣＧ法」（Ｃｏａｒｓｅ  Ｇｒｉｄ  ＣＧ  ｍｅｔｈｏｄ，ＣＧＣＧ  ｍｅｔｈｏｄ）と称する。Ｋを当該構造問題の剛性行列とし、領域分割による部分領域単位の自由度の空間をｃｏａｒｓｅ空間、それにＫ−直交な補空間をｆｉｎｅ空間と称する。ＣＧＣＧ法は、射影としてｆｉｎｅ空間へのＫ−直交射影を採用した共役射影勾配（Ｃｏｎｊｕｇａｔｅ  Ｐｒｏｊｅｃｔｅｄ  Ｇｒａｄｉｅｎｔ  ａｌｇｏｒｉｔｈｍ：以下射影ＣＧ法と称する）を並列処理で行なうものである。
背景技術
構造問題は、通常次のように扱われる。まず、対象となる構造物を連続体として定式化し、連続体の運動方程式（静的問題の場合は釣合い方程式）を定める。この方程式はほとんどの場合厳密に解くことが不可能であり、数値解析的に解かざるを得ない。そのための連続体問題の離散近似的再定式法が必要となる。その方法の一つとして有限要素法（日本機械学会  計算力学ハンドブック（Ｉ  有限要素法  構造編）日本機械学会１９９８）がある。有限要素法では、まず連続体の占める空間領域を複数の要素（有限要素法の要素）に分割し、非ゼロの値をとる領域が各要素に局在した関数（形状関数）を導入して、連続体変位場をそれらの関数の重ね合わせで表現されるものに限定する近似により、連続体の変位場およびその運動方程式（あるいは釣合い方程式）を離散化する。この離散化により、方程式は一つ（静的線形問題の場合）あるいは複数（例えば非線形問題の場合増分ステップ毎、動的問題の場合時間ステップ毎）の線形方程式に帰着される。有限要素法による連続体領域の要素分割の分割数により線形方程式が扱う自由度（次元）数が異なり、一般に近似をよくするために分割数を増やせば自由度数は増加し、対応する線形方程式の求解の難易度が高まる。
微小変形の静的構造問題は有限要素法においては式１のような、有限次元内積ベクトル空間Ｖ上の線形問題に帰着される。
【式１】
2002093412.tif 000003
ここで、Ｖは許容変位場（変位場の境界条件を満たすベクトル場）の空間、ＫはＶの次元ｄｉｍＶ次剛性マトリクス（正定値対称）、
2002093412.tif 000004
である。Ｖの次元ｄｉｍＶは問題の自由度数と等しい。Ｖの中から式１を満たすベクトルｕを探し出すことがソルバの課題になる。
先ず、本発明の基盤となる共役射影勾配法（ＣＰＧ：Ｃｏｎｊｕｇａｔｅ  Ｐｒｏｊｅｃｔｅｄ  Ｇｒａｄｉｅｎｔ  Ａｌｇｏｒｉｔｈｍ：Ｃ．Ｆａｒｈａｔ  Ｆ．−Ｘ．Ｒｏｕｘ：Ｉｍｐｌｉｃｉｔ  Ｐａｒａｌｌｅｌ  Ｐｒｏｃｅｓｓｉｎｇ  ｉｎ  Ｓｔｒｕｃｔｕｒａｌ  Ｍｅｃｈａｎｉｃｓ，Ｃｏｍｐｕｔａｔｉｏｎａｌ  Ｍｅｃｈａｎｉｃｓ  Ａｄｖａｎｃｅｓ  ２，１−１２４，１９９４）について説明する。一般に有限次元内積ベクトル空間Ｖ上の線型方程式式２
【式２】
2002093412.tif 000005
を考える。Ｋは正定値対称線形変換である。Ｖの部分空間Ｙを一つ選ぶ。ＶをＹへＫ−直交射影するＫ−直交射影子Ｐ（Ｙ），
2002093412.tif 000006
に条件Ｐ（Ｙ）＋Ｐ（ａ）＝１を満たすＫ−直交射影子Ｐ（ａ）（Ｐ（Ｙ）の補射影子），
2002093412.tif 000007
Ｙ，Ｖ（ａ）へ式３とＫ−直交直和分解される。
【式３】
2002093412.tif 000008
また、Ｐ（Ｙ）とＰ（ａ）のＫ−直交性をまとめて書けば式４なので、
【式４】
2002093412.tif 000009
線形方程式式２は
【式５】
2002093412.tif 000010
2002093412.tif 000011
いは  ｕのＹ−成分Ｐ（Ｙ）ｕとＶ（ａ）−成分Ｐ（ａ）ｕは独立に扱えることになる。つまり方程式式２は独立な２つの方程式式６、式７に式８とＫ−直交直和分離される。便宜のため、以下Ｙを直接法空間、ＹのＫ−直交補空間となる  Ｖ（ａ）を反復法空間と称する。
【式６】
2002093412.tif 000012
【式７】
2002093412.tif 000013
【式８】
2002093412.tif 000014
直接法空間Ｙの部分空間Ｗを特別に設定することを考える。このＷを以下ｃｏａｒｓｅ空間と称する。ＷのＫ−直交補空間をｆｉｎｅ空間と称する。ＹをＷとＷのＫ−直交補空間ＷｃにＫ−直交直和分解すれば、方程式式６と同様にＷ上の方程式式９
【式９】
2002093412.tif 000015
が設定される。ここでＰ（Ｗ）はＷへのＫ−直交射影子である。こ
2002093412.tif 000016
設定により定義されている場合、式９は具体的に式１０と定式化される。
【式１０】
2002093412.tif 000017
ここで、式１１と定義する。
【式１１】
2002093412.tif 000018
2002093412.tif 000019
ることからくる。
Ｋ−直交直和分解した方程式のうち式６を直接法で、式７を反復法で解くことを考える。方程式式７に反復法を適用すると第ｎステップ反復式が式１２、式１３と書ける。
【式１２】
2002093412.tif 000020
【式１３】
2002093412.tif 000021
2002093412.tif 000022
【式１４】
2002093412.tif 000023
【式１５】
2002093412.tif 000024
【式１６】
2002093412.tif 000025
2002093412.tif 000026
∈ＫＶ（ａ）に注意。このように、Ｋ−直交直和分解した方程式のうち片方について前処理付きＣＧ法で解く方法を共役射影勾配法（ＣＰＧ法）と称する。もう片方の方程式を直接法で解くことも含めて、便宜上以下、これを射影ＣＧ法と称する。
2002093412.tif 000027
のものを直接計算することは不可能だが、Ｐ（Ｗ＋ａ）は計算可能でＷの
2002093412.tif 000028
式１７と書けて、
【式１７】
2002093412.tif 000029
2002093412.tif 000030
ｇｒｉｄ問題式１８に帰着させる。
【式１８】
2002093412.tif 000031
ｃｏａｒｓｅ  ｇｒｉｄ問題は式１０の手順にしたがって解くことができる。
ＣＧＣＧ法の特徴づけのため、既存の有限要素法並列ソルバアルゴリズムＤＤＭ、ＢＤＤ法、並列ＣＧ法を説明する。以下に述べるとおり、領域分割に基づいた既存の有限要素法並列ソルバアルゴリズムＤＤＭ、ＢＤＤ法は、Ｋ−直交直和分解法に基づいた射影ＣＧ法の一種とみなせる。以下、有限要素法にしたがい定式化された構造問題の全自由度の空間をＶとし、解くべき線形方程式を式１９とする。
【式１９】
2002093412.tif 000032
ここで、Ｋは剛性マトリクスと呼ばれる。
超大規模な（自由度が１００万程度以上の）問題を効率的に（特に並列処理で）解くことを目的に考案された方法の一つとして領域分割法（ＤＤＭ：Ｄｏｍａｉｎ  Ｄｅｃｏｍｐｏｓｉｔｉｏｎ  Ｍｅｔｈｏｄ）がある。有限要素法により分割された各要素を互いに隣接するもの同士で適当にグループ分けする。各グループの占める空間領域を部分領域（ｓｕｂｄｏｍａｉｎ）と称する（第１図参照）。これにより、全領域がまず部分領域へ分割され、次に各部分領域の有限要素へ分割されるという具合に空間が階層的に分割離散化されることになる。これにしたがい、領域全体に関する処理が各部分領域ごとの処理と部分領域間の処理の２段階に分けられる。特に各部分領域ごとの処理は並列に行なうことができる。上記の処理の２段階分離は具体的には、領域全体の変位場を、部分領域内部に働く外力に応答する部分領域内部の変位場と、部分領域間の境界（内部境界と呼ぶ）上に働く外力に応答する残りの変位場に分けて別々に処理することにより実現している。後者の変位場はさらに部分領域内部の変位場と内部境界上の変位場に分けられ、前者は後者の従属変数になる。独立変数である内部境界の自由度についてはＣＧ法により解く。内部境界について解けば、部分領域内部の変位場はその境界上の変位を境界条件とし負荷ゼロの変位場として定まる。
2002093412.tif 000033
界Γ，部分領域内部全部の合併領域Ωｉおよび内部境界Γｓに分け
2002093412.tif 000034
散化の方法により変わる計算の便宜上設定されるものであり、物理
2002093412.tif 000035
の許容変位場の自由度の空間Ｖｉと内部境界Γｓ上の許容変位場の自由度の空間Ｖｓ直交し、全自由度空間Ｖは式２０のようにＶｉとＶｓに直交直和分解する。
【式２０】
2002093412.tif 000036
2002093412.tif 000037
とブロック分解されるが、これに基づき式２１とブロック対角化する。
【式２１】
2002093412.tif 000038
ＳはＳｃｈｕｒの補元と呼ばれる。射影子Ｐ（ｉ），Ｐ（ｓ）を式２２と定義する。
【式２２】
2002093412.tif 000039
Ｐ（ｓ）は内部境界Γｓ上の変位場を、それを境界条件とし外力を受け
2002093412.tif 000040
るから、
【式２３】
2002093412.tif 000041
【式２４】
2002093412.tif 000042
式２５となり、
【式２５】
2002093412.tif 000043
Ｐ（ｉ），Ｐ（ｓ）は条件Ｐ（ｉ）＋Ｐ（ｓ）＝１を満たすＫ−直交射影子である。特に、
2002093412.tif 000044
式２６と
【式２６】
2002093412.tif 000045
Ｋ−直交直和分解される。（ＶｉＶｓ）と（Ｖ（ｉ）Ｖ（ｓ））の関係は式２７である。
【式２７】
2002093412.tif 000046
2002093412.tif 000047
上の反力がゼロであるような変位場の空間であることを示している。一方第３式は、Ｖ（ｓ）がΓｓ上の変位場を幾何的境界条件として持つ
2002093412.tif 000048
ての２つの特徴づけは等価である。また、射影子Ｐ（ｓ）の固有値０に属する固有空間すなわちｋｅｒＰ（ｓ）についてｋｅｒＰ（ｓ）＝Ｖ（ｉ）＝ＶｉだからＶｓとＶ（ｓ）は線形同型である。したがってＰ（ｓ）のＶｓへの制限Ｖｓ→Ｖ（ｓ）は線形同型であることがわかる。このことは、Ｖ（ｓ）が変数空間である場合、代わりにＶｓを変数空間にとることができることを示している。
空間及び方程式のＫ−直交直和分解の直接法空間Ｙと反復法空間Ｖ（ａ）を式２８とおき、
【式２８】
2002093412.tif 000049
式１４〜式１６の前処理行列Ｇを単位行列Ｇ＝１ととって射影ＣＧ法（Ｊ．Ｍａｎｄｅｌ，Ｍ．Ｂｒｅｚｉｎａ：Ｂａｌａｎｃｉｎｇ  Ｄｏｍａｉｎ  Ｄｅｃｏｍｐｏｓｉｔｉｏｎ：Ｔｈｅｏｒｙ  ａｎｄ  Ｐｅｒｆｏｒｍａｎｃｅ  ｉｎ  ｔｗｏ  ａｎｄ  Ｔｈｒｅｅ  Ｄｉｍｅｎｓｉｏｎｓ，ＭＧＮｅｔ，ｈｔｔｐ：／／ｃａｓｐｅｒ．ｃｓ．ｙａｌｅ．ｅｄｕ／ｍｇｎｅｔ／ｗｗｗ／ｍｇｎｅｔ−ｐａｐｅｒｓ．ｈｔｍｌ）を適用したも
2002093412.tif 000050
自由度の空間、Ｖ（ｓ）はΓｓ上の許容変位場とそれを境界条件（正確
2002093412.tif 000051
由度の空間を表わす。ＫＶ（ｓ）はＶ（ｓ）の変位場に対応する内部境界上の反力の自由度の空間である。通常のＤＤＭアルゴリズムは上記にしたがいＶ（ｓ）の代わりにＶｓを変数空間にとる。
単位行列以外のＧを設定することにより前処理つきＤＤＭも可能
2002093412.tif 000052
とき前処理として有効であるといえる。
ＤＤＭのＣＧ法に用いる前処理として、Ｎｅｕｍａｎｎ前処理（Ｐ．Ｌｅ．Ｔａｌｌｅｃ：Ｄｏｍａｉｎ  ｄｅｃｏｍｐｏｓｉｔｉｏｎ  ｍｅｔｈｏｄｓ  ｉｎ  ｃｏｍｐｕｔａｔｉｏｎａｌｍｅｃｈａｎｉｃｓ，Ｃｏｍｐｕｔａｔｉｏｎａｌ  Ｍｅｃｈａｎｉｃｓ  Ａｄｖａｎｃｅｓ  １（２）（１９９４）１２１−２２０）がある。これは上記の前処理行列Ｇを単位行列１ではなくＮｅｕｍａｎｎ前処理行列にとったものである。Ｎｅｕｍａｎｎ前処理行列は各部分領域Ｉの局所剛性マトリクスＫＩのＳｃｈｕｒ補元ＳＩ（局
2002093412.tif 000053
するブロック分解表示で式２９で定義される。
【式２９】
2002093412.tif 000054
ここで、Ｉは領域分割による部分領域のインデクス、ＮＩは部分領域Ｉの自由度を領域全体の自由度へ写像する０−１成分行列、｛ＤＩ｝Ｉは
2002093412.tif 000055
い可能性がある。
2002093412.tif 000056
よって部分領域ごとの不定の剛体変位が反復処理ごとに混入する。これは部分領域ごとにばらばらの浮遊運動を生じさせ、反復収束効
2002093412.tif 000057
原因でもある。ＢＤＤ法（Ｂａｌａｎｃｉｎｇ  Ｄｏｍａｉｎ  Ｄｅｃｏｍｐｏｓｉｔｉｏｎ  ｍｅｔｈｏｄ，Ｊ．Ｍａｎｄｅｌ：Ｂａｌａｎｃｉｎｇ  Ｄｏｍａｉｎ  Ｄｅｃｏｍｐｏｓｉｔｉｏｎ，Ｃｏｍｍｕｎｉｃａｔｉｏｎｓ  ｏｎ  Ｎｕｍｅｒｉｃａｌ  Ｍｅｔｈｏｄｓ  ｉｎ  Ｅｎｇｉｎｅｅｒｉｎｇ  ９（１９９３）２３３−３４１．，Ｊ．Ｍａｎｄｅｌ，Ｍ．Ｂｒｅｚｉｎａ：Ｂａｌａｎｃｉｎｇ  Ｄｏｍａｉｎ  Ｄｅｃｏｍｐｏｓｉｔｉｏｎ：Ｔｈｅｏｒｙ  ａｎｄ  Ｐｅｒｆｏｒｍａｎｃｅ  ｉｎ  Ｔｗｏ  ａｎｄ  Ｔｈｒｅｅ  Ｄｉｍｅｎｓｉｏｎｓ，ＭＧＮｅｔ，ｈｔｔｐ：／／ｃａｓｐｅｒ．ｃｓ．ｙａｌｅ．ｅｄｕ／ｍｇｎｅｔ／ｗｗｗ／ｍｇｎｅｔ−ｐａｐｅｒｓ．ｈｔｍｌ，ＡＲＡＳＯＬ  Ａｎ  Ｉｎｔｅｇｒａｔｅｄ  Ｐｒｏｇｒａｍｍｉｎｇ  Ｅｎｖｉｒｏｎｍｅｎｔ  ｆｏｒ  Ｐａｒａｌｌｅｌ  Ｓｐａｒｓｅ  Ｍａｔｒｉｘ  Ｓｏｌｖｅｒｓ（Ｐｒｏｊｅｃｔ  Ｎｏ．２０１６０），Ｄｅｌｉｖｅｒａｂｌｅ  Ｄ  ２．４ｅ  Ｆｉｎａｌ  ｒｅｐｏｒｔ  Ｄｏｍａｉｎ  Ｄｅｃｏｍｐｏｓｉｔｉｏｎ  Ａｌｇｏｒｉｔｈｍｓ  ｆｏｒ  Ｌａｒｇｅ  Ｓｃａｌｅ  Ｉｎｄｕｓｔｒｉａｌ  Ｆｉｎｉｔｅ  Ｅｌｅｍｅｎｔ  Ｐｒｏｂｌｅｍｓ，Ｊｕｌｙ  ３０，１９９９．）はＤＤＭに基づき、内部境界上の変位について前処理つきＣＧ法を用いて解く解法であり、特に線形材料微小変形の静的問題に関しては、内部境界上に働く外力に応答する変位場を部分領域ごとの剛体変位と残りのひずみ変位とに分け、前者の自由度だけ先に直接法で解き、後者の自由度を次の前処理付きＣＧ法で解くものである。つまり、Ｎｅｕｍａｎｎ前処理に、先に解いた部分領域ごとの剛体変位の自由度の混入を排する射影を反復ごとに追加する。もう少し詳しく述べると、この射影は部分領域にでたらめな浮遊運動を生じさせるような外力を削除するという形で行なわれ、Ｋ−直交直和分解された方程式式７、式８でいうとＦ→Ｆ（ａ）の対応、射影ＣＧ法アルゴリズムでいえば残差の射影ｒ→ｒ（ａ）である。このように、Ｎｅｕｍａｎｎ前処理つきＤＤＭを改善して部分領域ごとの浮遊運動を押さえる上記の射影処理をｂａｌａｎｃｉｎｇと称する。
具体的には、ＤＤＭにおける部分空間Ｖ（ｓ）＝Ｐ（ｓ）Ｖｓをさらに次のようにＫ−直交直和分解する。ｋｅｒＳＩを含む部分空間のＩにわたる和空間、例えば部分領域ごとの剛体変位を部分領域全部について重ね合わせた運動を考え、その自由度の空間をｃｏａｒｓｅ空間Ｗとする。ＷはＶ（ｓ）の部分空間である。Ｖ（ｓ）をＷとそのＶ（ｓ）でのＫ−直交補空間Ｖ（ｔ）にＫ−直交直和分解する。対応するＫ−直交射影子のペアをＰ（Ｗ），Ｐ（ｔ）とする。式３０が成立する。
【式３０】
2002093412.tif 000058
全自由度空間Ｖは式３１、式３２とＫ−直交直和分解される。
【式３１】
2002093412.tif 000059
【式３２】
2002093412.tif 000060
空間及び方程式のＫ−直交直和分解の直接法空間と反復法空間を式３３とし、
【式３３】
2002093412.tif 000061
式１４〜式１６の前処理行列ＧをＮｅｕｍａｎｎ前処理行列にとって射影ＣＧ法（Ｊ．Ｍａｎｄｅｌ，Ｍ．Ｂｒｅｚｉｎａ：Ｂａｌａｎｃｉｎｇ  Ｄｏｍａｉｎ  Ｄｅｃｏｍｐｏｓｉｔｉｏｎ：Ｔｈｅｏｒｙ  ａｎｄ  Ｐｅｒｆｏｒｍａｎｃｅ  ｉｎ  Ｔｗｏ  ａｎｄ  Ｔｈｒｅｅ  Ｄｉｍｅｎｓｉｏｎｓ，ＭＧＮｅｔ，ｈｔｔｐ：／／ｃａｓｐｅｒ．ｃｓ．ｙａｌｅ．ｅｄｕ／ｍｇｎｅｔ／ｗｗｗ／ｍｇｎｅｔ−ｐａｐｅｒｓ．ｈｔｍｌ）を適用したものがＢＤＤ法である。特に、式１７、式１８に相当する前処理はＰ（Ｗ＋ａ）＝Ｐ（ｓ）であることから式３４、式３５となる。
【式３４】
2002093412.tif 000062
【式３５】
2002093412.tif 000063
2002093412.tif 000064
ＢＤＤ法の反復法空間Ｖ（ｔ）は同時に、後に局所ｃｏａｒｓｅ空間のＷＩの定義のところで一緒に定義されるｂａｌａｎｃｅｄ空間でもある。
2002093412.tif 000065
Ｄｅｃｏｍｐｏｓｉｔｉｏｎ，Ｃｏｍｍｕｎｉｃａｔｉｏｎｓ  ｏｎ  Ｎｕｍｅｒｉｃａｌ  Ｍｅｔｈｏｄｓ  ｉｎ  Ｅｎｇｉｎｅｅｒｉｎｇ  ９（１９９３）２３３−３４１．，Ｊ．Ｍａｎｄｅｌ，Ｍ．Ｂｒｅｚｉｎａ：Ｂａｌａｎｃｉｎｇ  Ｄｏｍａｉｎ  Ｄｅｃｏｍｐｏｓｉｔｉｏｎ：Ｔｈｅｏｒｙ  ａｎｄ  Ｐｅｒｆｏｒｍａｎｃｅ  ｉｎ  Ｔｗｏ  ａｎｄ  Ｔｈｒｅｅ  Ｄｉｍｅｎｓｉｏｎｓ，ＭＧＮｅｔ，ｈｔｔｐ：／／ｃａｓｐｅｒ．ｃｓ．ｙａｌｅ．ｅｄｕ／ｍｇｎｅｔ／ｗｗｗ／ｍｇｎｅｔ−ｐａｐｅｒｓ．ｈｔｍｌ」。ＢａｌａｎｃｉｎｇはＢＤＤ法に適用した射影ＣＧ法のアルゴリズム式１４〜式１６には直接現われず、Ｋ−
2002093412.tif 000066
る。Ｎｅｕｍａｎｎ前処理行列式２９はこの条件を満たし、本節冒頭で提示したＮｅｕｍａｎｎ前処理つきＤＤＭの問題点が解決する。
並列ＣＧ法は領域分割を行なわず、全自由度空間ＶをそのままＣＧ法で処理する。ＶのＫ−直交直和分解も行なわない。直接法空間をＹ＝｛０｝ととることに相当する。問題が大規模で、ベクトル空間Ｖの次元が大きい場合は、解析対象の空間領域を幾つかに分割し（これを部分（ｐａｒｔ）と称する）、それに伴なうＶの分解（境界の重複した分解）にしたがって各々の部分空間についての自由度を別々のプロセッサが処理する。ベクトルの内積や行列−ベクトル積のように部分空間相互の情報交換が必要な計算を行なうときに限りプロセス間通信が必要になる。
以上、超大規模構造問題を解く方法として、領域分割法（ＤＤＭ：Ｄｏｍａｉｎ  Ｄｅｃｏｍｐｏｓｉｔｉｏｎ  Ｍｅｔｈｏｄ）、Ｎｅｕｍａｎｎ前処理法つきＤＤＭ、ＢＤＤ法（Ｂａｌａｎｃｉｎｇ  Ｄｏｍａｉｎ  Ｄｅｃｏｍｐｏｓｉｔｉｏｎ  Ｍｅｔｈｏｄ）、並列ＣＧ法の解法について説明した。しかし、これらの方法では、超大規模構造問題を解く場合、発散して解が求まらない、計算時間が極めて長くなるという問題が指摘されている。本発明は、このような状況に鑑みてなされたもので、自由度が１００万以上の超大規模構造問題を解く際、解が発散することなく求まり、反復計算の回数が少なく、短時間の計算時間で解が求まる、超大規模構造問題を解くシステムＣＧＣＧ法を提供しようとするものである。
発明の開示
本発明は、超大規模構造問題を解くシステム、該システムを機能させるためのプログラム及び該プログラムを記録したコンピュータ読み取り可能な記録媒体を含むものである。ＣＧＣＧ法は、領域分割を行ない、それに基づいて部分領域の内部と境界の区別をすることなしにｃｏａｒｓｅ空間を設定し、前処理として簡便な対角スケーリングを採用した前処理つき射影ＣＧ法を用いる、超大規模構造問題のための有限要素法ソルバアルゴリズムである。
本発明は、自由度１００万以上の超大規模構造問題を解く並列有限要素法計算システムであって、領域分割を行う手段、部分領域を各プロセッサーの担当する部分へ分配する手段、剛性マトリックスを作成する手段、部分領域全体重ね合わせ運動を設定する手段、全自由度の前処理つき射影ＣＧ法の初期設定を行う手段、全自由度の前処理つき射影ＣＧ法の反復計算をする手段及び変位解を出力する手段を備えたことを特徴とする並列有限要素法計算システムである。
また、前記部分領域全体重ね合わせ運動を設定する手段が、全自由度表示の射影子を作成する手段、部分領域全体重ね合わせ運動マトリックスを作成する手段及び部分領域全体重ね合わせ運動マトリックスをＬＵ分解する手段を備えることができる。
前記全自由度の前処理つき射影ＣＧ法の初期設定を行う手段が、全自由度の初期変位を設定する手段、全自由度の初期残差を計算する手段、対角スケーリング前処理の計算をする手段、全自由度のＣＯＡＲＳＥ  ＧＲＩＤ前処理の計算をする手段及び全自由度のＣＧ法探索方向のベクトル初期値を設定する手段を備えることができる。
そして、前記全自由度の前処理つき射影ＣＧ法の反復計算をする手段が、全自由度の変位を更新する手段、全自由度の残差を更新する手段、対角スケーリング前処理の計算をする手段、全自由度のＣＯＡＲＳＥ  ＧＲＩＤ前処理の計算をする手段、全自由度のＣＧ法探索方向のベクトルを更新する手段及び収束を判定する手段を備えることができる。
本発明は、上記システムを機能させるプログラムでもある。即ち、自由度１００万以上の超大規模構造問題を解く並列有限要素法計算システムであって、領域分割を行う手段、部分領域を各プロセッサーの担当する部分へ分配する手段、剛性マトリックスを作成する手段、部分領域全体重ね合わせ運動を設定する手段、全自由度の前処理つき射影ＣＧ法の初期設定を行う手段、全自由度の前処理つき射影ＣＧ法の反復計算をする手段及び変位解を出力する手段として機能させるための並列有限要素法計算プログラムである。そして、前記部分領域全体重ね合わせ運動を設定する手段が、全自由度表示の射影子を作成する手段、部分領域全体重ね合わせ運動マトリックスを作成する手段及び部分領域全体重ね合わせ運動マトリックスをＬＵ分解する手段として機能させることができ、前記全自由度の前処理つき射影ＣＧ法の初期設定を行う手段が、全自由度の初期変位を設定する手段、全自由度の初期残差を計算する手段、対角スケーリング前処理の計算をする手段、全自由度のＣＯＡＲＳＥ  ＧＲＩＤ前処理の計算をする手段及び全自由度のＣＧ法探索方向のベクトル初期値を設定する手段として機能させることができ、前記全自由度の前処理つき射影ＣＧ法の反復計算をする手段が、全自由度の変位を更新する手段、全自由度の残差を更新する手段、対角スケーリング前処理の計算をする手段、全自由度のＣＯＡＲＳＥ  ＧＲＩＤ前処理の計算をする手段、全自由度のＣＧ法探索方向のベクトルを更新する手段及び収束を判定する手段として機能させることができる。
また、本発明は、上記システムを機能させるプログラムを記録したコンピュータ読み取り可能な記録媒体でもある。即ち、
自由度１００万以上の超大規模構造問題を解く並列有限要素法計算システムであって、領域分割を行う手段、部分領域を各プロセッサーの担当する部分へ分配する手段、剛性マトリックスを作成する手段、部分領域全体重ね合わせ運動を設定する手段、全自由度の前処理つき射影ＣＧ法の初期設定を行う手段、全自由度の前処理つき射影ＣＧ法の反復計算をする手段及び変位解を出力する手段として機能させるための並列有限要素法計算プログラムを記録したコンピュータ読み取り可能な記録媒体であり、前記部分領域全体重ね合わせ運動を設定する手段が、全自由度表示の射影子を作成する手段、部分領域全体重ね合わせ運動マトリックスを作成する手段及び部分領域全体重ね合わせ運動マトリックスをＬＵ分解する手段として機能させることができ、前記全自由度の前処理つき射影ＣＧ法の初期設定を行う手段が、全自由度の初期変位を設定する手段、全自由度の初期残差を計算する手段、対角スケーリング前処理の計算をする手段、全自由度のＣＯＡＲＳＥ  ＧＲＩＤ前処理の計算をする手段及び全自由度のＣＧ法探索方向のベクトル初期値を設定する手段として機能させることができ、前記全自由度の前処理つき射影ＣＧ法の反復計算をする手段が、全自由度の変位を更新する手段、全自由度の残差を更新する手段、対角スケーリング前処理の計算をする手段、全自由度のＣＯＡＲＳＥ  ＧＲＩＤ前処理の計算をする手段、全自由度のＣＧ法探索方向のベクトルを更新する手段及び収束を判定する手段として機能させることができる。
発明を実施するための最良の形態
ＣＧＣＧ法は、領域分割を行ない、それに基づいて部分領域の内部と境界の区別をすることなしにｃｏａｒｓｅ空間を設定し、前処理として簡便な対角スケーリングを採用した前処理つき射影ＣＧ法を用いる、超大規模構造問題のための有限要素法ソルバアルゴリズムである。ＤＤＭ、ＢＤＤ法と同様、ＣＧＣＧ法は領域分割に基づいた射影ＣＧ法を基盤とする。まずＹ＝Ｗととり、射影ＣＧ法で式３に倣って全自由度空間Ｖを式３６
【式３６】
2002093412.tif 000067
とＫ−直交分解する。この際、ｃｏａｒｓｅ空間ＷとしてＶのどの部分空間を選ぶかが問題になる。ＣＧＣＧ法では、Ｗとして領域分割における部分領域ごとのある限定された運動（例えば剛体運動）を表わすＶの部分空間を選ぶ。
例えばＷしてＢＤＤ法と同様に、部分領域ごとの剛体変位を部分領域全部について重ね合わせた運動の自由度の空間を採用した場合、Ｖ（ａ）は部分領域内部のひずみ運動を表わす変位場の空間を成す。つまりＢＤＤ法の設定式３３と比較して直接法空間と反復法空間を式３７ととることになる。
【式３７】
2002093412.tif 000068
この場合、ＢＤＤ法がＶ（ｉ）の自由度については直接法で解くのに対し、ＣＧＣＧ法ではそれも含めてＣＧ法で解く。この意味で、ＣＧＣＧ法の方がＢＤＤ法よりＣＧ法の負担が重く直接法の負担が軽い。射影ＣＧ法前処理でもＣＧＣＧ法ではＢＤＤ法における式３４の射影
2002093412.tif 000069
Ｃｏａｒｓｅ空間Ｗの実装面での取扱いもＤＤＭ法と並列ＣＧ法の相違を反映して次のように異なってくる。ＤＤＭ法では、Ｖｓを変数空間にとるため、領域全体の変位場を部分領域内の変位と内部境界上の変位に分けて扱う。並列ＣＧ法では、部分領域をその内部と内部境界とに分けることをせず、部分領域はその内部と境界を一体として扱う。ＢＤＤ法でのｃｏａｒｓｅ空間Ｗは次のように設定される。各部分領域の境界（内部境界の一部）の剛体変位を先に定式化し、その部分領域全部にわたる重ね合わせの空間Ｗｓを設定する。それを式２２のＰ（ｓ）でＫ−直交射影して部分領域内の剛体変位に拡張してＷとする（Ｗ＝Ｐ（ｓ）Ｗｓ）。一方、ＣＧＣＧ法でのｃｏａｒｓｅ空間Ｗは、各部分領域ごとにその内部と境界全体の剛体変位を直接定式化してその部分領域全部にわたる重ね合わせで設定する。これにしたがい、射影ＣＧ法の前処理では、空間Ｗｓの自由度の求解のほかに部分領域ごとの直接法処理の手間が掛かっている（Ｓｈｕｒ補元の計算のため）ＢＤＤ法に対し、ＣＧＣＧ法では、その内部も含めた部分領域ごとの剛体変位の自由度について解くだけでよい。
方程式式１を式６〜式８に倣って式３８、式３９、式４０
【式３８】
2002093412.tif 000070
【式３９】
2002093412.tif 000071
【式４０】
2002093412.tif 000072
と分解する。式３８は変位場Ｗ−成分ｕ（Ｗ）を決定するｃｏａｒｓｅ  ｇｒｉｄ問題、式３９は変位場Ｖ（ａ）−成分ｕ（ａ）を決定する方程式である。ＣＧＣＧ法では射影ＣＧ法に倣って式３８を直接法である修正コレスキー法、式３９を射影ＣＧ法の前処理つきＣＧ反復処理で解く。前処理つき
2002093412.tif 000073
４１
【式４１】
2002093412.tif 000074
2002093412.tif 000075
Ｋの対角スケーリングと反復法空間Ｖ（ａ）へのＫ−直交射影Ｐ（ａ）の合成である。
ＢＤＤ法ではＮｅｕｍａｎｎ前処理という、各部分領域で局所Ｓｃｈｕｒ
2002093412.tif 000076
に対し、ＣＧＣＧ法では対角スケーリングという、簡便な前処理で計算コストとメモリ使用量を大幅に削減している。
ＣＧＣＧ法でのＫ−直交分解式３６におけるｃｏａｒｓｅ空間Ｗの定義を述べる。Ｉを領域分割による部分領域のインデクスとする。部分領域Ｉ内の自由度が成すベクトル空間をＶＩとする。ＶＩの部分空間ＷＩを定める。これを局所ｃｏａｒｓｅ空間と呼ぶことにする。ｃｏａｒｓｅ空間Ｗを｛ＷＩ｝Ｉの重み付き重ね合わせで式４２を定義する。
【式４２】
2002093412.tif 000077
ここで、｛ＤＩ｝ＩはＮｅｕｍａｎｎ前処理と同様１の各部分領域への分解
2002093412.tif 000078
体運動）を表わす変位場の空間である。ｍＩ≡ｄｉｍＷＩをＷＩの次元、
2002093412.tif 000079
【式４３】
2002093412.tif 000080
2002093412.tif 000081
特に、ＷＩが条件ｋｅｒＳＩ⊂ＷＩを満たすとき、対応するｆｉｎｅ空間の部分空間をｂａｌａｎｃｅｄ空間と称する。ＢＤＤ法では、Ｎｅｕｍａｎｎ前処理を行なうため、ＷＩはｋｅｒＳＩ⊂ＷＩを満たす必要があるがＣＧＣＧ法ではそのような制約がない。すなわちＣＧＣＧ法の反復法空間Ｖ（ａ）はｂａｌａｎｃｅｄ空間である必要がない。
部分領域Ｉの規定された運動として剛体運動を考える。このときＷＩは具体的には次のように構成できる。部分領域Ｉの剛体運動は次の
2002093412.tif 000082
変形後の座標として、式４４、式１５、式４５となる。
【式４４】
2002093412.tif 000083
【式４５】
2002093412.tif 000084
【式４６】
2002093412.tif 000085
微小変形を考えると、式４７、式４８である。
【式４７】
2002093412.tif 000086
【式４８】
2002093412.tif 000087
そこでＷＩを、次の条件で定義される部分領域Ｉ上の６つの変位場
2002093412.tif 000088
４９、式５０である。
【式４９】
2002093412.tif 000089
【式５０】
2002093412.tif 000090
ただし、式５０では次の事情を考慮している。部分領域によっては、
2002093412.tif 000091
にする。番号ｊの振り直しにはＧｒａｍ−Ｓｃｈｍｉｄｔの直交化法を用いる。
2002093412.tif 000092
ではなく、その内部境界上の６つの変位場を定めるところがＣＧＣＧ法と異なる。
上記の規定された運動としての剛体運動の定式を微小変形問題ソリッドモデルの剛体運動へ適用することは容易である。ソリッドモデルの場合、式４９で定義される部分領域Ｉ上の６つの変位場
2002093412.tif 000093
【式５１】
2002093412.tif 000094
2002093412.tif 000095
つぎに、微小変形問題シェルモデルへの適用を述べる。節点αに
2002093412.tif 000096
点の座標ｘ＝ｘ（ξ１，ξ２，ξ３）は式５２と書ける。
【式５２】
2002093412.tif 000097
ここで（ξ１，ξ２，ξ３）は当該シェル要素の局所座標である。またＮα，Ｘα，Ｕα，φαは節点αのそれぞれシェル中立面上で定義された形状関数、
2002093412.tif 000098
当該シェル要素に含まれる節点にわたる和を示す。
まず、シェルの微小平行移動ｘ→ｘ＋ａは式５３である。
【式５３】
2002093412.tif 000099
つまり節点変位Ｕαと現ディレクタの回転量φαの変化でまとめると式５４である。
【式５４】
2002093412.tif 000100
次にシェルの微小剛体回転を考える。変位Ｕαおよび現ディレクタの回転φαは微小とする。座標ｘ＝ｘ（ξ１，ξ２，ξ３）に微小剛体回転１＋Ｏｉθｉ＝１＋θ×を作用させる。ここで、｛Ｏｉ｝ｉ＝１，２，３は式４６で定義されたマトリクス、θ＝（θ１θ２θ３）Ｔは回転方向ベクトルである。微小剛体
2002093412.tif 000101
は微小量一次まで考えて式５５である。
【式５５】
2002093412.tif 000102
2002093412.tif 000103
張られる部分空間）に射影したベクトル式５６である。
【式５６】
2002093412.tif 000104
従って、微小剛体回転による座標ｘ全体の変化量は式５７となる。
【式５７】
2002093412.tif 000105
従って、シェルの微小剛体回転は節点変位Ｕαと現ディレクタの回転量φαの変化でまとめると式５８である。
【式５８】
2002093412.tif 000106
微小平行移動ｘ→ｘ＋ａと微小剛体回転ｘ→（１＋θ×）ｘを合わせると式５９である。
【式５９】
2002093412.tif 000107
2002093412.tif 000108
【式６０】
2002093412.tif 000109
となる。すなわち、式４９で定義される部分領域Ｉ上の６つの変位
2002093412.tif 000110
【式６１】
2002093412.tif 000111
2002093412.tif 000112
部分領域Ｉの規定された運動としてアフィン変換（直線を直線に写し、その上の有向線分比を変えないと特徴づけられる変換）で表現される運動（以下、アフィン変換運動と呼ぶことにする）を考える。つまり、前節で扱った剛体運動を一般化したものであり、並進運動と一般線形変換運動の合成である。この場合ＷＩは具体的に次のように構成できる。
2002093412.tif 000113
【式６２】
2002093412.tif 000114
【式６３】
2002093412.tif 000115
【式６４】
2002093412.tif 000116
あるいは、式６５、式６６、式６７、式６８である。
【式６５】
2002093412.tif 000117
【式６６】
2002093412.tif 000118
【式６７】
2002093412.tif 000119
【式６８】
2002093412.tif 000120
微小変形を考えて、式６９、式７０である。
【式６９】
2002093412.tif 000121
【式７０】
2002093412.tif 000122
そこでＷＩを、次の条件で定義される部分領域Ｉ上の１２の変位場
2002093412.tif 000123
式７１で
【式７１】
2002093412.tif 000124
式７２である。
【式７２】
2002093412.tif 000125
剛体運動の場合と同様、式７２では次の事情を考慮している。部分領域によっては、そのアフィン変換運動の自由度が１１以下の場合
2002093412.tif 000126
ｍＩ≡ｄｉｍＷＩ＜１２となるようにする。番号ｊの振り直しにはＧｒａｍ−Ｓｃｈｍｉｄｔの直交化法を用いる。
次に、ＣＧＣＧ法による並列有限要素法ソルバの具体的な計算手順を述べる。
１．領域分割
解析領域全体を複数の部分領域に分割する。
２．部分領域の各プロセッサの担当する部分への分配
分割した部分領域の各計算ノード（ＣＰＵ）への割り当てを決める。すなわち、各計算ノード（ＣＰＵ）は、１つまたは複数の部分領域を担当する。各ノードが担当する部分領域全体を部分と称する（第１図参照）。
３．剛性マトリックスの作成
各ノードにおいて、受け持ちの部分領域に対する剛性マトリックスを作成する。
４．部分領域全体重ね合わせ運動設定
４．１  全自由度表示の射影子作成
部分領域全体重ね合わせ運動を抜き出すための全自由度表示のｃｏａｒｓｅ空間の基底（部分空間の基底はその部分空間への射影子とみ
2002093412.tif 000127
する。
４．２  部分領域全体重ね合わせ運動マトリックスＫｒｇｄの作成
2002093412.tif 000128
を作成する。この計算は、各部分領域の剛体マトリクスＫＩを保持するノードにてその部分領域からの寄与を求め、それを全体にわたって足しあわせることで求められる。
４．３  部分領域全体重ね合わせ運動マトリックスＫｒｇｄのＬＵ分解
Ｋｒｇｄに対し、ＬＵ分解Ｋｒｇｄ＝ＬｒｇｄＵｒｇｄを行なう。このＬＵ分解されたＬｒｇｄ，Ｕｒｇｄをすべてのノードで保持しておく。
５．全自由度の前処理つき射影ＣＧ法初期設定
５．１  全自由度の初期変位設定
５．１．１  部分領域全体重ね合わせ運動μ（Ｗ）
次の方程式式７３をμ（Ｗ）について解く。
【式７３】
2002093412.tif 000129
2002093412.tif 000130
に関して保持しており、各ノードで（ＮＤＺ）ＴＦを求めた後通信して領域全体にわたる右辺ベクトルを作成し、すべてのノードでＬＵ分解したＫｒｇｄを用いてμ（Ｗ）を求める。式１０のｕ（Ｗ）に相当するものはｕ（Ｗ）＝ＮＤＺμ（Ｗ）である。
５．１．２  初期変位ｕ０の設定
初期変位ｕ０を式７４で設定する。
【式７４】
2002093412.tif 000131
このｕ０は各ノードの担当領域部分に対してのみ保持する。
５．２  全自由度の初期残差ｇ０の計算
初期残差ｇ０を式７５で求める。
【式７５】
2002093412.tif 000132
この残差ベクトルは各ノードの担当分のみ保持する。式１４〜式
2002093412.tif 000133
５．３  対角スケーリング前処理計算
2002093412.tif 000134
分と等しい対角行列）である。各ノードにて担当部分について求める。
５．４  全自由度のＣｏａｒｓｅ  Ｇｒｉｄ前処理計算
2002093412.tif 000135
５．４．１  部分領域全体重ね合わせ運動変数μ（Ｗ）
次の方程式式７６をμ（Ｗ）について解く。
【式７６】
2002093412.tif 000136
右辺ベクトルは、通信してすべてのノードで全体領域にわたって作成し、すべてのノードでμ（Ｗ）を求める。
2002093412.tif 000137
2002093412.tif 000138
【式７７】
2002093412.tif 000139
2002093412.tif 000140
2002093412.tif 000141
【式７８】
2002093412.tif 000142
担当部分に関して求める。
５．５  全自由度のＣＧ法探索方向ベクトル初期値ｗ０の設定
全自由度のＣＧ法探索方向ベクトル初期値ｗ０を式７９で設定する。
【式７９】
2002093412.tif 000143
担当部分に関して求める。
６．全自由度の前処理つき射影ＣＧ法反復計算
６．１  全自由度の変位更新
ｎ＞１に対して、以下の手順にて全自由度の変位ｕｎ−１を更新する。
６．１．１  Ｋｗｎ−１の計算
Ｋｗｎ−１を担当部分に関して求める。
６．１．２  αｎの計算
通信を行ない、全ノードにてαｎを式８０のとおり求める。
【式８０】
2002093412.tif 000144
６．１．３  ｕｎ−１の担当部分の更新
式８１によりｕｎ−１の担当部分を更新する。
【式８１】
2002093412.tif 000145
６．２  全自由度の残差更新
ｎ＞１に対して、全自由度の残差ｇｎ−１の担当部分を式８２に従って更新する。
【式８２】
2002093412.tif 000146
2002093412.tif 000147
して求める。（対角スケーリング前処理）。
2002093412.tif 000148
ち２つの方程式式８３、式８４を解き式８５とする。
【式８３】
2002093412.tif 000149
【式８４】
2002093412.tif 000150
【式８５】
2002093412.tif 000151
６．５  全自由度のＣＧ法探索方向ベクトルｗｎ−１の更新
次のとおりｗｎ−１を更新する。
６．５．１  βｎの計算
βｎを式８６のとおり、通信して全ノードにて求める。
【式８６】
2002093412.tif 000152
６．５．２  ｗｎ−１の更新
ｗｎ−１を式８７のとおり担当部分ごとに更新する。
【式８７】
2002093412.tif 000153
６．６  収束判定
更新した残差ｇｎより収束を判定する。収束していない場合は、６．１に戻りそれ以降を再度実行する。
７．変位解ｕの出力
前処理つきＣＧ法収束時点でのｕｎを変位解ｕとする。必要ならｕより歪みや応力を求める。この計算手順は、第２図に示した。
実施例
以下実施例に基づいて、ＣＧＣＧ法が、ＤＤＭ、並列ＣＧ法、ＢＤＤと比較して、計算パフォーマンスに優れていることを説明する。鉄道台車をモデルとする有限要素法における計算パフォーマンスを比較する。このモデルでは、要素が４面体１次要素、節点数が３２３，６３９、要素数が１，１２３，８３６、自由度数が９７０，９１１（６自由度のみ拘束）であり、計算環境は、計算機Ｄｕａｌ  Ｐｅｎｔｉｕｍ  ＩＩＩ  ６００ＭＨｚ４台、ＰＥ数４（各マシンに１プロセス）、計算条件はＣＧ法トレランス１．０×１０−６であった。計算パーフォーマンスの比較を第１表に示した。第１表から、次のことがわかる。即ちＣＧＣＧ法はＤＤＭ、並列ＣＧ法に比較して、反復ステップ数、計算時間共に極めて少ない。ＣＧＣＧ法は反復ステップ数に関してはＢＤＤ法に劣るが、計算速度に関してはＢＤＤ法の３倍近く早い。
2002093412.tif 000154
比較例１
第１表には、Ｎｅｕｍａｎｎ前処理つきＤＤＭとの比較は示されていないが、ＤＤＭ、及びＮｅｕｍａｎｎ前処理つきＤＤＭの二者の計算パーフォーマンスを比較した計算比較例がある。比較例の計算パーフォンマンスを第２表に示した。この比較例は、複数の直方体を組み合わせて作った簡易モデルにおいて、要素が４面体２次要素、節点数が１，０２９、要素数が５０４、自由度数が３，０８７（底面を完全固定、上面に強制変位）であり、計算環境は計算機Ａｌｐｈａ２１１６４  ６００ＭＨｚ１台、ＰＥ数は１、計算条件はＣＧ法でトレランスが１．０×１０−７であった。第２表は、ＤＤＭとＮｅｕｍａｎｎ前処理つきＤＤＭの計算時間と反復ステップ数を示している。Ｎｅｕｍａｎｎ前処理つきＤＤＭで短縮された時間は１３％程度、反復ステップ数もせいぜい１／２程度である。この比較例１と実施例から、本発明のＣＧＣＧ法は、超大規模構造問題の解法として従来知られているＤＤＭ、並列ＣＧ法、Ｎｅｕｍａｎｎ前処理つきＤＤＭに比較して、計算時間が極めて短く、高速化及びロバスト性向上を目的として開発されたＢＤＤ法と比べてもより速い。反復ステップ数に関しては、ＢＤＤ法には劣るが、ＤＤＭ、並列ＣＧ法及びＮｅｕｍａｎｎ前処理付きＤＤＭとの比較では極めて少なくなることがわかる。このように、ＣＧＣＧ法は、超構造問題の解法として計算パーフォーマンスが極めて高いことがわかる。
2002093412.tif 000155
産業上の利用可能性
ＣＧＣＧ法は、大規模な問題を有効に解く有限要素法ソルバアルゴリズムであり、有限要素法に基づく解析一般に適用可能である。たとえば、有限要素法によって定式化される連続体力学全般、固体力学、構造力学、流体力学、熱伝導問題、電磁場問題に利用できる。また、有限要素法は微分方程式を支配方程式とする境界値問題の一般的解法とみなされるので、その意味で、ＣＧＣＧ法は、微分方程式を支配方程式とする境界値問題全般の解法として適用可能である。どの問題においても、ＣＧＣＧ法の適用は大規模な問題を効率よく解くことに有効である。ＣＧＣＧ法は、大規模な問題を解く有限要素法ソルバアルゴリズムとして従来から知られている、ＤＤＭ、並列ＣＧ法及びＮｅｕｍａｎｎ前処理付きＤＤＭに比較して計算時間が極めて短くなり、計算反復ステップが極めて少なくなるという特徴を有する。
更に、ＣＧＣＧ法は計算時間に関して、高速化を目的として開発されたＢＤＤ法に比較してもより高い性能を示す。
【図面の簡単な説明】
第１図は領域分割を説明する図であり、第２図はＣＧＣＧ法の計算フローチャートである。
【図１】
2002093412.tif 000156
【図２】
2002093412.tif 000157
【国際調査報告】
2002093412.tif 000158
2002093412.tif 000159
2002093412.tif 000160
2002093412.tif 000161