본문 바로가기

OpenFOAM/소스코드 파해치기

[OpenFOAM 소스코드 파해치기]대수방정식솔버

일본어 원문링크


작성일 : 2009년 12월  17일
번역일 : 2016년   2월  14일


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


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

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

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

 ----- OpenFOAM 소스코드를 다루는 문법 기본 -----

OpenFOAM 소스코드를 다루는 문법 기본으로 이동

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



 

  들어가기 


대수방정식솔버에 대해 살펴보자.



  사용 버전 


OpenFOAM 1.6

 


   파일


이번에 사용파는 파일이다.(11-solver.tar.gz).


  • Make/files
  • Make/options
  • system/controlDict
  • system/fvSchemes
  • system/fvSolution
  • constant/polyMesh/blockMeshDict
  • 0/T
  • heat.C
  • myPCG/PCG.C
  • myPCG/PCG.H

이번에도 열전도 프로그램을 기반으로 하고 있다.









  대수방정식솔버의 작성

 

자신만의 대수방정식 솔버를 만들어 보자. 열전도 프로그램의 system/fvSolution을 보면

solvers
{
    T
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-06;
        relTol          0;
    }
}

이다. 이것은 T의 방정식을 풀 때에 PCG(전처리를 포함하는 병역구배법)을 사용한다는 것을 의미한다. 여기서 이 PCG의 소스코드를 복사해 자신만의 PCG "myPCG"를 작성해 보자.


먼저,  myPCG라는 폴더를 작성한다. 그 안에 PCG의 소스코드인 src/OpenFOAM/matrices/lduMatrix/solvers/PCG의 내용 (PCG.C, PCG.H)를 복사하다. 복사한 파일의 파일명을 각각 myPCG.C, myPCG.H로 바꾸고 각각의 내용을 에디터의 치환기능을 사용하여 "PCG"->"myPCG"로 모두 바꾼다.


컴파일하기위해 Make/files에 myPCG/myPCG.C를 추가한다.


계산에 myPCG를 사용하려면, system/fvSolution의 "PCG"를 "myPCG"로 바꾼다.




  실행


컴파일하여 실행해 본다.

$ wmake $ blockMesh $ ./heat

계산중에 "DICmyPCG"라는 문자가 부이면 OK .

DICmyPCG: Solving for T, Initial residual = 1, Final residual = 6.43341e-07, No Iterations 13




  병역구배법


병역구배법(CG법)의 알고리즘은 아래와 같다.


대수방정식을 Ax = b 로 설정한다. x에 초기값을 설정하고, 초기잔사(Residual) r = b - Ax 를 계산한다. p = r 로 설정하며 아래의 조작을 반복한다.

1.  alpha = (r * r)/(p * Ap)

2.  x = x + alpha p

3.  r0 = r

4.  r = r - alpha A p

5.  수렴판정

6.  beta =  (r*r)/(r0*r0)

7.  p = r + beta p

전처리를 포함할 경우 대수방정식 Ax = b 대신에 전처리행렬 M을 사용해 MAx = Mb 를 푼다. M의 역행렬을 M'으로 하여 아래와 같은 알고리즘을 따른다.

x 에 초기값을 설정하고 초기잔사 r = b - Ax를 계산한다. p = M'r 로 하며 아래와 같이 반복한다.

1.  alpha = (r * M' r)/(p * Ap)

2.  x = x + alpha p

3.  r0 = r

4.  r = r - alpha A p

5.  수렴판정

6.  beta = (r * M' r)/(r0 * M' r0)

7.  p = M' r + beta p








  프로그램

 

myPCG의 멤버함수 solve()를 살펴보도록 하자.

myPCG.C (일부)

... register label nCells = psi.size(); scalar* __restrict__ psiPtr = psi.begin(); scalarField pA(nCells); scalar* __restrict__ pAPtr = pA.begin(); scalarField wA(nCells); scalar* __restrict__ wAPtr = wA.begin(); scalar wArA = matrix_.great_; scalar wArAold = wArA; // --- Calculate A.psi matrix_.Amul(wA, psi, interfaceBouCoeffs_, interfaces_, cmpt); // --- Calculate initial residual field scalarField rA(source - wA); scalar* __restrict__ rAPtr = rA.begin(); ... // --- Select and construct the preconditioner autoPtr preconPtr = lduMatrix::preconditioner::New ( *this, controlDict_ ); // --- Solver iteration do { // --- Store previous wArA wArAold = wArA; // --- Precondition residual preconPtr->precondition(wA, rA, cmpt); // --- Update search directions: wArA = gSumProd(wA, rA); if (solverPerf.nIterations() == 0) { for (register label cell=0; cell<nCells; cell+)

{

pAPtr[cell]= wAPtr[cell];

}

}

