OpenFOAM プログラミングメモ

2015年1月18日

バージョン

OpenFOAM 2.2 〜 3.0 ぐらい

基本型

  • label 整数
  • scalar スカラー (実数)
  • vector ベクトル
  • word 文字列
  • bool ブーリアン

単位付き型

  • dimensionedScalar スカラー
  • dimensionedVector ベクトル

dimensionSet で単位を定義できる。

// m/s
dimensionSet dim(0, 1, -1, 0, 0, 0, 0); // [kg m s K mol A Cd]
dimensionedVector vel("vel", dimensionSet(0, 1, -1, 0, 0, 0, 0), vector(0., 0., 0.));
dimensionedVector acc("acc", dimensionSet(0, 1, -2, 0, 0, 0, 0), vector(0., 0., 0.));

vector v = vel.value(); // 値の取り出し
dimensionSet dim(vel.dimensions()); // 単位の取り出し

単位のプリセット

  • dimless (無次元)
  • dimMass
  • dimLength
  • dimTemperature
  • dimArea
  • dimVolume
  • dimVelocity
  • dimAcceleration
  • dimDensity
  • dimForce
  • dimEnergy
  • dimPressure
dimensionedVector vel("vel", dimLength/dimTime, vector(0., 0., 0.));
dimensionedVector acc("acc", dimLength/sqr(dimTime), vector(0., 0., 0.));

dimensionedVector vel("vel", dimVelocity, vector(0., 0., 0.));
dimensionedVector acc("acc", dimAcceleration, vector(0., 0., 0.));

リスト型

List<label> a; // label のリスト

a.append(1);
a.append(2);
a.append(3);

forAll(a, i)
{
    Info<< a[i] << endl;
}

基本型リスト

  • labelList
  • scalarList
  • vectorList
  • wordList
  • boolList

単位付き型リスト

List<dimensionedScalar> a;

フィールド型

  • volScalarField セルのスカラー場
  • volVectorField セルのベクトル場
  • surfaceScalarField フェイスのスカラー場
  • surfaceVectorField フェイスのベクトル場

スカラー

絶対値

scalar a = mag(x);

累乗

scalar a;
a = sqr(x);    // 2 乗
a = pow2(x);   // 2 乗
a = pow3(x);   // 3 乗
a = pow4(x);   // 4 乗
a = pow5(x);   // 5 乗
a = pow6(x);   // 6 乗
a = pow025(x); // 0.25 乗
a = pow(x, r); // r 乗

0 割り防止

scalar a = y/stabilise(x, SMALL);

ベクトル

定義

vector v(0., 0., 0.); // Vector<scalar> と同じ

ベクトルの成分

scalar vx = v.x();
scalar vy = v.y();
scalar vz = v.z();

大きさ

scalar a = mag(v1);
scalar b = magSqr(v1); // sqr(mag(v1))

スカラー倍

vector v2 = a*v1;

和・差

vector v3 = v1 + v2;
vector v4 = v1 - v2;

内積

scalar a = v1 & v2;

外積

vector v3 = v1 ^ v2;

テンソル (2 階)

定義

Tensor<scalar> t(
	1., 2., 3.,
	4., 5., 6.,
	7., 8., 9.
);

0 テンソル

Tensor<scalar> t = Tensor<scalar>::zero

単位テンソル

Tensor<scalar> t = Tensor<scalar>::I

テンソルの成分

scalar txx = t.xx();
scalar txy = t.xy();
scalar txz = t.xz();
scalar tyx = t.yx();
...
scalar tzz = t.zz();

大きさ

scalar a = mag(t);    // ::sqrt(magSqr(t))
scalar b = magSqr(t); // 各成分の magSqrt() の和

転置

Tensor<scalar> t2 = t.T();

対称テンソル

SymmTensor<scalar> t2 = symm(t);  // 1./2.*(t + t.T())
SymmTensor<scalar> t3 = twoSymm(t); // 2.*sym(t)

