= Ex18.c = A parallel version of ex17.c {{{#!C /* Large Scale Computing 2D Heat/Mass Transfer ex18.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 Parallel Version of ex17.c Red-Black Gauss-Seidel */ #include #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,k,n,ie,iw,in,is; double x,y,dx,w,err; double cnc[NUM * NUM]; /* define 1D array of length NUM * NUM */ int rb; 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; #pragma omp parallel private(w,rb,i,j,ie,iw,in,is) for (rb = 0 ; rb <= 1 ; rb++){ /* rb = 0 or 1 */ #pragma omp for reduction (+:err) 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*/ if ((i + j) % 2 == rb){ /* rb=0 (red) rb=1 (black) */ 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; } } #pragma omp barrier /* wait until all come here*/ } /* end of parallel scope */ 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"); } }}}