OpenFOAM の adjoint Shape Optimization Foam を考える #03実装

ようやくOpenFOAMの話です。
今回は、adjoint法を用いてトポロジー最適化を行う際の実装について考えてみます。

コードの行数と内容はOpenFOAMv1812を参考にしています。

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

シリーズ

OpenFoamでの実装

ここまでの回で見てきたように、adjoint法を用いて感度を計算するためには、主方程式（速度と圧力についての式）とadjoint方程式（adjoint velocityとadjoint pressureの式）を解く必要があります。

この2つの方程式には似たような項が含まれているので、同じようなアルゴリズムを使用して解くことができます。
ソルバーは、部分的に収束した値を使用して感度を計算します。
ここでは圧力と速度を結合するSIMPLE法を用いて、adjoint方程式を主方程式と同じ方法で解きました。

ソルバー

adjointShapeOptimizationFoam ソルバーは $FOAM_SOLVERS/incompressible/adjointShapeOptimizationFoamの中にあります。
フォルダの中にはこのように色々なファイルが含まれています。

これを順番に見ていこうと思います。

adjointShapeOptimizationFoam.Cで呼ばれるヘッダファイル

まずは、adjointShapeOptimizationFoamの本体ともいえる「adjointShapeOptimizationFoam.C」です。

int main(int argc, char *argv[])
{
    #include "setRootCaseLists.H"
    #include "createTime.H"
    #include "createMesh.H"
    #include "createFields.H"
    #include "initContinuityErrs.H"
    #include "initAdjointContinuityErrs.H"

基本的に最初に読み込むヘッダファイルは、 src/OpenFoam/lnInclude/ と、 src/finiteVolume/lnInclude/ に含まれています。
createFields.H と initAdjointContinuityErrs.H だけは個別のものが、 adjointShapeOptimizationFoam/ フォルダにあります。

setRootCaseLists.H

このヘッダファイルでは、argList の設定をします。
argListはmain関数の引数を管理しています。
実際にソルバーを適用するcaseフォルダやrootフォルダの位置などを引数として持っているかのチェックをします。
caseに関する情報としては、プロセッサ数やプロセッサの含まれるフォルダについてのものなども含まれています。

createTime.H

TimeクラスとrunTimeクラスを作ります。
runTime オブジェクトは、時間についてのデータと、caseフォルダの中の system/controlDict で使うカウンタや値を持ちます。

createMesh.H

fvMesh クラスと mesh クラスオブジェクトを構成します。
mesh オブジェクトでは、ノードの位置や計算対象領域、セル体積、使用される補間スキームといった有限体積メッシュについての情報を持ちます。

createFields.H

adjointShapeOptimizationFoam で用いる変数の初期化をします。
速度/圧力とadjoint velocity / adjoint pressure、粘性や乱流モデルの設定などが含まれています。

adjoint velocityとadjoint pressureの設定には、通常の速度や圧力を設定するのと同じように IOobject クラスを用います。

Info<< "Reading field pa\n" << endl;
volScalarField pa
(
    IOobject
    (
        "pa",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

これによって、case フォルダの中の各タイムステップの名前がついたフォルダ内に pa ファイルを生成します。

同様に、 volVectorField Ua として、adjoint velocityのファイルも作ります。

この時点で、境界条件や境界部分のセルの値などの初期値を読み込みます。

#include "createPhia.H"


label paRefCell = 0;
scalar paRefValue = 0.0;
setRefCell
(
    pa,
    simple.dict(),
    paRefCell,
    paRefValue
);
mesh.setFluxRequired(pa.name());


singlePhaseTransportModel laminarTransport(U, phi);

autoPtr<incompressible::turbulenceModel> turbulence
(
    incompressible::turbulenceModel::New(U, phi, laminarTransport)
);


dimensionedScalar zeroSensitivity("0", dimVelocity*dimVelocity, 0.0);
dimensionedScalar zeroAlpha("0", dimless/dimTime, 0.0);

dimensionedScalar lambda(laminarTransport.lookup("lambda"));
dimensionedScalar alphaMax(laminarTransport.lookup("alphaMax"));

const labelList& inletCells = mesh.boundary()["inlet"].faceCells();
//const labelList& outletCells = mesh.boundary()["outlet"].faceCells();

volScalarField alpha
(
    IOobject
    (
        "alpha",
        runTime.timeName(),
        mesh,
        IOobject::READ_IF_PRESENT,
        IOobject::AUTO_WRITE
    ),
    lambda*max(Ua & U, zeroSensitivity)
);
zeroCells(alpha, inletCells);
//zeroCells(alpha, outletCells);

この部分では、ソルバーで使用するスカラーの定義とfvSolution で設定できる乱流モデルへのスマートポインタ(autoPtr)などを参照しています。

66行目では、createPhia.H を呼んで、surfaceScalarField phia の初期化をしています。
これは、createPhia.H で、

linearInterpolate(Ua) & mesh.Sf()

のように定義されています。
隣接する節点から、ある面における Ua の値が補間され、面領域ベクトルの内積を計算します。
93行目のalphaMaxのようなスカラーは、ケースフォルダ内のDictで設定しておくことができます。
それ以外の変数は、93行目にある、流入部の面ノードへの参照であるinletCells のように、このファイルで定義しています。

また、98−109行ではporosity場の定義をしています。

READ_IF_PRESENTというオプションを設定すると、値が存在する場合(if present) に読み込んでalphaを初期化します。

initAdjointContinuityErrs.H

通常のcontinuityErrs.H と同じように、adjoin velocityやadjoint pressureの収束判定に使う残差を初期化しています。

adjointShapeOptimizationFoam.C内の実装

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

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

//alpha +=
//    mesh.relaxationFactor("alpha")
//   *(lambda*max(Ua & U, zeroSensitivity) - alpha);
alpha +=
    mesh.fieldRelaxationFactor("alpha")
   *(min(max(alpha + lambda*(Ua & U), zeroAlpha), alphaMax) - alpha);

zeroCells(alpha, inletCells);
//zeroCells(alpha, outletCells);

// Pressure-velocity SIMPLE corrector
{
    // Momentum predictor

    tmp<fvVectorMatrix> tUEqn
    (
        fvm::div(phi, U)
      + turbulence->divDevReff(U)
      + fvm::Sp(alpha, U)
     ==
        fvOptions(U)
    );
    fvVectorMatrix& UEqn = tUEqn.ref();

    UEqn.relax();

    fvOptions.constrain(UEqn);

    solve(UEqn == -fvc::grad(p));

    fvOptions.correct(U);

106-108行目では、最急降下法を用いてporosity場が更新されています。

caseフォルダ内のfvSolutionで緩和係数を定義することができます。
定義された緩和係数はfieldRelaxationFactor()関数から取得しています。

createFields.H で以下のように定義されている、

dimensionedScalarlambda(laminarTransport.lookup("lambda"));

“lambda”は、セルの体積とステップの長さの積です。
これはcase フォルダの constant/transportProperties で与えることができます。

[1]　$\alpha_{n+1} = \alpha_n (1-\gamma) + \gamma min( max(( \alpha_n – \textbf{u}_{i} \cdot \textbf{v}_{i} V_{i} \delta), 0), \alpha_{max})$

元になっている式[1]と比較すると、alphaの新しい値が定義されている108行目の符号が間違っているようです。
これを修正すると以下のようになります。

alpha +=
    mesh.fieldRelaxationFactor("alpha")
   *(min(max(alpha - lambda*(Ua & U), zeroAlpha), alphaMax) - alpha);

110行目では、adjointShapeOptimizationFoam.C のmain関数の前に関数zeroCells を定義して、流入部のセルのporosityを0に設定しています。
zeroCellsは指定した変数を0に設定することができる関数で、スカラー、ベクトル、テンソルなど、複数の型の変数を使用することができます。

次に、117行目以降と以下の運動方程式を比較します。

[2]　$(\textbf{v}\cdot\nabla)\textbf{v}+\nabla p-\nabla(2\nu D(\textbf{v}))+\alpha\textbf{v}=0$

120行目の turbulence は createFields.H で

autoPtr<incompressible::turbulenceModel> turbulence
(
    incompressible::turbulenceModel::New(U, phi, laminarTransport)
);

のように設定されているRASModelへのポインタです。
これを見ると、incompressible::turbulenceModel型を使っているようなので、/src/TurbulenceModels/incompressible/IncompressibleTurbulenceModel/IncompressibleTurbulenceModel.Cを見てみます。
偏差応力項であるdivDevReff()は以下のように定義されています。

template<class TransportModel>
Foam::tmp<Foam::fvVectorMatrix>
Foam::IncompressibleTurbulenceModel<TransportModel>::divDevReff
(
    volVectorField& U
) const
{
    return divDevRhoReff(U);
}

divDevRhoReff(U) を使っているようです。

標準k-εモデルを使用する場合について考えると、src/TurbulenceModels/turbulenceModels/RAS/kEpsilon/kEpsilon.C に、

// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

template<class BasicTurbulenceModel>
kEpsilon<BasicTurbulenceModel>::kEpsilon
(
    const alphaField& alpha,
    const rhoField& rho,
    const volVectorField& U,
    const surfaceScalarField& alphaRhoPhi,
    const surfaceScalarField& phi,
    const transportModel& transport,
    const word& propertiesName,
    const word& type
)
:
    eddyViscosity<RASModel<BasicTurbulenceModel>>
    (
        type,
        alpha,
        rho,
        U,
        alphaRhoPhi,
        phi,
        transport,
        propertiesName
    ),

とあるので、eddyViscosity<RASModel<BasicTurbulenceModel>>を継承していることが分かります。
src/TurbulenceModels/turbulenceModels/eddyViscosity/eddyViscosity.C では、

// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

template<class BasicTurbulenceModel>
Foam::eddyViscosity<BasicTurbulenceModel>::eddyViscosity
(
    const word& type,
    const alphaField& alpha,
    const rhoField& rho,
    const volVectorField& U,
    const surfaceScalarField& alphaRhoPhi,
    const surfaceScalarField& phi,
    const transportModel& transport,
    const word& propertiesName
)
:
    linearViscousStress<BasicTurbulenceModel>
    (
        type,
        alpha,
        rho,
        U,
        alphaRhoPhi,
        phi,
        transport,
        propertiesName
    ),

linearViscousStress を継承しているようです。

template<class BasicTurbulenceModel>
Foam::tmp<Foam::fvVectorMatrix>
Foam::linearViscousStress<BasicTurbulenceModel>::divDevRhoReff
(
    volVectorField& U
) const
{
    return
    (
      - fvc::div((this->alpha_*this->rho_*this->nuEff())*dev2(T(fvc::grad(U))))
      - fvm::laplacian(this->alpha_*this->rho_*this->nuEff(), U)
    );
}

ということで、ようやくたどり着きましたが、これが

+ turbulence->divDevReff(U)

の指し示している内容のようです。
粘性応力の発散を計算しているようです。

src/OpenFOAM/primitives/Tensor/TensorI.H によると、dev2() は、

$$dev2(\textbf{A})=\textbf{A}-\frac{2}{3}(tr\textbf{A})\textbf{I}$$

です。

式[2]中の $2\nu D(\textbf{v})$ にあたる部分です。
一般的に非圧縮性流れの場合、連続性によってこの項は$0$になります。

porosity によるシンク項は他のソース項と同じように扱われ、109行目のfvm::Sp(alpha,U) によって陰的に実装されます。

次に、運動量方程式を解いている部分(Momentum Predictor)で予測した補正速度の値を用いて圧力を求めます。
そのようにして得た圧力pを使用して速度を計算するというSIMPLE法を用いています。

// Adjoint Pressure-velocity SIMPLE corrector
{
    // Adjoint Momentum predictor
 
    volVectorField adjointTransposeConvection((fvc::grad(Ua) & U));
    //volVectorField adjointTransposeConvection
    //(
    //    fvc::reconstruct
    //    (
    //        mesh.magSf()*fvc::dotInterpolate(fvc::snGrad(Ua), U)
    //    )
    //);
 
    zeroCells(adjointTransposeConvection, inletCells);
 
    tmp<fvVectorMatrix> tUaEqn
    (
        fvm::div(-phi, Ua)
        - adjointTransposeConvection
        + turbulence->divDevReff(Ua)
        + fvm::Sp(alpha, Ua)
        ==
        fvOptions(Ua)
    );
    fvVectorMatrix& UaEqn = tUaEqn.ref();
 
    UaEqn.relax();
 
    fvOptions.constrain(UaEqn);
 
    solve(UaEqn == -fvc::grad(pa));

adjoint運動量方程式の解き方は通常の運動方程式の解法にとてもよく似ています。
違いはいくつかの符号と追加の項 -adjointTransposeConvectionです。

adjoint方程式における非圧縮性の定常状態のNavier-Stokes方程式は以下のようになりますが、

[4]　$-2D(\textbf{u})\textbf{v}=-\nabla\textbf{u}\cdot\textbf{v}-(v\cdot\nabla)\textbf{u}$

その右辺第1項が-adjointTransposeConvection で、第2項は180行目で計算されている部分です。

176行目を見ると、adjointTransposeConvection項の流入部のセルにおける値が0に設定されているのがわかります。

まとめ

長くなりましたが、OpenFOAM の adjoint Shape Optimization Foam を考える #01手法に対応するコードの基本的な部分はここまでです。
adjoint方程式の解析には、通常の速度と圧力を解くのととても似た方法を使います。

次回は、ソルバー内でどのように境界条件の設定をしているのかを見ます。

- OpenFOAM の adjoint Shape Optimization Foam を考える #04境界条件の実装 – 白黒羊
2020/07/09 - 返信

[…] 前回は「OpenFOAM の adjoint Shape Optimization Foam を考える #03実装」でした。 […]
- OpenFOAM の adjoint Shape Optimization Foam を考える #05ソルバーの改造 – 白黒羊
2021/02/10 - 返信

[…] OpenFOAM の adjoint Shape Optimization Foam を考える #03実装 […]