Ex22.c

/*
  Large Scale Computing
  2D Heat/Mass Transfer
  Ex22.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 */
#define FAC 1.0e+5 /* Large factor */
 
int main(void ){
  int i,j,k,n;
  int iw,ie,is,in;
  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 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 (n = 0 ; n < NUM * NUM ; n++){
    phi[n] = 0.0;
  }
 
  /* Setup Matrix A and RHS */
  nnz = 0;
  for (n = 0 ; n < NUM * NUM ; n++){
    /* set general RHS */
    rhs[n] = -dx * dx;
 
    /* set first nonzero column of row index(i,j) */
    rowptr[n] = nnz;
 
    ie = n + 1;
    iw = n - 1;
    in = n + NUM;
    is = n - NUM;
 
    /* sounth */
    if (is >= 0){
      nzval[nnz] = 1.0;
      colind[nnz] = is;
      nnz++;
    }     
 
    /* west */
    if (iw >= 0){
      nzval[nnz] = 1.0;
      colind[nnz] = iw;
      nnz++;
    } 
 
    /* diagonal Element */
    nzval[nnz] = -4.0;
    colind[nnz] = n;
    nnz++;
 
    /* east */
    if (ie < NUM * NUM){
      nzval[nnz] = 1.0;
      colind[nnz] = ie;
      nnz++;
    } 
 
    /* north */
    if (in < NUM * NUM){
      nzval[nnz] = 1.0;
      colind[nnz] = in;
      nnz++;
    } 
  }
  rowptr[NUM * NUM] = nnz;
 
 
  /* Setting Boundary Conditions */
  /* Add Big Number to Both Sides */
  for (i = 0 ; i < NUM ; i++){
    /* Bottom */
    n = i;
    for (k = rowptr[n] ; k < rowptr[n+1] ; k++){
      if (n == colind[k]){
        nzval[k] += FAC;
        rhs[n] += FAC * phi[n];
      }
    }
    /* Top */
    n = i + NUM * (NUM - 1);
    for (k = rowptr[n] ; k < rowptr[n+1] ; k++){
      if (n == colind[k]){
        nzval[k] += FAC;
        rhs[n] += FAC * phi[n];
      }
    }
  }
  for (j = 0 ; j < NUM ; j++){
    /* Left */
    n = j * NUM;
    for (k = rowptr[n] ; k < rowptr[n+1] ; k++){
      if (n == colind[k]){
        nzval[k] += FAC;
        rhs[n] += FAC * phi[n];
      }
    }
    /* Right */
    n = NUM - 1 + j * NUM;
    for (k = rowptr[n] ; k < rowptr[n+1] ; k++){
      if (n == colind[k]){
        nzval[k] += FAC;
        rhs[n] += FAC * phi[n];
      }
    }
  }
 
 
  /* solve with CG */
  k = solve_cg(NUM * NUM, nnz, nzval, colind, rowptr, phi, rhs, TOL);
  if (k < 0){
    printf("calculation failed\n");
    exit(-1);
  }
  printf("Number of Interation=%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,phi[i + NUM * j]);
    }
  }
  fclose(fp);
 
  printf("Done\n");
}