はじめに
縁あって流体解析なんぞやっていて，オープンソースのOpenFOAMをよく使うので，それについてメモ書き。ややこしいコードなので間違ってるかもしれませんが，そのときはご指摘ください。
ブログの本筋からは大分外れててすいません。
概要
OpenFOAM2.3.xのpisoFOAM中にある非直交補正のループ内に
pisoFoam.C
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()へのリンクができています。リンクをたどると以下のような候補が表示されます。
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
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行目からの以下の部分
955 if (faceFluxCorrectionPtr_)
956 {
957 fieldFlux += *faceFluxCorrectionPtr_;
958 }
faceFluxCorrectionPtr_の示す値が追加されています。これが何をやっているかを調べてみます。
faceFluxCorrectionPtr_について
faceFluxCorrecitonPtr_のゲッターとしてfaceFluxCorrectionPtr()関数があります。fvMatrix.Hより
fvMatrix.H
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
(setq gdb-many-windows t)
EmacsをGUIで起動して，M-x gdbでGDBを起動します。起動時にはpisoFOAMが立ち上がるようにします。OpenFOAMのコードを走らせるときにはLocal bufferを閉じておきます。こうしないとデバッグが途中で停止します。
cavity flowのチュートリアルをもとにしてデバッグします。fvSchemesはそのままです。非直交補正のについて調べるので，laplacianスキーム，snGradスキームがcorrectedになっているか確認して下さい。
fvSchemes
/*--------------------------------*- 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
/*--------------------------------*- 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
/*--------------------------------*- 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で計算を走らせます。コールスタックを抜き出したのが以下。
#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
declareFvmLaplacianScalarGamma(scalar);
となっています。これはマクロでその定義はすぐ上にあります。
gaussLapacianSchemes.C
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で見てみます。
71 gammaMagSf*this->tsnGradScheme_().correction(vf)
GlobalでtsnGradScheme_().correction()を調べる
correction()関数の候補はたくさん表示されます。ここではcorrectedSnGradを使っているので，correctedSnGrad.Cのcorrection()を見てみます。
実装が以下で，correction()関数の計算には102行目にあるfullGradCorrection()関数が関わっています。
correctedSnGrad.C
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
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でいろいろたどれるので，使ってみてください。