【OpenFOAM TIPS】objectRegistryにオブジェクトを登録する

はじめに

　よく忘れるので備忘録としてまとめます。
　OpenFOAMの計算実行時には、objectRegistryオブジェクトがメモリ内のオブジェクトを制御しています。ここに登録してあるものには色んな場所からアクセスることができます。この記事では、coded境界条件の初回呼び出し時にフィールドを登録し、以降は同一のそれを呼び出せるようにする、という事をやっていきます。

計算環境

• OS : Ubuntu 20.04 (WSL / Windows 10 Home)
• OpenFOAM : v2012

メリット

　以下のようなメリットが考えられます。

• 一つのオブジェクトを色々な場所から共有できる
• 出力のタイミングを他のオブジェクトと合わせることができる
• 場として出力できればParaViewで他のフィールドと同様に可視化できる

　このやり方はcoded等の書き捨てコードに限らず、例えば独自の境界条件を作る場合などにも有効です。例えば速度と圧力の境界条件間で何かオブジェクトを渡したい場合に、計算を制御するobjectRegistryにそのオブジェクトを登録することでメモリ内で実現できたりします。

例として用いたケース

　ここでは例として、inletの質量流量をvolScalarFieldクラス"massInlet"の境界値として各時刻で保存する場合を考えます。

　以下のような脈動する水と油の二相流を例に説明します。（※本記事ではこの計算自体には触れません。）

　なおここで使用したケースはGithub(github.com/inabower/objectRegistryTest)で公開しています。

脈動を表す境界条件

　Uでは以下のようにcodedFixedValueで脈動を設定しているとします。内容としては時間によってx方向の速度が脈動するような条件です。上下の位置で分布を持たせたりもしています。

boundaryField
{
    inlet
    {
        type            codedFixedValue;
        name           wave;
        value           uniform (0 0 0);
        
        code
        #{
            const fvMesh& mesh = patch().boundaryMesh().mesh();
            const Time& t = mesh.time();

            vectorField ans(patch().size(), Zero);
            forAll(ans, fi)
            {
                ans[fi].x() = 5.0* (1.0-cos(t.value()*3.141592*2.0))*(0.25-pow(patch().Cf()[fi].y(), 2.0));
            }
            operator==(ans);
        #};
    }

    ...(他の境界)
}

　ここのcodedの中にいろいろ追加していきます。具体的にはmassInletというフィールドを作成し、その境界値に密度と速度を掛けた値を追加することを考えます。

ダメな例

　通常のcreateFields.Hのようにこの中でvolScalarFieldを定義することをやってみます。

source

boundaryField
{
    inlet
    {
        type            codedFixedValue;
        name           wave;
        value           uniform (0 0 0);
        
        code
        #{
            const fvMesh& mesh = patch().boundaryMesh().mesh();
            const Time& t = mesh.time();

            vectorField ans(patch().size(), Zero);
            forAll(ans, fi)
            {
                ans[fi].x() = 5.0* (1.0-cos(t.value()*3.141592*2.0))*(0.25-pow(patch().Cf()[fi].y(), 2.0));
            }
            operator==(ans);

            const word massInletName = "massInlet";
            volScalarField massInlet
            (
                IOobject
                (
                    massInletName,
                    mesh.time().timeName(),
                    mesh,
                    IOobject::NO_READ,
                    IOobject::AUTO_WRITE
                ),
                mesh,
                dimensionedScalar(dimless, Zero)
            );

            scalarField& massInletRef = massInlet.boundaryFieldRef()[patch().index()];

            const scalarField& alpha = patch().lookupPatchField<volScalarField, scalar>("alpha.water");;

            forAll(massInletRef, fi)
            {
                massInletRef[fi] = (alpha[fi] * 1000 + (1.0 - alpha[fi]) * 800) * ans[fi].x();
            }
        #};
    }

    ...(他の境界)
}

　これだと一応計算は回りますが、massInletは結果としては出力されません。これはvolScalarField オブジェクトの寿命がこの関数内で終了していることに由来します。オブジェクトを保存して後のイタレーションでも使用するためにはobjectRegistryに登録する必要があります。

よさげな例

source

　以下のように、ポインタをobjectRegistryに保存すると、以降の計算でも同一のオブジェクトを呼び出すことができ、さらにcontrolDictなどで指定したwriteFileのルールに従い他のフィールドと同じタイミングで出力してくれます。

