はじめに
縁あって流体解析なんぞやっていて,オープンソースのOpenFOAMをよく使うので,それについてメモ書き。ややこしいコードなので間違ってるかもしれませんが,そのときはご指摘ください。
ブログの本筋からは大分外れててすいません。
概要
OpenFOAM2.3.xのpisoFOAM中にある非直交補正のループ内に
pisoFoam.C
1 2 3 4 5 |
120 if (nonOrth == nNonOrthCorr) 121 { 122 phi = phiHbyA - pEqn.flux(); 123 } |
pEqn.flux()関数の呼び出しがあります。(行頭は行番号です)flux()関数がどうなっているかをGDBによるデバッグとGNU Globalによるソースコードの読解で調べてみました。
GDBとGlobal
GDB
GDBはGCCのデバッガです。OpenFOAMでGDBを使えるようにするには,コンパイルオプションWM_COMPILE_OPTIONをDebugにしてAllwmakeします。詳細は下記なんかがわかりやすいです。
http://www.tfd.chalmers.se/~hani/kurser/OS_CFD_2008/debugging.pdf
Global
GNU Globalはソースコードにタグをつけてくれるソフトです。関数のソースがどこにあるかなどの情報が見やすくなります。EclipseのF3を押しても同様のことができますが,候補が一覧で表示されるので,個人的にはEclipseよりも使いやすいです。
参考
http://www.machu.jp/diary/20090307.html#p01
http://uguisu.skr.jp/Windows/gtags.html
Globalはソースコードにタグ付けを行うので,OpenFOAMで使う場合にはtutorials, etc, bin, docを除外します。除外の仕方は以下を参考にしました。
http://d.hatena.ne.jp/ohtorii/20110219/1298092604
pisoFOAMのソースコードを見る
はじめにGlobalでflux()関数を見る
GlobalでpisoFoam.Cを見ると,flux()へのリンクができています。リンクをたどると以下のような候補が表示されます。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 |
flux 58 src/OpenFOAM/primitives/Vector/vector/vector.H class flux flux 65 src/OpenFOAM/primitives/Vector/vector/vector.H class flux<scalar> flux 57 src/finiteVolume/finiteVolume/convectionSchemes/boundedConvectionScheme/boundedConvectionScheme.C boundedConvectionScheme<Type>::flux flux 64 src/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionScheme.C gaussConvectionScheme<Type>::flux flux 61 src/finiteVolume/finiteVolume/convectionSchemes/multivariateGaussConvectionScheme/multivariateGaussConvectionScheme.C multivariateGaussConvectionScheme<Type>::flux flux 44 src/finiteVolume/finiteVolume/fvc/fvcFlux.C flux flux 62 src/finiteVolume/finiteVolume/fvc/fvcFlux.C flux flux 80 src/finiteVolume/finiteVolume/fvc/fvcFlux.C flux flux 98 src/finiteVolume/finiteVolume/fvc/fvcFlux.C flux flux 117 src/finiteVolume/finiteVolume/fvc/fvcFlux.C flux flux 132 src/finiteVolume/finiteVolume/fvc/fvcFlux.C flux flux 149 src/finiteVolume/finiteVolume/fvc/fvcFlux.C flux flux 166 src/finiteVolume/finiteVolume/fvc/fvcFlux.C flux flux 883 src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C flux() const flux 213 src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/limitedSurfaceInterpolationScheme/limitedSurfaceInterpolationScheme.C limitedSurfaceInterpolationScheme<Type>::flux |
ごちゃごちゃしてますが,pEqnがfvMatrixのオブジェクトなので,fvMatrix.Cのflux()を見てみます。
fvMatrix.C
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 |
880 template<class Type> 881 Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh> > 882 Foam::fvMatrix<Type>:: 883 flux() const 884 { 885 if (!psi_.mesh().fluxRequired(psi_.name())) 886 { 887 FatalErrorIn("fvMatrix<Type>::flux()") 888 << "flux requested but " << psi_.name() 889 << " not specified in the fluxRequired sub-dictionary" 890 " of fvSchemes." 891 << abort(FatalError); 892 } 893 894 // construct GeometricField<Type, fvsPatchField, surfaceMesh> 895 tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tfieldFlux 896 ( 897 new GeometricField<Type, fvsPatchField, surfaceMesh> 898 ( 899 IOobject 900 ( 901 "flux("+psi_.name()+')', 902 psi_.instance(), 903 psi_.mesh(), 904 IOobject::NO_READ, 905 IOobject::NO_WRITE 906 ), 907 psi_.mesh(), 908 dimensions() 909 ) 910 ); 911 GeometricField<Type, fvsPatchField, surfaceMesh>& fieldFlux = tfieldFlux(); 912 913 for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++) 914 { 915 fieldFlux.internalField().replace 916 ( 917 cmpt, 918 lduMatrix::faceH(psi_.internalField().component(cmpt)) 919 ); 920 } 921 922 FieldField<Field, Type> InternalContrib = internalCoeffs_; 923 924 forAll(InternalContrib, patchI) 925 { 926 InternalContrib[patchI] = 927 cmptMultiply 928 ( 929 InternalContrib[patchI], 930 psi_.boundaryField()[patchI].patchInternalField() 931 ); 932 } 933 934 FieldField<Field, Type> NeighbourContrib = boundaryCoeffs_; 935 936 forAll(NeighbourContrib, patchI) 937 { 938 if (psi_.boundaryField()[patchI].coupled()) 939 { 940 NeighbourContrib[patchI] = 941 cmptMultiply 942 ( 943 NeighbourContrib[patchI], 944 psi_.boundaryField()[patchI].patchNeighbourField() 945 ); 946 } 947 } 948 949 forAll(fieldFlux.boundaryField(), patchI) 950 { 951 fieldFlux.boundaryField()[patchI] = 952 InternalContrib[patchI] - NeighbourContrib[patchI]; 953 } 954 955 if (faceFluxCorrectionPtr_) 956 { 957 fieldFlux += *faceFluxCorrectionPtr_; 958 } 959 960 return tfieldFlux; 961 } |
途中boundaryFieldの量を各パッチで計算したりしています。これは
の資料にあるように,flux()は圧力勾配から来る流束の修正を行っているように見えます。
(上のNozakiさんの資料は大変参考になります)
気になるのが955行目からの以下の部分
1 2 3 4 5 |
955 if (faceFluxCorrectionPtr_) 956 { 957 fieldFlux += *faceFluxCorrectionPtr_; 958 } |
faceFluxCorrectionPtr_の示す値が追加されています。これが何をやっているかを調べてみます。
faceFluxCorrectionPtr_について
faceFluxCorrecitonPtr_のゲッターとしてfaceFluxCorrectionPtr()関数があります。fvMatrix.Hより
fvMatrix.H
1 2 3 4 5 6 |
321 //- Return pointer to face-flux non-orthogonal correction field 322 surfaceTypeFieldPtr& faceFluxCorrectionPtr() 323 { 324 return faceFluxCorrectionPtr_; 325 } |
pisoFoamのfaceFluxCorrectionPtr()関数の詳細について次で調べてみます。
GDBでデバッグしてみる
ここからGDBによるデバッグでfaceFluxCorrectionPtr()関数について調べてみます。
デバッグのためのセットアップ
GDBの結果を見るのにはEmacsが便利だと思います。以下の設定を.emacs.d/init.elや.emacsにします。Emacs+GDBは以下参照。 http://d.hatena.ne.jp/higepon/20090505/p1
1 2 |
(setq gdb-many-windows t) |
EmacsをGUIで起動して,M-x gdbでGDBを起動します。起動時にはpisoFOAMが立ち上がるようにします。OpenFOAMのコードを走らせるときにはLocal bufferを閉じておきます。こうしないとデバッグが途中で停止します。
cavity flowのチュートリアルをもとにしてデバッグします。fvSchemesはそのままです。非直交補正のについて調べるので,laplacianスキーム,snGradスキームがcorrectedになっているか確認して下さい。
fvSchemes
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 |
/*--------------------------------*- C++ -*----------------------------------* | ========= | | | / F ield | OpenFOAM: The Open Source CFD Toolbox | | / O peration | Version: 2.3.0 | | / 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 limitedLinearV 1; div(phi,k) Gauss limitedLinear 1; div(phi,epsilon) Gauss limitedLinear 1; div(phi,R) Gauss limitedLinear 1; div(R) Gauss linear; div(phi,nuTilda) Gauss limitedLinear 1; div((nuEff*dev(T(grad(U))))) Gauss linear; } laplacianSchemes { default Gauss linear corrected; } interpolationSchemes { default linear; } snGradSchemes { default corrected; } fluxRequired { default no; p ; } // ************************************************************************* // |
fvSolutionは以下のようにしました。
fvSolution
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 |
/*--------------------------------*- C++ -*----------------------------------* | ========= | | | / F ield | OpenFOAM: The Open Source CFD Toolbox | | / O peration | Version: 2.3.0 | | / A nd | Web: www.OpenFOAM.org | | / M anipulation | | *---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; location "system"; object fvSolution; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // solvers { p { solver PCG; preconditioner DIC; tolerance 1e-06; relTol 0; } pFinal { solver PCG; preconditioner DIC; tolerance 1e-06; relTol 0; } U { solver smoothSolver; smoother GaussSeidel; tolerance 1e-05; relTol 0; } } PISO { nCorrectors 2; nNonOrthogonalCorrectors 1; pRefCell 0; pRefValue 0; } |
メッシュ数が多いとデバッグで見ていくときに何かと大変なので,セル数5個のメッシュを作成しました。以下がそのメッシュを生成するためのblockMeshdictです。blockMeshでメッシュを生成します。初期条件,境界条件は適当に設定します。
blockMeshDict
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 |
/*--------------------------------*- C++ -*----------------------------------* | ========= | | | / F ield | OpenFOAM: The Open Source CFD Toolbox | | / O peration | Version: 2.3.0 | | / A nd | Web: www.OpenFOAM.org | | / M anipulation | | *---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; object blockMeshDict; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // convertToMeters 1; vertices ( (0 0 0) (1 0 0) (5 5 0) (1 5 0) (0 0 0.1) (1 0 0.1) (5 5 0.1) (1 5 0.1) ); blocks ( hex (0 1 2 3 4 5 6 7) (1 5 1) simpleGrading (1 1 1) ); edges ( ); boundary ( wall { type wall; faces ( (1 2 6 5) (0 3 7 4) ); } inlet { type patch; faces ( (0 1 5 4) ); } outlet { type patch; faces ( (2 3 7 6) ); } frontAndBack { type empty; faces ( (0 1 2 3) (4 5 6 7) ); } ); mergePatchPairs ( ); // ************************************************************************* // |
GDBでブレークポイントを設定する
実際にブレイクポイントを設定して,pisoFoamのfaceFluxCorrectionPtr()関数の挙動を見てみます。
fvMatrix.Hはsrc/finiteVolume/fvMatrices/fvMatrix/fvMatrix.Hにありますが,計算時に参照されるのはsrc/finiteVolume/lnInclude/fvMatrix.Hです。ここにブレイクポイントを置きます。
runで計算を走らせます。コールスタックを抜き出したのが以下。
1 2 3 4 5 6 7 |
#0 Foam::fvMatrix<double>::faceFluxCorrectionPtr (this=0x828fe0) at lnInclude/fvMatrix.H:324 #1 0x00007ffff66f5404 in Foam::fv::gaussLaplacianScheme<double, double>::fvmLaplacian (this=0x827ba0, gamma=..., vf=...) at finiteVolume/laplacianSchemes/gaussLaplacianScheme/gaussLaplacianSchemes.C:117 #2 0x00007ffff6713278 in Foam::fv::laplacianScheme<double, double>::fvmLaplacian (this=0x827ba0, gamma=..., vf=...) at lnInclude/laplacianScheme.C:107 #3 0x000000000045d2c1 in Foam::fvm::laplacian<double, double> (gamma=..., vf=..., name=...) at /home/ishigaki/OpenFOAM/Debug/OpenFOAM-2.3.xdbg/src/finiteVolume/lnInclude/fvmLaplacian.C:219 #4 0x0000000000451497 in Foam::fvm::laplacian<double, double> (gamma=..., vf=...) at /home/ishigaki/OpenFOAM/Debug/OpenFOAM-2.3.xdbg/src/finiteVolume/lnInclude/fvmLaplacian.C:251 #5 0x0000000000448efc in main (argc=1, argv=0x7fffffffbfd8) at pisoFoam.C:102 |
#1をクリックするとその呼出が表示されます。gaussLapalacianSchemes.Cの該当箇所を見ると,
gaussLapacianSchemes.C
1 2 |
declareFvmLaplacianScalarGamma(scalar); |
となっています。これはマクロでその定義はすぐ上にあります。
gaussLapacianSchemes.C
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 |
39 #define declareFvmLaplacianScalarGamma(Type) \ 40 \ 41 template<> \ 42 Foam::tmp<Foam::fvMatrix<Foam::Type> > \ 43 Foam::fv::gaussLaplacianScheme<Foam::Type, Foam::scalar>::fvmLaplacian \ 44 ( \ 45 const GeometricField<scalar, fvsPatchField, surfaceMesh>& gamma, \ 46 const GeometricField<Type, fvPatchField, volMesh>& vf \ 47 ) \ 48 { \ 49 const fvMesh& mesh = this->mesh(); \ 50 \ 51 GeometricField<scalar, fvsPatchField, surfaceMesh> gammaMagSf \ 52 ( \ 53 gamma*mesh.magSf() \ 54 ); \ 55 \ 56 tmp<fvMatrix<Type> > tfvm = fvmLaplacianUncorrected \ 57 ( \ 58 gammaMagSf, \ 59 this->tsnGradScheme_().deltaCoeffs(vf), \ 60 vf \ 61 ); \ 62 fvMatrix<Type>& fvm = tfvm(); \ 63 \ 64 if (this->tsnGradScheme_().corrected()) \ 65 { \ 66 if (mesh.fluxRequired(vf.name())) \ 67 { \ 68 fvm.faceFluxCorrectionPtr() = new \ 69 GeometricField<Type, fvsPatchField, surfaceMesh> \ 70 ( \ 71 gammaMagSf*this->tsnGradScheme_().correction(vf) \ 72 ); \ 73 \ 74 fvm.source() -= \ 75 mesh.V()* \ 76 fvc::div \ 77 ( \ 78 *fvm.faceFluxCorrectionPtr() \ 79 )().internalField(); \ 80 } \ 81 else \ 82 { \ 83 fvm.source() -= \ 84 mesh.V()* \ 85 fvc::div \ 86 ( \ 87 gammaMagSf*this->tsnGradScheme_().correction(vf) \ 88 )().internalField(); \ 89 } \ 90 } \ 91 \ 92 return tfvm; \ 93 } \ 94 \ 95 \ 96 template<> \ 97 Foam::tmp<Foam::GeometricField<Foam::Type, Foam::fvPatchField, Foam::volMesh> >\ 98 Foam::fv::gaussLaplacianScheme<Foam::Type, Foam::scalar>::fvcLaplacian \ 99 ( \ 100 const GeometricField<scalar, fvsPatchField, surfaceMesh>& gamma, \ 101 const GeometricField<Type, fvPatchField, volMesh>& vf \ 102 ) \ 103 { \ 104 const fvMesh& mesh = this->mesh(); \ 105 \ 106 tmp<GeometricField<Type, fvPatchField, volMesh> > tLaplacian \ 107 ( \ 108 fvc::div(gamma*this->tsnGradScheme_().snGrad(vf)*mesh.magSf()) \ 109 ); \ 110 \ 111 tLaplacian().rename("laplacian(" + gamma.name() + ',' + vf.name() + ')');\ 112 \ 113 return tLaplacian; \ 114 } |
64行目以降が非直交補正に関係有るところです。fvSchemesでfluxRequiredを指定してあるので,68行目でfaceFluxCorrectionPtr()関数に対して,snGradの非直交補正による補正量が代入されています。74行目でsource項に非直交補正による補正を代入しています。
このようにGDBを使うことで,OpenFOAMの深い関数呼び出しの関係を調べられます。途中マクロ展開で関数呼び出しを実装することがOpenFOAMの場合には多くあります。このため,関数にあたりをつけて,適切にブレイクポイントを設定する必要があります。
次に71行目のcorrection()関数についてGlobalで見てみます。
1 2 |
71 gammaMagSf*this->tsnGradScheme_().correction(vf) |
GlobalでtsnGradScheme_().correction()を調べる
correction()関数の候補はたくさん表示されます。ここではcorrectedSnGradを使っているので,correctedSnGrad.Cのcorrection()を見てみます。
実装が以下で,correction()関数の計算には102行目にあるfullGradCorrection()関数が関わっています。
correctedSnGrad.C
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 |
68 template<class Type> 69 Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh> > 70 Foam::fv::correctedSnGrad<Type>::correction 71 ( 72 const GeometricField<Type, fvPatchField, volMesh>& vf 73 ) const 74 { 75 const fvMesh& mesh = this->mesh(); 76 77 // construct GeometricField<Type, fvsPatchField, surfaceMesh> 78 tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tssf 79 ( 80 new GeometricField<Type, fvsPatchField, surfaceMesh> 81 ( 82 IOobject 83 ( 84 "snGradCorr("+vf.name()+')', 85 vf.instance(), 86 mesh, 87 IOobject::NO_READ, 88 IOobject::NO_WRITE 89 ), 90 mesh, 91 vf.dimensions()*mesh.nonOrthDeltaCoeffs().dimensions() 92 ) 93 ); 94 GeometricField<Type, fvsPatchField, surfaceMesh>& ssf = tssf(); 95 96 for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++) 97 { 98 ssf.replace 99 ( 100 cmpt, 101 correctedSnGrad<typename pTraits<Type>::cmptType>(mesh) 102 .fullGradCorrection(vf.component(cmpt)) 103 ); 104 } 105 106 return tssf; 107 } |
fullGradCorrection()はその上に定義されています。
correctedSnGrad.C
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 |
42 template<class Type> 43 Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh> > 44 Foam::fv::correctedSnGrad<Type>::fullGradCorrection 45 ( 46 const GeometricField<Type, fvPatchField, volMesh>& vf 47 ) const 48 { 49 const fvMesh& mesh = this->mesh(); 50 51 // construct GeometricField<Type, fvsPatchField, surfaceMesh> 52 tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tssf = 53 mesh.nonOrthCorrectionVectors() 54 & linear<typename outerProduct<vector, Type>::type>(mesh).interpolate 55 ( 56 gradScheme<Type>::New 57 ( 58 mesh, 59 mesh.gradScheme("grad(" + vf.name() + ')') 60 )().grad(vf, "grad(" + vf.name() + ')') 61 ); 62 tssf().rename("snGradCorr(" + vf.name() + ')'); 63 64 return tssf; 65 } |
mesh.nonOrthCorrectionVectors()とlinearで補間されたgradの値の内積を計算しています。これはJasakさんのD論 ( http://powerlab.fsb.hr/ped/kturbo/OpenFOAM/docs/HrvojeJasakPhD.pdf )にある,非直交補正の式と対応します。
以上の内容はNozakiさんの資料にも説明がありますので,そちらを参考にするとわかりやすいです。
まとめ
flux()関数を使った流束の補正とfaceFluxCorrectionPtr()関数から始まる非直交補正の詳細について調べてみました。GDBとGlobalでいろいろたどれるので,使ってみてください。