見出し画像

【OpenFOAM用のC++勉強(4)】速度場と圧力場の取得

物理・プログラミング・Python、Web制作・自己啓発

こんにちは。
OpenFOAM用のC++の勉強のメモを自分のために残しています。

OpenFOAM v2112

前回はメッシュのセルや面、境界に対しての情報を取得するというのをしました。
今回はセルが持つ速度や圧力の情報を取得します。
前回の記事で書いたコードを長くなってきたので必要な部分に絞って書きます。

mainOF.C

#include "argList.H"
#include "fvMesh.H"
#include "volFields.H"
#include "surfaceFields.H"

using namespace Foam;

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Main program:

int main(int argc, char *argv[])
{
    Info << "================== Create mesh ==================" << endl;
    argList args(argc, argv);

    #include "createTime.H"
    
    // mesh
    fvMesh mesh
    (
        IOobject
        (
            fvMesh::defaultRegion,
            runTime.timeName(),
            runTime,
            IOobject::MUST_READ
        )
    );
    
    Info << "============= volScalarField p ===============" << endl;
    // volScalarField typedef=> GeometricField => DimensionedField => Field
    // https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1DimensionedField.html

    word t_ = "0.5";
    
    volScalarField p
    (
        IOobject
        (
            "p",
            // runTime.timeName(), // word
            t_,
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

    Info << t_ << endl;
    vector P_pos(0, 0, 0);
    Info << "get p value at P_pos = " << P_pos << endl;

    // mesh is fvMesh class.
    // fvMesh => polyMesh
    /*
    Foam::label Foam::polyMesh::findCell
    (
        const point& p,
        const cellDecomposition decompMode
    )
    */
    label cellPNo = mesh.findCell(P_pos);
    const labelList& cellFaces = mesh.cells()[cellPNo]; // face no. list of cell
    vector cellCentor = mesh.C()[cellPNo]; // centor cord. of cell　fvMesh => polyMesh => primitiveMesh 
    scalar cellVolume = mesh.V()[cellPNo]; // volume of cell

    Info << "cellNo1 = " << cellPNo << endl;
    Info << "cellFaces = " << cellFaces << endl;
    Info << "cellCentor = " << cellCentor << endl;
    Info << "cellVolume = " << cellVolume << endl;

    Info << "p" << P_pos << " = " << p[cellPNo] << endl;
    Info << "p" << P_pos << " = " << p.ref()[cellPNo] << endl; //https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1GeometricField.html#a77a3ea1ce7e2adc04d27301292b095ae

    const dimensionedScalar& p1 = p.internalField()[cellPNo]; // Return a const-reference to the dimensioned internal field. 
    const scalar& p2 = p.primitiveField()[cellPNo]; // don't have dimension

    Info << "p1 = " << p1 << endl;
    Info << "p2 = " << p2 << endl;

    Info<< "End\n" << endl;

    return 0;
}

Make/files

meshOF.C

EXE = $(FOAM_USER_APPBIN)/mesh005

Make/options

EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude

EXE_LIBS = \
    -lfiniteVolume

以下のコマンドでコンパイル。

 wmake 2>&1 | tee log

※チュートリアルは以下からコピーしてくる。

cp -r /opt/OpenFOAM/OpenFOAM-v2012/tutorials/incompressible/icoFoam/cavity/cavity mycavity

以下のコマンドでmycavityのケースファイルに移動します。

cd mycavity

そしてblockMeshを実行しておきます。

blockMesh

これはメッシュ情報には以下のファイルが必要だからです。

  • constant/polyMesh/boundary

  • constant/polyMesh/faces

  • constant/polyMesh/neighbour

  • constant/polyMesh/owner

  • constant/polyMesh/points

そして流速Uや圧力pの情報を取得するには計算を実行しておく必要があります。

icoFoam

では、mycavityの中で以下のmesh005を実行してます。
実行結果は長いのでlog.mesh005に掃き出します。

mesh005 > log.meesh005

log.mesh005

================== Create mesh ==================
/*---------------------------------------------------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
Build  : _7bdb509494-20201222 OPENFOAM=2012
Arch   : "LSB;label=32;scalar=64"
Exec   : mesh005
Date   : Dec 06 2022
Time   : 21:09:37
Host   : DESKTOP-CT961AV
PID    : 2259
I/O    : uncollated
Case   : /mnt/d/openfoam/20221011_wolf_benchmark/OpenFOAM_Training/OF9/101programming/myWork/mesh005/mycavity
nProcs : 1
trapFpe: Floating point exception trapping enabled (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster (fileModificationSkew 5, maxFileModificationPolls 20)
allowSystemOperations : Allowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

============= volScalarField p ===============
0.5
get p value at P_pos = (0 0 0)
cellNo1 = 0
cellFaces = 6(0 1 780 820 840 1240)
cellCentor = (0.0025 0.0025 0.005)
cellVolume = 2.5e-07
p(0 0 0) = 4.29931e-06
p(0 0 0) = 4.29931e-06
p1 = 4.29931e-06 [0 0 0 0 0 0 0] 4.29931e-06
p2 = 4.29931e-06

End

プログラムの引数情報を渡す

argList のオブジェクトは、mesh005の実行ファイルの実行時に引数を付けて実行した情報をargc, argvを渡して作成するものです。
argList を持つヘッダーファイルは、プログラムの冒頭でインクルードしています。

#include "argList.H"

argListクラスのインスタンス化を行っています。
変数名はargsとしていますが、これはcreateTime.Hで必要なのでそのような名前にしています。

次にメッシュ情報のための記述をしています。

    // mesh
    fvMesh mesh
    (
        IOobject
        (
            fvMesh::defaultRegion,
            runTime.timeName(),
            runTime,
            IOobject::MUST_READ
        )
    );

メッシュ情報の取得は前回の記事で解説した内容ですね。
ちなみにrunTime.timeName()はwordクラスで文字列を意味しています。
はじめはrunTime.timeName()は0なので初期のメッシュを取得しているという意味です。
これはcreateTime.H内のrunTime関数の返りがFoam::Timeなのでそれを追っていくとtimeName()関数があるので、返ってくる型を確認すればわかります。runTime.timeName()は、wordとなっているので文字列が返ってきます。

圧力場の取得

runTime.timeName()がどの時刻を取得するか（どのファイルを取得するか）を意味するので、word t_ = "0.5";として0.5のファイルから圧力を取得するようにします。
runTime.timeName()をコメントアウトしてword型（文字列）で定義したt_をIOobjectの引数にしました。

    word t_ = "0.5";
    
    volScalarField p
    (
        IOobject
        (
            "p",
            // runTime.timeName(), // word
            t_,
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

圧力の取得はまず、座標を指定します。
それにはvector型を使って定義をすれば良いです。

    vector P_pos(0, 0, 0);
    Info << "get p value at P_pos = " << P_pos << endl;

出力結果は以下のようになっています。

get p value at P_pos = (0 0 0)

座標を指定して、近傍のセルNo.を取得します。
セルNo.は整数なのでOpenFOAM用の整数型であるlabel型で定義します。

label cellPNo = mesh.findCell(P_pos);

ちなみに、findCell()はfvMeshが継承しているpolyMeshのメンバ関数です。

https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1polyMesh.html

セルNo.が取得できれば圧力はp[セルNo.]とすることで取得できます。
「座標に対応したセルNo.→セルNo.に対応した圧力」という流れですね。

ついでにセル中心の座標や体積も取得しておきます。

    const labelList& cellFaces = mesh.cells()[cellPNo]; // face no. list of cell
    vector cellCentor = mesh.C()[cellPNo]; // centor cord. of cell　fvMesh => polyMesh => primitiveMesh 
    scalar cellVolume = mesh.V()[cellPNo]; // volume of cell

    Info << "cellNo1 = " << cellPNo << endl;
    Info << "cellFaces = " << cellFaces << endl;
    Info << "cellCentor = " << cellCentor << endl;
    Info << "cellVolume = " << cellVolume << endl;

圧力は次元を持った型であるvolScalarFieldです。
volScalarFieldはGemmetricFieldを定義しなおしたものなのでGemmetricFieldクラスのメンバ関数を扱うことができます。
以下のようにp[cellPNo]や p.ref()[cellPNo]とすると値を取得できます。

    Info << "p" << P_pos << " = " << p[cellPNo] << endl;
    Info << "p" << P_pos << " = " << p.ref()[cellPNo] << endl; //https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1GeometricField.html#a77a3ea1ce7e2adc04d27301292b095ae

またセルでの圧力を取得したい場合は以下のようにします。

    const dimensionedScalar& p1 = p.internalField()[cellPNo]; // Return a const-reference to the dimensioned internal field. 
    const scalar& p2 = p.primitiveField()[cellPNo]; // don't have dimension

    Info << "p1 = " << p1 << endl;
    Info << "p2 = " << p2 << endl;

internalField()とprimitiveField()はGemmetricFieldが持つメンバ関数で違いは以下のようになります。

  • internalField()：次元を持つ。constなので値を変えれない。

  • primitiveField()：値のみ参照する。値を変更できる。

出力結果は以下のようになります。

p1 = 4.29931e-06 [0 0 0 0 0 0 0] 4.29931e-06
p2 = 4.29931e-06

p1は次元を持っていることがわかります。
でも圧力のdimensionは[0 2 -2 0 0 0 0]のはずなのにおかしいな(*_*;
後で調べるとしてひとまず置いておこう・・・

流速Uの値も出力

圧力の取得方法とほぼ同じです。
流速はベクトルなのでvolVectorField型として定義しておく必要があります。
volVectorFieldもGemmetricFieldなので、GemmetricFieldが持つメンバ関数を使うことができます。
取得する座標は vector U_pos(0.05, 0.05, 0);としておきました。

#include "argList.H"
#include "fvMesh.H"
#include "volFields.H"
#include "surfaceFields.H"

using namespace Foam;

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Main program:

int main(int argc, char *argv[])
{
    Info << "================== Create mesh ==================" << endl;
    argList args(argc, argv);

    #include "createTime.H"
    
    // mesh
    fvMesh mesh
    (
        IOobject
        (
            fvMesh::defaultRegion,
            runTime.timeName(),
            runTime,
            IOobject::MUST_READ
        )
    );
    Info << "runTime.timeName() = " << runTime.timeName() << endl;
    
    Info << "============= volScalarField p ===============" << endl;
    // volScalarField typedef=> GeometricField => DimensionedField => Field
    // https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1DimensionedField.html

    word t_ = "0.5";
    
    volScalarField p
    (
        IOobject
        (
            "p",
            // runTime.timeName(), // word
            t_,
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

    Info << t_ << endl;
    vector P_pos(0.05, 0.05, 0);
    Info << "get p value at P_pos = " << P_pos << endl;

    // mesh is fvMesh class.
    // fvMesh => polyMesh
    /*
    Foam::label Foam::polyMesh::findCell
    (
        const point& p,
        const cellDecomposition decompMode
    )
    */
    label cellPNo = mesh.findCell(P_pos);
    const labelList& cellFaces = mesh.cells()[cellPNo]; // face no. list of cell
    vector cellCentor = mesh.C()[cellPNo]; // centor cord. of cell　fvMesh => polyMesh => primitiveMesh 
    scalar cellVolume = mesh.V()[cellPNo]; // volume of cell

    Info << "cellNo1 = " << cellPNo << endl;
    Info << "cellFaces = " << cellFaces << endl;
    Info << "cellCentor = " << cellCentor << endl;
    Info << "cellVolume = " << cellVolume << endl;

    Info << "p" << P_pos << " = " << p[cellPNo] << endl;
    Info << "p" << P_pos << " = " << p.ref()[cellPNo] << endl; //https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1GeometricField.html#a77a3ea1ce7e2adc04d27301292b095ae

    const dimensionedScalar& p1 = p.internalField()[cellPNo]; // Return a const-reference to the dimensioned internal field. 
    const scalar& p2 = p.primitiveField()[cellPNo]; // don't have dimension

    Info << "p1 = " << p1 << endl;
    Info << "p2 = " << p2 << endl;

    Info << "============= volVectorField U ===============" << endl;

    volVectorField U
    (
        IOobject
        (
            "U",
            // runTime.timeName(), // word
            t_,
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

    vector U_pos(0.05, 0.05, 0);
    Info << "get U value at U_pos = " << U_pos << endl;

    label cellUNo = mesh.findCell(U_pos);
    Info << "U" << U_pos << " = " << U[cellUNo] << endl;
    Info << "U" << U_pos << " = " << U.ref()[cellUNo] << endl;

    return 0;
}

出力結果は以下のようになります。

================== Create mesh ==================
/*---------------------------------------------------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
Build  : _7bdb509494-20201222 OPENFOAM=2012
Arch   : "LSB;label=32;scalar=64"
Exec   : mesh005
Date   : Dec 06 2022
Time   : 22:05:37
Host   : DESKTOP-CT961AV
PID    : 3076
I/O    : uncollated
Case   : /mnt/d/openfoam/20221011_wolf_benchmark/OpenFOAM_Training/OF9/101programming/myWork/mesh005/mycavity
nProcs : 1
trapFpe: Floating point exception trapping enabled (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster (fileModificationSkew 5, maxFileModificationPolls 20)
allowSystemOperations : Allowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

runTime.timeName() = 0
============= volScalarField p ===============
0.5
get p value at P_pos = (0.05 0.05 0)
cellNo1 = 189
cellFaces = 6(369 370 1029 1429 331 367)
cellCentor = (0.0475 0.0475 0.005)
cellVolume = 2.5e-07
p(0.05 0.05 0) = -0.00133453
p(0.05 0.05 0) = -0.00133453
p1 = -0.00133453 [0 0 0 0 0 0 0] -0.00133453
p2 = -0.00133453
============= volVectorField U ===============
get U value at U_pos = (0.05 0.05 0)
U(0.05 0.05 0) = (-0.196068 0.0263898 0)
U(0.05 0.05 0) = (-0.196068 0.0263898 0)

以上のように座標からセルNo.を取得すれば、セルNo.からpやUの値を取得することができます。

ちなみにOpenFOAMのコードの中で出てくるforAllはC++の文法というわけではなく、OpenFOAMが独自で定義したforループです。

#define forAll(list, i) \
    for (Foam::label i=0; i<(list).size(); ++i)

と定義しています。

まだまだわかっていない部分が多いですが、引き続き勉強します。

参考にしているOpenFOAMの辞書です。
Foundation版をベースにしていますがわからないことはこちらで調べています。

引き続きC++の勉強を続けます。

Twitter➡@t_kun_kamakiri
ブログ➡宇宙に入ったカマキリ（物理ブログ）

いいなと思ったら応援しよう！