OpenFOAM マルチリージョンの取り扱い（境界条件と並列処理）

使用環境

Ubuntu 16.04
OpenFOAM 3.0.1

概要

マルチリージョンソルバーを自作/カスタマイズする際に、隣接regionとのやりとりはとても重要になる。
その辺についてメモしていく。

標準の境界条件

chtMultiRegion(SIMPLE)Foamでは温度について隣接regionとやりとりする境界条件が備わっている。
それが↓

• turbulentTemperatureCoupledBaffleMixed (GitHub)
• turbulentTemperatureRadCoupledMixed (GitHub)

2者の違いは輻射熱の無/有。
前者を基準に説明していく。

なおこれらは
temperatureCoupledBase (GitHub)
mixed (GitHub)
を継承したものである。

境界条件内で行っている計算

以下に前者境界条件の中のupdateCoeffs()関数を示す。
（スペース省略のためにコメントやデバッグ操作を省略）

ここでこの境界での

• refValue　：　参照値
• refGrad　：　参照勾配
• valueFraction　：　値と勾配の固定するバランス（0 - 1）

を計算する。
詳しくはmixed(documentation)もしくはロビン境界条件(wikipedia)を参照。

$FOAM_SRC/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/turbulentTemperatureCoupledBaffleMixed/turbulentTemperatureCoupledBaffleMixedFvPatchScalarField.C

void turbulentTemperatureCoupledBaffleMixedFvPatchScalarField::updateCoeffs()
{
    if (updated()){ return; }

    // 並列処理用のtagを一時的に逃がす。
    int oldTag = UPstream::msgType();
    UPstream::msgType() = oldTag+1;

    // 自patchをmappedPatchBaseにrefCast。これにより隣patchの情報を引き出せる。
    const mappedPatchBase& mpp = refCast<const mappedPatchBase>(patch().patch());      // refCastで型変換
    const polyMesh& nbrMesh = mpp.sampleMesh();                                 // 隣patchが所属するpolyMesh
    const label samplePatchI = mpp.samplePolyPatch().index();                   // 隣Mesh内での隣patchの番号
    const fvPatch& nbrPatch = refCast<const fvMesh>(nbrMesh).boundary()[samplePatchI]; // 隣fvPatch

    // 隣patchから自patch内でのrobin境界条件(mixed)を計算していく。

    // 隣patchの中の(名前==TnbrName_)のフィールドを探し出してここと同じ型にrefCastで変換
    const turbulentTemperatureCoupledBaffleMixedFvPatchScalarField& nbrField =
    refCast<const turbulentTemperatureCoupledBaffleMixedFvPatchScalarField> (
        nbrPatch.lookupPatchField<volScalarField, scalar>( TnbrName_ )
    );

    // 変換されたpatchから隣の値を取り出す(一度しか使わないのでtmp<>なのかな)
    tmp<scalarField> nbrIntFld(new scalarField(nbrField.size(), 0.0));
    tmp<scalarField> nbrKDelta(new scalarField(nbrField.size(), 0.0)); // まずはサイズを決めて中身0で生成

    if (contactRes_ == 0.0) {
        nbrIntFld() = nbrField.patchInternalField(); // 隣のセル値
        nbrKDelta() = nbrField.kappa(nbrField)*nbrPatch.deltaCoeffs(); // 隣の熱伝達率
    }
    else { // 境界熱抵抗がある場合には
        nbrIntFld() = nbrField;    // 隣の境界値
        nbrKDelta() = contactRes_; // 熱伝達率
    }

    mpp.distribute(nbrIntFld());
    mpp.distribute(nbrKDelta()); // 並列処理をする場合には必須

    tmp<scalarField> myKDelta = kappa(*this)*patch().deltaCoeffs(); // 自patchの熱伝達率

    this->refValue() = nbrIntFld(); // fixedValueに相当
    this->refGrad() = 0.0;          // zeroGradientに相当
    this->valueFraction() = nbrKDelta()/(nbrKDelta() + myKDelta()); // 上二つの割合を伝津伝達率から計算

    mixedFvPatchScalarField::updateCoeffs(); // 3つの値をmixedに渡す。

    // 並列処理用のtagを戻す
    UPstream::msgType() = oldTag;
}