反対称テンソル

Tensor<scalar> t2 = skew(t); // 1./2.*(t - t.T())

トレース

scalar a = tr(t) // t.xx() + t.yy() + t.zz()

偏差テンソル

Tensor<scalar> t2 = dev(t);  // t - 1./3.*tr(t)
Tensor<scalar> t3 = dev2(t); // t - 2./3.*tr(t)

行列式 (determinant)

scalar a = det(t);

逆行列

Tensor<scalar> t = inv(t);

不変量

scalar a = invariantI(t);   // 第 1 不変量
scalar b = invariantII(t);  // 第 2 不変量
scalar c = invariantIII(t); // 第 3 不変量

固有値・固有ベクトル

vector v = eigenValues(t);
Tensor<scalar> t2 = eigenVectors(t);

フィールド

平均・最大・最小

scalar pavg = average(p).value();
scalar pmax = max(p).value();
scalar pmin = min(p).value();

値の制限

// 0 - 1 の範囲に制限
Yi.max(0.);
Yi.min(1.);

テンソルフィールドの転置

T(fvc::grad(U))

入出力

出力

Info<< "hello" << endl;

Info<< "x = " << x << endl;

Info<< "v = ("
    << v[0] << ", " << v[1] << ", " << v[2]
    << ")" << endl;

エラー停止

FatalErrorIn
(
    "main()" // 場所
)    << "Error" // エラーメッセージ
     << exit(FatalError); // 終了

入力

既存の辞書を使う

