🐺

【OpenFOAM】手法を変えて処理を追加③ : 関数やクラスとして処理を追加する

2020/12/08に公開

C++

オープンCAEアドベントカレンダー 8日目のエントリーです。

はじめに

『【OpenFOAM】手法を変えて処理を追加② : ソルバーに処理を追加する』の続きです。
　OpenFOAMの計算に手法を変えて処理を追加していきます。前回はソルバーのmain()関数に処理を記述しました。
今回はその処理を関数もしくはクラスとして作成して、main()関数ではそれを呼び出すという方法を紹介します。

手法一覧

この記事は様々な手法を紹介するシリーズの一つであり、全体としては以下の方法で行う予定です。実行確認済みのファイルはGitHubで公開しています。

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

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

動作環境

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

関数で処理を追加する

まずは処理の関数化を行ってみます。これにより、処理を複数の箇所で呼び出す際等にコードが大きく簡略化されます。

下準備

以下のようなスクリプトを作成して実行します。これにより、元となるpimpleFoamとそのチュートリアルをコピーしています。以降ではこのスクリプトを実行したという前提で進めていきます。

prepareFuncCase.sh

#!/bin/bash

# (任意) 今回作成するディレクトリとソルバーとケースに名前をつける
dirName=$FOAM_RUN/003.class
solverName=scoreFuncPimpleFoam
caseName=funcCase

mkdir -p $dirName

# ソルバーのコピー
cp -r $FOAM_SOLVERS/incompressible/pimpleFoam $dirName/$solverName
rm -r $dirName/$solverName/overPimpleDyMFoam $dirName/$solverName/SRFPimpleFoam 
sed -i -e "s/APPBIN)\/pimpleFoam/USER\_APPBIN)\/${solverName}/g" $dirName/$solverName/Make/files
cp $FOAM_SRC/finiteVolume/cfdTools/general/CorrectPhi/CorrectPhi.H $dirName/$solverName/cpCorrectPhi.H
sed -i -e "s/CorrectPhi.H/cpCorrectPhi.H/" $dirName/$solverName/pimpleFoam.C
echo "Solver : $dirName/$solverName"

# チュートリアルのコピー
cp -r $FOAM_TUTORIALS/incompressible/pimpleFoam/RAS/TJunction $dirName/$caseName
sed -i -e "s/pimpleFoam/${solverName}/" ${dirName}/${caseName}/system/controlDict
# 元々の1.5秒を計算する必要は無いので減らす
sed -i -e "s/^endTime.*/endTime         0.3;/" ${dirName}/${caseName}/system/controlDict
# functionObjectの記述を削除
startLine=`cat ${dirName}/${caseName}/system/controlDict | grep -n "^functions" | sed -e 's/:.*//g'`
endLine=`cat ${dirName}/${caseName}/system/controlDict | wc -l`
sed -i -e "${startLine},${endLine}d" ${dirName}/${caseName}/system/controlDict
echo "Case   : $dirName/$caseName"

関数の定義

関数の定義はmain()関数の前に記述します。以下のようにpimpleFoam.Cの#includeとmain()の間に関数がどのように動作するかを記述します。
　以下ではgetScoreという名前の関数を定義しています。fvMeshクラスを引数とし、出入り口の流量の総和をscalarクラスとして返すように定義しています。

#include "fvCFD.H"
#include "dynamicFvMesh.H"
#include "singlePhaseTransportModel.H"
#include "turbulentTransportModel.H"
#include "pimpleControl.H"
#include "cpCorrectPhi.H"
#include "fvOptions.H"

// 85行目付近 ここを追加
scalar getScore(const fvMesh& mesh_)
{
    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);
    }

    return score_;
}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{

関数の呼び出し

この定義した関数を呼び出すために、同じくpimpleFoam.Cのmain()関数の中で以下のように記述します。

            if (pimple.turbCorr())
            {
                laminarTransport.correct();
                turbulence->correct();
            }
        }

        // 187行目付近。この2行を追加する。
        const scalar score = getScore(mesh);
        Info << nl << "score: " << runTime.value() << tab << score << nl << endl;

        runTime.write();

これにより、scoreという変数にgetScore()関数で計算された流量総和が格納され、そのあとにInfoで出力しています。

コンパイルとテスト

以下のようにコンパイルを行います。

$ cd $FOAM_RUN/003.class/scoreFuncPimpleFoam
$ wclean
$ wmake

