【OpenFOAM】手法を変えて処理を追加① : codedFunctionObject

この記事はオープンCAEアドベントカレンダー2020 3日目の記事です。

はじめに

僕はOpenFOAMの最大の利点はそのカスタマイズ性にあると思っています。今回はこのカスタマイズの有用性をなんとなく感じてもらうために、ある一つの簡単な処理を複数の方法で追加してみたいと思います。

動作環境

• OS : Ubuntu 20.04 LTS (WSL2 on Windows 10)
• OpenFOAM : ESI版 v2006 (from source)
• gcc : 9.3.0

手法一覧

具体的には以下の方法で行う予定です。実行確認済みのファイルはGitHubで公開しています。

• 例題を用意する ★今ここ (①)（キーワード：pimpleFoam）
• 既存機能でやるとしたら ★今ここ (①)（キーワード：functionObjects）

1. 計算ファイルに処理を書き込む ★今ここ (①)（キーワード：codedFunctionObject ）
2. ソルバーに処理を追加する (②)
3. 関数とクラスとして処理を追加する
4. 継承クラスに処理を追加
5. 設定項目を実行時に読み込む（キーワード：IOdictionary）
6. 継承クラスを実行時に選択する（キーワード：runTimeSelectionTable）
7. クラスを別でコンパイルする（キーワード：wmake libso）
8. 継承クラスの一つとしてcodedを追加（キーワード：codedTemplete）

● 例題を用意する (pimpleFoam)

追加する処理の動作確認を行うためのケースを用意します。今回は非圧縮流体の非定常計算を行うpimpleFoamのTJunctionというチュートリアルを例題として使用します。

$ mkdir -p $FOAM_RUN/001.codedFunctionObject
$ cp -r $FOAM_TUTORIALS/incompressible/pimpleFoam/RAS/TJunction $FOAM_RUN/001.codedFunctionObject/orgCase
$ cd $FOAM_RUN/001.codedFunctionObject/orgCase
$ foamRunTutorials
$ paraFoam

追加する処理

追加する処理として、今回は『出入口の流量の総和』をスコアとして出力する処理を作成することとします。これはつまり、上の図の中のinlet, outlet1, outlet2の三つの面の流量phiの和を求めることとなります。

● 既存機能でやるとしたら (functionObject)

この各面の流量の和はカスタマイズを行わずとも求めることができます。まず今回の例題を具体的に示すために、OpenFOAMの標準機能であるfunctionObjectsを使用して上記の総和をみましょう。
　functionObjectsはsystem/controlDictにfunctionsという項目を追加して色々な処理を追加する機能です。今回の処理の場合は、system/controlDictの最後の方を以下のように変更します。

$FOAM_RUN/001.codedFunctionObject/orgCase/system/controlDict

//（上略）
adjustTimeStep  yes;
maxCo           5; // ここまでは変更なし

functions
{
    inletFlow
    {
        type                surfaceFieldValue;
        libs                ("libfieldFunctionObjects.so");
        log                 true;
        executeControl      writeTime;
        executeInterval     1;
        writeFields         false;
        regionType          patch;
        operation           sum;
        surfaceFormat       none;
        fields              (phi);
        name                inlet;
    }
    outlet1Flow
    {
        $inletFlow;
        name                outlet1;
    }
    outlet2Flow
    {
        $inletFlow;
        name                outlet2;
    }
}

今回はsurfaceFieldValueという種類のfunctionObjectを使用しています。これを用いて流量場phiのpatchの合計を求めています。
　これにより計算後に出力されるファイルlog.pimpleFoamの一回のイタレーションの内容は以下のようになっています。

$ cd $FOAM_RUN/001.codedFunctionObject/orgCase
$ foamCleanTutorials
$ foamRunTutorials

log.pimpleFoam(抜粋)

deltaT = 0.00120482
Time = 0.00120482

PIMPLE: iteration 1
smoothSolver:  Solving for Ux, Initial residual = 1, Final residual = 5.02829e-10, No Iterations 2
smoothSolver:  Solving for Uy, Initial residual = 1, Final residual = 5.02829e-10, No Iterations 2
smoothSolver:  Solving for Uz, Initial residual = 0, Final residual = 0, No Iterations 0
GAMG:  Solving for p, Initial residual = 1, Final residual = 0.00324045, No Iterations 2
time step continuity errors : sum local = 7.52981e-05, global = -1.61311e-05, cumulative = -1.61311e-05
GAMG:  Solving for p, Initial residual = 0.00482057, Final residual = 7.62938e-07, No Iterations 10
time step continuity errors : sum local = 2.40386e-08, global = 3.40153e-09, cumulative = -1.61277e-05
smoothSolver:  Solving for epsilon, Initial residual = 0.620822, Final residual = 1.2033e-06, No Iterations 1
smoothSolver:  Solving for k, Initial residual = 1, Final residual = 2.61872e-10, No Iterations 2
ExecutionTime = 0.1 s  ClockTime = 0 s

surfaceFieldValue inletFlow write:
    total faces   = 25
    total area    = 0.0004

    sum(inlet) of phi = -7.86843e-06

surfaceFieldValue outlet1Flow write:
    total faces   = 25
    total area    = 0.0004

    sum(outlet1) of phi = -7.56386e-06

surfaceFieldValue outlet2Flow write:
    total faces   = 25
    total area    = 0.0004

    sum(outlet2) of phi = 1.5433e-05

Courant Number mean: 0.00787256 max: 0.0212701

ちなみに各timeStepの値などはposeProcessingに抽出して出力されています。
あとは目的に応じてposeProcessingやlog.pimpleFoamをGnuplotやJupyterや表計算ソフトなどで処理することになります。
　以上がこれから行う内容をOpenFOAMの既存機能でやるとしたらという内容です。これと同じ処理をいろいろな方法で行うのが今回の趣旨です。

