[OpenFOAM 기본] OpenFOAM 문법 기본

일본어 원문링크

작성일 : 2015년 1월 8일
번역일 : 2016년 2월 3일

크롬 브라우저로 보시는 것을 권장해 드립니다.

------ OpenFOAM 배경이론 간단정리 시리즈 ------

OpenFOAM 배경이론 간단정리 목차로 이동

------------------------------------------------------------

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

단위의 프리셋

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

대수방정식

// 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

------ OpenFOAM 소스코드 파해치기 시리즈 ------

OpenFOAM 소스코드 파해치기 목차로 이동

------------------------------------------------------------

저작자표시 비영리 변경금지

'OpenFOAM > OpenFOAM 기본' 카테고리의 다른 글

[OpenFOAM 기본] OpenFOAM 설정파일 형식 (2)	2016.02.01
[OpenFOAM 기본] OpenFOAM 개요 및 기본 (4)	2016.02.01

memorandus

[OpenFOAM 기본] OpenFOAM 문법 기본

'OpenFOAM > OpenFOAM 기본' 카테고리의 다른 글

티스토리툴바

[OpenFOAM 기본] OpenFOAM 문법 기본

'OpenFOAM > OpenFOAM 기본' 카테고리의 다른 글

'OpenFOAM/OpenFOAM 기본' Related Articles

티스토리툴바