= Ex19.c = {{{#!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 Conjugate Gradient Method */ #include #include #include #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*/ if (j == NUM-1) continue; /*skip top 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*/ if (j == NUM-1) continue; /*skip top 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"); } }}}