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


Ignore:
Timestamp:
May 14, 2015 3:12:59 PM (4 years ago)
Author:
cmaggio
Comment:

migrated content from ccs wiki

Legend:

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

    v1 v1  
     1= Ex13.c =
     2a parallel version of ex3.c
     3{{{#!C
     4/*
     5  Large Scale Computing
     6  Heat/Mass Transfer
     7  ex13.c
     8 */
     9#include<stdio.h>
     10#include<stdlib.h>
     11#include<math.h>
     12#include<omp.h>
     13 
     14int main(int argc, char **argv){
     15  double a,b,d;
     16  double *phi,*phi0;
     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  phi0 = (double *)calloc(num, sizeof(double));
     46  if (phi0 == (double *)NULL){
     47    printf("Memory Allocation Failed\n");
     48    exit(-1);   
     49  }
     50 
     51  rhs = (double *)calloc(num, sizeof(double));
     52  if (rhs == (double *)NULL){
     53    printf("Memory Allocation Failed\n");
     54    exit(-1);   
     55  } 
     56 
     57  /** Setup **/
     58  dx = 1.0 / (double)(num - 1);
     59  for (i = 0 ; i < num ; i++){
     60    rhs[i] = -dx * dx / d * (dx * (double)i);
     61    phi[i] = phi0[i] = 0.0;
     62  }
     63 
     64  /* Boundary Condition*/
     65  phi[0] = phi0[0] = a;
     66  phi[num-1] = phi0[num-1] = b;
     67 
     68  /* Solve with Jacobi Method*/
     69  itc = 0;
     70 
     71  while(1){
     72    share_err = 0.0;
     73 
     74#pragma omp parallel private(merr,err)
     75    {
     76      nthreads = omp_get_num_threads(); 
     77 
     78      /* update phi */
     79#pragma omp for
     80      for (i = 1 ; i < (num - 1) ; i++){
     81        phi[i] = -0.5 * (rhs[i] - phi0[i + 1] - phi0[i - 1]);
     82      }
     83 
     84      /* Check Convergence */
     85      merr = 0.0;
     86#pragma omp for
     87      for (i = 0 ; i < num ; i++){
     88        err = fabs(phi[i] - phi0[i]);
     89        if (merr < err) merr = err;
     90        phi0[i] = phi[i];
     91      }
     92 
     93#pragma omp critical
     94      {
     95        /* find max in whole threads */
     96        if (share_err < merr) share_err = merr;
     97      }
     98    }
     99    itc++;
     100 
     101    if (share_err < tol) break;
     102    if ((itc % 10000) == 0) printf("itc=%d err=%e\n",itc,share_err);
     103  }
     104  printf("Number of Iteration=%d\n",itc);
     105  printf("Number of Threads=%d\n",nthreads);
     106 
     107 
     108  /* Output Result */
     109  fp = fopen("res.dat","w");
     110  if (fp == NULL){
     111    printf("File could not create\n");
     112    exit(-1);
     113  }
     114 
     115  for (i = 0 ; i < num ; i++){
     116    fprintf(fp,"%e %e\n",dx * (double)i,phi[i]);
     117  }
     118  fclose(fp);
     119 
     120  /* END : Deallocate Memory */
     121  free(phi);
     122  free(phi0);
     123  free(rhs);
     124  printf("Done\n");
     125}
     126}}}