= Ex22.c = {{{#!C /* Large Scale Computing 2D Heat/Mass Transfer Ex22.c : Numerical Analysis Solve for d2c/dx2 + d2c/dy2 + 1 = 0 With the boundary conditions of c = 0 along lines of x=1,-1, and y=1, and -1 */ #include #include #include #include"cg.h" #define NUM 100 /* number of points for x and y*/ #define TOL 1.0e-12 /* convergence tolerance */ #define FAC 1.0e+5 /* Large factor */ int main(void ){ int i,j,k,n; int iw,ie,is,in; double x,y,dx; 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; /*assuming dx = dy : domain is 2x2 */ dx = 2.0 / (double)(NUM - 1); /* initial & boundary conditions */ for (n = 0 ; n < NUM * NUM ; n++){ phi[n] = 0.0; } /* Setup Matrix A and RHS */ nnz = 0; for (n = 0 ; n < NUM * NUM ; n++){ /* set general RHS */ rhs[n] = -dx * dx; /* set first nonzero column of row index(i,j) */ rowptr[n] = nnz; ie = n + 1; iw = n - 1; in = n + NUM; is = n - NUM; /* sounth */ if (is >= 0){ nzval[nnz] = 1.0; colind[nnz] = is; nnz++; } /* west */ if (iw >= 0){ nzval[nnz] = 1.0; colind[nnz] = iw; nnz++; } /* diagonal Element */ nzval[nnz] = -4.0; colind[nnz] = n; nnz++; /* east */ if (ie < NUM * NUM){ nzval[nnz] = 1.0; colind[nnz] = ie; nnz++; } /* north */ if (in < NUM * NUM){ nzval[nnz] = 1.0; colind[nnz] = in; nnz++; } } rowptr[NUM * NUM] = nnz; /* Setting Boundary Conditions */ /* Add Big Number to Both Sides */ for (i = 0 ; i < NUM ; i++){ /* Bottom */ n = i; for (k = rowptr[n] ; k < rowptr[n+1] ; k++){ if (n == colind[k]){ nzval[k] += FAC; rhs[n] += FAC * phi[n]; } } /* Top */ n = i + NUM * (NUM - 1); for (k = rowptr[n] ; k < rowptr[n+1] ; k++){ if (n == colind[k]){ nzval[k] += FAC; rhs[n] += FAC * phi[n]; } } } for (j = 0 ; j < NUM ; j++){ /* Left */ n = j * NUM; for (k = rowptr[n] ; k < rowptr[n+1] ; k++){ if (n == colind[k]){ nzval[k] += FAC; rhs[n] += FAC * phi[n]; } } /* Right */ n = NUM - 1 + j * NUM; for (k = rowptr[n] ; k < rowptr[n+1] ; k++){ if (n == colind[k]){ nzval[k] += FAC; rhs[n] += FAC * phi[n]; } } } /* solve with CG */ k = solve_cg(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 = -1.0 + 2.0 * (double)i / (double)(NUM - 1); y = -1.0 + 2.0 * (double)j / (double)(NUM - 1); fprintf(fp,"%e %e %e\n",x,y,phi[i + NUM * j]); } } fclose(fp); printf("Done\n"); } }}}