fix main file

3 files changed, 166 insertions, 175 deletions

diff --git a/main.c b/main.c
index fb9548a..98708c8 100644
--- a/main.c
+++ b/main.c
@@ -38,10 +38,18 @@ int main (void)
 	double kf1 = 1.0;
 	double kf2 = 1.0;
 
-	// initialize permeability matrix
+	int ids[5] = { - nx, -1, 0, 1, nx }; // system matrix
+	int pre_ids[1] = { 0 };              // preconditioner
 
+	// allocate memory
 	double* kf = (double*) calloc(dim, sizeof(double));
+	double* a  = (double*) calloc(5 * dim, sizeof(double));
+	double* c  = (double*) calloc(dim, sizeof(double));
+	double* r  = (double*) calloc(dim, sizeof(double));
+	double* ar = (double*) calloc(dim, sizeof(double));
+	double* x  = (double*) calloc(dim, sizeof(double));
 
+	// initialize permeability matrix
 	for (int i = 0; i < ny; i ++) {
 		for (int j = 0; j < nx; j ++) {
 			double xij = xL + j * dx + 0.5 * dx;
@@ -56,11 +64,6 @@ int main (void)
 
 	// initialize system matrix
 	//
-	// store diagonal ids for each non-zero diagonal
-	int ids[5] = { - nx, -1, 0, 1, nx };
-	// giving a matrix that with shape (5, dim)
-	double* a = (double*) calloc(5 * dim, sizeof(double));
-
 	struct diag m;
 	m.values = a;
 	m.ids    = ids;
@@ -70,74 +73,70 @@ int main (void)
 	m.len_id = 5;
 
 	system_matrix(&m, nx, ny, kf);
-	printf("... \n");
+
+	// initialize boundary conditions
+	//
+
+	// neumann boundary
+	r[7] = -10.0;
+
+	// dirichlet boundary
+	for (int i = 0; i < nx; i ++) {
+		x[i] = 40.0;
+	}
+
+	multiply(&m, x, ar);
+
 	for (int i = 0; i < dim; i ++) {
+		r[i] = r[i] - ar[i];
+		if (i < nx) {
+			r[i] = x[i];
+		}
+	}
+
+
+	for (int i = 0; i < nx; i ++) {
 		for (int j = 0; j < dim; j ++) {
-			printf(" %g ", get(&m, i, j));
+			if (i != j) {
+				put(&m, j, i, 0.0);
+				put(&m, i, j, 0.0);
+			} else {
+				put(&m, j, i, 1.0);
+			}
 		}
-		printf("\n");
 	}
 
-	// initialize boundary conditions
+	// solve the system A*x = b
 	//
 
-//	double* amod = (double*) calloc(5 * dim, sizeof(double));
-//	amod = memcpy(amod, m.values, 5 * dim * sizeof(double));
-//
-//	double* r = (double*) calloc(dim, sizeof(double));
-//
-//	// dirichlet boundary
-//	for (int i = 0; i < nx; i ++) {
-//		r[i] = 10.0 * kf[i];
-//		for (int j = 0; j < dim; j ++) {
-//			if (i != j) {
-//				put(m, j, i, 0.0);
-//			} else {
-//				put(m, j, i, kf[i]);
-//			}
-//		}
-//	}
-//
-// m.values = amod;
-//
-//	// neumann boundary
-//	for (int i = dim - nx; i < dim; i ++) {
-//		r[i] = -1.8e-5;
-//	}
-//
-//	// solve the system A*x = b
-//	//
-//
-//	double* x = (double*) calloc(dim, sizeof(double));
-//	for (int i = 0; i < dim; i ++) {
-//		x[i] = 1.0;
-//	}
-//
-
-
-//	// define preconditioner
-//	double* c = (double*) calloc(dim, sizeof(double));
-//	for (int i = 0; i < dim; i ++) {
-//		c[i] = 1.0 / get(m, i, i);
-//	}
-//	int pre_ids[1] = { 0 };
-//
-//	struct mg pre = { c, pre_ids, dim, dim, 1 };
-//
-//	int flag = pcg(m, r, x, pre, 1.0e-6);
-//	printf("\n\niteration finished w/ flag : %d\n", flag);
-//
-//	printf("\nsolution:\n");
-//	for (int i = 0; i < ny; i ++) {
-//		for (int j = 0; j < nx; j ++) {
-//			printf("%g; ", x[i * nx + j]);
-//		}
-//		printf("\n");
-//	}
-//
+	for (int i = 0; i < dim; i ++) {
+		x[i] = 1.0;
+	}
+
+	// define preconditioner
+	for (int i = 0; i < dim; i ++) {
+		c[i] = 1.0 / get(&m, i, i);
+	}
+
+	struct diag pre = { c, pre_ids, dim, dim, dim, 1 };
+
+	pcg(&m, r, x, &pre, 1.0e-12);
+
+	printf("\n>>> (app) :: solution:\n");
+	for (int i = 0; i < ny; i ++) {
+		printf(">>>");
+		for (int j = 0; j < nx; j ++) {
+			printf(" %6g ", x[i * nx + j]);
+		}
+		printf("\n");
+	}
+
 	// clean up
 	free(kf);
 	free(a);
+	free(c);
+	free(r);
+	free(x);
 
 	return 0;
 
diff --git a/solver.c b/solver.c
index 4f5d683..3e21b75 100644
--- a/solver.c
+++ b/solver.c
@@ -71,12 +71,12 @@ int pcg(struct diag* a, double* b, double* x0, struct diag* c, double tol)
 	}
 
 	multiply(c, r0, z0);
-	
+
 
 	memcpy(p, z0, len_c * sizeof(double));
 
 	// loop
-	while (count < 100) {
+	while (count < 1000) {
 
 		multiply(a, p, t);
 		alph = dot(r0, z0, len_c) / dot(p, t, len_c);
diff --git a/sysmat.c b/sysmat.c
index 2c422f7..b1626f4 100644
--- a/sysmat.c
+++ b/sysmat.c
@@ -57,115 +57,107 @@ void system_matrix(struct diag* m, int nx, int ny, double* val)
 		}
 	}
 
-	printf("%g\n", get(m, 4, 4));
-	printf("%g\n", get(m, 4, 3));
-	printf("%g\n", get(m, 4, 5));
-	printf("%g\n", get(m, 4, 1));
-	printf("%g\n", get(m, 4, 7));
-
-//	// boundary cells with 2 neighbors
-//
-//	// upper left corner : r = 0, c = 0;
-//	double kb1 = harmonic(val[0], val[1]);
-//	double kb2 = harmonic(val[0], val[nx]);
-//
-//	put(m, 0, 0,    kb1 + kb2);
-//	put(m, 0, 1,  - kb1);
-//	put(m, 0, nx, - kb2);
-//
-//	// upper right corner : r = 0, c = nx
-//	kb1 = harmonic(val[nx - 1], val[nx - 2]);
-//	kb2 = harmonic(val[nx - 1], val[2 * nx - 1]);
-//
-//	put(m, nx - 1, nx - 1,       kb1 + kb2);
-//	put(m, nx - 1, nx - 2,     - kb1);
-//	put(m, nx - 1, 2 * nx - 1, - kb2);
-//
-//
-//	// lower left corner : r = ny, c = 0
-//	kb1 = harmonic(val[dim - nx], val[dim - nx + 1]);
-//	kb2 = harmonic(val[dim - nx], val[dim - 2 * nx]);
-//
-//	put(m, dim - nx, dim - nx,       kb1 + kb2);
-//	put(m, dim - nx, dim - nx + 1, - kb1);
-//	put(m, dim - nx, dim - 2 * nx, - kb2);
-//
-//	// lower right corner
-//	kb1 = harmonic(val[dim - 1], val[dim - 2]);
-//	kb2 = harmonic(val[dim - 1], val[dim - 1 - nx]);
-//
-//	put(m, dim - 1, dim - 1,        kb1 + kb2);
-//	put(m, dim - 1, dim - 2,      - kb1);
-//	put(m, dim - 1, dim - 1 - nx, - kb2);
-//
-//
-//	// boundary cells with 3 neighbors
-//
-//	// additional value required
-//	double kb3 = 0.0;
-//
-//	// upper boundary (i = 0)
-//	for (int j = 1; j < nx - 1; j ++) {
-//		kb1 = harmonic(val[j], val[j - 1]);
-//		kb2 = harmonic(val[j], val[j + 1]);
-//		kb3 = harmonic(val[j], val[nx + j]);
-//
-//		put(m, j, j,        kb1 + kb2 + kb3);
-//		put(m, j, j - 1,  - kb1);
-//		put(m, j, j + 1,  - kb2);
-//		put(m, j, nx + j, - kb3);
-//
-//	}
-//
-//	// lower boundary (i = n_max)
-//	for (int j = 1; j < nx - 1; j ++) {
-//
-//		int c = dim - nx;
-//
-//		kb1 = harmonic(val[c + j], val[c + j - 1]);
-//		kb2 = harmonic(val[c + j], val[c + j + 1]);
-//		kb3 = harmonic(val[c + j], val[dim - 2 * nx + j]);
-//
-//		put(m, c + j, c + j,              kb1 + kb2 + kb3);
-//		put(m, c + j, c + j - 1,        - kb1);
-//		put(m, c + j, c + j + 1,        - kb2);
-//		put(m, c + j, dim - 2 * nx + j, - kb3);
-//
-//	}
-//
-//	// left boundary (j = 0)
-//	for (int i = 1; i < ny - 1; i ++) {
-//
-//		int c = i * nx;
-//
-//		kb1 = harmonic(val[c], val[(i - 1) * nx]);
-//		kb2 = harmonic(val[c], val[(i + 1) * nx]);
-//		kb3 = harmonic(val[c], val[i * nx + 1]);
-//
-//		put(m, c, c,              kb1 + kb2 + kb3);
-//		put(m, c, (i - 1) * nx, - kb1);
-//		put(m, c, (i + 1) * nx, - kb2);
-//		put(m, c, i * nx + 1,   - kb3);
-//
-//	}
-//
-//	// right boundary (j = n_max)
-//	for (int i = 1; i < ny - 1; i ++) {
-//
-//		int c = (i + 1) * nx - 1;
-//
-//		kb1 = harmonic(val[c], val[i * nx - 1]);
-//		kb2 = harmonic(val[c], val[(i + 2) * nx - 1]);
-//		kb3 = harmonic(val[c], val[i * nx + nx - 2]);
-//
-//		put(m, c, c,                  kb1 + kb2 + kb3);
-//		put(m, c, i * nx - 1,       - kb1);
-//		put(m, c, (i + 2) * nx - 1, - kb2);
-//		put(m, c, i * nx + nx - 2,  - kb3);
-//
-//	}
-
-	//return m;
+	// boundary cells with 2 neighbors
+
+	// upper left corner : r = 0, c = 0;
+	double kb1 = harmonic(val[0], val[1]);
+	double kb2 = harmonic(val[0], val[nx]);
+
+	put(m, 0, 0,    kb1 + kb2);
+	put(m, 0, 1,  - kb1);
+	put(m, 0, nx, - kb2);
+
+	// upper right corner : r = 0, c = nx
+	kb1 = harmonic(val[nx - 1], val[nx - 2]);
+	kb2 = harmonic(val[nx - 1], val[2 * nx - 1]);
+
+	put(m, nx - 1, nx - 1,       kb1 + kb2);
+	put(m, nx - 1, nx - 2,     - kb1);
+	put(m, nx - 1, 2 * nx - 1, - kb2);
+
+
+	// lower left corner : r = ny, c = 0
+	kb1 = harmonic(val[dim - nx], val[dim - nx + 1]);
+	kb2 = harmonic(val[dim - nx], val[dim - 2 * nx]);
+
+	put(m, dim - nx, dim - nx,       kb1 + kb2);
+	put(m, dim - nx, dim - nx + 1, - kb1);
+	put(m, dim - nx, dim - 2 * nx, - kb2);
+
+	// lower right corner
+	kb1 = harmonic(val[dim - 1], val[dim - 2]);
+	kb2 = harmonic(val[dim - 1], val[dim - 1 - nx]);
+
+	put(m, dim - 1, dim - 1,        kb1 + kb2);
+	put(m, dim - 1, dim - 2,      - kb1);
+	put(m, dim - 1, dim - 1 - nx, - kb2);
+
+
+	// boundary cells with 3 neighbors
+
+	// additional value required
+	double kb3 = 0.0;
+
+	// upper boundary (i = 0)
+	for (int j = 1; j < nx - 1; j ++) {
+		kb1 = harmonic(val[j], val[j - 1]);
+		kb2 = harmonic(val[j], val[j + 1]);
+		kb3 = harmonic(val[j], val[nx + j]);
+
+		put(m, j, j,        kb1 + kb2 + kb3);
+		put(m, j, j - 1,  - kb1);
+		put(m, j, j + 1,  - kb2);
+		put(m, j, nx + j, - kb3);
+
+	}
+
+	// lower boundary (i = n_max)
+	for (int j = 1; j < nx - 1; j ++) {
+
+		int c = dim - nx;
+
+		kb1 = harmonic(val[c + j], val[c + j - 1]);
+		kb2 = harmonic(val[c + j], val[c + j + 1]);
+		kb3 = harmonic(val[c + j], val[dim - 2 * nx + j]);
+
+		put(m, c + j, c + j,              kb1 + kb2 + kb3);
+		put(m, c + j, c + j - 1,        - kb1);
+		put(m, c + j, c + j + 1,        - kb2);
+		put(m, c + j, dim - 2 * nx + j, - kb3);
+
+	}
+
+	// left boundary (j = 0)
+	for (int i = 1; i < ny - 1; i ++) {
+
+		int c = i * nx;
+
+		kb1 = harmonic(val[c], val[(i - 1) * nx]);
+		kb2 = harmonic(val[c], val[(i + 1) * nx]);
+		kb3 = harmonic(val[c], val[i * nx + 1]);
+
+		put(m, c, c,              kb1 + kb2 + kb3);
+		put(m, c, (i - 1) * nx, - kb1);
+		put(m, c, (i + 1) * nx, - kb2);
+		put(m, c, i * nx + 1,   - kb3);
+
+	}
+
+	// right boundary (j = n_max)
+	for (int i = 1; i < ny - 1; i ++) {
+
+		int c = (i + 1) * nx - 1;
+
+		kb1 = harmonic(val[c], val[i * nx - 1]);
+		kb2 = harmonic(val[c], val[(i + 2) * nx - 1]);
+		kb3 = harmonic(val[c], val[i * nx + nx - 2]);
+
+		put(m, c, c,                  kb1 + kb2 + kb3);
+		put(m, c, i * nx - 1,       - kb1);
+		put(m, c, (i + 2) * nx - 1, - kb2);
+		put(m, c, i * nx + nx - 2,  - kb3);
+
+	}
 
 }