From b9bf4d7b63a1c9eed90c0a199e34449fd0b14a98 Mon Sep 17 00:00:00 2001 From: Jan Eitzinger Date: Wed, 6 Mar 2024 07:27:28 +0100 Subject: [PATCH] Split red-balack and standard sor solvers --- BasicSolver/2D-seq/config.mk | 4 +- BasicSolver/2D-seq/dcavity.par | 4 +- BasicSolver/2D-seq/src/solver-rb.c | 76 +++++++++++++++++++++++++++++ BasicSolver/2D-seq/src/solver-sor.c | 62 +---------------------- 4 files changed, 81 insertions(+), 65 deletions(-) create mode 100644 BasicSolver/2D-seq/src/solver-rb.c diff --git a/BasicSolver/2D-seq/config.mk b/BasicSolver/2D-seq/config.mk index a782baf..dc874ee 100644 --- a/BasicSolver/2D-seq/config.mk +++ b/BasicSolver/2D-seq/config.mk @@ -1,8 +1,8 @@ # Supported: GCC, CLANG, ICC TAG ?= CLANG ENABLE_OPENMP ?= false -# Supported: sor, mg -SOLVER ?= sor +# Supported: sor, rb, mg +SOLVER ?= rb # Run in debug settings DEBUG ?= false diff --git a/BasicSolver/2D-seq/dcavity.par b/BasicSolver/2D-seq/dcavity.par index b1c238d..798819f 100644 --- a/BasicSolver/2D-seq/dcavity.par +++ b/BasicSolver/2D-seq/dcavity.par @@ -26,8 +26,8 @@ p_init 0.0 # initial value for pressure xlength 1.0 # domain size in x-direction ylength 1.0 # domain size in y-direction -imax 100 # number of interior cells in x-direction -jmax 100 # number of interior cells in y-direction +imax 128 # number of interior cells in x-direction +jmax 128 # number of interior cells in y-direction # Time Data: # --------- diff --git a/BasicSolver/2D-seq/src/solver-rb.c b/BasicSolver/2D-seq/src/solver-rb.c new file mode 100644 index 0000000..265456a --- /dev/null +++ b/BasicSolver/2D-seq/src/solver-rb.c @@ -0,0 +1,76 @@ +/* + * 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 "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; +} + +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; + 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) { + + 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 +} diff --git a/BasicSolver/2D-seq/src/solver-sor.c b/BasicSolver/2D-seq/src/solver-sor.c index f817b5f..3d21ea5 100644 --- a/BasicSolver/2D-seq/src/solver-sor.c +++ b/BasicSolver/2D-seq/src/solver-sor.c @@ -15,7 +15,7 @@ void initSolver(Solver* s, Discretization* d, Parameter* p) s->omega = p->omg; } -void solveSOR(Solver* solver, double* p, double* rhs) +void solve(Solver* solver, double* p, double* rhs) { int imax = solver->grid->imax; int jmax = solver->grid->jmax; @@ -66,63 +66,3 @@ void solveSOR(Solver* solver, 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; - 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) { - - 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 -}