본문 바로가기

OpenFOAM/소스코드 파해치기

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

일본어 원문링크


작성일 : 2011년    1월   9일
번역일 : 2016년   3월  28일


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


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

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

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

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

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

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



 

  들어가기 


이전의 대수방정식솔버 포스팅에 추가해, 대수방정식의 데이터 구조에 대해 알아보자



  사용 버전 


OpenFOAM 1.7.1

 


   파일


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


  • Make/files
  • Make/options
  • system/controlDict
  • system/fvSchemes
  • system/fvSolution
  • constant/transportProperties
  • constant/polyMesh/blockMeshDict
  • OutputMatrix/OutputAsymMatrix.C
  • OutputMatrix/OutputAsymMatrix.H
  • OutputMatrix/OutputSymMatrix.C
  • OutputMatrix/OutputSymMatrix.H
  • createFields.H
  • icoFoam.C









  대수방정식의 출력

 

대수방정식솔버를 작성하는 요령으로 대수방정식 솔버의 계수행렬과 우변벡터를 출력해보도록 하자.

솔버로는  icoFoam을 사용한다. icoFoam 솔버의 소스코드를 복사하여 그 안에 OutputMatrix라는 폴더를 추가한다. 추가한 폴더 안에 $FOAM_SRC/OpenFOAM/matricies/lduMatrix/solvers/PCG에 있는 PCG.C, PCG.H를 복사해 파일명을 OutputSymMatrix.C, OuputSymMatrix.H로 변경한다. 파일 내 문자열 "PCG"를 OutputSymMatrix"로 치환한다. 같은 방법으로 OutputAsymMatrix.C, OutputAsymMatrix.H를 작성한다. OutputSymMatrix는  p (대칭행렬) 을 출력하기위해,  OutputAsymMatrix는 U (비대칭행렬)을 다룬다.

솔버의 함수 solve() 로 계수행렬을 출력한다. 아래에서는 OutputAsymMatrix의 코드내용을 나타내고 있으나, OutputSymMatrix도 같은 방식이다.



OutputAsymMatrix.C

...
Foam::lduMatrix::solverPerformance Foam::OutputAsymMatrix::solve
(
    scalarField& psi,
    const scalarField& source,
    const direction cmpt
) const
{
    // --- Setup class containing solver performance data
    lduMatrix::solverPerformance solverPerf
    (
        lduMatrix::preconditioner::getName(controlDict_) + typeName,
        fieldName_
    );

    int n = psi.size();
    scalar *a = (scalar *)malloc(sizeof(scalar)*n*n);

    for(int i = 0; i < n; i++){
	    for(int j = 0; j < n; j++){
		    a[i*n + j] = 0.;
	    }
    }

    ...

    for(int i = 0; i < matrix_.lower().size(); i++){
	    int r = matrix_.lduAddr().upperAddr()[i];
	    int c = matrix_.lduAddr().lowerAddr()[i];
	    a[r*n + c] = matrix_.lower()[i];
    }

    for(int i = 0; i < matrix_.diag().size(); i++){
	    a[i*n + i] = matrix_.diag()[i];
    }

    for(int i = 0; i < matrix_.upper().size(); i++){
	    int r = matrix_.lduAddr().lowerAddr()[i];
	    int c = matrix_.lduAddr().upperAddr()[i];
	    a[r*n + c] = matrix_.upper()[i];
    }

    char filename[32];
    sprintf(filename, "matrix-%s.txt", fieldName_.c_str());

    FILE *fp = fopen(filename, "w");
    if(fp == NULL){
        printf("error: Can't open file: OutputAsymMatrix.txt\n");
	std::exit(0);
    }

    fprintf(fp, "%d\n", n);

    for(int i = 0; i < n; i++){
	    for(int j = 0; j < n; j++){
		    fprintf(fp, "%15.7e", a[i*n + j]);
	    }
	    fprintf(fp, "\n");
    }

    for(int i = 0; i < n; i++) fprintf(fp, "%15.7e\n", source[i]);

    fclose(fp);

    free(a);

    return solverPerf;
}

a 라는 포인터에 메모리를 할당하여, 계수행렬을 입력한다. matrix_가 계수행렬의 정보를 가지고 있다. 여기서 laduMatrix(LDU 행렬)은 sparse matrix이므로 0이 아닌 값만을 저장하고 있다.

matrix_.lduAddr(),upperAddr(), matrix_.lduAddr(),lowerAddr()에 의해 각각 위삼각요소, 아래삼각요소의 주소를 얻어낸다. matrix_.upper(), matrix_.lower(), matrix_.diag() 에 의해 각각 위삼각요소, 아래삼각요소, 대칭요소의 값을 얻어낸다. 각각의 정보를 통해 전체행렬 a를 구성하고 있다.

