Version 1 (modified by 9 years ago) ( diff ) | ,
---|
Ex25.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"); }
Note:
See TracWiki
for help on using the wiki.