白黒羊

OpenFOAM の adjoint Shape Optimization Foam を考える #05ソルバーの改造

2020/07/19

これまで既存の実装を色々と見てきましたが、ここからはソルバーを書き換えて、境界条件を変更することで求めるコスト関数を設定するフェーズに移ろうと思います。

前回は「OpenFOAM の adjoint Shape Optimization Foam を考える #04境界条件の実装 」でした。

  1. OpenFOAM の adjoint Shape Optimization Foam を考える #01手法
  2. OpenFOAM の adjoint Shape Optimization Foam を考える #02境界条件
  3. OpenFOAM の adjoint Shape Optimization Foam を考える #03実装
  4. OpenFOAM の adjoint Shape Optimization Foam を考える #04境界条件の実装
  5. OpenFOAM の adjoint Shape Optimization Foam を考える #05ソルバーの改造 （今回）

ソルバーの作成

今回adjointShapeOptimizationFoamを書き換えて自作ソルバーを作るにあたって、元のデータをいじらないように、$WMPROJECTUSERDIR 内にadjointShapeOptimizationFoam をコピーしておきましょう。
ここではコピーしたソルバーに、myAdjointShapeOptimizationFoamという名前を付けて呼びます。
諸々書き換え終わったらコンパイルなどもする必要があるのですが本記事では割愛します。

境界条件の変更

流出部分での均質流

全圧損失を低減させる目的で境界条件が設定されているのを書き換えて流出部分が均質流になるような設定に変更してみます。
これまで何度か見てきたように、その場合の流出部分の境界条件は以下のようになります。

[1]　$q = \textbf{u}\cdot\textbf{v} + u_n v_n + \nu (\textbf{n}\cdot\nabla ) u_n + \frac{\partial J_\Gamma}{\partial v_n}$

[2]　$0 = v_n \textbf{u}_t + \nu (\textbf{n}\cdot\nabla ) \textbf{u}_t + \frac{\partial J\Gamma}{\partial \textbf{v}_t}$

流出部分(patch)であることを示す添字 ${}_p$ と、その隣接ノードであることを示す添字 ${}_neigh$ 、さらにそれぞれ接線要素と法線要素であることを示す添字 ${}_t$、${}_n$ を用いて表現すると、

[3]　$q = \textbf{u}_p\cdot\textbf{v}_p + u_{p,n} v_{p,n} + \frac{\nu}{\Delta} (u_{p,n}-u_{neigh,n}) + c(v_{p,n}-v_{p,n}^d)$

[4]　$\textbf{u}_t =\frac{\frac{\nu}{\Delta} \textbf{u}{neigh,t}-c(\textbf{v}_{p,n}-\textbf{v}{p,n}^d)}{v_{p,n}+\frac{\nu}{\Delta}}$

ここで、元々の $\nu (\textbf{n}\cdot\nabla ) u_n$ 項が、対象パッチである流出部分とその隣接ノードの速度差を領域間の距離 $\Delta$ で割ったものとして近似されていることに注意します。
また、以下のような近似も行うことができます。

[5]　$\nu (\textbf{n}\cdot\nabla ) \textbf{u}_t \sim \nu\frac{\textbf{u}_{p,t}-\textbf{u}_{neigh,t}}{\Delta}$

adjoint圧力

ここからは、「applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointOutletPressure/adjointOutletPressureFvPatchScalarField.C」からコピーしているadjointOutletPressureFvPatchScalarField.Cを編集して境界条件を書き換えます。
なお、以下では同様に、単純にコピーしたファイルを指してadjointOutletPressureFvPatchScalarField.C と呼びます。
元のソルバーは編集しないでください。
まずは、turbulentTransportModel.Hを呼んでおきます。

#include "adjointOutletPressureUniFvPatchScalarField.H"
#include "addToRunTimeSelectionTable.H"
#include "fvPatchMapper.H"
#include "volFields.H"
#include "surfaceFields.H"
#include "turbulentTransportModel.H"

次に、メンバ関数がまとめられている箇所を見てupdateCoeffs()関数を書き換えます。

// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

