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

Timestamp:

May 14, 2015 3:21:20 PM (7 years ago)

Author:

cmaggio

Comment:

migrated content from ccs wiki

cypress/Programming/HeatMassTransfer/ex19.c

@@ +1 @@
+= Ex19.c =
+{{{#!C
+/*
+  Large Scale Computing
+  2D Heat/Mass Transfer
+  ex19.c : Numerical Analysis
+  Solve for
+  d2c/dx2 + d2c/dy2 + 1 = 0
+  With the boundary conditions of c = 0 
+  along lines of x=1,-1, and y=1, and -1
+ 
+  Conjugate Gradient Method
+ */
+#include<stdio.h>
+#include<stdlib.h>
+#include<math.h>
+ 
+#define NUM 100  /* number of points for x and y*/
+#define TOL 1.0e-12  /* convergence tolerance */
+ 
+int main(void ){
+  int i,j,n;
+  double x,y,dx;
+  double cnc[NUM * NUM]; /* define 1D array of length NUM * NUM */
+  double r[NUM * NUM];
+  double p[NUM * NUM];
+  double Ap[NUM * NUM];
+  double rr,pAp,rr1,alpha,beta;
+  int k;
+  FILE *fp;
+ 
+  /*assuming dx = dy : domain is 2x2 */
+  dx = 2.0 / (double)(NUM - 1);
+ 
+  /* Initialize */
+  for (n = 0 ; n < NUM * NUM ; n++){
+    cnc[n] = 0.0;
+    r[n] = 0.0;
+  }
+ 
+  /* computing for C with Conjugate Gradient Method */
+  /* r0 = b - Ax0 */
+  for (n = 0 ; n < NUM * NUM ; n++){
+    i = n % NUM; /* i is residual */
+    j = n / NUM; /* j is division */
+    if (i == 0) continue; /*skip left end*/
+    if (i == NUM-1) continue; /*skip right end*/
+    if (j == 0) continue; /*skip bottom end*/
+    if (j == NUM-1) continue; /*skip top end*/
+    r[n] = -dx * dx;
+    r[n] -= -4.0 * cnc[n];
+    r[n] -= cnc[n + 1];
+    r[n] -= cnc[n + NUM];
+    r[n] -= cnc[n - 1];
+    r[n] -= cnc[n - NUM];
+  }
+ 
+  rr = 0.0;
+  for (n = 0 ; n < NUM * NUM ; n++){
+    p[n] = r[n];    /* p = r */
+    rr += r[n] * r[n];  /* rr = r.r */
+  }
+ 
+  k = 0;
+  while(rr > TOL){    
+    // Ap = A*p
+    for (n = 0 ; n < NUM * NUM ; n++){
+      Ap[n] = 0.0;
+      i = n % NUM; /* i is residual */
+      j = n / NUM; /* j is division */
+      if (i == 0) continue; /*skip left end*/
+      if (i == NUM-1) continue; /*skip right end*/
+      if (j == 0) continue; /*skip bottom end*/
+      if (j == NUM-1) continue; /*skip top end*/
+ 
+      /* computing A*p */
+      Ap[n] = -4.0 * p[n];
+      Ap[n] += p[n + 1];
+      Ap[n] += p[n + NUM];
+      Ap[n] += p[n - 1];
+      Ap[n] += p[n - NUM];
+    }    
+ 
+    // alpha = r.r / p.Ap
+    pAp = 0.0;
+    for (n = 0 ; n < NUM * NUM ; n++){
+      pAp += p[n] * Ap[n];
+    }
+    alpha = rr / pAp;
+ 
+    //Beta
+    rr1 = 0.0;
+    for (n = 0 ; n < NUM * NUM ; n++){
+      cnc[n] += alpha * p[n];
+      r[n] -= alpha * Ap[n];
+      rr1 += r[n] * r[n];
+    }  
+ 
+    beta = rr1 / rr;
+ 
+    for (n = 0 ; n < NUM * NUM ; n++){
+      p[n] = r[n] + beta * p[n];
+    }
+ 
+    rr = rr1;
+    k++;
+  }
+ 
+  printf("# of Iteration=%d\n",k);
+ 
+  /* Output Result */
+  fp = fopen("res.dat","w");
+  if (fp == NULL){
+    printf("File could not create\n");
+    exit(-1);
+  }
+ 
+  for (i = 0 ;i < NUM ; i++){
+    for (j = 0 ; j < NUM ; j++){
+      x = -1.0 + 2.0 * (double)i / (double)(NUM - 1);
+      y = -1.0 + 2.0 * (double)j / (double)(NUM - 1);
+      fprintf(fp,"%e %e %e\n",x,y,cnc[i + NUM * j]);
+    }
+  }
+  fclose(fp);
+ 
+  printf("Done\n");
+}
+}}}