Ex23.c

/*
  Large Scale Computing
  Convective Heat/Mass Transfer
  ex23.c : Gauss-Seidel 
 
  duphi/dx = d2phi/dx2
 */
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
 
/* min and max macros */
#define max(a,b)        (((a) > (b)) ? (a) : (b))
#define min(a,b)        (((a) < (b)) ? (a) : (b))
 
int main(int argc, char **argv){
  double a,b,d,u;
  double *phi;
  double w;
  double dx,x,sol;
  double err;
  double ae,aw,ap;
  double tol = 1.0e-10;
  int i,num,itc;
  FILE *fp;
 
  if (argc < 6){
    printf("Usage:%s [NUM] [A] [B] [D] [U]\n",argv[0]);
    exit(-1);
  }
 
  /* set parameters */
  num = atoi(argv[1]);
  a = atof(argv[2]);
  b = atof(argv[3]);
  d = atof(argv[4]);
  u = atof(argv[5]);
 
  printf("num=%d A=%e B=%e D=%e U=%e\n",num,a,b,d,u);
 
  /* Memory Allocation and Check*/
  phi = (double *)calloc(num, sizeof(double));
  if (phi == (double *)NULL){
    printf("Memory Allocation Failed\n");
    exit(-1);    
  }
 
  /** Setup **/
  dx = 1.0 / (double)(num - 1);
  for (i = 0 ; i < num ; i++){
    phi[i] = 0.0;
  }
  ae = d / dx + max(-u,0.0);
  aw = d / dx + max(u,0.0);
  ap = ae + aw;
 
  /* Boundary Condition*/
  phi[0] = a;
  phi[num-1] = b;
 
  /* Solve with Gauss-Seidel Method*/
  itc = 0;
  while(1){
    /* update phi   &  Check Convergence */
    err = 0.0;    
    for (i = 1 ; i < (num - 1) ; i++){
      w =  (ae *  phi[i + 1] + aw *  phi[i - 1]) / ap;
      err += (w - phi[i]) * (w - phi[i]);      
      phi[i] = w;
    }
 
    itc++;
    if (err < tol) break;
    if ((itc % 100) == 0) printf("itc=%d err=%e\n",itc,err);
  }
 
  printf("Number of Iteration=%d\n",itc);
 
  /* Output Result */
  fp = fopen("res.dat","w");
  if (fp == NULL){
    printf("File could not create\n");
    exit(-1);
  }
 
  for (i = 0 ; i < num ; i++){
    x = dx * (double)i;
    sol = a + (b - a) * (exp(x * u / d) - 1.0) / (exp(u / d) - 1.0);
    fprintf(fp,"%e %e %e\n",x,phi[i],sol);
  }
  fclose(fp);
 
  /* END : Deallocate Memory */
  free(phi);
  printf("Done\n");
}