エラーが出力されなければ以下のようにソルバーを実行してみます。

$ cd $FOAM_RUN/003.class/funcCase
$ foamCleanTutorials
$ foamRunTutorials

以下のようにlog.scoreFuncPimpleFoamにscoreが出力されていることを確認します。

$ cat $FOAM_RUN/003.class/funcCase/log.scoreFuncPimpleFoam | grep -e "score:" | sed "s/score: //"
0.00120482      7.00171e-10
0.00265769      1.88353e-10
0.00439595      -7.65986e-10
0.00647429      2.55232e-09
0.0089355       3.27205e-09
0.0118731       -1.48406e-09
0.0153981       2.49265e-09
0.0196282       2.54628e-09
0.0246515       -9.17793e-09
0.0304475       -6.57753e-09
0.0374028       -4.7443e-09
0.0452274       5.67256e-09
0.0543562       -5.79236e-09
0.0619635       -3.46103e-09
0.0695708       -4.28805e-09
0.0756566       -5.23626e-09
0.0817425       2.75405e-09
...

上のように出力されていればOKです。

クラスで処理を追加する

次は処理を行うためのクラスを作成してみましょう。

下準備

以下のスクリプトを実行して先ほどとは別の名前のソルバー(scoreClassPimpleFoam)とケース(funcCase)を作成します。先ほどのスクリプトと違いのは最初の方だけです。

prepareClassCase.sh

#!/bin/bash

# (任意) 今回作成するディレクトリとソルバーとケースに名前をつける
dirName=$FOAM_RUN/003.class
solverName=scoreClassPimpleFoam #scoreFuncPimpleFoam
caseName=classCase #funcCase

mkdir -p $dirName

# ソルバーのコピー
cp -r $FOAM_SOLVERS/incompressible/pimpleFoam $dirName/$solverName
rm -r $dirName/$solverName/overPimpleDyMFoam $dirName/$solverName/SRFPimpleFoam 
sed -i -e "s/APPBIN)\/pimpleFoam/USER\_APPBIN)\/${solverName}/g" $dirName/$solverName/Make/files
cp $FOAM_SRC/finiteVolume/cfdTools/general/CorrectPhi/CorrectPhi.H $dirName/$solverName/cpCorrectPhi.H
sed -i -e "s/CorrectPhi.H/cpCorrectPhi.H/" $dirName/$solverName/pimpleFoam.C
echo "Solver : $dirName/$solverName"

# チュートリアルのコピー
cp -r $FOAM_TUTORIALS/incompressible/pimpleFoam/RAS/TJunction $dirName/$caseName
sed -i -e "s/pimpleFoam/${solverName}/" ${dirName}/${caseName}/system/controlDict
# 元々の1.5秒を計算する必要は無いので減らす
sed -i -e "s/^endTime.*/endTime         0.3;/" ${dirName}/${caseName}/system/controlDict
# functionObjectの記述を削除
startLine=`cat ${dirName}/${caseName}/system/controlDict | grep -n "^functions" | sed -e 's/:.*//g'`
endLine=`cat ${dirName}/${caseName}/system/controlDict | wc -l`
sed -i -e "${startLine},${endLine}d" ${dirName}/${caseName}/system/controlDict
echo "Case   : $dirName/$caseName"

".C"と".H"

OpenFOAMではクラスを作成するときに、宣言と定義を分けることが慣習化されているため、今回もそれに則って作成します。つまり、以下のような構成になっています。

同名の.Cと.Hを作成 (例：getScore.CとgetScore.H)
宣言：.Hでは『どんな関数変数引数があるか』を記述する。
定義：.Cでは『関数等がどのように動作するか』を記述する。
ソルバーでは.Hをincludeする。
.Cはソルバーとは別にコンパイルされてソルバーとリンクされる

これらがどのようにコンパイルされるのかなどについてはQiita(【OpenFOAM】wmakeって何してるの)でまとめたことがありますのでこちらを参照ください。
　それでは、この構成を踏襲して、まずはgetScore.HとgetScore.Cを作成していきます。

ディレクトリ構造

今回は最終的に以下のようなディレクトリ構成を目指します。

├── getScore
│   ├── getScore.C
│   └── getScore.H
├── Make
│   ├── files
│   └── options
├── UEqn.H
├── correctPhi.H
├── cpCorrectPhi.H
├── createFields.H
├── pEqn.H
└── pimpleFoam.C

このため、以下のようにソルバーのディレクトリにgetScoreディレクトリを作成します。

