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


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

migrated content from ccs wiki

Legend:

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

    v1 v1  
     1= Ex19.c =
     2{{{#!C
     3/*
     4  Large Scale Computing
     5  2D Heat/Mass Transfer
     6  ex19.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  Conjugate Gradient Method
     13 */
     14#include<stdio.h>
     15#include<stdlib.h>
     16#include<math.h>
     17 
     18#define NUM 100  /* number of points for x and y*/
     19#define TOL 1.0e-12  /* convergence tolerance */
     20 
     21int main(void ){
     22  int i,j,n;
     23  double x,y,dx;
     24  double cnc[NUM * NUM]; /* define 1D array of length NUM * NUM */
     25  double r[NUM * NUM];
     26  double p[NUM * NUM];
     27  double Ap[NUM * NUM];
     28  double rr,pAp,rr1,alpha,beta;
     29  int k;
     30  FILE *fp;
     31 
     32  /*assuming dx = dy : domain is 2x2 */
     33  dx = 2.0 / (double)(NUM - 1);
     34 
     35  /* Initialize */
     36  for (n = 0 ; n < NUM * NUM ; n++){
     37    cnc[n] = 0.0;
     38    r[n] = 0.0;
     39  }
     40 
     41  /* computing for C with Conjugate Gradient Method */
     42  /* r0 = b - Ax0 */
     43  for (n = 0 ; n < NUM * NUM ; n++){
     44    i = n % NUM; /* i is residual */
     45    j = n / NUM; /* j is division */
     46    if (i == 0) continue; /*skip left end*/
     47    if (i == NUM-1) continue; /*skip right end*/
     48    if (j == 0) continue; /*skip bottom end*/
     49    if (j == NUM-1) continue; /*skip top end*/
     50    r[n] = -dx * dx;
     51    r[n] -= -4.0 * cnc[n];
     52    r[n] -= cnc[n + 1];
     53    r[n] -= cnc[n + NUM];
     54    r[n] -= cnc[n - 1];
     55    r[n] -= cnc[n - NUM];
     56  }
     57 
     58  rr = 0.0;
     59  for (n = 0 ; n < NUM * NUM ; n++){
     60    p[n] = r[n];    /* p = r */
     61    rr += r[n] * r[n];  /* rr = r.r */
     62  }
     63 
     64  k = 0;
     65  while(rr > TOL){   
     66    // Ap = A*p
     67    for (n = 0 ; n < NUM * NUM ; n++){
     68      Ap[n] = 0.0;
     69      i = n % NUM; /* i is residual */
     70      j = n / NUM; /* j is division */
     71      if (i == 0) continue; /*skip left end*/
     72      if (i == NUM-1) continue; /*skip right end*/
     73      if (j == 0) continue; /*skip bottom end*/
     74      if (j == NUM-1) continue; /*skip top end*/
     75 
     76      /* computing A*p */
     77      Ap[n] = -4.0 * p[n];
     78      Ap[n] += p[n + 1];
     79      Ap[n] += p[n + NUM];
     80      Ap[n] += p[n - 1];
     81      Ap[n] += p[n - NUM];
     82    }   
     83 
     84    // alpha = r.r / p.Ap
     85    pAp = 0.0;
     86    for (n = 0 ; n < NUM * NUM ; n++){
     87      pAp += p[n] * Ap[n];
     88    }
     89    alpha = rr / pAp;
     90 
     91    //Beta
     92    rr1 = 0.0;
     93    for (n = 0 ; n < NUM * NUM ; n++){
     94      cnc[n] += alpha * p[n];
     95      r[n] -= alpha * Ap[n];
     96      rr1 += r[n] * r[n];
     97    } 
     98 
     99    beta = rr1 / rr;
     100 
     101    for (n = 0 ; n < NUM * NUM ; n++){
     102      p[n] = r[n] + beta * p[n];
     103    }
     104 
     105    rr = rr1;
     106    k++;
     107  }
     108 
     109  printf("# of Iteration=%d\n",k);
     110 
     111  /* Output Result */
     112  fp = fopen("res.dat","w");
     113  if (fp == NULL){
     114    printf("File could not create\n");
     115    exit(-1);
     116  }
     117 
     118  for (i = 0 ;i < NUM ; i++){
     119    for (j = 0 ; j < NUM ; j++){
     120      x = -1.0 + 2.0 * (double)i / (double)(NUM - 1);
     121      y = -1.0 + 2.0 * (double)j / (double)(NUM - 1);
     122      fprintf(fp,"%e %e %e\n",x,y,cnc[i + NUM * j]);
     123    }
     124  }
     125  fclose(fp);
     126 
     127  printf("Done\n");
     128}
     129}}}