티스토리 뷰

OpenFOAM/OpenFOAM 기본

[OpenFOAM 기본] OpenFOAM 문법 기본

하루생각 하루생각 2016.02.03 14:50

일본어 원문링크


작성일 : 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>

단위의 프리셋

  • 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
직접 딕셔너리를 작성
기존의 딕셔너리를 사용
// 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 소스코드 파해치기 목차로 이동

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


저작자 표시 비영리 변경 금지
신고
댓글
  • 프로필사진 비밀댓글입니다 2016.09.01 11:30
  • 프로필사진 BlogIcon 하루생각 하루생각 안녕하세요. 애틀란타에 계시는군요^^
    정확히 어떤 매시간의 결과들을 원하시는지요?

    postprocessing의 개념이면, controlDict에서 저장하기로 설정한 시각의 값들이 폴더로 저장되어있어야 읽어들일 수 있습니다.

    solver가 계산하는 도중의 매번의 값을 이용한다면, solver자체를 변경하거나, controlDict에 Cp나 힘을 계산하는 라이브러리코드를 수정하는 방법이 있습니다. (키워드 : Run-time Postprocessing)
    https://www.hpc.ntnu.no/display/hpc/OpenFOAM+-+Run-time+Postprocessing
    가 참고가 될지 모르겠습니다. (자세히 살펴보지는 않았습니다)

    보통 솔버에서 변수를 건드리는 경우는, 새로운 변수를 정의하는 식 뿐만아니라, 기존의 변수들과의 관계식들을 추가하여야하는 경우가 많은데, 이럴때 코드가 매우 불안정해집니다... (실패경험, 일전에 Multiphase문제에서, 비압축성 icoFoam에 온도항만 추가해서(온도항 추가는 크게 코드 안정성에 무리 없음), 이것을 가지고 상변화를 보정시키려 했었습니다만, 뭔가 잘 되지 않았습니다.)

    먼저 변수간의 관계식이 직접적으로 솔버에 식을 추가하지않는 방향으로 가정을 더해 변형시켜 Post 로 풀어내는 것을 추천드립니다.
    2016.09.01 14:39 신고
공지사항
글 보관함
링크
«   2017/12   »
          1 2
3 4 5 6 7 8 9
10 11 12 13 14 15 16
17 18 19 20 21 22 23
24 25 26 27 28 29 30
31            
Total
80,139
Today
17
Yesterday
144