$ mkdir -p $FOAM_RUN/003.class/scoreClassPimpleFoam/getScore

この下にgetScore.CとgetScore.Hを作成していきます。

getScore/getScore.H

まずは以下のようなgetScore/getScore.Hを作成します。

#ifndef getScore_H
#define getScore_H

#include "fvMesh.H"

namespace Foam
{
class getScore
{
protected:
    // Protected data
        const fvMesh& mesh_;
        scalar score_;
public:
    // Constructors
        getScore(const fvMesh& mesh);
    //- Destructor
        ~getScore()
        {}
    // Member Functions
        void calculate();

        const scalar& value() const
        {
            return score_;
        }
};
}

#endif

まずこのクラスの骨格としては以下のようになっています。

namespace Foam
{
    class getScore
    {
    public:
        // Constructors
            getScore(const fvMesh& mesh);
        //- Destructor
            ~getScore()
            {}
    };
}

FoamというOpenFOAMのnamespaceの中で、getScoreという名前のクラスが宣言されています。またコンストラクタが宣言されており、デストラクタについては{}によって『何もしない』という定義もされています。

また、以下の部分ではこのクラスの中で保持される変数が定義されています。

protected:
    // Protected data
        const fvMesh& mesh_;
        scalar score_;

protectedであるため、このクラスの中でのみmesh_とscore_を呼び出すことができます。

以下の部分では関数が宣言されています。

public:
    // Member Functions
        void calculate();

        const scalar& value() const
        {
            return score_;
        }

これらの関数はpublicであるため、クラスの外、つまりソルバーなどから呼び出すことができます。今回はcalculate()でスコアの計算を行いvalue()でその値を返すという構成にしています。なおvalue()は単純な関数であるためここで定義まで行っています。

また、このクラスではfvMeshクラスを変数として使用するため、fvMeshクラスの宣言を読み込むために以下のように最初に記述されています。

#include "fvMesh.H"

以下の部分は、他の場所でgetScoreクラスを宣言する記述があった場合に重複を避けるための記述です。この場合先に読み込まれたgetScoreが使われることになります。ただしこれは念のためのものであり、本来は名前の重複は無い方がいいので、名前の重複は無いようにしましょう。

#ifndef getScore_H
#define getScore_H

...

#endif

getScore/getScore.C

次に以下のように定義部分を作成します。

#include "getScore.H"
#include "surfaceFields.H"

Foam::getScore::getScore
(
    const fvMesh& mesh
)
:
    mesh_(mesh)
{
}

void Foam::getScore::calculate()
{

    Info << "Running calculate()" << endl;

    const surfaceScalarField& phi = mesh_.lookupObject<surfaceScalarField>("phi");
    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);
    }
}

このクラスにはコンストラクタ、デストラクタ、関数calculate()、関数value()がありますが、デストラクタと関数value()についてはgetScore.Hで既に定義されていますので、残りのコンストラクタと**関数calculate()**についてgetScore.Cで定義しています。

まず以下のようにコンストラクタを定義しています。

Foam::getScore::getScore
(
    const fvMesh& mesh
)
:
    mesh_(mesh)
{
}

fvMeshクラスのポインタを受け取り、クラス内変数のmesh_の定義としています。これによりソルバー側のメッシュにアクセスできるようになります。OpenFOAMの場合はfvMeshにアクセスできれば大抵何でもできますが、もし他にアクセスしたいソルバー側の変数や情報等があれば、同じように引数と変数を作成することになります。

以下のcalculate()では、今回の目的である『出入り口の流量の総和を求める』という処理を行っています。返し値はvoidにしていますので、この関数自体は何も返さず、『クラス内変数のscore_を更新する』という事を行っています。

void Foam::getScore::calculate()
{

    Info << "Running calculate()" << endl;

    const surfaceScalarField& phi = mesh_.lookupObject<surfaceScalarField>("phi");
    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);
    }
}

最初の#includeではこのgetScoreクラスの宣言であるgetScore.Hと、calculate()内で登場するsurfaceScalarFieldの宣言がされているsurfaceFields.Hを読み込んでいます。calculate()で他に使用したいクラスがある場合は同様に宣言を読み込む必要があります。

#include "getScore.H"
#include "surfaceFields.H"

pimpleFoam.C

以上のように作成したgetScore.CとgetScore.Hを呼び出すためにmain()関数の中で呼び出します。そのためには以下のようにpimpleFoam.Cを3か所変更します。

