OpenFOAM ３次元解析における計算コストに関して
723 views
Skip to first unread message
takemitsu
Jun 10, 2020, 4:04:29 AM6/10/20
to OpenFOAM
いつもお世話になっていますtakemitsuです．
本日も皆さんのアドバイスを頂きたく投稿させてもらいました．
現在，buoyantPimpleFoamを用いて熱流動解析を行っているのですが，大型計算機センターを利用しても計算の進行が非常に遅いと感じています．
Re=100で行っています．
blockMeshは以下に通りです
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 0.01;
vertices
(
(-4 -0.11 1) //0
(10 -0.11 1) //1
(10 0.79 1) //2
(-4 0.79 1) //3
(-4 -0.11 6) //4
(10 -0.11 6) //5
(10 0.79 6) //6
(-4 0.79 6) //7
);
blocks
(
hex (0 1 2 3 4 5 6 7) (1400 90 40) simpleGrading (1 1 1)
);
edges
(
);
boundary
(
bottom
{
type cyclic;
neighbourPatch top;
faces
(
(1 5 4 0)
);
}
top
{
type cyclic;
neighbourPatch bottom;
faces
(
(2 3 7 6)
);
}
inlet
{
type patch;
faces
(
(0 4 7 3)
);
}
outlet
{
type patch;
faces
(
(2 6 5 1)
);
}
symFront
{
type zeroGradient;
faces
(
(4 5 6 7)
);
}
symBack
{
type zeroGradient;
faces
(
(0 1 2 3)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
また，snappyHexMeshで形状付近を細かくした結果，格子点数は7310259点です．
メッシュ自体はおかしいないとcheckMeshによって確認しました．
離散化は，以下のようにしています
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 4.1 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
ddtSchemes
{
default Euler;
}
gradSchemes
{
default Gauss linear;
}
divSchemes
{
default none;
div(phi,U) Gauss QUICK;
div(phi,h) Gauss QUICK;
div(phi,e) Gauss QUICK;
div(phi,k) Gauss QUICK;
div(phi,epsilon) Gauss QUICK;
div(phi,R) Gauss QUICK;
div(phi,K) Gauss QUICK;
div(phi,Ekp) Gauss QUICK;
div(R) Gauss QUICK;
div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
}
laplacianSchemes
{
default Gauss linear corrected;
}
interpolationSchemes
{
default linear;
}
snGradSchemes
{
default corrected;
}
// ************************************************************************* //
代数方程式のソルバーは以下のようにしています．
*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 4.1 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
solvers
{
"rho.*"
{
solver PCG;
preconditioner DIC;
tolerance 0;
relTol 0;
}
p_rgh
{
solver PCG;
preconditioner DIC;
tolerance 1e-8;
relTol 0.01;
}
p_rghFinal
{
$p_rgh;
relTol 0;
}
"(U|h|e|k|epsilon|R)"
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-6;
relTol 0.1;
}
"(U|h|e|k|epsilon|R)Final"
{
$U;
relTol 0;
}
}
PIMPLE
{
momentumPredictor yes;
nOuterCorrectors 2;
nCorrectors 2;
nNonOrthogonalCorrectors 0;
}
// ************************************************************************* //
このように設定した場合，実時間で24時間計算を行っても書き出し時間は0.1ｓ程です．これはOpenFOAMの計算コストとして正当なのでしょうか？
教えて下さい．
takemitsu
Jun 12, 2020, 12:21:43 AM6/12/20
to OpenFOAM
上記の情報に追加します．
並列計算は40分割で行っています．その際，大型計算機センターのコアを40コア用いて行っている状況です．
その他，情報が足りないのであれば言って下さい．
宜しくおねがいします．
ishigaki
Jun 12, 2020, 7:29:54 PM6/12/20
to OpenFOAM
最小のメッシュ幅と時間刻みと流速からクーラン数計算してみて、時間刻みが想定より小さいですか？時間刻みが小さすぎるなら、メッシュが潰れてるところがあって、最小メッシュ幅が小さすぎるかもしれません。
takemitsu
Jun 13, 2020, 5:09:36 AM6/13/20
to OpenFOAM
ishigakiさんへ
質問に答えて下さりありがとうございます．最小メッシュ幅はベースメッシュ幅0.01cmを形状表面でlevel3まで細分化しているので1.25e-5 mだと考えられます．
入り口流速は0.15m/sで，logより時間刻みは1.14e-5程度でした．したがってこれよりクーラん数C=0.136程度になります．
ここで2つ疑問があります．数値解析の経験がほとんどないことやOpenFOAMに未だよく理解していないのですごく初歩的な質問になるのですが宜しくおねがいします．
1つ目は，ishigakiさんがおっしゃる時間刻みの想定とは，どのようにされていますか？
2つ目は，controlDictにて最大クーラン数maxCoを指定しており，これを私は0.1と指定しました．これは上記の計算によって導かれた値と不適合だと思いますがこれが原因で計算負荷がかかっているのでしょうか？
長々とすみませんが宜しくおねがいします．
2020年6月13日土曜日 8時29分54秒 UTC+9 ishigaki:
ishigaki
Jun 13, 2020, 9:06:55 AM6/13/20
to OpenFOAM
>1つ目は，ishigakiさんがおっしゃる時間刻みの想定とは，どのようにされていますか？
時間刻みはクーラン数が最大で0.3程度くらいに設定して計算することが多いです。今回の計算もそこそこ妥当じゃないですか？
>2つ目は，controlDictにて最大クーラン数maxCoを指定しており，これを私は0.1と指定しました．
controlDictの内容がわからないので、はっきりしませんが、buoyantPimpleFoamでは時間刻み固定の計算と最大クーラン数で時間刻みを調整する
計算の両方が使えたはずです。
最大Coを0.1にしてて、0.136なら、妥当なところだとおもいます。maxCo=0.3〜0.5くらいなら、計算できることも
結構ありますが、結果を見て判断というところだと思います。
回答の順序がおかしくなって申し訳ないですが、24時間で0.1秒分しか計算進まない件ですが、特別遅いってこともないとおもいます。
730万メッシュを40並列で、1ステップあたり10秒弱かかるくらいなら仕方ないかという気もします。
早くしたいなら
・メッシュ数を見直す（メッシュサイズは必要最低限になってます？）
・並列数を増やす
・macCoを上げる
・反復解法を変えてみる
くらいじゃないでしょうか。
2020年6月13日土曜日 18時09分36秒 UTC+9 takemitsu:
takemitsu
Jun 14, 2020, 3:37:49 AM6/14/20
to OpenFOAM
ishigakiさんへ
ありがとうございます．色々勉強になります．
現在の自分の計算はスパコンを使ってもこの程度の計算しか進行しないことがよくわかりました．提案していただいたように，メッシュ数を減らすなど見直しをし，クーラン数を少し上げたりなどしてみます．
現在の自分の計算はスパコンを使ってもこの程度の計算しか進行しないことがよくわかりました．提案していただいたように，メッシュ数を減らすなど見直しをし，クーラン数を少し上げたりなどしてみます．
それでも，こちらが思うほど計算コストに成果が出なければ反復解法を変えていきたいと思います．
2020年6月13日土曜日 22時06分55秒 UTC+9 ishigaki:
小南秀彰
Jun 15, 2020, 10:05:04 AM6/15/20
to OpenFOAM
takemitsu ishigakiさん、ishigakiさん （おふたりへ）
fvShhmesの中のddtSchemesを Euler;（オイラー陽解法）としているようです。、
このスキームの場合、通常はcontrolDictの中の maxCo(最大クーラン数)は1です。(数値計算理論での安定領域の上限値が1のためです)
流れが不安定な場合は繰返収束計算途中での局所クーラン数が急に大きくなって計算が発散することがあします。それを予防するたために、maxCoを1以下にする場合もあります。しかし局所クーラン数が1以上となっても収束計算が安定的に進行する場合もあります。
buoyantPimpleFoamというソルバーで、必ずmaxCo=0.3〜0.5としないといけないというわけではありません。
さて、
ddtSchemesを CrankNicolson 1; 変えてみてください。
このスキームの場合、maxCo(最大クーラン数)の数値計算理論上の安定領域の上限は2ですが、もっと大きくしても収束計算が安定的に進行する場合が多いです。
（OpenFoam v1806では、buoyantPimpleFoamソルバーでCrankNicolsonスキームが使用できました。）
2020年6月14日日曜日 16時37分49秒 UTC+9 takemitsu:
Masashi Imano
Jun 15, 2020, 8:54:53 PM6/15/20
to OpenFOAM
今野です．
system/fvSolution
```
p_rgh
{
solver PCG;
preconditioner DIC;
tolerance 1e-8;
relTol 0.01;
}
PIMPLE
{
momentumPredictor yes;
nOuterCorrectors 2;
nCorrectors 2;
nNonOrthogonalCorrectors 0;
}
```
p_rghのtoleranceが1e-8になっていますが，
これを解(簡単に調べるなら，ログの平均・最大クーラン数，もしくは，probeによる重要地点の時系列モニター値をグラフで比較)
が僅かにしか変わらない程度に緩和する手法もあります．
また，同様な手法で nOuterCorrectors を 1 にする事も考えられます．
なお，40並列と並列数が少ないのであれば，
恐らく p_rgh の solver もしくは，preconditioner として GAMG を用いたほうが速いので，
例えば以下のようにして，チュートリアルでGAMGを使用している例を調べて，
設定を変更してみてください．
```
find $FOAM_TUTORIALS -name fvSolution | xargs grep "solver.*GAMG"
find $FOAM_TUTORIALS -name fvSolution | xargs grep "precon.*GAMG"
```
system/fvSchemes
```
divSchemes
{
default none;
div(phi,U) Gauss QUICK;
div(phi,h) Gauss QUICK;
div(phi,e) Gauss QUICK;
div(phi,k) Gauss QUICK;
div(phi,epsilon) Gauss QUICK;
div(phi,R) Gauss QUICK;
div(phi,K) Gauss QUICK;
div(phi,Ekp) Gauss QUICK;
div(R) Gauss QUICK;
div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
}
```
解法の妥当性や高速性を模索している段階では，
移流項スキームの不安定性の影響を受けないように，
hotRoomチュートリアル同様，
移流項の離散化スキームは upwind としていほうが良いと思います．
また，linear で良いものは linear のままにしたほうが良いと思います．
tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/system/fvSchemes
```
divSchemes
{
default none;
div(phi,U) Gauss upwind;
div(phi,h) Gauss upwind;
div(phi,e) Gauss upwind;
div(phi,k) Gauss upwind;
div(phi,epsilon) Gauss upwind;
div(phi,R) Gauss upwind;
div(phi,K) Gauss linear;
div(phi,Ekp) Gauss linear;
div(R) Gauss linear;
div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
}
```
解法が定まった後，U, hや乱流統計量等の移流拡散方程式で解いている場の
移流項の離散化スキームを高精度化すると思いますが，
QUICKはアンダーシュートや(工学的な意味での)オーバーシュートを伴なうので，
直交格子でない場合には，hや0以下にならない乱流統計量に用いると，
一般に安定性に欠けます．
ベンチマークテスト目的で，離散化スキームが変更できない等の拘束が無いのであれば，
hや乱流統計量アンダーシュートやオーバーシュートが少ない高次制度のスキームをお勧めします．
なお，チュートリアルで h の移流項スキームとして何が使われているかは，
例えば以下のようにして調べられます．
```
grep -r "div(phi,h)" $FOAM_TUTORIALS
grep -rh "div(phi,h)" $FOAM_TUTORIALS | tr -s "[:blank:]" | sort | uniq
```
東京大学情報基盤センターお試しアカウント付き並列プログラミング講習会「OpenFOAM講習会」
「OpenFOAM中級 演習資料」 の "19.2 速度Uの移流項スキームを変更"も参考にしてください．
以上です．
takemitsu
Jun 16, 2020, 9:45:54 AM6/16/20
to OpenFOAM
小南秀彰さん
夜分に失礼します．質問に対して回答していただきありがとうございます．
アドバイス通り，fvShhmesの中のddtSchemesを Euler法をddtSchemesを CrankNicolson 1; 変え，maxCo(最大クーラン数)を1にしてみたところ，時間刻み7e-4あたりで温度場が発散しました．エラーは以下のとおりです．
[33]
[33]
[33] --> FOAM FATAL ERROR:
[33] Negative initial temperature T0: -358.854
[33]
[33] From function Foam::scalar Foam::species::thermo<Thermo, Type>::T(Foam::scalar, Foam::scalar, Foam::scalar, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar, Foam::scalar)const, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar, Foam::scalar)const, Foam::scalar (Foam::species::thermo<Thermo, Type>::*)(Foam::scalar)const) const [with Thermo = Foam::hConstThermo<Foam::perfectGas<Foam::specie> >; Type = Foam::sensibleEnthalpy; Foam::scalar = double; Foam::species::thermo<Thermo, Type> = Foam::species::thermo<Foam::hConstThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleEnthalpy>]
[33] in file /home/w10207/source/openfoam/v1806/OpenFOAM-v1806/src/thermophysicalModels/specie/lnInclude/thermoI.H at line 54.
[33]
FOAM parallel run aborting
[33]
[33] #0 Foam::error::printStack(Foam::Ostream&) at ??:?
[33] #1 Foam::error::abort() at ??:?
[33] #2 Foam::species::thermo<Foam::hConstThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleEnthalpy>::THs(double, double, double) const at ??:?
[33] #3 Foam::heRhoThermo<Foam::rhoThermo, Foam::pureMixture<Foam::constTransport<Foam::species::thermo<Foam::hConstThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleEnthalpy> > > >::calculate(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&, bool) at ??:?
[33] #4 Foam::heRhoThermo<Foam::rhoThermo, Foam::pureMixture<Foam::constTransport<Foam::species::thermo<Foam::hConstThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleEnthalpy> > > >::correct() at ??:?
[33] #5 ? at ??:?
[33] #6 __libc_start_main in /lib64/libc.so.6
[33] #7 ? at ??:?
Abort(1) on node 33 (rank 33 in comm 0): application called MPI_Abort(MPI_COMM_WORLD, 1) - process 33
この原因がわかりません．最後のクーラン数はこちらになります．
Courant Number mean: 0.000464019 max: 0.732827
deltaT = 1.13861e-07
最大クーラン数は1以下になっているのですが…，何かご存知であれば教えて下さい．
また，先程の文章でEuler;（オイラー陽解法）とありましたが，OpenFOAMにおけるEulerは陰的とどこかで見たことがあるのですがこれは私の認識違いなのでしょうか？まだまだ未熟者で申し訳ないのですが宜しくおねがいします．
2020年6月15日月曜日 23時05分04秒 UTC+9 小南秀彰:
Message has been deleted
ishigaki
Jun 16, 2020, 10:06:32 AM6/16/20
to OpenFOAM
小南様
> buoyantPimpleFoamというソルバーで、必ずmaxCo=0.3〜0.5としないといけないというわけではありません。
コメントありがとうございました。おっしゃるとおりです。言葉足らずでした。経験的にこれくらいだと
うまくいくことが多いくらいの意味合いでした。
> Euler;（オイラー陽解法）
オイラー陰解法だと思いますが，いかがでしょうか。OpenFOAMのユーザーガイドにも陰解法となってますが。。。
https://www.openfoam.com/documentation/guides/latest/doc/guide-schemes-time-euler.html2020年6月15日月曜日 23時05分04秒 UTC+9 小南秀彰:
takemitsu
Jun 16, 2020, 10:10:09 AM6/16/20
to OpenFOAM
今野さん
夜分に失礼します．
いつも返信してくださりありがとうございます．OpenFOAM内の代数方程式ソルバーはいままで変更したことなく，また，ほとんど勉強してきていなかったので勉強した後頂いたアドバイスを参考の基で計算を行いたいと思います．計算した後また問題点が生じた際は宜しくおねがいします．
また，参考になるサイトのURLを添付してくださりありがとうございます．合わせて勉強させていただきます．
2020年6月16日火曜日 9時54分53秒 UTC+9 Masashi Imano:
小南秀彰
Jun 19, 2020, 8:22:20 AM6/19/20
to OpenFOAM
FATAL ERROR:の内容は、温度が負値とのことです。
最大クーラン数を理論上の安定上限以下にしていて起きるということは、現象そのものの不安定が高いためだと思います。
もともとの質問は、非定常計算での計算時間が掛かり過ぎることについてでしたの。現象そのものの不安定性が強いため、オイラー陰解法を使っても最大クーラン数を小さくしないといけないという状況に陥っているようです。
今の計算の目的あるいは意図を把握できていませんが、(a) メッシュ数を少なくする。(b) TVDでないQUICKスキームをやめて upwindスキームにする というようなことをして、とりあえず、そこそこの実時間で計算が終わる条件を見つけたほうが良いと思います。
2020年6月16日火曜日 22時45分54秒 UTC+9 takemitsu:
takemitsu
Jun 20, 2020, 11:22:49 AM6/20/20
to OpenFOAM
小南秀彰さん
2020年6月19日金曜日 21時22分20秒 UTC+9 小南秀彰:
夜分に失礼します．
クーラン数が安定上限以下でこのようなことが起こる理由が，現象の不安定性によるものなのですね．
頂いたアドバイスを基に試行錯誤を続けて解決に努めたいと思います．
ありがとうございます．
2020年6月19日金曜日 21時22分20秒 UTC+9 小南秀彰:
Reply all
Reply to author
Forward
0 new messages