else

{

scalar beta = wArA/wArAold;


for (register label cell=0; cell < nCells; cell++)

{

pAPtr[cell] = wAPtr[cell] + beta*pAPtr[cell];  

}

}  


// --- Update preconditioned residual matrix_.Amul(wA, pA, interfaceBouCoeffs_, interfaces_, cmpt); scalar wApA = gSumProd(wA, pA); ... // --- Update solution and residual: scalar alpha = wArA/wApA;

for (register label cell=0; cell < nCells; cell++)

{

psiPtr[cell] += alpha*pAPtr[cell];

rAPtr[cell] -= alpha*wAPtr[cell];

}

...


}while

(

solvePerf.nIterations()++ < maxIter_

&& !(solverPerf.checkConvergence(tolerance_, relTol_))

);

...

알고리즘과 프로그램의 변수의 대응은 다음과 같다.

  • psi, psiPtr : x
  • pA, pAPtr : p
  • wA, wAPtr : A x, M' r, A p
  • wArA : r * M' r
  • wArAold : r0 * M' r0
  • rA : r
  • wApA : p * A p

wA는 복수의 목적으로 사용되고 있다.

자세히 살펴보면, Ax -> wA

// --- Calculate A.psi

matrix_.Amul(wA, psi, interfaceBouCoeffs_, interfaces_, cmpt);


r = b - Ax

// --- Calculate initial residual field

scalarField rA(source - wA);


루프에 넣어 r * M' r 를 r0 * M' r0 로 하여 저장

// --- Store previous wArA

wArAold = wArA;


M'r -> wA

// --- Precondition residual

preconPtr -> precondition(wA, rA, cmpt);


r * M' r

// --- Update search directions :

wArA = gSumProd(wA, rA);


루프의 처음이라면 p = M' r

if (solverPerf.nIterations() == 0)

{

for (register label cell=0; cell<nCells; cell+)

         {

   pAPtr[cell]= wAPtr[cell];

}

}


이외에서는

else

{

scalar beta = wArA/wArAold;

 

for (register label cell=0; cell < nCells; cell++)

{

pAPtr[cell] = wAPtr[cell] + beta*pAPtr[cell];

 }

}


A p -> wA

// --- Update preconditioned residual

matrix_.Amul(wA, pA, interfaceBouCoeffs_, interfaces_, cmpt);


p * Ap

scalar wApA = gSumProd(wA, pA);


alpha = (r * M' r)/(p * A p)

scalar alpha = wArA/wApA;


x = x + alpha p 와 r = r - alpha A p

for (register label cell=0; cell<nCells; cell++)

{

psiPtr[cell] += alpha*pAPtr[cell];

rAPtr[cell] -= alpha*wAPtr[cell];

}


그리고 수렴판정.

while

 ( 

  solvePerf.nIterations()++ < maxIter_

  && !(solverPerf.checkConvergence(tolerance_, relTol_))

  );





  추가정보


프로그램을 전처리가 없는 CG법으로 수정해보았다(11-solver2.tar.gz).

CG.C (일부)

... // --- Solver iteration do { // --- Store previous rArA rArAold = rArA; // --- Update search directions: rArA = gSumProd(rA, rA); if (solverPerf.nIterations() == 0) {

for (register label cell=0; cell<nCells; cell++)

{

pAPtr[cell] = rAPtr[cell];

}

}

else

{

scalar beta = rArA/rArAold;


for (register label cell=0; cell<nCells; cell++)

{

pAPtr[cell] = rAPtr[cell] + beta*pAPtr[cell];

}

}

matrix_.Amul(wA, pA, interfaceBouCoeffs_, interfaces_, cmpt); scalar wApA = gSumProd(wA, pA); // --- Test for singularity if (solverPerf.checkSingularity(mag(wApA)/normFactor)) break; // --- Update solution and residual: scalar alpha = rArA/wApA;


for (register label cell=0; cell<nCells; cell++)

{

psiPtr[cell] += alpha*pAPtr[cell];

rAPtr[cell[ -= alpha*wAPtr[cell];

}


solverPerf.finalResidual() = gSumMag(rA)/normFactor;


}while

(

solvePerf.nIterations()++ < maxIter_

&& !(solverPerf.checkConvergence(tolerance_, relTol_))

);

...

실행해보자.

CG:  Solving for T, Initial residual = 1, Final residual = 8.75773e-10, No Iterations 10


전처리를 하지 않으므로 이터레이션 횟수가... 줄었다?? PCG에서 preconditioner를 none으로 해도 같은 결과를 얻으므로, 이렇게 해도 괜찮은 듯 하다.




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