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


Ignore:
Timestamp:
05/14/15 15:11:29 (9 years ago)
Author:
cmaggio
Comment:

migrated content from ccs wiki

Legend:

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

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