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


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

migrated content from ccs wiki

Legend:

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

    v1 v1  
     1= Ex4.c =
     2{{{#!C
     3/*
     4  Large Scale Computing
     5  Heat/Mass Transfer
     6  ex4.c : Gauss-Seidel
     7 */
     8#include<stdio.h>
     9#include<stdlib.h>
     10#include<math.h>
     11 
     12int main(int argc, char **argv){
     13  double a,b,d;
     14  double *phi;
     15  double w;
     16  double *rhs;
     17  double dx;
     18  double err,merr;
     19  double tol = 1.0e-8;
     20  int i,num,itc;
     21  FILE *fp;
     22 
     23  if (argc < 5){
     24    printf("Usage:%s [NUM] [A] [B] [D]\n",argv[0]);
     25    exit(-1);
     26  }
     27 
     28  /* set parameters */
     29  num = atoi(argv[1]);
     30  a = atof(argv[2]);
     31  b = atof(argv[3]);
     32  d = atof(argv[4]);
     33 
     34  printf("num=%d A=%e B=%e D=%e\n",num,a,b,d);
     35 
     36  /* Memory Allocation and Check*/
     37  phi = (double *)calloc(num, sizeof(double));
     38  if (phi == (double *)NULL){
     39    printf("Memory Allocation Failed\n");
     40    exit(-1);   
     41  }
     42 
     43  rhs = (double *)calloc(num, sizeof(double));
     44  if (rhs == (double *)NULL){
     45    printf("Memory Allocation Failed\n");
     46    exit(-1);   
     47  } 
     48 
     49  /** Setup **/
     50  dx = 1.0 / (double)(num - 1);
     51  for (i = 0 ; i < num ; i++){
     52    rhs[i] = -dx * dx / d * (dx * (double)i);
     53    phi[i] = 0.0;
     54  }
     55 
     56  /* Boundary Condition*/
     57  phi[0] = a;
     58  phi[num-1] = b;
     59 
     60  /* Solve with Gauss-Seidel Method*/
     61  itc = 0;
     62  while(1){
     63    /* update phi   &  Check Convergence */
     64    merr = 0.0;   
     65    for (i = 1 ; i < (num - 1) ; i++){
     66      w = -0.5 * (rhs[i] - phi[i + 1] - phi[i - 1]);
     67      err = fabs(w - phi[i]);     
     68      phi[i] = w;
     69      if (merr < err) merr = err;
     70    }
     71 
     72    itc++;
     73    if (merr < tol) break;
     74    if ((itc % 10000) == 0) printf("itc=%d err=%e\n",itc,merr);
     75  }
     76 
     77  printf("Number of Iteration=%d\n",itc);
     78 
     79  /* Output Result */
     80  fp = fopen("res.dat","w");
     81  if (fp == NULL){
     82    printf("File could not create\n");
     83    exit(-1);
     84  }
     85 
     86  for (i = 0 ; i < num ; i++){
     87    fprintf(fp,"%e %e\n",dx * (double)i,phi[i]);
     88  }
     89  fclose(fp);
     90 
     91  /* END : Deallocate Memory */
     92  free(phi);
     93  free(rhs);
     94  printf("Done\n");
     95}
     96}}}