boundaryField
{
    inlet
    {
        type            codedFixedValue;
        name           wave;
        value           uniform (0 0 0);
        
        code
        #{
            const fvMesh& mesh = patch().boundaryMesh().mesh();
            const Time& t = mesh.time();

            vectorField ans(patch().size(), Zero);
            forAll(ans, fi)
            {
                ans[fi].x() = 5.0* (1.0-cos(t.value()*3.141592*2.0))*(0.25-pow(patch().Cf()[fi].y(), 2.0));
            }
            operator==(ans);

            fvMesh& meshRef = const_cast<fvMesh&>(mesh);
            const word massInletName = "massInlet";
            if(!meshRef.foundObject<volScalarField>(massInletName))
            {
                autoPtr<volScalarField> massInletPtr(
                    new volScalarField
                    (
                        IOobject
                        (
                            massInletName,
                            meshRef.time().timeName(),
                            meshRef,
                            IOobject::NO_READ,
                            IOobject::AUTO_WRITE
                        ),
                        meshRef,
                        dimensionedScalar(dimless, Zero)
                    )
                );
                meshRef.objectRegistry::store(massInletPtr);
            }
            const scalarField& massInletFld = patch().lookupPatchField<volScalarField, scalar>(massInletName);
            scalarField& massInletRef = const_cast<scalarField&>(massInletFld);

            const scalarField& alpha = patch().lookupPatchField<volScalarField, scalar>("alpha.water");;

            forAll(massInletRef, fi)
            {
                massInletRef[fi] = (alpha[fi] * 1000 + (1.0 - alpha[fi]) * 800) * ans[fi].x();
            }
        #};
    }

    ...(他の境界)
}

ここでは以下の動作を行っています。

• objectRegistryに"massInlet"が登録されていない場合
  • volScalarFieldクラスの"massInlet"を新しく生成しそのポインタをmassInletPtrとして取り出す
  • massInletPtrポインタを引数にしてobjectRegistryに登録する
• patch().lookupPatchField関数を用いて、objectRegistryに登録されているmassInletのpatchの値を呼び出す
• 内容を変更する

注意すべき点をいくつか書いておきます。

① 発動するタイミング

   if(!meshRef.foundObject<volScalarField>(massInletName))

　objectRegistryに登録するのは必ず一回にしなくてはいけません。そのためここではfound関数を用いて、すでに登録されていない場合のみstoreするようにしています。

② store関数

      autoPtr<volScalarField> massInletPtr( new volScalarField(io) );
      meshRef.objectRegistry::store(massInletPtr);

　登録するためのobjectRegistryクラスのstore関数を用います。ここではfvMeshがobjectRegistryクラスを継承しており、OpenFOAMの計算の根幹のobjectRegistryを呼び出すことができるため、meshRef.objectRegistry::store(massInletPtr);のようにしています。
　ちなみにautoPtr<>はOpenFOAMで用意しているスマートポインタです。可能な限りC++ネイティブのものではなくこれを使うようにした方が個人的には思っています。（OpenFOAMの作者の意図に反したエラーを回避できるため）

③ const_cast

    const fvMesh& mesh = patch().boundaryMesh().mesh();

    fvMesh& meshRef = const_cast<fvMesh&>(mesh);

ここ(fvPatchField)からは以下のようにしてfvMeshクラスのポインタにアクセスできます。

• const fvPatch& patch()
• const fvBoundaryMesh& patch().boundaryMesh()
• const fvMesh& patch().boundaryMesh().mesh()

ただこれらはconstであるため、const_castで内容を変更可能にします。

そのほか応用例

　空のオブジェクトを登録する操作を入れることで、最初の一回かどうかのフラグとして使用することもできます。 ここでは最初の１回目のイタレーションのみ速度0、それ以外は速度1にしています。

source

boundaryField
{
    inlet
    {
        type            codedFixedValue;
        name           wave;
        value           uniform (0 0 0);
        
        code
        #{
            vector ans(Zero);
            
            const fvMesh& mesh = patch().boundaryMesh().mesh();
            fvMesh& meshRef = const_cast<fvMesh&>(mesh);

            bool isFirstTime = false; //ここ
            if(!mesh.foundObject<IOdictionary>("firstIterFlag_"+patch().name()))
            {
                autoPtr<IOdictionary> flagPtr(
                    new IOdictionary
                    (
                        IOobject
                        (
                            "firstIterFlag_"+patch().name(),
                            meshRef.time().constant(),
                            meshRef,
                            IOobject::NO_READ,
                            IOobject::NO_WRITE
                        )
                    )
                );
                meshRef.objectRegistry::store(flagPtr);
                isFirstTime = true; //ここ
            }

            // 最初のイタレーションでない場合はx方向の速度を1.0
            if(!isFirstTime) ans.x() = 1.0;
            operator==(ans);
        #};
    }

    ...(他の境界)
}

　計算の安定化などで役立つこともあるかもしれません。

最後に

　もしも何かの役に立ったら幸いです。

参考

• 公式Doxygen
  • fvPatchField
  • fvMesh
  • objectRegistry