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