작성일 : 2011년 1월 9일
번역일 : 2016년 3월 28일
크롬 브라우저로 보시는 것을 권장해 드립니다.
------------------------------------------------------------
----- 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 > 소스코드 파해치기' 카테고리의 다른 글
| [OpenFOAM 소스코드 파해치기] 열물성 다루기 (0) | 2016.03.28 |
|---|---|
| [OpenFOAM 소스코드 파해치기] simpleFoam (2) | 2016.03.28 |
| [OpenFOAM 소스코드 파해치기] 데이터읽기 (4) | 2016.02.22 |
| [OpenFOAM 소스코드 파해치기]대수방정식솔버 (0) | 2016.02.14 |
| [OpenFOAM 소스코드 파해치기] 시간미분 (0) | 2016.02.12 |