OpenFOAM プログラミングメモ2020年6月27日 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
バージョンOpenFOAM 2.2 〜 型基本型
単位付き型
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()); // 単位の取り出し
単位のプリセット
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;
}
基本型リスト
単位付き型リスト List<dimensionedScalar> a; フィールド型
スカラー絶対値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); 入出力出力
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"));
List<word> names(transportProperties.lookup("names"));
// 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())
{
...
}
メッシュメッシュに関するクラスとメンバ関数の早見表
セルループ
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");
平均・最大・最小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)) 前回の値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, NO_WRITE を同時に指定する場合、入出力設定の引数を省略してもよい。 初期化
C = dimensionedScalar("C", C.dimensions(), 0.);
初期境界条件の設定
volScalarField C
(
IOobject
(
"C",
runTime.timeName(),
mesh,
IOobject::NO_READ
),
mesh,
dimensionedScalar("C", dimless, 0.),
"zeroGradient"
);
最後の引数で紹介条件を指定。デフォルトは "calculated"。 フィールドの出力設定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() を呼んでおかないと値が変更されたと認識されない。 セルの探索点を含むセルを探す。 const label celli = mesh.findCell(p); 補間点の位置で値を補間する。
#include "interpolation.H"
autoPtr<interpolation<scalar>> interpolator
(
interpolation<scalar>::New("cellPoint", T)
);
const label celli = mesh.findCell(p);
const scalar Tp = interpolator().interpolate(p, celli, -1);
スキームは "cell", "cellPoint", "cellPointFace" などを指定できる。 代数方程式代数方程式
// 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() 以外に以下のものがある。
リストの統合各プロセッサでデータを append() したリストを統合するにはつぎのようにする。 List<label> ids; ... Pstream::combineGather(ids, ListOps::appendEqOp 出力プロセッサごとに出力させたい場合は、Info の代わりに cout を用いる。 ファイルの出力値をファイルに出力する。
#include "writer.H"
autoPtr<writer<scalar>> scalarFormatter
(
writer<scalar>::New("raw")
);
coordSet cs
(
"",
"xyz",
points,
mag(points)
);
List<const scalarField*> valueList
(
1,
&values
);
OFstream os(filename);
scalarFormatter().write
(
cs,
wordList(0),
valueList,
os
);
点列 points に対応した値の列 values を raw フォーマットで出力する。ライブラリのリンクが必要かも。 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| PENGUINITIS | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||