収束性悪化の原因箇所の推定方法

[OKUYAMA]

chtMultiRegionSimpleFoamを使って円管内を流体が流れるモデルの定常解析をしていますが、収束性が思わしくありません。

条件やメッシュの切り方により、イタレーション数が多いため時間がかかったり発散したりするケースと、イタレーション数は多くないが初期残差がなかなか減少しないケースがあります。

メッシュ品質(サイズ・アスペクト比)を疑っていますが、特におかしなところはなく、全体的に細かくしても改善するとは限らない(むしろ悪化することもある)ため、原因箇所の特定に苦慮しています。

そこで、残差や連続の式の誤差の分布をみて、収束性の悪化の原因となっている場所の見当をつけたいと考えているのですが、このようなことは可能でしょうか。

(現在使っている構造解析ソフトにはこのような機能があるため、OpenFOAMにもないかと思い、質問させていただきました。)

また、他に原因箇所の推定方法があれば、ご教示いただけますでしょうか。

[Inaba]

こんにちは。
稲葉と申します。


残念ながら残渣の分布出力機能は現時点では無いと思います。（3.0.1と4.1で私の見た範囲では）

ただ参考までに、並列計算を行っている場合に限りますが
私が良く分からない発散で躓いたときはよく停止したプロセッサの番号を見ています。
例えばおなじみのゼロ割りエラーの場合には、
[2] #0 Foam::error::printStack(Foam::Ostream&) at ??:?
[2] #1 Foam::sigFpe::sigHandler(int) at ??:?
...
の[2]の部分でprocessor2で該当エラーが起きたことが分かります。

一方でdecomposePar -celldistを行うことで、どの部位をどのプロセッサが担当しているかが分かります。
分割の仕方を何回か変えて計算してみることで、エラーが起きる場所が大体どの辺りかを推測できるかと思います。
最終的にエラーになってしまうのであればこのような確認の仕方もありますよということで。


ただこれだとエラーで計算終了しない限り分かりません。
おっしゃるようにやっぱりparaViewで経時的に残渣を見てみたいものです。
私も前々からそんな機能が欲しいと思っていたので少し調べてみました。

残渣の計算は行列計算の中でそれぞれ行っているようです。

例えばsmoothSolverの場合に計算中に出力される残渣は
src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.C　の中のsolve関数の中で
solverPerf.initialResidual() = cmptDivide(gSumCmptMag(this->matrix_.source() - Apsi),normFactor);
みたいに計算しているようです。
"cmptDivide(gSumCmptMag"はベクトルの軸ごとに絶対値を積分してそれぞれの積分値を出力、みたいな感じでしょうか。

つまり、
this->matrix_.source() - Apsi
の値をフィールドに書き出せば残渣の分布を出力できそうです。

