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


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

migrated content from ccs wiki

Legend:

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

    v1 v1  
     1= Ex22.c =
     2{{{#!C
     3/*
     4  Large Scale Computing
     5  2D Heat/Mass Transfer
     6  Ex22.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 */
     13#include<stdio.h>
     14#include<stdlib.h>
     15#include<math.h>
     16#include"cg.h"
     17 
     18#define NUM 100  /* number of points for x and y*/
     19#define TOL 1.0e-12  /* convergence tolerance */
     20#define FAC 1.0e+5 /* Large factor */
     21 
     22int main(void ){
     23  int i,j,k,n;
     24  int iw,ie,is,in;
     25  double x,y,dx;
     26  double phi[NUM * NUM]; /* define 1D array of length NUM * NUM */
     27  double rhs[NUM * NUM]; /* define 1D array of length NUM * NUM */
     28  double nzval[5 * NUM * NUM]; /* non-zero elements of A */
     29  int colind[5 * NUM * NUM]; /* column indices of nonzeros */
     30  int rowptr[NUM * NUM + 1]; /**/
     31  int nnz;
     32  FILE *fp;
     33 
     34  /*assuming dx = dy : domain is 2x2 */
     35  dx = 2.0 / (double)(NUM - 1);
     36 
     37  /* initial & boundary conditions */
     38  for (n = 0 ; n < NUM * NUM ; n++){
     39    phi[n] = 0.0;
     40  }
     41 
     42  /* Setup Matrix A and RHS */
     43  nnz = 0;
     44  for (n = 0 ; n < NUM * NUM ; n++){
     45    /* set general RHS */
     46    rhs[n] = -dx * dx;
     47 
     48    /* set first nonzero column of row index(i,j) */
     49    rowptr[n] = nnz;
     50 
     51    ie = n + 1;
     52    iw = n - 1;
     53    in = n + NUM;
     54    is = n - NUM;
     55 
     56    /* sounth */
     57    if (is >= 0){
     58      nzval[nnz] = 1.0;
     59      colind[nnz] = is;
     60      nnz++;
     61    }     
     62 
     63    /* west */
     64    if (iw >= 0){
     65      nzval[nnz] = 1.0;
     66      colind[nnz] = iw;
     67      nnz++;
     68    }
     69 
     70    /* diagonal Element */
     71    nzval[nnz] = -4.0;
     72    colind[nnz] = n;
     73    nnz++;
     74 
     75    /* east */
     76    if (ie < NUM * NUM){
     77      nzval[nnz] = 1.0;
     78      colind[nnz] = ie;
     79      nnz++;
     80    }
     81 
     82    /* north */
     83    if (in < NUM * NUM){
     84      nzval[nnz] = 1.0;
     85      colind[nnz] = in;
     86      nnz++;
     87    }
     88  }
     89  rowptr[NUM * NUM] = nnz;
     90 
     91 
     92  /* Setting Boundary Conditions */
     93  /* Add Big Number to Both Sides */
     94  for (i = 0 ; i < NUM ; i++){
     95    /* Bottom */
     96    n = i;
     97    for (k = rowptr[n] ; k < rowptr[n+1] ; k++){
     98      if (n == colind[k]){
     99        nzval[k] += FAC;
     100        rhs[n] += FAC * phi[n];
     101      }
     102    }
     103    /* Top */
     104    n = i + NUM * (NUM - 1);
     105    for (k = rowptr[n] ; k < rowptr[n+1] ; k++){
     106      if (n == colind[k]){
     107        nzval[k] += FAC;
     108        rhs[n] += FAC * phi[n];
     109      }
     110    }
     111  }
     112  for (j = 0 ; j < NUM ; j++){
     113    /* Left */
     114    n = j * NUM;
     115    for (k = rowptr[n] ; k < rowptr[n+1] ; k++){
     116      if (n == colind[k]){
     117        nzval[k] += FAC;
     118        rhs[n] += FAC * phi[n];
     119      }
     120    }
     121    /* Right */
     122    n = NUM - 1 + j * NUM;
     123    for (k = rowptr[n] ; k < rowptr[n+1] ; k++){
     124      if (n == colind[k]){
     125        nzval[k] += FAC;
     126        rhs[n] += FAC * phi[n];
     127      }
     128    }
     129  }
     130 
     131 
     132  /* solve with CG */
     133  k = solve_cg(NUM * NUM, nnz, nzval, colind, rowptr, phi, rhs, TOL);
     134  if (k < 0){
     135    printf("calculation failed\n");
     136    exit(-1);
     137  }
     138  printf("Number of Interation=%d\n",k);
     139 
     140 
     141  /* Output Result */
     142  fp = fopen("res.dat","w");
     143  if (fp == NULL){
     144    printf("File could not create\n");
     145    exit(-1);
     146  }
     147 
     148  for (i = 0 ;i < NUM ; i++){
     149    for (j = 0 ; j < NUM ; j++){
     150      x = -1.0 + 2.0 * (double)i / (double)(NUM - 1);
     151      y = -1.0 + 2.0 * (double)j / (double)(NUM - 1);
     152      fprintf(fp,"%e %e %e\n",x,y,phi[i + NUM * j]);
     153    }
     154  }
     155  fclose(fp);
     156 
     157  printf("Done\n");
     158}
     159}}}