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

Timestamp:

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

Author:

cmaggio

Comment:

migrated content from ccs wiki

cypress/Programming/HeatMassTransfer/ex20.c

@@ +1 @@
+= Ex20.c =
+{{{#!C
+/*
+  Large Scale Computing
+  2D Heat/Mass Transfer
+  Ex20.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
+ 
+ */
+#include<stdio.h>
+#include<stdlib.h>
+#include<math.h>
+#include"cg.h"
+ 
+#define NUM 100  /* number of points for x and y*/
+#define TOL 1.0e-12  /* convergence tolerance */
+ 
+int index(int ii,int jj);
+int main(void ){
+  int i,j,k,n;
+  double x,y,dx;
+  double phi[NUM][NUM]; /* define 1D array of length NUM * NUM */
+  double rhs[NUM * NUM]; /* define 1D array of length NUM * NUM */
+  double sol[NUM * NUM]; /* define 1D array of length NUM * NUM */
+  double nzval[5 * NUM * NUM]; /* non-zero elements of A */
+  int colind[5 * NUM * NUM]; /* column indices of nonzeros */
+  int rowptr[NUM * NUM + 1]; /**/
+  int nnz;
+  FILE *fp;
+ 
+  /*assuming dx = dy : domain is 2x2 */
+  dx = 2.0 / (double)(NUM - 1);
+ 
+  /* initial & boundary conditions */
+  for (j = 0 ; j < NUM  ; j++){
+    for (i = 0 ; i < NUM ; i++){
+      phi[i][j] = 0.0;
+    }
+  }
+ 
+  /* Setup Matrix A and RHS */
+  nnz = 0;
+  for (j = 1 ; j < NUM - 1 ; j++){
+    for (i = 1 ; i < NUM - 1 ; i++){
+      /* initial guess */
+      sol[index(i,j)] = phi[i][j];
+ 
+      /* set general RHS */
+      rhs[index(i,j)] = -dx * dx;
+ 
+      /* set first nonzero column of row index(i,j) */
+      rowptr[index(i,j)] = nnz;
+ 
+      /* sounth */
+      if (j == 1){
+        rhs[index(i,j)] -= phi[i][0]; 
+      }else{
+        nzval[nnz] = 1.0;
+        colind[nnz] = index(i,j-1);
+        nnz++;
+      } 
+ 
+      /* west */
+      if (i == 1){
+        rhs[index(i,j)] -= phi[0][j]; 
+      }else{
+        nzval[nnz] = 1.0;
+        colind[nnz] = index(i-1,j);
+        nnz++;
+      } 
+ 
+      /* diagonal Element */
+      nzval[nnz] = -4.0;
+      colind[nnz] = index(i,j);
+      nnz++;
+ 
+      /* east */
+      if (i == NUM - 2){
+        rhs[index(i,j)] -= phi[NUM-1][j]; 
+      }else{
+        nzval[nnz] = 1.0;
+        colind[nnz] = index(i+1,j);
+        nnz++;
+      } 
+ 
+      /* north */
+      if (j == NUM - 2){
+        rhs[index(i,j)] -= phi[i][NUM - 1]; 
+      }else{
+        nzval[nnz] = 1.0;
+        colind[nnz] = index(i,j+1);
+        nnz++;
+      } 
+    }
+  }
+ 
+  /* last element of rowptr si nnz */
+  rowptr[(NUM-2) * (NUM-2)] = nnz;
+ 
+  /* solve with CG */
+  k = solve_cg((NUM-2) * (NUM-2), nnz, nzval, colind, rowptr, sol, rhs, TOL);
+  if (k < 0){
+    printf("calculation failed\n");
+    exit(-1);
+  }
+  printf("Number of Interation=%d\n",k);
+ 
+  /* store solution to main aray */
+  for (j = 1 ; j < NUM - 1 ; j++){
+    for (i = 1 ; i < NUM - 1 ; i++){
+      phi[i][j] = sol[index(i,j)];
+    }
+  }
+ 
+  /* 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,phi[i][j]);
+    }
+  }
+  fclose(fp);
+ 
+  printf("Done\n");
+}
+ 
+/* index for interior Points */
+int index(int ii,int jj){
+  return ii - 1 + (NUM - 2) * (jj - 1);
+}
+}}}
+
+[[cypress/Programming/Cg|cg.c]]