diff --git a/EnhancedSolver/2D-seq/config.mk b/EnhancedSolver/2D-seq/config.mk index ccdfa2c..ceed6bf 100644 --- a/EnhancedSolver/2D-seq/config.mk +++ b/EnhancedSolver/2D-seq/config.mk @@ -1,7 +1,7 @@ # Supported: GCC, CLANG, ICC TAG ?= CLANG ENABLE_OPENMP ?= false -# Supported: sor, mg +# Supported: sor, rb, mg SOLVER ?= mg # Run in debug settings DEBUG ?= false diff --git a/EnhancedSolver/2D-seq/src/solver-rb.c b/EnhancedSolver/2D-seq/src/solver-rb.c new file mode 100644 index 0000000..13d253b --- /dev/null +++ b/EnhancedSolver/2D-seq/src/solver-rb.c @@ -0,0 +1,93 @@ +/* + * Copyright (C) NHR@FAU, University Erlangen-Nuremberg. + * All rights reserved. This file is part of nusif-solver. + * Use of this source code is governed by a MIT style + * license that can be found in the LICENSE file. + */ +#include "grid.h" +#include "solver.h" +#include "util.h" + +void initSolver(Solver* s, Discretization* d, Parameter* p) +{ + s->grid = &d->grid; + s->itermax = p->itermax; + s->eps = p->eps; + s->omega = p->omg; + + Grid* g = s->grid; + int imax = s->grid->imax; + int jmax = s->grid->jmax; + + s->totalFluidCells = 0; + for (int j = 0; j < jmax + 2; j++) { + for (int i = 0; i < imax + 2; i++) { + if (gridIsFluid(g, i, j)) { + s->totalFluidCells++; + } + } + } +} + +void solve(Solver* s, double* p, double* rhs) +{ + int imax = s->grid->imax; + int jmax = s->grid->jmax; + double eps = s->eps; + int itermax = s->itermax; + double dx2 = s->grid->dx * s->grid->dx; + double dy2 = s->grid->dy * s->grid->dy; + double idx2 = 1.0 / dx2; + double idy2 = 1.0 / dy2; + double factor = s->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2); + double epssq = eps * eps; + int it = 0; + Grid* g = s->grid; + double res = 1.0; + int pass, jsw, isw; + + while ((res >= epssq) && (it < itermax)) { + res = 0.0; + jsw = 1; + + for (pass = 0; pass < 2; pass++) { + isw = jsw; + + for (int j = 1; j < jmax + 1; j++) { + for (int i = isw; i < imax + 1; i += 2) { + if (gridIsFluid(g, i, j)) { + double r = RHS(i, j) - + ((P(i + 1, j) - 2.0 * P(i, j) + P(i - 1, j)) * idx2 + + (P(i, j + 1) - 2.0 * P(i, j) + P(i, j - 1)) * + idy2); + + P(i, j) -= (factor * r); + res += (r * r); + } + } + isw = 3 - isw; + } + jsw = 3 - jsw; + } + + for (int i = 1; i < imax + 1; i++) { + P(i, 0) = P(i, 1); + P(i, jmax + 1) = P(i, jmax); + } + + for (int j = 1; j < jmax + 1; j++) { + P(0, j) = P(1, j); + P(imax + 1, j) = P(imax, j); + } + + res = res / (double)(s->totalFluidCells); +#ifdef DEBUG + printf("%d Residuum: %e\n", it, res); +#endif + it++; + } + +#ifdef VERBOSE + printf("Solver took %d iterations to reach %f\n", it, sqrt(res)); +#endif +} diff --git a/EnhancedSolver/2D-seq/src/solver-sor.c b/EnhancedSolver/2D-seq/src/solver-sor.c index 165e103..de5ceb9 100644 --- a/EnhancedSolver/2D-seq/src/solver-sor.c +++ b/EnhancedSolver/2D-seq/src/solver-sor.c @@ -29,7 +29,7 @@ void initSolver(Solver* s, Discretization* d, Parameter* p) } } -void solveSOR(Solver* s, double* p, double* rhs) +void solve(Solver* s, double* p, double* rhs) { int imax = s->grid->imax; int jmax = s->grid->jmax; @@ -42,6 +42,7 @@ void solveSOR(Solver* s, double* p, double* rhs) double factor = s->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2); double epssq = eps * eps; int it = 0; + Grid* g = s->grid; double res = 1.0; while ((res >= epssq) && (it < itermax)) { @@ -49,13 +50,14 @@ void solveSOR(Solver* s, double* p, double* rhs) for (int j = 1; j < jmax + 1; j++) { for (int i = 1; i < imax + 1; i++) { + if (gridIsFluid(g, i, j)) { + double r = RHS(i, j) - + ((P(i + 1, j) - 2.0 * P(i, j) + P(i - 1, j)) * idx2 + + (P(i, j + 1) - 2.0 * P(i, j) + P(i, j - 1)) * idy2); - double r = RHS(i, j) - - ((P(i + 1, j) - 2.0 * P(i, j) + P(i - 1, j)) * idx2 + - (P(i, j + 1) - 2.0 * P(i, j) + P(i, j - 1)) * idy2); - - P(i, j) -= (factor * r); - res += (r * r); + P(i, j) -= (factor * r); + res += (r * r); + } } } @@ -80,66 +82,3 @@ void solveSOR(Solver* s, double* p, double* rhs) printf("Solver took %d iterations to reach %f\n", it, sqrt(res)); #endif } - -void solve(Solver* solver, double* p, double* rhs) -{ - int imax = solver->grid->imax; - int jmax = solver->grid->jmax; - double eps = solver->eps; - int itermax = solver->itermax; - double dx2 = solver->grid->dx * solver->grid->dx; - double dy2 = solver->grid->dy * solver->grid->dy; - double idx2 = 1.0 / dx2; - double idy2 = 1.0 / dy2; - double factor = solver->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2); - double epssq = eps * eps; - int it = 0; - Grid* g = solver->grid; - double res = 1.0; - int pass, jsw, isw; - - while ((res >= epssq) && (it < itermax)) { - res = 0.0; - jsw = 1; - - for (pass = 0; pass < 2; pass++) { - isw = jsw; - - for (int j = 1; j < jmax + 1; j++) { - for (int i = isw; i < imax + 1; i += 2) { - if (gridIsFluid(g, i, j)) { - double r = RHS(i, j) - - ((P(i + 1, j) - 2.0 * P(i, j) + P(i - 1, j)) * idx2 + - (P(i, j + 1) - 2.0 * P(i, j) + P(i, j - 1)) * - idy2); - - P(i, j) -= (factor * r); - res += (r * r); - } - } - isw = 3 - isw; - } - jsw = 3 - jsw; - } - - for (int i = 1; i < imax + 1; i++) { - P(i, 0) = P(i, 1); - P(i, jmax + 1) = P(i, jmax); - } - - for (int j = 1; j < jmax + 1; j++) { - P(0, j) = P(1, j); - P(imax + 1, j) = P(imax, j); - } - - res = res / (double)(imax * jmax); -#ifdef DEBUG - printf("%d Residuum: %e\n", it, res); -#endif - it++; - } - -#ifdef VERBOSE - printf("Solver took %d iterations to reach %f\n", it, sqrt(res)); -#endif -}