Ex18.c

A parallel version of ex17.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<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<omp.h>
 
#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");
}