2016

2016. 02. 03 / 하루생각 Blog


[OpenFOAM 기본] OpenFOAM 문법 기본

일본어 원문링크


들어가기

OpenFOAM 을 이용한 코드 작성 및 분석을 위한 기본 문법



기준 OpenFOAM 버전

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;
}


키워드 부분은, 다음과 같은 형식이다.


<키워드> <항목 1> <항목 2> ... <항목 N>;



단위의 프리셋

  • labelList
  • scalarList
  • vectorList
  • wordList
  • boolList


단위를 포함하는 리스트


List<label> List 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;


텐서

정의

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("b", true);
label n = piso.lookupOrDefault("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())
{
    ...
}


격자

격자와 관련된 클래스와 멤버함수 정리표

격자

 

polyMesh

fvMesh

(polyMesh 계승)

cells

-

faces

-

points

-

면 소유 셀

faceOwner

owner (내부 면)

면 접근 셀

faceNeighbour

neighbour

셀 중심

cellCentres

C

면 중심

faceCenteres

Cf (내부 면)

셀 체적

cellVolumes

V

면 면적

faceAreas

Sf, magSf (내부 면)

경계 격자

boundaryMesh

boundary

셀 존

cellZones

-

면 존

faceZones

-


경계 격자

 

polyBoundaryMesh

fvBoundaryMesh

polyPatch

operator[]

patch

fvPatch

-

operator[]


패치

 

polyMesh

fvPatch

이름

 name

 name

타입

 type

 type

물리타입

 physicalType

 -

 operator[]

 -

셀 ID

 faceCells

 faceCells

면 중심

 faceCentres

 Cf

셀 중심

 faceCellCentres

 Cn

면 면적

 faceAreas

Sf, mag Sf 

polyPatch

 

 patch


셀루프

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 &pb = patch.lookupPatchField("p");
const fvPatchField &Ub = patch.lookupPatchField("U");

이전 루프의 값

const volVectorField &U0 = U.oldTime(); // 이전 타임스텝에서의 값
const volVectorField &U0 = U.prevIter(); // 이전 반복문에서의 값

각각 storeOldTimes(), storePrevIter()로 데이터를 저장해 놓아야 한다. storeOldTimes()은 이미 코드 중에 사용되어있는 경우가 있다(internalFiled() 명령어 등에서)


격자

const fvMesh &mesh = U.mesh();

오브젝트 취득

const volVectorField &U = mesh.lookupObject("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_READ : 출력한다

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()를 이용해 불러와 놓지 않으면 값의 변경을 인식하지 않는다.



대수방정식

f

// 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이 되고 그 이외의 값이 입력된다.


Residual

r = m.solve().initialResidual(); // 풀기 이전의 Residual
r = m.solve().finalResidual(); // 풀고 난 후의 Residual

행렬의 계수 주소

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());

reduce()의 첫번째 인수는 처리하고싶은 변수, 두번째 변수는 연산의 지정으로 sumOp()이외에 아래와 같은 것이 있다.


  • sumOp
  • minOp
  • maxOp
  • andOp
  • orOp

프로필사진

하루생각

OneDay by Eric

에릭의 티스토리 http://onedayof.tistory.com


[원문출처 : http://onedayof.tistory.com/entry/OpenFOAM-기본-OpenFOAM-문법-기본]

Creative Commons License
크리에이티브 커먼즈 라이선스
맨 위로
맨 위로