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


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

Legend:

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

    v1 v1  
     1= Ex24.c =
     2{{{#!C
     3/*
     4  Large Scale Computing
     5  2D Convective Heat/Mass Transfer
     6  Ex24.c : Numerical Analysis
     7  Solve for
     8  duphi/dx dvphi/dy = d2phi/dx2 + d2phi/dy2
     9  the boundary conditions:
     10  phi = 0 along y=0, phi = 1 along y = 1
     11  dphi/dx=0 along y=+-0.5
     12 
     13 */
     14#include<stdio.h>
     15#include<stdlib.h>
     16#include<math.h>
     17#include"bcg.h"
     18 
     19/* min and max macros */
     20#define max(a,b)        (((a) > (b)) ? (a) : (b))
     21#define min(a,b)        (((a) < (b)) ? (a) : (b))
     22 
     23#define NUM 100  /* number of points for x and y*/
     24#define TOL 1.0e-14  /* convergence tolerance */
     25 
     26int main(int argc, char **argv){
     27  int i,j,k,n;
     28  int iw,ie,is,in;
     29  double d,x,y,u0,u,v;
     30  double dl;
     31  double ae,aw,an,as,ap;
     32  double phi[NUM * NUM]; /* define 1D array of length NUM * NUM */
     33  double rhs[NUM * NUM]; /* define 1D array of length NUM * NUM */
     34  double nzval[5 * NUM * NUM]; /* non-zero elements of A */
     35  int colind[5 * NUM * NUM]; /* column indices of nonzeros */
     36  int rowptr[NUM * NUM + 1]; /**/
     37  int nnz;
     38  FILE *fp;
     39 
     40  if (argc < 3){
     41    printf("Usage:%s [D] [U0]\n",argv[0]);
     42    exit(-1);
     43  }
     44 
     45  /* set parameters */
     46  d = atof(argv[1]);
     47  u0 = atof(argv[2]);
     48  printf("D=%e U0=%e\n",d,u0);
     49 
     50  /*assuming dx = dy : domain is 1x1 */
     51  dl = 1.0 / (double)(NUM - 1);
     52 
     53  /* initial & boundary conditions */
     54  for (n = 0 ; n < NUM * NUM ; n++){
     55    phi[n] = 0.0;
     56  }
     57  for (i = 0 ; i < NUM ; i++){
     58    n = i + NUM * (NUM - 1);
     59    phi[n] = 1.0; /* fixed phi=1 */
     60  }
     61 
     62  /* Setup Matrix A and RHS */
     63  nnz = 0;
     64  for (n = 0 ; n < NUM * NUM ; n++){
     65    /* compute i,j indices */
     66    i = n % NUM; /* i is residual */
     67    j = n / NUM; /* j is division */
     68    ie = n + 1;
     69    iw = n - 1;
     70    in = n + NUM;
     71    is = n - NUM;
     72 
     73    /* (x,y) at p-pont */
     74    x = -0.5 + (double)i / (double)(NUM - 1);
     75    y = (double)j / (double)(NUM - 1);
     76 
     77    /* set general RHS */
     78    rhs[n] = 0.0;
     79 
     80    /* set first nonzero column of row index(i,j) */
     81    rowptr[n] = nnz;
     82 
     83    /* general coefficients */
     84    /* (u,v) = u0(x,-y) */
     85    u = u0 * (x + dl * 0.5); /* u at e-point*/
     86    ae = d / dl + max(-u,0.0);
     87    u = u0 * (x - dl * 0.5); /* u at w-point*/
     88    aw = d / dl + max(u,0.0);
     89    v = -u0 * (y + dl * 0.5); /* v at n-point*/
     90    an = d / dl + max(-v,0.0);
     91    v = -u0 * (y - dl * 0.5); /* v at s-point*/
     92    as = d / dl + max(v,0.0);
     93    ap = ae + aw + an + as;
     94 
     95    /* sounth */
     96    if (is >= 0){
     97      nzval[nnz] = as;
     98      colind[nnz] = is;
     99      nnz++;
     100    }     
     101 
     102    /* west */
     103    if (iw >= 0){
     104      nzval[nnz] = aw;
     105      colind[nnz] = iw;
     106      nnz++;
     107    }
     108 
     109    /* diagonal Element */
     110    nzval[nnz] = -ap;
     111    colind[nnz] = n;
     112    nnz++;
     113 
     114    /* east */
     115    if (ie < NUM * NUM){
     116      nzval[nnz] = ae;
     117      colind[nnz] = ie;
     118      nnz++;
     119    }
     120 
     121    /* north */
     122    if (in < NUM * NUM){
     123      nzval[nnz] = an;
     124      colind[nnz] = in;
     125      nnz++;
     126    }
     127  }
     128  rowptr[NUM * NUM] = nnz;
     129 
     130 
     131  /* Setting Boundary Conditions */
     132  for (i = 0 ; i < NUM ; i++){
     133    /* Bottom fixed phi=0 */
     134    n = i;
     135    for (k = rowptr[n] ; k < rowptr[n+1] ; k++){
     136      if (n == colind[k]){
     137        nzval[k] = -1.0; /* diag=-1 */
     138        rhs[n] = nzval[k] * phi[n];
     139      }else{
     140        nzval[k] = 0.0; /* off-diag=0 */
     141      }
     142    }
     143    /* Top fixed phi=1 */
     144    n = i + NUM * (NUM - 1);
     145    for (k = rowptr[n] ; k < rowptr[n+1] ; k++){
     146      if (n == colind[k]){
     147        nzval[k] = -1.0; /* diag=-1 */
     148        rhs[n] = nzval[k] * phi[n];
     149      }else{
     150        nzval[k] = 0.0; /* off-diag=0 */
     151      }
     152    }
     153  }
     154  for (j = 0 ; j < NUM ; j++){
     155    /* Left  dphi/dx=0 */
     156    n = j * NUM;
     157    ie = n + 1;
     158    rhs[n] = 0.0;
     159    for (k = rowptr[n] ; k < rowptr[n+1] ; k++){
     160      if (n == colind[k]){
     161        nzval[k] = -1.0; /* diag=-1 */
     162      }else{
     163        if (ie == colind[k]){
     164          nzval[k] = 1.0; /* off-diag_e=1 */
     165        }else{
     166          nzval[k] = 0.0; /* off-diag_else=0 */
     167        }
     168      }
     169    }
     170    /* Right dphi/dx=0 */
     171    n = NUM - 1 + j * NUM;
     172    iw = n - 1;
     173    rhs[n] = 0.0;
     174    for (k = rowptr[n] ; k < rowptr[n+1] ; k++){
     175      if (n == colind[k]){
     176        nzval[k] = -1.0; /* diag=-1 */
     177      }else{
     178        if (iw == colind[k]){
     179          nzval[k] = 1.0; /* off-diag_w=1 */
     180        }else{
     181          nzval[k] = 0.0; /* off-diag_else=0 */
     182        }
     183      }
     184    }
     185  }
     186 
     187  /* solve with BiCG */
     188  k = solve_bcg(NUM * NUM, nnz, nzval, colind, rowptr, phi, rhs, TOL);
     189  if (k < 0){
     190    printf("calculation failed\n");
     191    exit(-1);
     192  }
     193  printf("Number of Interation=%d\n",k);
     194 
     195 
     196  /* Output Result */
     197  fp = fopen("res.dat","w");
     198  if (fp == NULL){
     199    printf("File could not create\n");
     200    exit(-1);
     201  }
     202 
     203  for (i = 0 ;i < NUM ; i++){
     204    for (j = 0 ; j < NUM ; j++){
     205      x = -0.5 + (double)i / (double)(NUM - 1);
     206      y = (double)j / (double)(NUM - 1);
     207      fprintf(fp,"%e %e %e\n",x,y,phi[i + NUM * j]);
     208    }
     209  }
     210  fclose(fp);
     211 
     212  printf("Done\n");
     213}
     214}}}
     215[[cypress/Programming/BiConjugateGradient|bcg.c]]