pimpleFoam.C:110行目付近

（前略）
#include "fvCFD.H"
#include "dynamicFvMesh.H"
#include "singlePhaseTransportModel.H"
#include "turbulentTransportModel.H"
#include "pimpleControl.H"
#include "cpCorrectPhi.H"
#include "fvOptions.H"

// 85行目付近
#include "getScore.H" // ★これを追加

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    #incl...（略）
    #include "CourantNo.H"
    #include "setInitialDeltaT.H"

    // 110行目付近
    getScore myScore(mesh); // ★これを追加

    turbulence->validate();

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    Info<< "\nStarting time loop\n" << endl;

    while (runTime.run())
    {
（中略）

            if (pimple.turbCorr())
            {
                laminarTransport.correct();
                turbulence->correct();
            }
        }

        // ★この2行を追加
        myScore.calculate();
        Info << nl << "score: " << runTime.value() << tab << myScore.value() << nl << endl;

        runTime.write();

        runTime.printExecutionTime(Info);
    }
（後略）

以下のように最初にgetScoreクラスのオブジェクトをmyScoreという名前で作成しています。

        getScore myScore(mesh);

そしてcalculate()関数でmyScoreの中のscore_変数を更新して、

        myScore.calculate();

value()関数でそのscore_変数を呼び出しています。

        Info << nl << "score: " << runTime.value() << tab << myScore.value() << nl << endl;

というわけで以下のような構成になりました。

getScore.Cのコンパイル

このままではgetScore.Cに関する情報は読み込まれません。getScore.Cをコンパイルするためには以下のようにMake/ディレクトリの内容を変更します。

Make/files

getScore/getScore.C # ★ここを追加
pimpleFoam.C

EXE = $(FOAM_USER_APPBIN)/scoreClassPimpleFoam

Make/options

EXE_INC = \
    -IgetScore \ # ★ここを追加
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude \
    -I$(LIB_SRC)/sampling/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
    -I$(LIB_SRC)/transportModels \
    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
    -I$(LIB_SRC)/dynamicMesh/lnInclude \
    -I$(LIB_SRC)/dynamicFvMesh/lnInclude

EXE_LIBS = \
    -lfiniteVolume \
    -lfvOptions \
    -lmeshTools \
    -lsampling \
    -lturbulenceModels \
    -lincompressibleTurbulenceModels \
    -lincompressibleTransportModels \
    -ldynamicMesh \
    -ldynamicFvMesh \
    -ltopoChangerFvMesh \
    -latmosphericModels

コンパイルとテスト

以下のようにコンパイルを行います。

$ cd $FOAM_RUN/003.class/scoreClassPimpleFoam
$ wclean
$ wmake

エラーが出力されなければ以下のようにソルバーを実行してみます。

$ cd $FOAM_RUN/003.class/classCase
$ foamCleanTutorials
$ foamRunTutorials

以下のようにlog.scoreFuncPimpleFoamにscoreが出力されていることを確認します。

$ cat $FOAM_RUN/003.class/classCase/log.scoreClassPimpleFoam | grep -e "score:" | sed "s/score: //"
0.00120482      7.00171e-10
0.00265769      1.88353e-10
0.00439595      -7.65986e-10
0.00647429      2.55232e-09
0.0089355       3.27205e-09
0.0118731       -1.48406e-09
0.0153981       2.49265e-09
0.0196282       2.54628e-09
0.0246515       -9.17793e-09
0.0304475       -6.57753e-09
0.0374028       -4.7443e-09
0.0452274       5.67256e-09
0.0543562       -5.79236e-09
0.0619635       -3.46103e-09
0.0695708       -4.28805e-09
0.0756566       -5.23626e-09
0.0817425       2.75405e-09
...

上のように出力されていればOKです。

最後に

今回は関数とクラスを追加する方法を紹介しました。

処理が複雑になって来た場合はこのように関数やクラスとして作成することで、管理や変更などが行いやすくなってきます。
　次回はこのクラスについて、継承という機能を紹介します。

はじめに

手法一覧

動作環境

関数で処理を追加する

下準備

関数の定義

関数の呼び出し

コンパイルとテスト

クラスで処理を追加する

下準備

".C"と".H"

ディレクトリ構造

getScore/getScore.H

getScore/getScore.C

pimpleFoam.C

getScore.Cのコンパイル

コンパイルとテスト

最後に

Discussion