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

migrated content from ccs wiki

# Ex19.c

```/*
Large Scale Computing
2D Heat/Mass Transfer
ex19.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 100  /* number of points for x and y*/
#define TOL 1.0e-12  /* convergence tolerance */

int main(void ){
int i,j,n;
double x,y,dx;
double cnc[NUM * NUM]; /* define 1D array of length NUM * NUM */
double r[NUM * NUM];
double p[NUM * NUM];
double Ap[NUM * NUM];
double rr,pAp,rr1,alpha,beta;
int k;
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;
r[n] = 0.0;
}

/* computing for C with Conjugate Gradient Method */
/* r0 = b - Ax0 */
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*/
r[n] = -dx * dx;
r[n] -= -4.0 * cnc[n];
r[n] -= cnc[n + 1];
r[n] -= cnc[n + NUM];
r[n] -= cnc[n - 1];
r[n] -= cnc[n - NUM];
}

rr = 0.0;
for (n = 0 ; n < NUM * NUM ; n++){
p[n] = r[n];    /* p = r */
rr += r[n] * r[n];  /* rr = r.r */
}

k = 0;
while(rr > TOL){
// Ap = A*p
for (n = 0 ; n < NUM * NUM ; n++){
Ap[n] = 0.0;
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*/

/* computing A*p */
Ap[n] = -4.0 * p[n];
Ap[n] += p[n + 1];
Ap[n] += p[n + NUM];
Ap[n] += p[n - 1];
Ap[n] += p[n - NUM];
}

// alpha = r.r / p.Ap
pAp = 0.0;
for (n = 0 ; n < NUM * NUM ; n++){
pAp += p[n] * Ap[n];
}
alpha = rr / pAp;

//Beta
rr1 = 0.0;
for (n = 0 ; n < NUM * NUM ; n++){
cnc[n] += alpha * p[n];
r[n] -= alpha * Ap[n];
rr1 += r[n] * r[n];
}

beta = rr1 / rr;

for (n = 0 ; n < NUM * NUM ; n++){
p[n] = r[n] + beta * p[n];
}

rr = rr1;
k++;
}

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");
}
```
Note: See TracWiki for help on using the wiki.