というわけで、例えばsmoothSolverで計算するUの残渣を見たい場合には
まずソルバーの中で
volVectorField ridualFieldU
(
    IOobject
    (
        "residual_U",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    mesh,
    dimensionedVector("residualVectorU", dimless, 0.0)
);
みたいに書き出せるフィールドを作成し、

SmoothSolver.Cの該当関数内の残渣計算式の直前か直後あたりで
const fvMesh& mesh = this->matrix_.mesh(); // 戻り値はlduMesh型ですがfvMeshを継承しているのできっと大丈夫
if(mesh.foundObject<volVectorField>("residual_U")) {
    volVectorField& resU = const_cast<volVectorField&> (
        mesh.lookupObject<volVectorField>("residual_U")
    );
    vectorField& priResU = resU.primitiveFieldRef(); //3.0.1以前は確かinternalField()
    priResU = this->matrix_.source() - Apsi;
}

みたいに値をフィールドに移してあげればよさそうです。
ただ全くの未確認です。（多分どこかでProtectedになりそう）
時間があったら実装してみてまた報告します。

以上です。
直接的な解決にはなっていませんが、ご参考になれば幸いです。
（こんなに長々と書いて出力機能ありますよとかだったら恥ずかしいですね。）

[OKUYAMA]

稲葉さん、

　ご回答いただき、ありがとうございました。

2017年3月14日火曜日 23時02分34秒 UTC+9 Inaba:

ただ参考までに、並列計算を行っている場合に限りますが
私が良く分からない発散で躓いたときはよく停止したプロセッサの番号を見ています。

　なるほど! こんな方法があるのですね。目から鱗でした。

 

ただこれだとエラーで計算終了しない限り分かりません。
おっしゃるようにやっぱりparaViewで経時的に残渣を見てみたいものです。
私も前々からそんな機能が欲しいと思っていたので少し調べてみました。

　わざわざお調べいただき、さらに実装方法までご教示いただき、ありがとうございました。

　私のほうでも試してみたいと思います。

[OKUYAMA]

奥山です。

　教えていただいた方法を試そうとしたのですが、残念ながらだめでした。

2017年3月14日火曜日 23時02分34秒 UTC+9 Inaba:

SmoothSolver.Cの該当関数内の残渣計算式の直前か直後あたりで
const fvMesh& mesh = this->matrix_.mesh(); // 戻り値はlduMesh型ですがfvMeshを継承しているのできっと大丈夫

　のところで、fvMeshの方がlduMeshを継承しているので、コンパイルを通りませんでした。

　その後自分でも調べてみたところ、fvMeshクラスがresidualという関数を持っていることが分かったのですが、これから要素ごとの残差をとれないものでしょうか。

　実際に値をとってみると、負になっているところが多かったので、そのまま残差として使えるものではなさそうなのですが。。。

　私の力では、ソースを見てもこの関数が何をやるものかわかりませんでしたので、ご存じであれば教えていただけると助かります。

　

　

　

　

[Inaba]

奥山さん

　fvMeshの方がlduMeshを継承しているので、コンパイルを通りませんでした。

逆でしたか。。
適当なことを言ってごめんなさい。

今もう一回Doxygenを見てたら思いついたのですが一回フィールドを通してメッシュを呼び出すのはどうでしょうか。
const fvMesh& mesh = this->matrix_.psi().mesh();

 

　その後自分でも調べてみたところ、fvMeshクラスがresidualという関数を持っていることが分かったのですが、これから要素ごとの残差をとれないものでしょうか。　

ｆｖMatrixのresidual()ですね。確かにありますね。というかこれが正解なのではないでしょうか。
ちゃんとフィールドで返ってくるみたいですし、これならソルバーの中だけで呼び出せますね。
計算内容もちゃんと行列式 A x = b の右辺-左辺みたいです。

この値に違和感を覚える理由として、まず負の値に関しては問題ないかと思います。
出力する際には絶対値にして合計するみたいなので、このresidual()で呼び出した値をmagしてしまえばいいのではないでしょうか。

ただここまで書いておいてあれなのですが、この 右辺-左辺 を normFactor のフィールドで割る必要があるかもしれません。
出力される残渣は sum(mag(右辺-左辺)) / sum(normFactor)　となっているようで、セル毎の残渣を得るためにはセル毎のnormFactorのフィールドも使用する必要がありそうです。
normFactorは
https://www.cfd-online.com/Forums/openfoam-solving/57903-residuals-convergence-segregated-solvers.html
を見たところ、
xの平均値のフィールドxrefとして | A x - A xref | + | b - A xref |がnormFactorにあたるようです。

セル毎に値の違うものみたいなのでどうにかnormFactorをフィールドとして呼び出すか計算をする必要があるかもしれません。
もうちょっと調べてみますね。
ただあまりに大きな残渣であればこれ無しでも検知できるのではないでしょうか。

不確かなことばかりで申し訳ありません。
詳しい人が現れることを祈るばかりです。

[Inaba]

奥山さん

ちなみになんですが、最初の回答の際のプロセッサの番号でエラー特定という方法に関して

下記のような実行ファイルを前に作ったのですがもしよかったら使ってみますか？
https://github.com/inabower/OpenFOAM4.1-utilities/tree/master/functionProcessor

実行ファイル自体はただcellDistとdecomposeParのmanual分割用のファイルを書き出すだけのものなのですが、
functionObjectsのcodedFunctionObjectを使ってゴリゴリと場合分けでプロセッサを好き勝手に分配していくことができます。
エラー箇所の特定にも使えますし、うまく分配するとたまに速度が上がったりエラーが減ったりします。

ご参考まで。

[Inaba]

連投で失礼します。

残渣の表示に関して、fvMatrixのresidual関数を試してみたら良さそうでしたので報告します。
https://github.com/inabower/OpenFOAM4.1-utilities/tree/master/residualSimpleFoam

ソルバーの内部のpEqnなどを生成した直後、createFieldで作ったフィールドに残渣を出力してみました。
下の図のようにinitialResidualが出力できていると思います。

pとUとの残渣の比較はできなさそうですが、どの場所がウィークポイントかは確認できるかと思います。

[OKUYAMA]

稲葉さん、

　奥山です。何から何まで丁寧に教えていただき、ありがとうございました。

　私の方でも試してみて、うまくいったら報告します。

　ところで、少し話が変わるのですが、稲葉さんはOF-4.1をお使いのようですが、やはり新しいバージョンに追随する方がいいものでしょうか。

　私が使っているのはOF-2.4.xで、chtMultiRegionSimpleFoamで熱伝達を解析しているのですが、とくだん不都合がないので1年以上そのまま使っています。

　(バージョンアップ直後は不具合がありそうで怖かったのと、設定を変えたりするのが面倒というのが理由です。)

[OKUYAMA]

奥山です。

　だいぶ時間がたってしまいましたが、私の方でも試してみましたので報告します。

　なお、私が使っているのはOF-2.4.0です。

　じつは、結局 fvMatrix::residualを使わずに、fvMatrix::solveSegregated を改造して残差を出力する処理を入れるようにしました。

　理由は以下の2点です。

　(1) fvMatrix::residual だと、なぜか並列実行時の残差がおかしくなった。

　　この理由は分かっていません。直列だと正しい値になります。

　(2) Uの残差について、実際のソルバー実行時だとUx→Uy→Uzの順に解き、それぞれでinitialResidualを計算してログに出力するが、

　　fvMatrix::residual ではまとめて計算するのでUy, Uzの残差がログに出力されるものと違う。

　　(pyFoamPlotWatcherで残差を表示しているので、ログの残差と違うと少し困るのです。細かいことではありますが。)

　githubが使えないので、修正したソースの断片を添付します。

　ソルバーの方でres_Uxとかres_p_rghという名前で残差を出力したい変数に対して出力用のフィールドを作っておくと、そこに出力されるようになっています。

　一応、これで残差を出力できるようになりました。

　いろいろと教えていただき、ありがとうございました。

　たいへん勉強になりました。

[Inaba]

ご無沙汰しております。稲葉です。

今年度からしばらくOpenFOAMを離れることになりましたのでかなり間隔が空いてしまいました。

まずバージョンの追随については、私はできるだけ最新版を使ったほうがいいと考えております。

個人的な意見としては、

　①まず当然ですがバージョンを重ねていく毎に多くのバグが解消されていっています。

この場合のバグはプログラミング的なものと科学計算的なものの両方を含みます。

もしOF2.4で原因不明のバグに悩まされていたとしてもそれは3.0.1や4.1では解消されている可能性も少なくありません。

その場合に例えばここやCFD Onlineで質問をしても他者は再現できないため解決しない恐れがあります。

　②またバージョンの変更作業は特に自分のソルバーを持っている場合には結構大変な作業になりますが、個人的にはこの作業の中で新たなエラーを見つけることが何度かありました。

ケースフォルダの修正に関しても同様です。

変更作業がリファクタリング作業を兼ねることにもなるという点でもメリットがあるかと思います。

　③最後に最新版でバグを見つけた場合に開発者に報告するだけでもOFの開発に貢献できると思います。

　一方でどのバージョンに追随するかという問題があります。

GPUを使えるRapidCFDやchtのリージョン間の行列を結合して陰的に計算が出来るextendなど、OF4.1以外にも別系統で開発されているものがあります。（それぞれ他にもいろいろな機能があります。）

今上げた二つは確か3.0.1から派生したものですので、逆に4.1から変更するのは大変かもしれないので今のうちに2.4のものから修正していくほうが良いかと思います。

というわけでバージョンをあげる前に一度色々な派生OpenFOAMを調べてみることをお勧めします。

また添付していただいたファイルにつきまして、

残渣の出力ができたようでなによりです。

やりとりの中で私のほうも大変勉強になりました。

ありがとうございました。