中身の簡単な説明

熱計算

・Fld ：　この場合は温度。
・kappa()　：　熱伝導率。temperatureCoupledBaseから継承した関数。指定した熱モデルとkappaNameから熱伝導率を返す。
・KDelta　：　熱伝達率。熱伝導率と、境界とセルとの距離から算出される。
・contactRes　：　境界熱抵抗がある場合の熱伝達率。以下のようにコンストラクタの時点で計算される。

コンストラクタ内

if (dict.found("thicknessLayers"))
{
    dict.lookup("thicknessLayers") >> thicknessLayers_;
    dict.lookup("kappaLayers") >> kappaLayers_;

    if (thicknessLayers_.size() > 0) {
        // Calculate effective thermal resistance by harmonic averaging
        forAll (thicknessLayers_, iLayer) {
            contactRes_ += thicknessLayers_[iLayer]/kappaLayers_[iLayer];
        }
        contactRes_ = 1.0/contactRes_;
    }
}

境界層は複数指定できるようだ。
まず各層の（厚さ/熱伝導率）を計算して、それを合計して、最後に逆数にする。

並列処理用

・ UPstream::msgType()
勉強不足のためよく分からないです。
多分標準出力（Infoとか）の設定だと思うのですが、値を変えてみてもよく分かりませんでした。
きっと後日追記します。
　

・ mappedPatchBase::distribute(Field)

これを怠ったがために一時期ハマった。
単コア計算ではうまくいくけど並列になるとエラーが起きる、という厄介な事態になる。
具体的にどんな処理かというと、

regionを分割した際に、境界面の表裏で計算するprocessorが違うことがよくある。（自動的に分割するとこうなる。）
「processor0の中ではregionAのpatchXの面数は10だが、それに対応するregionBのpatchYの面数が20である」
こんな表裏の面数が違うという事態が起きるので、例えばforAllループやフィールド同士の計算や代入を行うとエラーとなる。（コンパイルは通る）
ここでdistribute()を使うと、その名の通りいい感じに分配してくれる。

distributeを使ってみる

例：自patchのAというフィールドの値だけジャンプする

myFvPatchScalarField.C

myFvPatchScalarField::updateCoeffs()
{
    ...

    const mappedPatchBase& mpp = refCast<const mappedPatchBase>(patch().patch());      // refCastで型変換
    const polyMesh& nbrMesh = mpp.sampleMesh();                                 // 隣patchが所属するpolyMesh
    const label samplePatchI = mpp.samplePolyPatch().index();                   // 隣Mesh内での隣patchの番号
    const fvPatch& nbrPatch = refCast<const fvMesh>(nbrMesh).boundary()[samplePatchI]; // 隣fvPatch

    const myFvPatchScalarField& nbrField =
    refCast<const myFvPatchScalarField> (nbrPatch.lookupPatchField<volScalarField, scalar>( nbrName_ ) );

    scalarField nbrIntFld = nbrField.patchInternalField();
    scalarField nbrKDelta = nbrField.kappa(nbrField)*nbrPatch.deltaCoeffs();

    scalarField jump = patch().lookupPatchField<volScalarField, scalar>("A"); // jumpする値を読み込む

    mpp.distribute(nbrIntFld);
    mpp.distribute(nbrKDelta);

    scalarField myKDelta = kappa(*this)*patch().deltaCoeffs();

    this->refValue() = nbrIntFld + jump; // distribute()しないと(nbrIntFld.size() != jump.size())になりエラー
    this->refGrad() = 0.0;
    this->valueFraction() = nbrKDelta/(nbrKDelta + myKDelta);

    mixedFvPatchScalarField::updateCoeffs();
    ...
}

