pisoFOAMのpEqn.flux()について調べてみる

はじめに

縁あって流体解析なんぞやっていて,オープンソースの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でいろいろたどれるので,使ってみてください。

こちらもどうぞ

コメントを残す

メールアドレスが公開されることはありません。