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

Timestamp:

05/14/15 15:12:59 (10 years ago)

Author:

cmaggio

Comment:

migrated content from ccs wiki

cypress/Programming/HeatMassTransfer/ex13.c

@@ +1 @@
+= Ex13.c = 
+a parallel version of ex3.c
+{{{#!C
+/*
+  Large Scale Computing
+  Heat/Mass Transfer
+  ex13.c 
+ */
+#include<stdio.h>
+#include<stdlib.h>
+#include<math.h>
+#include<omp.h>
+ 
+int main(int argc, char **argv){
+  double a,b,d;
+  double *phi,*phi0;
+  double *rhs;
+  double dx;
+  double err,merr,share_err;
+  double tol = 1.0e-8;
+  int i,num,itc;
+  int nthreads;
+  FILE *fp;
+ 
+  if (argc < 5){
+    printf("Usage:%s [NUM] [A] [B] [D]\n",argv[0]);
+    exit(-1);
+  }
+ 
+  /* set parameters */
+  num = atoi(argv[1]);
+  a = atof(argv[2]);
+  b = atof(argv[3]);
+  d = atof(argv[4]);
+ 
+  printf("num=%d A=%e B=%e D=%e\n",num,a,b,d);
+ 
+  /* Memory Allocation and Check*/
+  phi = (double *)calloc(num, sizeof(double));
+  if (phi == (double *)NULL){
+    printf("Memory Allocation Failed\n");
+    exit(-1);    
+  }
+ 
+  phi0 = (double *)calloc(num, sizeof(double));
+  if (phi0 == (double *)NULL){
+    printf("Memory Allocation Failed\n");
+    exit(-1);    
+  }
+ 
+  rhs = (double *)calloc(num, sizeof(double));
+  if (rhs == (double *)NULL){
+    printf("Memory Allocation Failed\n");
+    exit(-1);    
+  }  
+ 
+  /** Setup **/
+  dx = 1.0 / (double)(num - 1);
+  for (i = 0 ; i < num ; i++){
+    rhs[i] = -dx * dx / d * (dx * (double)i); 
+    phi[i] = phi0[i] = 0.0;
+  }
+ 
+  /* Boundary Condition*/
+  phi[0] = phi0[0] = a;
+  phi[num-1] = phi0[num-1] = b;
+ 
+  /* Solve with Jacobi Method*/
+  itc = 0;
+ 
+  while(1){
+    share_err = 0.0;
+ 
+#pragma omp parallel private(merr,err)
+    {
+      nthreads = omp_get_num_threads();  
+ 
+      /* update phi */
+#pragma omp for
+      for (i = 1 ; i < (num - 1) ; i++){
+        phi[i] = -0.5 * (rhs[i] - phi0[i + 1] - phi0[i - 1]);
+      }
+ 
+      /* Check Convergence */
+      merr = 0.0;
+#pragma omp for 
+      for (i = 0 ; i < num ; i++){
+        err = fabs(phi[i] - phi0[i]);
+        if (merr < err) merr = err;
+        phi0[i] = phi[i];
+      }
+ 
+#pragma omp critical
+      {
+        /* find max in whole threads */
+        if (share_err < merr) share_err = merr;
+      }
+    }
+    itc++;
+ 
+    if (share_err < tol) break;
+    if ((itc % 10000) == 0) printf("itc=%d err=%e\n",itc,share_err);
+  }
+  printf("Number of Iteration=%d\n",itc);
+  printf("Number of Threads=%d\n",nthreads);
+ 
+ 
+  /* Output Result */
+  fp = fopen("res.dat","w");
+  if (fp == NULL){
+    printf("File could not create\n");
+    exit(-1);
+  }
+ 
+  for (i = 0 ; i < num ; i++){
+    fprintf(fp,"%e %e\n",dx * (double)i,phi[i]);
+  }
+  fclose(fp);
+ 
+  /* END : Deallocate Memory */
+  free(phi);
+  free(phi0);
+  free(rhs);
+  printf("Done\n");
+}
+}}}