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

Timestamp:

05/14/15 15:27:28 (10 years ago)

Author:

cmaggio

Comment:

migrated content from ccs wiki

cypress/Programming/HeatMassTransfer/ex25.c

@@ +1 @@
+= Ex25.c = 
+{{{#!C
+/*
+  Large Scale Computing
+  2D Convective Heat/Mass Transfer
+  ex25.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"slu.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*/
+ 
+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 SuperLU */
+  k = solve_slu(NUM * NUM, nnz, nzval, colind, rowptr, phi, rhs);
+  if (k != 0){
+    printf("calculation failed\n");
+    exit(-1);
+  }
+ 
+  /* 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");
+}
+}}}