Changes between Initial Version and Version 1 of cypress/Programming/HeatMassTransfer/ex18.c

Timestamp:

05/14/15 15:20:24 (10 years ago)

Author:

cmaggio

Comment:

—

cypress/Programming/HeatMassTransfer/ex18.c

@@ +1 @@
+= 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<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");
+}
+}}}