Ex24.c

/*
  Large Scale Computing
  2D Convective Heat/Mass Transfer
  Ex24.c : Numerical Analysis
  Solve for
  duphi/dx dvphi/dy = d2phi/dx2 + d2phi/dy2
  the boundary conditions:
  phi = 0 along y=0, phi = 1 along y = 1
  dphi/dx=0 along y=+-0.5
 
 */
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include"bcg.h"
 
/* min and max macros */
#define max(a,b)        (((a) > (b)) ? (a) : (b))
#define min(a,b)        (((a) < (b)) ? (a) : (b))
 
#define NUM 100  /* number of points for x and y*/
#define TOL 1.0e-14  /* convergence tolerance */
 
int main(int argc, char **argv){
  int i,j,k,n;
  int iw,ie,is,in;
  double d,x,y,u0,u,v;
  double dl;
  double ae,aw,an,as,ap;
  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;
 
  if (argc < 3){
    printf("Usage:%s [D] [U0]\n",argv[0]);
    exit(-1);
  }
 
  /* set parameters */
  d = atof(argv[1]);
  u0 = atof(argv[2]);
  printf("D=%e U0=%e\n",d,u0);
 
  /*assuming dx = dy : domain is 1x1 */
  dl = 1.0 / (double)(NUM - 1);
 
  /* initial & boundary conditions */
  for (n = 0 ; n < NUM * NUM ; n++){
    phi[n] = 0.0;
  }
  for (i = 0 ; i < NUM ; i++){
    n = i + NUM * (NUM - 1);
    phi[n] = 1.0; /* fixed phi=1 */
  }
 
  /* Setup Matrix A and RHS */
  nnz = 0;
  for (n = 0 ; n < NUM * NUM ; n++){
    /* compute i,j indices */
    i = n % NUM; /* i is residual */
    j = n / NUM; /* j is division */
    ie = n + 1;
    iw = n - 1;
    in = n + NUM;
    is = n - NUM;
 
    /* (x,y) at p-pont */
    x = -0.5 + (double)i / (double)(NUM - 1);
    y = (double)j / (double)(NUM - 1);
 
    /* set general RHS */
    rhs[n] = 0.0;
 
    /* set first nonzero column of row index(i,j) */
    rowptr[n] = nnz;
 
    /* general coefficients */
    /* (u,v) = u0(x,-y) */
    u = u0 * (x + dl * 0.5); /* u at e-point*/
    ae = d / dl + max(-u,0.0);
    u = u0 * (x - dl * 0.5); /* u at w-point*/
    aw = d / dl + max(u,0.0);
    v = -u0 * (y + dl * 0.5); /* v at n-point*/
    an = d / dl + max(-v,0.0);
    v = -u0 * (y - dl * 0.5); /* v at s-point*/
    as = d / dl + max(v,0.0);
    ap = ae + aw + an + as;
 
    /* sounth */
    if (is >= 0){
      nzval[nnz] = as;
      colind[nnz] = is;
      nnz++;
    }     
 
    /* west */
    if (iw >= 0){
      nzval[nnz] = aw;
      colind[nnz] = iw;
      nnz++;
    } 
 
    /* diagonal Element */
    nzval[nnz] = -ap;
    colind[nnz] = n;
    nnz++;
 
    /* east */
    if (ie < NUM * NUM){
      nzval[nnz] = ae;
      colind[nnz] = ie;
      nnz++;
    } 
 
    /* north */
    if (in < NUM * NUM){
      nzval[nnz] = an;
      colind[nnz] = in;
      nnz++;
    } 
  }
  rowptr[NUM * NUM] = nnz;
 
 
  /* Setting Boundary Conditions */
  for (i = 0 ; i < NUM ; i++){
    /* Bottom fixed phi=0 */
    n = i;
    for (k = rowptr[n] ; k < rowptr[n+1] ; k++){
      if (n == colind[k]){
        nzval[k] = -1.0; /* diag=-1 */
        rhs[n] = nzval[k] * phi[n];
      }else{
        nzval[k] = 0.0; /* off-diag=0 */
      }
    }
    /* Top fixed phi=1 */
    n = i + NUM * (NUM - 1);
    for (k = rowptr[n] ; k < rowptr[n+1] ; k++){
      if (n == colind[k]){
        nzval[k] = -1.0; /* diag=-1 */
        rhs[n] = nzval[k] * phi[n];
      }else{
        nzval[k] = 0.0; /* off-diag=0 */
      }
    }
  }
  for (j = 0 ; j < NUM ; j++){
    /* Left  dphi/dx=0 */
    n = j * NUM;
    ie = n + 1;
    rhs[n] = 0.0;
    for (k = rowptr[n] ; k < rowptr[n+1] ; k++){
      if (n == colind[k]){
        nzval[k] = -1.0; /* diag=-1 */
      }else{
        if (ie == colind[k]){
          nzval[k] = 1.0; /* off-diag_e=1 */
        }else{
          nzval[k] = 0.0; /* off-diag_else=0 */
        }
      }
    }
    /* Right dphi/dx=0 */
    n = NUM - 1 + j * NUM;
    iw = n - 1;
    rhs[n] = 0.0;
    for (k = rowptr[n] ; k < rowptr[n+1] ; k++){
      if (n == colind[k]){
        nzval[k] = -1.0; /* diag=-1 */
      }else{
        if (iw == colind[k]){
          nzval[k] = 1.0; /* off-diag_w=1 */
        }else{
          nzval[k] = 0.0; /* off-diag_else=0 */
        }
      }
    }
  }
 
  /* solve with BiCG */
  k = solve_bcg(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 = -0.5 + (double)i / (double)(NUM - 1);
      y = (double)j / (double)(NUM - 1);
      fprintf(fp,"%e %e %e\n",x,y,phi[i + NUM * j]);
    }
  }
  fclose(fp);
 
  printf("Done\n");
}

bcg.c