From 9c8585c7bc0919ce91a92fce6ca838a4c6969851 Mon Sep 17 00:00:00 2001 From: Jan Eitzinger Date: Sun, 18 Jun 2023 10:08:56 +0200 Subject: [PATCH] Add solver variants --- BasicSolver/2D-mpi-v3/src/solver.c | 1 + BasicSolver/2D-seq/src/solver.c | 62 ------- PoissonSolver/2D-seq/poisson.par | 3 +- PoissonSolver/2D-seq/src/main.c | 20 ++- PoissonSolver/2D-seq/src/solver.c | 278 ++++++++++++++++++++--------- PoissonSolver/2D-seq/src/solver.h | 10 +- 6 files changed, 217 insertions(+), 157 deletions(-) diff --git a/BasicSolver/2D-mpi-v3/src/solver.c b/BasicSolver/2D-mpi-v3/src/solver.c index cf9baf5..4721a95 100644 --- a/BasicSolver/2D-mpi-v3/src/solver.c +++ b/BasicSolver/2D-mpi-v3/src/solver.c @@ -428,6 +428,7 @@ int solve(Solver* solver) double dy2 = solver->dy * solver->dy; double idx2 = 1.0 / dx2; double idy2 = 1.0 / dy2; + // identical to 1/((2/dx2)+(2/dy2)) double factor = solver->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2); double* p = solver->p; double* rhs = solver->rhs; diff --git a/BasicSolver/2D-seq/src/solver.c b/BasicSolver/2D-seq/src/solver.c index 1005a93..7720be1 100644 --- a/BasicSolver/2D-seq/src/solver.c +++ b/BasicSolver/2D-seq/src/solver.c @@ -247,68 +247,6 @@ void solveRB(Solver* solver) #endif } -void solveRBA(Solver* solver) -{ - int imax = solver->imax; - int jmax = solver->jmax; - double eps = solver->eps; - int itermax = solver->itermax; - double dx2 = solver->dx * solver->dx; - double dy2 = solver->dy * solver->dy; - double idx2 = 1.0 / dx2; - double idy2 = 1.0 / dy2; - double factor = 0.5 * (dx2 * dy2) / (dx2 + dy2); - double omega = solver->omega; - double* p = solver->p; - double* rhs = solver->rhs; - double epssq = eps * eps; - int it = 0; - double anorm = 1.0, anormf = 0.0; - int pass, jsw, isw; - - while ((anorm >= eps * anormf) && (it < itermax)) { - anorm = 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); - - anorm += fabs(r); - P(i, j) -= (omega * factor * 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); - } - -#ifdef DEBUG - printf("%d Residuum: %e\n", it, res); -#endif - it++; - } - -#ifdef VERBOSE - printf("Solver took %d iterations to reach %f\n", it, anorm); -#endif -} - static double maxElement(Solver* solver, double* m) { int size = (solver->imax + 2) * (solver->jmax + 2); diff --git a/PoissonSolver/2D-seq/poisson.par b/PoissonSolver/2D-seq/poisson.par index 1f23554..08ce202 100644 --- a/PoissonSolver/2D-seq/poisson.par +++ b/PoissonSolver/2D-seq/poisson.par @@ -16,6 +16,7 @@ jmax 200 # number of interior cells in y-direction itermax 10000 # maximal number of pressure iteration in one time step eps 0.000001 # stopping tolerance for pressure iteration -omg 1.9 # relaxation parameter for SOR iteration +# omg 0.9999 # relaxation parameter for SOR iteration +omg 1.93 # relaxation parameter for SOR iteration #=============================================================================== diff --git a/PoissonSolver/2D-seq/src/main.c b/PoissonSolver/2D-seq/src/main.c index 588a848..2a4497b 100644 --- a/PoissonSolver/2D-seq/src/main.c +++ b/PoissonSolver/2D-seq/src/main.c @@ -4,24 +4,25 @@ * Use of this source code is governed by a MIT-style * license that can be found in the LICENSE file. */ -#include -#include -#include -#include #include +#include +#include +#include +#include #include "parameter.h" #include "solver.h" +#include "timing.h" - -int main (int argc, char** argv) +int main(int argc, char** argv) { + double startTime, endTime; Parameter params; Solver solver; initParameter(¶ms); - if ( argc != 2 ) { - printf("Usage: %s \n",argv[0]); + if (argc != 2) { + printf("Usage: %s \n", argv[0]); exit(EXIT_SUCCESS); } @@ -29,7 +30,10 @@ int main (int argc, char** argv) printParameter(¶ms); initSolver(&solver, ¶ms, 2); + startTime = getTimeStamp(); solve(&solver); + endTime = getTimeStamp(); + printf("Solution took %.2fs\n", endTime - startTime); writeResult(&solver); return EXIT_SUCCESS; diff --git a/PoissonSolver/2D-seq/src/solver.c b/PoissonSolver/2D-seq/src/solver.c index 3a3accd..c0395f8 100644 --- a/PoissonSolver/2D-seq/src/solver.c +++ b/PoissonSolver/2D-seq/src/solver.c @@ -4,140 +4,254 @@ * Use of this source code is governed by a MIT style * license that can be found in the LICENSE file. */ -#include "stdlib.h" -#include "stdio.h" #include "math.h" +#include "stdio.h" +#include "stdlib.h" -#include "solver.h" -#include "parameter.h" #include "allocate.h" +#include "parameter.h" +#include "solver.h" -#define PI 3.14159265358979323846 -#define P(i,j) p[(j)*(imax+2) + (i)] -#define RHS(i,j) rhs[(j)*(imax+2) + (i)] +#define PI 3.14159265358979323846 +#define P(i, j) p[(j) * (imax + 2) + (i)] +#define RHS(i, j) rhs[(j) * (imax + 2) + (i)] -void initSolver(Solver *solver, Parameter *params, int problem) +void initSolver(Solver* solver, Parameter* params, int problem) { - solver->imax = params->imax; - solver->jmax = params->jmax; - solver->dx = params->xlength/params->imax; - solver->dy = params->ylength/params->jmax; - solver->eps = params->eps; - solver->omega = params->omg; + solver->imax = params->imax; + solver->jmax = params->jmax; + solver->dx = params->xlength / params->imax; + solver->dy = params->ylength / params->jmax; + solver->eps = params->eps; + solver->omega = params->omg; solver->itermax = params->itermax; - int imax = solver->imax; - int jmax = solver->jmax; - size_t bytesize = (imax+2) * (jmax+2) * sizeof(double); - solver->p = allocate(64, bytesize); - solver->rhs = allocate(64, bytesize); + int imax = solver->imax; + int jmax = solver->jmax; + size_t bytesize = (imax + 2) * (jmax + 2) * sizeof(double); + solver->p = allocate(64, bytesize); + solver->rhs = allocate(64, bytesize); - double dx = solver->dx; - double dy = solver->dy; - double* p = solver->p; + double dx = solver->dx; + double dy = solver->dy; + double* p = solver->p; double* rhs = solver->rhs; - for( int j=0; jimax; + int jmax = solver->jmax; + double eps = solver->eps; + int itermax = solver->itermax; + double dx2 = solver->dx * solver->dx; + double dy2 = solver->dy * solver->dy; + double idx2 = 1.0 / dx2; + double idy2 = 1.0 / dy2; + double factor = solver->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2); + double* p = solver->p; + double* rhs = solver->rhs; + double epssq = eps * eps; + int it = 0; + double res = 1.0; - int imax = solver->imax; - int jmax = solver->jmax; - double eps= solver->eps; - double omega = solver->omega; - int itermax = solver->itermax; - - double dx2 = solver->dx * solver->dx; - double dy2 = solver->dy * solver->dy; - double idx2 = 1.0/dx2; - double idy2 = 1.0/dy2; - double factor = omega * 0.5 * (dx2*dy2) / (dx2+dy2); - double* p = solver->p; - double* rhs = solver->rhs; - - res = eps + 1.0; - - while((res >= eps) && (it < itermax)) { + while ((res >= epssq) && (it < itermax)) { res = 0.0; - for( int j=1; jimax; + int jmax = solver->jmax; + double eps = solver->eps; + int itermax = solver->itermax; + double dx2 = solver->dx * solver->dx; + double dy2 = solver->dy * solver->dy; + double idx2 = 1.0 / dx2; + double idy2 = 1.0 / dy2; + double factor = solver->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2); + double* p = solver->p; + double* rhs = solver->rhs; + 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++; } + + printf("Solver took %d iterations to reach %f\n", it, sqrt(res)); +} + +void solveRBA(Solver* solver) +{ + int imax = solver->imax; + int jmax = solver->jmax; + double eps = solver->eps; + int itermax = solver->itermax; + double dx2 = solver->dx * solver->dx; + double dy2 = solver->dy * solver->dy; + double idx2 = 1.0 / dx2; + double idy2 = 1.0 / dy2; + double factor = 0.5 * (dx2 * dy2) / (dx2 + dy2); + double rho = solver->omega; + double* p = solver->p; + double* rhs = solver->rhs; + double epssq = eps * eps; + int it = 0; + double res = 1.0; + int pass, jsw, isw; + double omega = 1.0; + + 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) -= (omega * factor * r); + res += (r * r); + } + isw = 3 - isw; + } + jsw = 3 - jsw; + omega = (it == 0 && pass == 0 ? 1.0 / (1.0 - 0.5 * rho * rho) + : 1.0 / (1.0 - 0.25 * rho * rho * omega)); + } + + 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 Omega: %e\n", it, res, omega); +#endif + it++; + } + + printf("Final omega: %f\n", omega); + printf("Solver took %d iterations to reach %f\n", it, sqrt(res)); } void writeResult(Solver* solver) { - int imax = solver->imax; - int jmax = solver->jmax; + int imax = solver->imax; + int jmax = solver->jmax; double* p = solver->p; - FILE *fp; - fp= fopen("p.dat", "w"); + FILE* fp; + fp = fopen("p.dat", "w"); - if (fp== NULL) { + if (fp == NULL) { printf("Error!\n"); exit(EXIT_FAILURE); } - for( int j=0; j