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

# 使用環境
Ubuntu 16.04
OpenFOAM 3.0.1

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

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

- **turbulentTemperatureCoupledBaffleMixed** ([GitHub](https://github.com/OpenFOAM/OpenFOAM-3.0.x/tree/master/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/turbulentTemperatureCoupledBaffleMixed))
- **turbulentTemperatureRadCoupledMixed** ([GitHub](https://github.com/OpenFOAM/OpenFOAM-3.0.x/tree/master/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/turbulentTemperatureRadCoupledMixed))

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

なおこれらは
**temperatureCoupledBase** ([GitHub](https://github.com/OpenFOAM/OpenFOAM-3.0.x/tree/master/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/temperatureCoupledBase))
**mixed** ([GitHub](https://github.com/OpenFOAM/OpenFOAM-3.0.x/tree/master/src/finiteVolume/fields/fvPatchFields/basic/mixed))
を継承したものである。

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

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

ここでこの境界での

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

を計算する。
詳しくはmixed([documentation](http://cpp.openfoam.org/v4/a01549.html))もしくはロビン境界条件([wikipedia](https://ja.wikipedia.org/wiki/%E3%83%AD%E3%83%93%E3%83%B3%E5%A2%83%E7%95%8C%E6%9D%A1%E4%BB%B6))を参照。

```c++:$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　：　境界熱抵抗がある場合の熱伝達率。以下のようにコンストラクタの時点で計算される。

```c++:コンストラクタ内
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というフィールドの値だけジャンプする

```c++: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"`を事前に設置するのを忘れずに。

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

```c++: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と元メッシュとの関係を取得　（後述）


```c++: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"の中身が

```c++:適当なケース
//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ではなかった）

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

```c++: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と一緒に読み込んで

```c++:ソルバーのどこか.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
```

こんな感じに使う。