[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 소스코드 파해치기 목차로 이동