自分はまだメモリ管理が下手なのでほとんどtmpは使わない。（パソコンさんごめんなさい）
普通のscalarFieldでも分配できる。

ソルバー内での境界間のやりとり

ソルバーの中で境界間のやり取りをする場合にもdistributeは必須になる。
#include "mappedPatchBase.H"を事前に設置するのを忘れずに。

例：境界間の温度差を取り出す。

mySolver.C

forAll(mesh.boundary(), patchI) {
    if(myPatch.type() != "mappedWall") continue; // このpatchが境界かどうかの判別

    const fvPatch& myPatch = mesh.boundary()[patchI];
    const scalarField& myT = T.boundaryField()[patchI]; // 自patchの境界値を取得

    const mappedPatchBase& mpp = refCast<const mappedPatchBase>(myPatch.patch());
    const polyMesh& nbrMesh = mpp.sampleMesh();
    const label samplePatchI = mpp.samplePolyPatch().index();
    const fvPatch& nbrPatch = refCast<const fvMesh>(nbrMesh).boundary()[samplePatchI];

    scalarField nbrT = nbrPatch.lookupPatchField<volScalarField, scalar>("T"); // 隣patchの値を取得（サイズ変更するので値渡し）
    mpp.distribute(nbrT);

    const scalarField deltaT = myT - nbrT; // 2者の演算ができるようになる
}

patchが境界かどうかの判別

上のコードの中の2行目の部分。
これをつけないと、普通のpatchをrefCastしようとしてエラーが起きる。
色々方法がある。

①　patchのcoupledを取得 (ぜんぜん違いました)
②　patchのtypeを取得　（上の方法）
③　patchと元メッシュとの関係を取得　（後述）

mySolver.C

forAll(mesh.boundary(), patchI) {
    const fvPatch& myPatch = mesh.boundary()[patchI];

    if(myPatch.type() != "mappedWall") continue; // (2) patchのtypeを取得
    if(grobalPatch[patchI] != -1) continue; // (3) patchと元メッシュとの関係を取得
}

おまけ：③について

これは正直メリットはないが、こんなこともできるよということで。

splitMeshRegionsでメッシュをマルチリージョンに分割した際に、constant/(regionName)/polyMeshに~RegionAddressingというファイルが生成される。
チュートリアルを実行してみると確認できると思う。
例えばregionAの"boundaryRegionAddressing"の中身が

適当なケース

//Header
...
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

5(2 3 4 -1 -1)

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

}

こんな感じになっているとする。
これはつまり、

• meshA.boundary()[0] <- 元mesh.boundary()[2]
• meshA.boundary()[1] <- 元mesh.boundary()[3]
• meshA.boundary()[2] <- 元mesh.boundary()[4]
• meshA.boundary()[3] は 境界（元メッシュではpatchではなかった）
• meshA.boundary()[4] は 境界（元メッシュではpatchではなかった）

という対応を示している。例えば

createFluidMesh.H

...
PtrList<labelIOList> fluidGrobalPatch(fluidNames.size());

forAll(fluidNames, i) {
    ...
    const fvMesh& mesh = fluidRegions[i];

    fluidGrobalPatch.set (
        i,
        new labelIOList (
            IOobject (
                "boundaryRegionAddressing",
                mesh.time().findInstance(mesh.meshDir(), "faces"),
                polyMesh::meshSubDir,
                mesh,
                IOobject::MUST_READ
            )
        )
     );
}

こんな感じでmeshと一緒に読み込んで

ソルバーのどこか.H

forAll(fluidNames, i) {
    const fvMesh& mesh = fluidRegions[i];

    forAll(mesh.boundary(), patchI) {
        const fvPatch& patch = mesh.boundary()[patchI];

        if(fluidGrobalPatch[i][patchI] == -1)
            Info << patch.name() << " is interface." << endl;
        }
        else {
            Info << patch.name() << " is surface." << endl;
        }
    } // patch loop end
} // fluid loop end

こんな感じに使う。