Version 1 (modified by cmaggio, 7 years ago) (diff)

migrated content from ccs wiki

# Ex20.c

```/*
Large Scale Computing
2D Heat/Mass Transfer
Ex20.c : Numerical Analysis
Solve for
d2c/dx2 + d2c/dy2 + 1 = 0
With the boundary conditions of c = 0
along lines of x=1,-1, and y=1, and -1

*/
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include"cg.h"

#define NUM 100  /* number of points for x and y*/
#define TOL 1.0e-12  /* convergence tolerance */

int index(int ii,int jj);
int main(void ){
int i,j,k,n;
double x,y,dx;
double phi[NUM][NUM]; /* define 1D array of length NUM * NUM */
double rhs[NUM * NUM]; /* define 1D array of length NUM * NUM */
double sol[NUM * NUM]; /* define 1D array of length NUM * NUM */
double nzval[5 * NUM * NUM]; /* non-zero elements of A */
int colind[5 * NUM * NUM]; /* column indices of nonzeros */
int rowptr[NUM * NUM + 1]; /**/
int nnz;
FILE *fp;

/*assuming dx = dy : domain is 2x2 */
dx = 2.0 / (double)(NUM - 1);

/* initial & boundary conditions */
for (j = 0 ; j < NUM  ; j++){
for (i = 0 ; i < NUM ; i++){
phi[i][j] = 0.0;
}
}

/* Setup Matrix A and RHS */
nnz = 0;
for (j = 1 ; j < NUM - 1 ; j++){
for (i = 1 ; i < NUM - 1 ; i++){
/* initial guess */
sol[index(i,j)] = phi[i][j];

/* set general RHS */
rhs[index(i,j)] = -dx * dx;

/* set first nonzero column of row index(i,j) */
rowptr[index(i,j)] = nnz;

/* sounth */
if (j == 1){
rhs[index(i,j)] -= phi[i];
}else{
nzval[nnz] = 1.0;
colind[nnz] = index(i,j-1);
nnz++;
}

/* west */
if (i == 1){
rhs[index(i,j)] -= phi[j];
}else{
nzval[nnz] = 1.0;
colind[nnz] = index(i-1,j);
nnz++;
}

/* diagonal Element */
nzval[nnz] = -4.0;
colind[nnz] = index(i,j);
nnz++;

/* east */
if (i == NUM - 2){
rhs[index(i,j)] -= phi[NUM-1][j];
}else{
nzval[nnz] = 1.0;
colind[nnz] = index(i+1,j);
nnz++;
}

/* north */
if (j == NUM - 2){
rhs[index(i,j)] -= phi[i][NUM - 1];
}else{
nzval[nnz] = 1.0;
colind[nnz] = index(i,j+1);
nnz++;
}
}
}

/* last element of rowptr si nnz */
rowptr[(NUM-2) * (NUM-2)] = nnz;

/* solve with CG */
k = solve_cg((NUM-2) * (NUM-2), nnz, nzval, colind, rowptr, sol, rhs, TOL);
if (k < 0){
printf("calculation failed\n");
exit(-1);
}
printf("Number of Interation=%d\n",k);

/* store solution to main aray */
for (j = 1 ; j < NUM - 1 ; j++){
for (i = 1 ; i < NUM - 1 ; i++){
phi[i][j] = sol[index(i,j)];
}
}

/* Output Result */
fp = fopen("res.dat","w");
if (fp == NULL){
printf("File could not create\n");
exit(-1);
}

for (i = 0 ;i < NUM ; i++){
for (j = 0 ; j < NUM ; j++){
x = -1.0 + 2.0 * (double)i / (double)(NUM - 1);
y = -1.0 + 2.0 * (double)j / (double)(NUM - 1);
fprintf(fp,"%e %e %e\n",x,y,phi[i][j]);
}
}
fclose(fp);

printf("Done\n");
}

/* index for interior Points */
int index(int ii,int jj){
return ii - 1 + (NUM - 2) * (jj - 1);
}
```