Ex17.c
/*
Large Scale Computing
2D Heat/Mass Transfer
ex17.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>
#define NUM 20 /* number of points for x and y*/
#define TOL 1.0e-12 /* convergence tolerance */
int main(void ){
int i,j,k,n,ie,iw,in,is;
double x,y,dx,w,err;
double cnc[NUM * NUM]; /* define 1D array of length NUM * NUM */
FILE *fp;
/*assuming dx = dy : domain is 2x2 */
dx = 2.0 / (double)(NUM - 1);
/* Initialize */
for (n = 0 ; n < NUM * NUM ; n++){
cnc[n] = 0.0;
}
/* computing for C with Gauss-Seidel */
k = 0;
while(1){
err = 0.0;
for (n = 0 ; n < NUM * NUM ; n++){
i = n % NUM; /* i is residual */
j = n / NUM; /* j is division */
if (i == 0) continue; /*skip left end*/
if (i == NUM-1) continue; /*skip right end*/
if (j == 0) continue; /*skip bottom end*/
if (j == NUM-1) continue; /*skip top end*/
ie = n + 1;
iw = n - 1;
in = n + NUM;
is = n - NUM;
w = 0.25 * (dx * dx + cnc[ie] + cnc[iw] + cnc[in] + cnc[is]);
err += (w - cnc[n]) * (w - cnc[n]);
cnc[n] = w;
}
k++;
if (err < TOL) break;
}
printf("# of Iteration=%d\n",k);
/* 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,cnc[i + NUM * j]);
}
}
fclose(fp);
printf("Done\n");
}