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


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

migrated content from ccs wiki

Legend:

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

    v1 v1  
     1= Ex17.c =
     2{{{#!C
     3/*
     4  Large Scale Computing
     5  2D Heat/Mass Transfer
     6  ex17.c : Numerical Analysis
     7  Solve for
     8  d2c/dx2 + d2c/dy2 + 1 = 0
     9  With the boundary conditions of c = 0
     10  along lines of x=1,-1, and y=1, and -1
     11 */
     12#include<stdio.h>
     13#include<stdlib.h>
     14#include<math.h>
     15 
     16#define NUM 20  /* number of points for x and y*/
     17#define TOL 1.0e-12  /* convergence tolerance */
     18 
     19int main(void ){
     20  int i,j,k,n,ie,iw,in,is;
     21  double x,y,dx,w,err;
     22  double cnc[NUM * NUM]; /* define 1D array of length NUM * NUM */
     23  FILE *fp;
     24 
     25  /*assuming dx = dy : domain is 2x2 */
     26  dx = 2.0 / (double)(NUM - 1);
     27 
     28  /* Initialize */
     29  for (n = 0 ; n < NUM * NUM ; n++){
     30    cnc[n] = 0.0;
     31  }
     32 
     33  /* computing for C with Gauss-Seidel */
     34  k = 0;
     35  while(1){
     36    err = 0.0;
     37    for (n = 0 ; n < NUM * NUM ; n++){
     38      i = n % NUM; /* i is residual */
     39      j = n / NUM; /* j is division */
     40      if (i == 0) continue; /*skip left end*/
     41      if (i == NUM-1) continue; /*skip right end*/
     42      if (j == 0) continue; /*skip bottom end*/
     43      if (j == NUM-1) continue; /*skip top end*/
     44 
     45      ie = n + 1;
     46      iw = n - 1;
     47      in = n + NUM;
     48      is = n - NUM;
     49      w = 0.25 * (dx * dx + cnc[ie] + cnc[iw] + cnc[in] + cnc[is]);
     50      err += (w - cnc[n]) * (w - cnc[n]);
     51      cnc[n] = w;
     52    }
     53    k++;
     54    if (err < TOL) break;
     55  }
     56  printf("# of Iteration=%d\n",k);
     57 
     58  /* Output Result */
     59  fp = fopen("res.dat","w");
     60  if (fp == NULL){
     61    printf("File could not create\n");
     62    exit(-1);
     63  }
     64 
     65  for (i = 0 ;i < NUM ; i++){
     66    for (j = 0 ; j < NUM ; j++){
     67      x = -1.0 + 2.0 * (double)i / (double)(NUM - 1);
     68      y = -1.0 + 2.0 * (double)j / (double)(NUM - 1);
     69      fprintf(fp,"%e %e %e\n",x,y,cnc[i + NUM * j]);
     70    }
     71  }
     72  fclose(fp);
     73 
     74  printf("Done\n");
     75}
     76}}}