Ex22.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<stdio.h>
#include<stdlib.h>
#include<math.h>
#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");
}