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


Ignore:
Timestamp:
May 14, 2015 3:20:24 PM (7 years ago)
Author:
cmaggio
Comment:

Legend:

Unmodified
Added
Removed
Modified
  • cypress/Programming/HeatMassTransfer/ex18.c

    v1 v1  
     1= Ex18.c =
     2
     3A parallel version of ex17.c
     4{{{#!C
     5/*
     6  Large Scale Computing
     7  2D Heat/Mass Transfer
     8  ex18.c : Numerical Analysis
     9  Solve for
     10  d2c/dx2 + d2c/dy2 + 1 = 0
     11  With the boundary conditions of c = 0
     12  along lines of x=1,-1, and y=1, and -1
     13 
     14  Parallel Version of ex17.c
     15  Red-Black Gauss-Seidel
     16 */
     17#include<stdio.h>
     18#include<stdlib.h>
     19#include<math.h>
     20#include<omp.h>
     21 
     22#define NUM 100  /* number of points for x and y*/
     23#define TOL 1.0e-12  /* convergence tolerance */
     24 
     25int main(void ){
     26  int i,j,k,n,ie,iw,in,is;
     27  double x,y,dx,w,err;
     28  double cnc[NUM * NUM]; /* define 1D array of length NUM * NUM */
     29  int rb;
     30  FILE *fp;
     31 
     32  /*assuming dx = dy : domain is 2x2 */
     33  dx = 2.0 / (double)(NUM - 1);
     34 
     35  /* Initialize */
     36  for (n = 0 ; n < NUM * NUM ; n++){
     37    cnc[n] = 0.0;
     38  }
     39 
     40  /* computing for C with Gauss-Seidel */
     41  k = 0;
     42  while(1){
     43    err = 0.0;
     44#pragma omp parallel private(w,rb,i,j,ie,iw,in,is)   
     45    for (rb = 0 ; rb <= 1 ; rb++){ /* rb = 0 or 1 */
     46#pragma omp for reduction (+:err)   
     47      for (n = 0 ; n < NUM * NUM ; n++){
     48        i = n % NUM; /* i is residual */
     49        j = n / NUM; /* j is division */
     50        if (i == 0) continue; /*skip left end*/
     51        if (i == NUM-1) continue; /*skip right end*/
     52        if (j == 0) continue; /*skip bottom end*/
     53        if (j == NUM-1) continue; /*skip top end*/
     54        if ((i + j) % 2 == rb){ /* rb=0 (red) rb=1 (black) */
     55          ie = n + 1;
     56          iw = n - 1;
     57          in = n + NUM;
     58          is = n - NUM;
     59          w = 0.25 * (dx * dx + cnc[ie] + cnc[iw] + cnc[in] + cnc[is]);
     60          err += (w - cnc[n]) * (w - cnc[n]);
     61          cnc[n] = w;
     62        }
     63      }
     64#pragma omp barrier /* wait until all come here*/
     65    } /* end of parallel scope */
     66 
     67    k++;
     68    if (err < TOL) break;
     69  }
     70  printf("# of Iteration=%d\n",k);
     71 
     72  /* Output Result */
     73  fp = fopen("res.dat","w");
     74  if (fp == NULL){
     75    printf("File could not create\n");
     76    exit(-1);
     77  }
     78 
     79  for (i = 0 ;i < NUM ; i++){
     80    for (j = 0 ; j < NUM ; j++){
     81      x = -1.0 + 2.0 * (double)i / (double)(NUM - 1);
     82      y = -1.0 + 2.0 * (double)j / (double)(NUM - 1);
     83      fprintf(fp,"%e %e %e\n",x,y,cnc[i + NUM * j]);
     84    }
     85  }
     86  fclose(fp);
     87 
     88  printf("Done\n");
     89}
     90}}}