// constant/transportProperties
label n(readLabel(transportProperties.lookup("n")); 
scalar a(readScalar(transportProperties.lookup("a")));
vector v(transportProperties.lookup("v"));
word w(transportProperties.lookup("w"));
bool b(readBool(transportProperties.lookup("b")));
dimensionedScalar x(transportProperties.lookup("x"));

// system/fvSolutions
dictionary piso(mesh.solutionDict().subDict("PISO"));
bool b = piso.lookupOrDefault<bool>("b", true);
label n = piso.lookupOrDefault<label>("n", 0);

自分で辞書を作る

// constant/myDict
IOdictionary myDict
(   
    IOobject(
        "myDict",
        mesh.time().constant(),
        mesh,
        IOobject::MUST_READ,
        IOobject::NO_WRITE
    )   
);

時間

時間オブジェクトの準備

#include "createTime.H"

時間

const word &t = runTime.timeName();
const word &t = mesh.time().timeName();

scalar t = runTime.timeOutputValue();
scalar t = mesh.time().timeOutputValue();

時間刻み幅

scalar deltaT = runTime.deltaTValue();
scalar deltaT = mesh.time().deltaTValue();

時間ループ

while(runTime.loop())
{
    ...
}

メッシュ

メッシュに関するクラスとメンバ関数の早見表

メッシュ
polyMeshfvMesh (polyMesh を継承)
セルcells-
フェイスfaces-
points-
フェイスオーナーセルfaceOwnerowner (内部フェイス)
フェイス隣接セルfaceNeighbourneighbour
セル中心cellCentresC
フェイス中心faceCentresCf (内部フェイス)
セル体積cellVolumesV
フェイス面積faceAreasSf, magSf (内部フェイス)
境界メッシュboundaryMeshboundary
セルゾーンcellZones-
フェイスゾーンfaceZones-

境界メッシュ
polyBoundaryMeshfvBoundaryMesh
polyPatchoperator[]patch
fvPatch-operator[]

パッチ
polyPatchfvPatch
名前namename
タイプtypetype
物理タイプphysicalType-
フェイスoperator[]-
セル IDfaceCellsfaceCells
フェイス中心faceCentresCf
セル中心faceCellCentresCn
フェイス面積faceAreasSf, magSf
polyPatchpatch

セルループ

forAll(mesh.cells(), cid)
{
    const cell &c = mesh.cells()[cid];
    ...
}

セルフェイスループ

forAll(mesh.cells(), cid)
{
    const cell &c = mesh.cells()[cid];

    forAll(c, cfid)
    {
        label fid = c[cfid];
        ...
    }
}

セル点ループ

forAll(mesh.cells(), cid)
{
    const cell &c = mesh.cells()[cid];

    forAll(mesh.cellPoints()[cid], cpid)
    {
        label pid = mesh.cellPoints()[cid][cpid];
        ...
    }
}

// フェイスごと
forAll(mesh.faces(), fid)
{
    const face &f = mesh.faces()[fid];

    forAll(f, fpid)
    {
        label pid = f[fpid];
        ...
    }
}

セル隣接セルループ

forAll(mesh.cells(), cid)
{
    const cell &c = mesh.cells()[cid];

    forAll(mesh.cellCells()[cid], ccid)
    {
        label ncid = mesh.cellCells()[cid][ccid];
        ...
    }
}

// フェイスごと
forAll(mesh.cells(), cid)
{
    const cell &c = mesh.cells()[cid];

    forAll(c, cfid)
    {
        label fid = c[cfid];
        
        if (!mesh.isInternalFace(fid)) continue;
        
        label ncid = mesh.faceNeighbour()[fid];
        ...
    }
}

フェイスループ

// 全フェイス
forAll(mesh.faces(), fid)
{
    const face &f = mesh.faces()[fid];
    ...
}

// 内部フェイス
forAll(mesh.faces(), fid)
{
    if (!mesh.isInternalFace(fid)) continue;

    const face &f = mesh.faces()[fid];
    ...
}

フェイス点ループ

forAll(mesh.faces(), fid)
{
    const face &f = mesh.faces()[fid];

    forAll(f, fpid)
    {
        label pid = f[fpid];
        ...
    }
}

点ループ

forAll(mesh.points(), pid)
{
    const point &p = mesh.points()[pid];
    ...
}

セル中心

const vector &C = mesh.C()[cid];

セル体積

scalar V = mesh.V()[cid];

フェイス中心座標

const vector &x = mesh.faceCentres()[fid]; // 全フェイス
const vector &x = mesh.Cf()[fid]; // 内部フェイス

フェイス面積

const vector &S = mesh.faceAreas()[fid]; // 全フェイス
const vector &S = mesh.Sf()[fid]; // 内部フェイス
scalar magS = mesh.magSf()[fid]; // 内部フェイス

フェイスオーナーセル

label cid = mesh.faceOwner()[fid]; // 全フェイス
label cid = mesh.owner()[fid]; // 内部フェイス

フェイス隣接セル

// 両者同じ
label cid = mesh.faceNeighbour()[fid];
label cid = mesh.neighbour()[fid];

点座標

const point &p = mesh.points()[pid];
const vector &x = p;

境界ループ

forAll(mesh.boundaryMesh(), bid)
{
    const polyPatch &patch = mesh.boundaryMesh()[bid];
    ...
}

forAll(mesh.boundary(), bid)
{
    const polyPatch &patch = mesh.boundary()[bid].patch();
    ...
}

境界フェイスループ

forAll(mesh.faces(), fid)
{
    if (mesh.isInternalFace(fid)) continue;

    const face &f = mesh.faces()[fid];
    ...
}

// パッチごと
forAll(mesh.boundaryMesh(), bid)
{
    const polyPatch &patch = mesh.boundaryMesh()[bid];

    forAll(patch, pfid)
    {
        const face &f = patch[pfid];
        ...
    }
}

// empty 境界を除く
forAll(mesh.boundary(), bid)
{
    const fvPatch &patch = mesh.boundary()[bid];
    
    forAll(patch, pfid)
    {
        ...
    }
}

境界セルループ

// フェイスごと
forAll(mesh.faces(), fid)
{
    if (mesh.isInternalFace(fid)) continue;

    label cid = mesh.owner()[fid];
    const cell &c = mesh.cells()[cid];
    ...
}

// パッチごと
forAll(mesh.boundaryMesh(), bid)
{
    const polyPatch &patch = mesh.boundaryMesh()[bid];

    forAll(patch.faceCells(), pcid)
    {
        label cid = patch.faceCells()[pcid];
        const cell &c = mesh.cells()[cid];
        ...
    }
}

// empty 境界を除く
forAll(mesh.boundary(), bid)
{
    const fvPatch &patch = mesh.boundary()[bid];
    
    forAll(patch.faceCells(), pcid)
    {
        label cid = patch.faceCells()[pcid];
        const cell &c = mesh.cells()[cid];
        ...
    }
}

パッチの取得

// polyPatch
// 境界 ID から
const polyPatch &patch = mesh.boundaryMesh()[bid];
const polyPatch &patch = mesh.boundary()[bid].patch();

// パッチ名から
const polyPatch &patch = mesh.boundaryMesh()[name];
const polyPatch &patch = mesh.boundary()[name].patch();

// 境界フェイス ID から
label bid = mesh.boundaryMesh().whichPatch(fid);
const polyPatch &patch = mesh.boundaryMesh()[bid];

// fvPatch
const fvPatch &patch = mesh.boundary()[bid];
const fvPatch &patch = mesh.boundary()[name];

パッチ ID の取得

label pid = mesh.boundaryMesh().findPatchID(name);
label pid = mesh.boundary().findPatchID(name);

境界情報

const polyPatch &patch = mesh.boundaryMesh()[bid];
const word &name = patch.name(); // 名前
const word &type = patch.type(); // タイプ
const word &physicalType = patch.physicalType(); // 物理タイプ

const fvPatch &patch = mesh.boundary()[bid];
const word &name = patch.name(); // 名前
const word &type = patch.type(); // タイプ

フェイス中心 (境界フェイス)

const vector &x = mesh.boundaryMesh()[bid].faceCentres()[fid];
const vector &x = mesh.boundary()[bid].Cf()[fid];

セル中心 (境界セル)

const vector &x = mesh.boundaryMesh()[bid].cellCentres()[cid];
const vector &x = mesh.boundary()[bid].Cn()[cid];

フェイス面積 (境界フェイス)

const vector &S = mesh.boundaryMesh()[bid].faceAreas()[fid];
const vector &S = mesh.boundary()[bid].Sf()[fid];
const scalar magS = mesh.boundary()[bid].magSf()[fid];

セルゾーン

label id = mesh.cellZones().findZoneID(name); // 見つからなければ -1
const labelList& cells = mesh.cellZones()[id];

フェイスゾーン

label id = mesh.faceZones().findZoneID(name); // 見つからなければ -1
const labelList& faces = mesh.faceZones()[id];

フィールド

セル値

scalar pc = p[cid]; // volScalarField
vector Uc = U[cid]; // volVectorField

フェイス値 (内部フェイス)

surfaceScalarField pf = fvc::interpolate(p);
scalar pfv = pf[fid];

surfaceVectorField Uf = fvc::interpolate(U);
vector Ufv = Uf[fid];

境界値

label bfid = fid - mesh.boundaryMesh()[bid].start();
scalar pfv = p.boundaryField()[bid][bfid];
vector Ufv = U.boundaryField()[bid][bfid];

境界タイプ

const word &type = U.boundaryField()[bid].type();

パッチの取得

const fvPatch &patch = U.boundaryField()[bid].patch();
const polyPatch &patch = U.boundaryField()[bid].patch().patch();

パッチから境界フィールドの取得

const fvPatch &patch = mesh.boundary()[bid];
const fvPatchField<scalar> &pb = patch.lookupPatchField<volScalarField, scalar>("p");
const fvPatchField<vector> &Ub = patch.lookupPatchField<volVectorField, vector>("U");

前回の値

const volVectorField &U0 = U.oldTime();  // 前回の時間ステップの値
const volVectorField &U0 = U.prevIter(); // 前回のイテレーションの値

それぞれ storeOldTimes(), storePrevIter() でデータを保存しておく必要がある。storeOldTimes() についてはなにかのタイミング (internalField() 呼び出しなど) で勝手に呼ばれていることがある。

メッシュ

const fvMesh &mesh = U.mesh();

オブジェクトの取得

const volVectorField &U = mesh.lookupObject<volVectorField>("U");

フィールドの作成

すでにあるフィールドをコピーして名前をつける

volScalarField my_p("my_p", p);
volVectorField my_U("my_U", U);

ちゃんと作る

volScalarField my_p
(
    IOobject
    (   
        "my_p",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),  
    mesh,
    dimensionedScalar("my_p", sqr(dimLength)/sqr(dimTime), 0.)
);

volVectorField my_U
(
    IOobject
    (   
        "my_U",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),  
    mesh,
    dimensionedVector("my_U", dimVelocity, vector(0., 0., 0.))
);

IOobject::NO_READ, IOobject::NO_WRITE はフィールドデータの入出力設定。

入力

  • NO_READ : 読まない
  • MUST_READ : 必ず読む
  • READ_IF_PRESENT : あったら読む

出力

  • NO_WRITE : 書き出さない
  • AUTO_WRITE : 書き出す

NO_READ, NO_WRITE を同時に指定する場合、入出力設定の引数を省略してもよい。

フィールドの出力設定

U.writeOpt() = IOobject::AUTO_WRITE; // 書き出す
U.writeOpt() = IOobject::NO_WRITE; // 書き出さない

フィールドの更新

    scalarField &fp = f.internalField();

    forAll(fp, celli)
    {
        fp[celli] = ... ;
    }

    forAll(f.boundaryField(), patchi)
    {
        fvPatchScalarField &fbp = f.boundaryField()[patchi];
        const fvPatch &patch = fbp.patch();

        forAll(patch, facei)
        {
            const label owner = patch.faceCells()[facei]; // フェイスセルの ID (必要であれば)

            fbp[facei] = ... ;
        }
    }

注意: フィールド内部の値は直接編集できるが、internalField() を呼んでおかないと値が変更されたと認識されない。

代数方程式

代数方程式

// 3 x 1
volScalarField x("x", p);
fvScalarMatrix m(x, dimless);

//  4   1   0
//  2   6   1
//  0   2   5
m.diag()[0] = 4;
m.diag()[1] = 6;
m.diag()[2] = 5;
m.upper()[0] = 1;
m.upper()[1] = 1;
m.lower()[0] = 2;
m.lower()[1] = 2;

// RHS
m.source()[0] = 2;
m.source()[1] = -8;
m.source()[2] = 6;

m.solve();

// x = [1 -2 2]'
Info<< "x =" << endl;
forAll(x, i)
{
    Info<< x[i] << endl;
}

3 x 1 のメッシュだと、セルが 0-1-2 と並び、セル 0 とセル 2 は関係がないので、行列の (0, 2) および (2, 0) の成分は 0 となり、それ以外に値が入れることになる。

残差

r = m.solve().initialResidual(); // 解く前の残差
r = m.solve().finalResidual(); // 解いた後の残差

行列の係数のアドレス

i = m.lduAddr().triIndex(cid, nid);
m.upper()[i] = a; // cid < nid
m.lower()[i] = a; // cid > nid

並列処理

フィールドの平均値や最大・最小値などを計算する場合、ふつうにやると並列計算時にパーティションごとに計算されて、全体の値を計算できない。全体の値を計算するためには reduce() を使う。

たとえば、全セル数を計算する場合、次のようにする。

label cells = mesh.cells().size();
reduce(cells, sumOp<label>());

reduce() の第 1 引数は処理したい変数、第 2 引数は演算の指定で、sumOp() 以外に以下のものがある。

  • sumOp
  • minOp
  • maxOp
  • andOp
  • orOp