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


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

Legend:

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

    v1 v1  
     1= Ex15.c =
     2A parallel version of ex4.c
     3{{{#!C
     4/*
     5  Large Scale Computing
     6  Heat/Mass Transfer
     7  ex15.c : Red-Black Gauss-Seidel
     8 */
     9#include<stdio.h>
     10#include<stdlib.h>
     11#include<math.h>
     12 
     13int main(int argc, char **argv){
     14  double a,b,d;
     15  double *phi;
     16  double w;
     17  double *rhs;
     18  double dx;
     19  double err,merr,share_err;
     20  double tol = 1.0e-8;
     21  int i,num,itc;
     22  int nthreads;
     23  FILE *fp;
     24 
     25  if (argc < 5){
     26    printf("Usage:%s [NUM] [A] [B] [D]\n",argv[0]);
     27    exit(-1);
     28  }
     29 
     30  /* set parameters */
     31  num = atoi(argv[1]);
     32  a = atof(argv[2]);
     33  b = atof(argv[3]);
     34  d = atof(argv[4]);
     35 
     36  printf("num=%d A=%e B=%e D=%e\n",num,a,b,d);
     37 
     38  /* Memory Allocation and Check*/
     39  phi = (double *)calloc(num, sizeof(double));
     40  if (phi == (double *)NULL){
     41    printf("Memory Allocation Failed\n");
     42    exit(-1);   
     43  }
     44 
     45  rhs = (double *)calloc(num, sizeof(double));
     46  if (rhs == (double *)NULL){
     47    printf("Memory Allocation Failed\n");
     48    exit(-1);   
     49  } 
     50 
     51  /** Setup **/
     52  dx = 1.0 / (double)(num - 1);
     53  for (i = 0 ; i < num ; i++){
     54    rhs[i] = -dx * dx / d * (dx * (double)i);
     55    phi[i] = 0.0;
     56  }
     57 
     58  /* Boundary Condition*/
     59  phi[0] = a;
     60  phi[num-1] = b;
     61 
     62  /* Solve with Gauss-Seidel Method*/
     63  itc = 0;
     64  while(1){
     65    share_err = 0.0;
     66 
     67#pragma omp parallel private(merr,err,w)
     68    {
     69      nthreads = omp_get_num_threads(); 
     70      /* update phi   &  Check Convergence */
     71      merr = 0.0;   
     72#pragma omp for /* Black (odd) */
     73      for (i = 1 ; i < (num - 1) ; i+=2){
     74        w = -0.5 * (rhs[i] - phi[i + 1] - phi[i - 1]);
     75        err = fabs(w - phi[i]);     
     76        phi[i] = w;
     77        if (merr < err) merr = err;
     78      }
     79 
     80#pragma omp barrier /* wait until all come here*/
     81 
     82#pragma omp for /* Red (even) */
     83      for (i = 2 ; i < (num - 1) ; i+=2){
     84        w = -0.5 * (rhs[i] - phi[i + 1] - phi[i - 1]);
     85        err = fabs(w - phi[i]);     
     86        phi[i] = w;
     87        if (merr < err) merr = err;
     88      }
     89 
     90#pragma omp critical
     91      {
     92        /* find max in whole threads */
     93        if (share_err < merr) share_err = merr;
     94      }
     95    }
     96 
     97    itc++;
     98    if (share_err < tol) break;
     99    if ((itc % 10000) == 0) printf("itc=%d err=%e\n",itc,share_err);
     100  }
     101 
     102  printf("Number of Iteration=%d\n",itc);
     103   printf("Number of Threads=%d\n",nthreads);
     104 
     105  /* Output Result */
     106  fp = fopen("res.dat","w");
     107  if (fp == NULL){
     108    printf("File could not create\n");
     109    exit(-1);
     110  }
     111 
     112  for (i = 0 ; i < num ; i++){
     113    fprintf(fp,"%e %e\n",dx * (double)i,phi[i]);
     114  }
     115  fclose(fp);
     116 
     117  /* END : Deallocate Memory */
     118  free(phi);
     119  free(rhs);
     120  printf("Done\n");
     121}
     122}}}