1. 計算ファイルに処理を書き込む (codedFunctionObject)

ここからが本題です。

まず一つ目の方法としてcodedFunctionObjectとう機能を使用します。これは先述のfunctionObjectの一種であり、計算実行時に読み込むsystem/controlDictにC++のコードを記述するというものです。コードは実行時にコンパイルされるため、PythonやJuliaといったインタープリタ言語のような使い方ができます。
　OpenFOAMで用意されていないちょっと複雑な処理を行いたい場合にとても有用な機能です。
また、計算結果と実行コードが同時に保存されるという点も大きなメリットです。
　それでは具体的に、先ほどのsurfaceFieldValueと同じメッセージを出力させてみましょう。

$ mkdir -p $FOAM_RUN/001.codedFunctionObject
$ cp -r $FOAM_RUN/001.codedFunctionObject/orgCase $FOAM_RUN/001.codedFunctionObject/codedCase
$ cd $FOAM_RUN/001.codedFunctionObject/codedCase
$ foamCleanTutorials

$FOAM_RUN/001.codedFunctionObject/codedCase/system/controlDict

上略

functions
{
    flows
    {
        type                coded;
        name                flows;
        libs                (utilityFunctionObjects);

        codeExecute
        #{
            Info << "Running code execute" << endl;

            const surfaceScalarField& phi = mesh().lookupObject<surfaceScalarField>("phi");
            scalar score = 0.0;

            forAll(phi.boundaryField(), patchI)
            {
                const fvPatch& patch = phi.boundaryField()[patchI].patch();
                const scalarField& pphi =  phi.boundaryField()[patchI];
                if (!(patch.name()=="inlet" || patch.name()=="outlet1" || patch.name()=="outlet2")) continue;
                Info << "    sum(" << patch.name() << ") of phi = " << gSum(pphi) << endl;
                score += gSum(pphi);
            }
            Info << nl << "score: " << mesh().time().value() << tab << score << nl << endl;
        #};
    }
}

このようにcodeExecuteの中にC++のコードが記載されています。
　codeExecuteは各イタレーションごとに実行されるコードですが、例えばcodeEndを設置することで全てのイタレーションが終わった後に実行したりもできます。
　これを実行すると以下のように出力されます。

$ cd $FOAM_RUN/001.codedFunctionObject/codedCase
$ foamRunTutorials

$FOAM_RUN/001.codedFunctionObject/codedCase/log.pimpleFoam

Starting time loop

Using dynamicCode for functionObject flows at line 56 in "???/system/controlDict.functions.flows"
Could not load "???/dynamicCode/platforms/linux64GccDPInt32Opt/lib/libflows_1c917a97787494bcbf87e1064c2b107276011b9c.so"
???/dynamicCode/platforms/linux64GccDPInt32Opt/lib/libflows_1c917a97787494bcbf87e1064c2b107276011b9c.so: cannot open shared object file: No such file or directory
Creating new library in "dynamicCode/flows/platforms/linux64GccDPInt32Opt/lib/libflows_1c917a97787494bcbf87e1064c2b107276011b9c.so"
Invoking wmake libso ???/dynamicCode/flows
wmake libso ???/dynamicCode/flows
    ln: ./lnInclude
    dep: functionObjectTemplate.C
    Ctoo: functionObjectTemplate.C
    ld: ???/dynamicCode/flows/../platforms/linux64GccDPInt32Opt/lib/libflows_1c917a97787494bcbf87e1064c2b107276011b9c.so

中略

PIMPLE: iteration 1
smoothSolver:  Solving for Ux, Initial residual = 0.320314, Final residual = 3.21164e-07, No Iterations 2
smoothSolver:  Solving for Uy, Initial residual = 0.485008, Final residual = 5.38689e-07, No Iterations 2
smoothSolver:  Solving for Uz, Initial residual = 0.900395, Final residual = 6.13833e-07, No Iterations 2
GAMG:  Solving for p, Initial residual = 0.00181342, Final residual = 1.74422e-05, No Iterations 2
time step continuity errors : sum local = 7.46994e-07, global = 7.68894e-08, cumulative = -1.60508e-05
GAMG:  Solving for p, Initial residual = 2.19072e-05, Final residual = 9.34896e-07, No Iterations 4
time step continuity errors : sum local = 4.0072e-08, global = 1.10344e-09, cumulative = -1.60497e-05
smoothSolver:  Solving for epsilon, Initial residual = 0.482851, Final residual = 6.94855e-10, No Iterations 2
smoothSolver:  Solving for k, Initial residual = 0.586274, Final residual = 4.68515e-07, No Iterations 2
ExecutionTime = 0.07 s  ClockTime = 6 s

Running code execute
    sum(inlet) of phi = -1.73889e-05
    sum(outlet1) of phi = -1.66022e-05
    sum(outlet2) of phi = 3.39913e-05

score: 0.00265769	1.88353e-10

Courant Number mean: 0.0209096 max: 0.0563732

出力の最初の方では、system/controlDictに記載したC++をコンパイルしています。
　そのあとの出力は前述のsurfaceFieldValueと同じ内容を出力していることが確認できました。なおついでにその下にはscore: 時間 スコアのように出力しています。

最後に

今回は導入編としてfunctionObjectで処理を追加してみました。
　functionObjectはコンパイルの手間なしでC++に触れることができる良い機能なので、個人的にはC++の入り口としてとても良いのではないかと思っています。ぜひ試してみてください。
　次回以降はこれと同じ処理を別の手段で行っていきます。次回はまずソルバーに追加してやってみます。