위삼각요소/아래삼각요소의 주소가 어떤것인지 출력해 확인해 보자.

upper 0 1 1 2 2 3 3 3 lower 0 0 1 0 2 1 3 2

위의 값은 icoFoam의 튜토리얼 케이스 cavity의 모델을 2x2 격자로 구성한 결과이다. 이 결과는 아래와 같이 데이터가 저장되어있음을 의미한다.

0 1 2 3 0 D U U . 1 L D . U 2 L . D U 3 . L L D

주소는 0이 아닌 요소의 열번호를 저장하고 있다. 행의 위치는 행렬이 대칭임을 전제로 표현되고 있다. 따라서 아래와 같은 방식으로 a 를 구성하게 된다.

for(int i = 0; i < matrix_.lower().size(); i++){ int r = matrix_.lduAddr().upperAddr()[i]; int c = matrix_.lduAddr().lowerAddr()[i]; a[r*n + c] = matrix_.lower()[i]; } for(int i = 0; i < matrix_.diag().size(); i++){ a[i*n + i] = matrix_.diag()[i]; } for(int i = 0; i < matrix_.upper().size(); i++){ int r = matrix_.lduAddr().lowerAddr()[i]; int c = matrix_.lduAddr().upperAddr()[i]; a[r*n + c] = matrix_.upper()[i]; }


위삼각요소의 주소는 u, 아래삼각요소의 주소는 l 이라고 한다면, 아래삼각요소를 구성할 때는 a(u,l), 위삼각요소를 구성할때는 a(l,u)으로 나타낼 수 있다.

우변벡터에 대해서는 source로부터 값을 추출하면 된다.

대수방정식 솔버의 작성을 완성하면, system/fvSolution을 설정하여 사용할 수 있다.

system/fvSolution

... solvers { p { //solver PCG; //preconditioner DIC; solver OutputSymMatrix; preconditioner none; tolerance 1e-06; relTol 0; } U { //solver PBiCG; //preconditioner DILU; solver OutputAsymMatrix; preconditioner none; tolerance 1e-05; relTol 0; } } ...




  실행


컴파일하여 실행해 본다.

$ wmake

$ blockMesh $ ./icoFoam



대수방정식을 출력하는것 뿐이므로, 실행은 1스텝만 수행해도 된다. matrix-*.txt 라는 형태로 얻어진다.







  추가정보

 

하는김에 자신만의 대수방정식 솔버를 작성해 보자.(13-solver2.tar.gz).

대수방정식을 출력하는 것과 같은 방법으로 정보를 얻어내고, 그것을 별도 모듈의 솔버(여기서는 자코비법 솔버)에 전달하여 받은 해(벡터)를 설정하고 있다.


AsymJacobi.C


...
Foam::lduMatrix::solverPerformance Foam::AsymJacobi::solve
(
    scalarField& psi,
    const scalarField& source,
    const direction cmpt
) const
{
    // --- Setup class containing solver performance data
    lduMatrix::solverPerformance solverPerf
    (
        lduMatrix::preconditioner::getName(controlDict_) + typeName,
        fieldName_
    );

    int n = psi.size();
    scalar *a = (scalar *)malloc(sizeof(scalar)*n*n);
    scalar *b = (scalar *)malloc(sizeof(scalar)*n);
    scalar *x = (scalar *)malloc(sizeof(scalar)*n);

    for(int i = 0; i < n; i++){
	    for(int j = 0; j < n; j++){
		    a[i*n + j] = 0.;
	    }
    }

    for(int i = 0; i < matrix_.lower().size(); i++){
	    int r = matrix_.lduAddr().upperAddr()[i];
	    int c = matrix_.lduAddr().lowerAddr()[i];
	    a[r*n + c] = matrix_.lower()[i];
    }

    for(int i = 0; i < matrix_.diag().size(); i++){
	    a[i*n + i] = matrix_.diag()[i];
    }

    for(int i = 0; i < matrix_.upper().size(); i++){
	    int r = matrix_.lduAddr().lowerAddr()[i];
	    int c = matrix_.lduAddr().upperAddr()[i];
	    a[r*n + c] = matrix_.upper()[i];
    }

    for(int i = 0; i < n; i++) b[i] = source[i];
    for(int i = 0; i < n; i++) x[i] = psi[i];

    scalar r0, r;
    int iter;
    jacobi(n, a, b, x, tolerance_, maxIter_, &r0, &r, &iter);
    
    for(int i = 0; i < n; i++) psi[i] = x[i];

    free(a);
    free(b);
    free(x);

    solverPerf.initialResidual() = r0;
    solverPerf.finalResidual() = r;
    solverPerf.nIterations() = iter;

    return solverPerf;
}


psi 가 미지변수 벡터에 해당한다. 대수방정식솔버에서 받은 해 x를 psi에 대입시키고 있다.


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