void Foam::adjointOutletPressureUniFvPatchScalarField::updateCoeffs()
{
    if (updated())
    {
        return;
    }

    const fvsPatchField<scalar>& phip =
        patch().lookupPatchField<surfaceScalarField, scalar>("phi"); // 削除

    const fvsPatchField<scalar>& phiap =
        patch().lookupPatchField<surfaceScalarField, scalar>("phia"); // 削除

    const fvPatchField<vector>& Up =
        patch().lookupPatchField<volVectorField, vector>("U");

    const fvPatchField<vector>& Uap =
        patch().lookupPatchField<volVectorField, vector>("Ua");

フラックスを使わないので、phipとphiapを削除します。

    scalarField Up_n = Up&patch().nf(); // 速度の法線要素の計算
    scalarField Uap_n = Uap&patch().nf(); // adjoint速度の法線要素の計算

また、速度とadjoint速度の法線ベクトルの絶対値を計算するためにpatch().nf() メンバ関数を用います。

    const incompressible::RASModel& rasModel =
        db().lookupObject<incompressible::RASModel>("turbulenceProperties");
    scalarField nueff = rasModel.nuEff()().boundaryField()[patch().index()];

式[1] の右辺の第1項と第2項を計算するためにRAS Modelを使います。
乱流粘性については、lookupObject()関数を使用してIOobjectである乱流モデルを検索して求めます。
有効粘度nueffは、上記のコードのようにして乱流モデルのクラスのメンバ関数を使用して求めることができます。

    const scalarField& deltainv = patch().deltaCoeffs(); // distance^(-1)

ノード間距離 $\Delta$ は、deltaCoeffs()関数を用いて求めます。
この関数は、ノード間の距離の逆数を返します。

    scalarField Uaneigh_n = (Uap.patchInternalField()&patch().nf());

patchInternalField()関数を使うことで、隣接ノードにおけるadjoint速度の法線要素の絶対値を求めることができます。
これはadjoint速度の流出部分の境界条件に用いることができますが、実際には使われることはありません。

残りの変数は、流出部分の速度です。
各caseの設定を決めるときにユーザがDictで与えることができるようにしておきます。
そのためには、速度、adjoint速度とは別のvolVectorFieldをcrateField.Hに設定する必要があります。
crateField.Hを開いて、次の行を追加します。

dictionary optFunc =
    mesh.solutionDict().subDict("objectiveFunctionDict");

vector Udlookup =
    optFunc.lookupOrDefault<vector>("Uduserdefnodim" , vector(1,0,0));
volVectorField Ud
(
    IOobject
    (
        "Ud",
        runTime.timeName(),
        mesh,
        IOobject::READ_IF_PRESENT,
        IOobject::AUTO_WRITE
    ),
    mesh,
    Udlookup
);

objectiveFunctionDictというsubDict内に Uduserdefnodim を定義しました。
保存して、adjointOutletPressureFvPatchScalarField.Cに戻ります。

    const fvPatchField<vector>& Udp =
        patch().lookupPatchField<volVectorField, vector>("Ud");
    scalarField Udp_n = (Udp&patch().nf());

今定義したベクトルの法線要素の絶対値は、上記のコードを追加することで確認できます。

    operator==((Uap&Up) + (Up_n*Uap_n) + nueff*deltainv*(Uap_n - Uaneigh_n) + (Up_n - Udp_n));
    fixedValueFvPatchScalarField::updateCoeffs();

演算子を定義します。
adjoint圧力の値を計算するために必要な変数をまとめています。

なお、ここでの定数 $c$ は1であり、単位は省略されていることに注意してください。

adjoint速度

adjoint速度の境界条件も同様に設定します。
adjoint圧力の境界条件の設定のときと全く同じように、adjointOutletVelocityFvPatchVectorField.Cを開いて、ヘッダーにturbulentTransportModel.Hを追加します。

// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

// Update the coefficients associated with the patch field
void Foam::adjointOutletVelocityUniFvPatchVectorField::updateCoeffs()
{
    if (updated())
    {
        return;
    }

    const fvsPatchField<scalar>& phiap =
        patch().lookupPatchField<surfaceScalarField, scalar>("phia"); // 削除

    const fvPatchField<vector>& Up =
        patch().lookupPatchField<volVectorField, vector>("U");

    scalarField Un(mag(patch().nf() & Up)); // 削除
    vectorField UtHat((Up - patch().nf()*Un)/(Un + small)); // 削除

    vectorField Uan(patch().nf()*(patch().nf() & patchInternalField())); // 削除

updateCoeffs()関数の内部のadjointフラックスphiap、スカラー値Ua、UtHat、Uanを削除します。

    const fvPatchField<vector>& Uap =
        patch().lookupPatchField<volVectorField, vector>("Ua");

adjoint速度としてUapを導入します。

    const incompressible::RASModel& rasModel =
        db().lookupObject<incompressible::RASModel>("turbulenceProperties");
    scalarField nueff = rasModel.nuEff()().boundaryField()[patch().index()];
    const scalarField& deltainv = patch().deltaCoeffs(); // dist^(-1)

乱流モデル（RAS Model）を設定し、乱流粘性nueffとノード間の距離の逆数deltainvを定義します。

    const fvPatchField<vector>& Udp =
        patch().lookupPatchField<volVectorField, vector>("Ud");
    vectorField Udp_n = (Udp& patch().nf())*patch().nf();
    vectorField Udp_t = Udp - Udp_n;

adjoint速度の設定には法線要素の絶対値ではなく接線要素を用います。
まず、$v_n^d$ に $\textbf{n}$を乗じて法線要素を計算し、それを$\textbf{v}_t^d=\textbf{v}^d-\textbf{v}_n^d$ のようにadjoint速度全体から引くことで接線要素$\textbf{v}_t^d$ を得ます。
この添字${}^d$ は、ユーザー定義(user defined)の「d」です。

    // Primal
    vectorField Up_n = (Up& patch().nf())*patch().nf(); // Normal
    scalarField Up_ns = (Up& patch().nf()); // Mag.of normal
    vectorField Up_t = Up - Up_n; // Tangential

同様に、速度の法線要素、接線要素を求めます。

    // Adjoint
    vectorField Uaneigh = Uap.patchInternalField();
    vectorField Uaneigh_n = (Uaneigh& patch().nf())*patch().nf(); // Normal
    vectorField Uaneigh_t = Uaneigh - Uaneigh_n; // Tangential

また、隣接ノードのadjoint速度とその両要素も求めます。

    vectorField Uap_n = (Uap& patch().nf())*patch().nf(); //Normal comp.
    vectorField Uap_t = (-(Up_t-Udp_t) + nueff*deltainv*Uaneigh_t) /
        (Up_ns+nueff*deltainv);
    operator==(Uap_t+Uap_n);

    fixedValueFvPatchVectorField::updateCoeffs();
}

これで式[2]の右辺にある、adjoint速度の接線要素を計算するために必要なすべての項の定義ができました。
演算子==の定義は、接線要素と法線要素を合わせた合計のadjoint速度の式である必要があります。
そのため、法線要素の前のステップの値を取得して、接線要素の現在のステップの値と足し合わせています。

各セルの感度の設定

各セルのコスト関数に対する感度について考えてみます。
今回設定したような流出部分の流れを均一にするコスト関数を考慮する場合では、流出部分での速度が過度に速いセルはporosityを増加させることによって影響力が抑えられます。
しかし、これは必ずしも設計上の目的からは望ましくない状況を引き起こすこともあります。
その場合には、porosityを変化させるのではなく感度を更新することができます。

各セルの感度を確認するためには、createFields.Hを修正し、myAdjointShapeOptimizationFoam.C内のループで値を更新します。

volScalarField sens
(
    IOobject
    (
    "sensitivity",
    runTime.timeName(),
    mesh,
    IOobject::READ_IF_PRESENT,
    IOobject::AUTO_WRITE
    ),
    Ua&U
);

createFields.H の末尾にvolScalarField sensを導入しました。

laminarTransport.correct();
turbulence->correct();

sens=Ua&U; // （追加）sensの更新

runTime.write();

そして、myAdjointShapeOptimizationFoam.Cのメインループ内に上記の1行を追加することで、毎ステップsensを更新します。
この値は、caseフォルダ内のタイムステップ名が付いた結果フォルダ内にsensitivityとして書き出されます。

コスト関数の表示

ここまではソルバーの計算部分の改造でしたが、ここからは分かりやすさのための表示の改良を行います。

コスト関数がどれだけ最小化されていくのかを可視化するためにmyAdjointShapeOptimizationFoam.CのSIMPLEループ内にコードを追加します。
一度だけ初期化する必要があるオブジェクトは、createFields.H ファイルで設定する必要があります。

解析対象となるcaseには、異なる名前の複数の流出部分が存在する場合があるため、コスト関数の計算時に使用する流出面をcaseフォルダで設定できるようにします。
また、計算の対象となる流出面の数もcaseフォルダのDictで定義します。

dictionary optFunc =
    mesh.solutionDict().subDict("objectiveFunctionDict");

int nObjPatch = optFunc.lookupOrDefault<scalar>("numberObjectivePatches", 0);
int objFunction = optFunc.lookupOrDefault<scalar>("objectiveFunction", 0);
wordList objPatchNames = optFunc.lookup("objectivePatchesNames");

流出部分の速度を決めるのと同じDict fvSolutionから値を持ってくることができるように、createFields.Hに書き加えました。

なお、整数objFunctionは、どの目的関数を使うかを選択するために使う値です。
0 のときは目的関数を設定していない状態、1のときは今回定義した出口での均一流を定義する状態、さらに2、3……と自分で作成したコスト関数に合わせて決めることができます。

Info << "Initializing objective function calculation:" << endl;
Info << "The objective function chosen is:" << objFunction << endl;
Info << "Name of the patches for which the cost function will be calculated:" 
    << objPatchNames << endl;
Info << "Number of patches:" << nObjPatch << endl;

label objPatchList[nObjPatch];
int iLoop;
for (iLoop=0; iLoop<nObjPatch; iLoop++)
{
    objPatchList[iLoop] =
        mesh.boundaryMesh().findPatchID(objPatchNames[iLoop]);
}

分かりやすさのために、どの目的関数を使ってどの流出面を何個対象にしているのかを書き出す設定をしました。
また、流出面の名前でアクセスできるようにobjPatchListというリストも作成しています。

コスト関数の計算は、costFunction.H という新しいファイルを作成して、その中で行います。

while (simple.loop(runTime))
{
    Info<< "Time = " << runTime.timeName() << nl << endl;
    #include "costFunction.H"

costFunction.Hを同じフォルダ内に作成し、myAdjointShapeOptimizationFoam.Cのメインループ内でcostFunction.Hを呼び出す設定を書いておきます。

scalar jFlowUni(0);

まず、目的関数の値となるスカラーを定義します。

for (iLoop=0; iLoop<nObjPatch; iLoop++)
{
    if (objFunction==1)
    {
        jFlowUni =
            jFlowUni + sum(0.5*magSqr(U.boundaryField()[objPatchList[iLoop]]
            -Ud.boundaryField()[objPatchList[iLoop]])
            );
    }
}

次に、既に定義されているnObjPatchを使用して計算を終了させる基準を定式化します。
どのコスト関数を使っているのかを決めるif節を使用して、流出面の数だけforループを実行します。
積分は、対象にしている流出面に存在する値の総和なのでsum関数で計算できます。

if (objFunction==1)
{
    Info << "Objective function (Flow uni at outlet) J: "<< jFlowUni << endl;
}

使用するコスト関数の計算結果を出力します。

まとめ

今回はコスト関数に応じたソルバーの書き換えについて見てみました。
基本的には境界条件の変更が全てのようです。
これを実行してみると確かにコスト関数の値が下がっていくのが分かりました。
ソルバーの改造についてはここまでにして、次回は実行する際のcaseフォルダの設定について書こうと思います。
このコスト関数の値をソルバーの実行時にログファイルに出力する設定をcaseフォルダ内のファイルに書いておくなどの設定の話です。