From 2fad29b9254dcfb97de3bfbdb2e722bf3cc4e05d Mon Sep 17 00:00:00 2001 From: Jan Eitzinger Date: Tue, 21 Nov 2023 05:27:11 +0100 Subject: [PATCH] Reformat. Merge improved solvers. --- .clangd | 4 +- PoissonSolver/2D-mpi/README.md | 4 +- PoissonSolver/2D-mpi/src/main.c | 17 +- PoissonSolver/2D-mpi/src/parameter.c | 43 +-- PoissonSolver/2D-mpi/src/parameter.h | 2 +- PoissonSolver/2D-mpi/src/solver.c | 426 +++++++++++++++++++-------- PoissonSolver/2D-mpi/src/solver.h | 4 +- PoissonSolver/2D-seq/p.dat | 302 +++++++++++++++++++ PoissonSolver/2D-seq/poisson.par | 8 +- PoissonSolver/2D-seq/src/main.c | 10 +- PoissonSolver/2D-seq/src/solver.c | 11 +- 11 files changed, 655 insertions(+), 176 deletions(-) create mode 100644 PoissonSolver/2D-seq/p.dat diff --git a/.clangd b/.clangd index 1ad7985..372227d 100644 --- a/.clangd +++ b/.clangd @@ -1,3 +1,3 @@ -CompileFlags: # Tweak the parse settings - Add: [-I/usr/local/include] # treat all files as C++, enable more warnings +CompileFlags: + Add: [-I/usr/local/include, -I/opt/homebrew/include] Compiler: clang diff --git a/PoissonSolver/2D-mpi/README.md b/PoissonSolver/2D-mpi/README.md index b0a80a6..2a898e3 100644 --- a/PoissonSolver/2D-mpi/README.md +++ b/PoissonSolver/2D-mpi/README.md @@ -2,7 +2,7 @@ ## Build -1. Configure the toolchain and additional options in `config.mk`: +1. Configure the tool chain and additional options in `config.mk`: ``` # Supported: GCC, CLANG, ICC TAG ?= GCC @@ -22,7 +22,7 @@ The verbosity options enable detailed output about affinity settings, allocation make ``` -You can build multiple toolchains in the same directory, but notice that the Makefile is only acting on the one currently set. +You can build multiple tool chains in the same directory, but notice that the Makefile is only acting on the one currently set. Intermediate build results are located in the `` directory. To output the executed commands use: diff --git a/PoissonSolver/2D-mpi/src/main.c b/PoissonSolver/2D-mpi/src/main.c index c797be3..981cac2 100644 --- a/PoissonSolver/2D-mpi/src/main.c +++ b/PoissonSolver/2D-mpi/src/main.c @@ -4,18 +4,17 @@ * 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 #include "parameter.h" #include "solver.h" - -int main (int argc, char** argv) +int main(int argc, char** argv) { int rank; Parameter params; @@ -25,13 +24,13 @@ int main (int argc, char** argv) MPI_Comm_rank(MPI_COMM_WORLD, &rank); initParameter(¶ms); - if ( argc != 2 ) { - printf("Usage: %s \n",argv[0]); + if (argc != 2) { + printf("Usage: %s \n", argv[0]); exit(EXIT_SUCCESS); } readParameter(¶ms, argv[1]); - if ( rank == 0 ) { + if (rank == 0) { printParameter(¶ms); } diff --git a/PoissonSolver/2D-mpi/src/parameter.c b/PoissonSolver/2D-mpi/src/parameter.c index 6c08af9..0a4b455 100644 --- a/PoissonSolver/2D-mpi/src/parameter.c +++ b/PoissonSolver/2D-mpi/src/parameter.c @@ -12,41 +12,47 @@ #include "util.h" #define MAXLINE 4096 -void initParameter(Parameter *param) { +void initParameter(Parameter* param) +{ param->xlength = 1.0; param->ylength = 1.0; - param->imax = 100; - param->jmax = 100; + param->imax = 100; + param->jmax = 100; param->itermax = 1000; - param->eps = 0.0001; - param->omg = 1.8; + param->eps = 0.0001; + param->omg = 1.8; } -void readParameter(Parameter *param, const char *filename) { - FILE *fp = fopen(filename, "r"); +void readParameter(Parameter* param, const char* filename) +{ + FILE* fp = fopen(filename, "r"); char line[MAXLINE]; int i; - if(!fp) { + if (!fp) { fprintf(stderr, "Could not open parameter file: %s\n", filename); exit(EXIT_FAILURE); } - while(!feof(fp)) { + while (!feof(fp)) { line[0] = '\0'; fgets(line, MAXLINE, fp); - for(i = 0; line[i] != '\0' && line[i] != '#'; i++); + for (i = 0; line[i] != '\0' && line[i] != '#'; i++) + ; line[i] = '\0'; - char *tok = strtok(line, " "); - char *val = strtok(NULL, " "); + char* tok = strtok(line, " "); + char* val = strtok(NULL, " "); - #define PARSE_PARAM(p,f) if(strncmp(tok, #p, sizeof(#p) / sizeof(#p[0]) - 1) == 0) { param->p = f(val); } - #define PARSE_STRING(p) PARSE_PARAM(p, strdup) - #define PARSE_INT(p) PARSE_PARAM(p, atoi) - #define PARSE_REAL(p) PARSE_PARAM(p, atof) +#define PARSE_PARAM(p, f) \ + if (strncmp(tok, #p, sizeof(#p) / sizeof(#p[0]) - 1) == 0) { \ + param->p = f(val); \ + } +#define PARSE_STRING(p) PARSE_PARAM(p, strdup) +#define PARSE_INT(p) PARSE_PARAM(p, atoi) +#define PARSE_REAL(p) PARSE_PARAM(p, atof) - if(tok != NULL && val != NULL) { + if (tok != NULL && val != NULL) { PARSE_REAL(xlength); PARSE_REAL(ylength); PARSE_INT(imax); @@ -60,7 +66,8 @@ void readParameter(Parameter *param, const char *filename) { fclose(fp); } -void printParameter(Parameter *param) { +void printParameter(Parameter* param) +{ printf("Parameters:\n"); printf("Geometry data:\n"); printf("\tDomain box size (x, y): %e, %e\n", param->xlength, param->ylength); diff --git a/PoissonSolver/2D-mpi/src/parameter.h b/PoissonSolver/2D-mpi/src/parameter.h index d57ce19..6344a31 100644 --- a/PoissonSolver/2D-mpi/src/parameter.h +++ b/PoissonSolver/2D-mpi/src/parameter.h @@ -8,7 +8,7 @@ #define __PARAMETER_H_ typedef struct { - double xlength, ylength; + double xlength, ylength; int imax, jmax; int itermax; double eps, omg, gamma; diff --git a/PoissonSolver/2D-mpi/src/solver.c b/PoissonSolver/2D-mpi/src/solver.c index 7264c44..62f135e 100644 --- a/PoissonSolver/2D-mpi/src/solver.c +++ b/PoissonSolver/2D-mpi/src/solver.c @@ -4,49 +4,53 @@ * 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 "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)] static int sizeOfRank(int rank, int size, int N) { - return N/size + ((N%size>rank) ? 1 : 0); + return N / size + ((N % size > rank) ? 1 : 0); } static void print(Solver* solver) { double* p = solver->p; - int imax = solver->imax; + int imax = solver->imax; - printf("### RANK %d #######################################################\n", solver->rank); - for( int j=0; j < solver->jmaxLocal+2; j++ ) { + printf("### RANK %d #######################################################\n", + solver->rank); + for (int j = 0; j < solver->jmaxLocal + 2; j++) { printf("%02d: ", j); - for( int i=0; i < solver->imax+2; i++ ) { + for (int i = 0; i < solver->imax + 2; i++) { printf("%12.8f ", P(i, j)); } printf("\n"); } - fflush( stdout ); + fflush(stdout); } static void exchange(Solver* solver) { - MPI_Request requests[4] = { MPI_REQUEST_NULL, MPI_REQUEST_NULL, MPI_REQUEST_NULL, MPI_REQUEST_NULL }; + MPI_Request requests[4] = { MPI_REQUEST_NULL, + MPI_REQUEST_NULL, + MPI_REQUEST_NULL, + MPI_REQUEST_NULL }; /* exchange ghost cells with top neighbor */ if (solver->rank + 1 < solver->size) { - int top = solver->rank + 1; - double* src = solver->p + (solver->jmaxLocal) * (solver->imax+2) + 1; - double* dst = solver->p + (solver->jmaxLocal+1) * (solver->imax+2) + 1; + int top = solver->rank + 1; + double* src = solver->p + (solver->jmaxLocal) * (solver->imax + 2) + 1; + double* dst = solver->p + (solver->jmaxLocal + 1) * (solver->imax + 2) + 1; MPI_Isend(src, solver->imax, MPI_DOUBLE, top, 1, MPI_COMM_WORLD, &requests[0]); MPI_Irecv(dst, solver->imax, MPI_DOUBLE, top, 2, MPI_COMM_WORLD, &requests[1]); @@ -54,222 +58,382 @@ static void exchange(Solver* solver) /* exchange ghost cells with bottom neighbor */ if (solver->rank > 0) { - int bottom = solver->rank - 1; - double* src = solver->p + (solver->imax+2) + 1; + int bottom = solver->rank - 1; + double* src = solver->p + (solver->imax + 2) + 1; double* dst = solver->p + 1; - MPI_Isend(src, solver->imax, MPI_DOUBLE, bottom, 2, MPI_COMM_WORLD, &requests[2]); - MPI_Irecv(dst, solver->imax, MPI_DOUBLE, bottom, 1, MPI_COMM_WORLD, &requests[3]); + MPI_Isend(src, solver->imax, MPI_DOUBLE, bottom, 2, MPI_COMM_WORLD, &requests[2]); + MPI_Irecv(dst, solver->imax, MPI_DOUBLE, bottom, 1, MPI_COMM_WORLD, &requests[3]); } MPI_Waitall(4, requests, MPI_STATUSES_IGNORE); } -void getResult(Solver *solver) +void getResult(Solver* solver) { double* Pall = NULL; int *rcvCounts, *displs; - if ( solver->rank == 0 ) { - Pall = allocate(64, (solver->imax+2) * (solver->jmax+2) * sizeof(double)); - rcvCounts = (int*) malloc(solver->size * sizeof(int)); - displs = (int*) malloc(solver->size * sizeof(int)); - rcvCounts[0] = solver->jmaxLocal * (solver->imax+2); - displs[0] = 0; - int cursor = rcvCounts[0]; + if (solver->rank == 0) { + Pall = allocate(64, (solver->imax + 2) * (solver->jmax + 2) * sizeof(double)); + rcvCounts = (int*)malloc(solver->size * sizeof(int)); + displs = (int*)malloc(solver->size * sizeof(int)); + rcvCounts[0] = solver->jmaxLocal * (solver->imax + 2); + displs[0] = 0; + int cursor = rcvCounts[0]; - for ( int i=1; i < solver->size; i++ ) { - rcvCounts[i] = sizeOfRank(i, solver->size, solver->jmax) * (solver->imax+2); - displs[i] = cursor; + for (int i = 1; i < solver->size; i++) { + rcvCounts[i] = sizeOfRank(i, solver->size, solver->jmax) * (solver->imax + 2); + displs[i] = cursor; cursor += rcvCounts[i]; } } - int cnt = solver->jmaxLocal*(solver->imax+2); - double* sendbuffer = solver->p + (solver->imax+2); - MPI_Gatherv(sendbuffer, cnt, MPI_DOUBLE, Pall, - rcvCounts, displs, MPI_DOUBLE, 0, MPI_COMM_WORLD); + int cnt = solver->jmaxLocal * (solver->imax + 2); + double* sendbuffer = solver->p + (solver->imax + 2); + MPI_Gatherv(sendbuffer, + cnt, + MPI_DOUBLE, + Pall, + rcvCounts, + displs, + MPI_DOUBLE, + 0, + MPI_COMM_WORLD); - if ( solver->rank == 0 ) { + if (solver->rank == 0) { writeResult(solver, Pall, "p.dat"); } } -void initSolver(Solver *solver, Parameter *params, int problem) +void initSolver(Solver* solver, Parameter* params, int problem) { MPI_Comm_rank(MPI_COMM_WORLD, &(solver->rank)); MPI_Comm_size(MPI_COMM_WORLD, &(solver->size)); - solver->imax = params->imax; - solver->jmax = params->jmax; + solver->imax = params->imax; + solver->jmax = params->jmax; solver->jmaxLocal = sizeOfRank(solver->rank, solver->size, solver->jmax); - printf("RANK %d: %d\n", solver->rank, solver->jmaxLocal); + printf("RANK %d: imaxLocal : %d, jmaxLocal : %d\n", + solver->rank, + solver->imax, + solver->jmaxLocal); - solver->dx = params->xlength/params->imax; - solver->dy = params->ylength/params->jmax; - solver->ys = solver->rank * solver->jmaxLocal * solver->dy; - solver->eps = params->eps; - solver->omega = params->omg; + solver->dx = params->xlength / params->imax; + solver->dy = params->ylength / params->jmax; + solver->ys = solver->rank * solver->jmaxLocal * solver->dy; + solver->eps = params->eps; + solver->omega = params->omg; solver->itermax = params->itermax; - int imax = solver->imax; - int jmax = solver->jmax; + int imax = solver->imax; + int jmax = solver->jmax; int jmaxLocal = solver->jmaxLocal; - solver->p = allocate(64, (imax+2) * (jmaxLocal+2) * sizeof(double)); - solver->rhs = allocate(64, (imax+2) * (jmax+2) * sizeof(double)); + solver->p = allocate(64, (imax + 2) * (jmaxLocal + 2) * sizeof(double)); + solver->rhs = allocate(64, (imax + 2) * (jmax + 2) * sizeof(double)); - 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; jys + j * dy; - for( int i=0; iimax; - int rank = solver->rank; + int imax = solver->imax; + int rank = solver->rank; double* p = solver->p; -/* for( int j=0; j < solver->jmaxLocal+2; j++ ) { */ -/* for( int i=0; i < solver->imax+2; i++ ) { */ -/* P(i, j) = (double) rank; */ -/* } */ -/* } */ + /* for( int j=0; j < solver->jmaxLocal+2; j++ ) { */ + /* for( int i=0; i < solver->imax+2; i++ ) { */ + /* P(i, j) = (double) rank; */ + /* } */ + /* } */ -/* for ( int i=0; i < solver->size; i++) { */ -/* if ( i == rank ) { */ -/* print(solver); */ -/* } */ -/* MPI_Barrier(MPI_COMM_WORLD); */ -/* } */ + /* for ( int i=0; i < solver->size; i++) { */ + /* if ( i == rank ) { */ + /* print(solver); */ + /* } */ + /* MPI_Barrier(MPI_COMM_WORLD); */ + /* } */ -/* if ( rank == 0 ) { */ -/* printf("##########################################################\n"); */ -/* printf("## Exchange ghost layers\n"); */ -/* printf("##########################################################\n"); */ -/* } */ -/* exchange(solver); */ + /* if ( rank == 0 ) { */ + /* printf("##########################################################\n"); */ + /* printf("## Exchange ghost layers\n"); */ + /* printf("##########################################################\n"); */ + /* } */ + /* exchange(solver); */ - for ( int i=0; i < solver->size; i++) { - if ( i == rank ) { - print(solver); + for (int i = 0; i < solver->size; i++) { + if (i == rank) { + print(solver); } MPI_Barrier(MPI_COMM_WORLD); } } -int solve(Solver *solver) +int solve(Solver* solver) { double r; int it = 0; - double res; + double res, res1; - int imax = solver->imax; - int jmax = solver->jmax; + int imax = solver->imax; + int jmax = solver->jmax; int jmaxLocal = solver->jmaxLocal; - double eps= solver->eps; - double omega = solver->omega; - int itermax = solver->itermax; + 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; + 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; + double epssq = eps * eps; res = eps + 1.0; - while((res >= eps) && (it < itermax)) { + while ((res >= epssq) && (it < itermax)) { res = 0.0; exchange(solver); - for( int j=1; jrank == 0 ) { - for( int i=1; irank == 0) { + for (int i = 1; i < imax + 1; i++) { + P(i, 0) = P(i, 1); } } - if ( solver->rank == (solver->size-1) ) { - for( int i=1; irank == (solver->size - 1)) { + for (int i = 1; i < imax + 1; i++) { + P(i, jmaxLocal + 1) = P(i, jmaxLocal); } } - for( int j=1; jrank == 0 ) { - printf("%d Residuum: %e\n",it, res); + if (solver->rank == 0) { + printf("%d Residuum: %e\n", it, res1); } #endif it++; } - if ( solver->rank == 0 ) { - printf("Solver took %d iterations\n",it); + if (solver->rank == 0) { + printf("Solver took %d iterations\n", it); } - if( res < eps ){ + if (res < eps) { return 1; - } else{ + } else { return 0; } } +int solveRB(Solver* solver) +{ + double r; + int it = 0; + double res, res1; + + int imax = solver->imax; + int jmax = solver->jmax; + int jmaxLocal = solver->jmaxLocal; + 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; + int pass, jsw, isw; + double epssq = eps * eps; + + res = eps + 1.0; + + while ((res >= epssq) && (it < itermax)) { + res = 0.0; + jsw = 1; + + for (pass = 0; pass < 2; pass++) { + isw = jsw; + exchange(solver); + + for (int j = 1; j < jmaxLocal + 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, jmaxLocal + 1) = P(i, jmaxLocal); + } + + for (int j = 1; j < jmaxLocal + 1; j++) { + P(0, j) = P(1, j); + P(imax + 1, j) = P(imax, j); + } + MPI_Allreduce(&res, &res1, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + res = res1; + res = res / (double)(imax * jmax); +#ifdef DEBUG + printf("%d Residuum: %e\n", it, res); +#endif + it++; + } + + if (solver->rank == 0) { + printf("Solver took %d iterations\n", it); + } + if (res < eps) { + return 1; + } else { + return 0; + } +} + +int solveRBA(Solver* solver) +{ + double r; + int it = 0; + double res; + + int imax = solver->imax; + int jmax = solver->jmax; + int jmaxLocal = solver->jmaxLocal; + 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; + int pass, jsw, isw; + double rho = solver->rho; + double epssq = eps * eps; + + res = eps + 1.0; + + while ((res >= epssq) && (it < itermax)) { + res = 0.0; + jsw = 1; + + for (pass = 0; pass < 2; pass++) { + isw = jsw; + exchange(solver); + + for (int j = 1; j < jmaxLocal + 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, jmaxLocal + 1) = P(i, jmaxLocal); + } + + for (int j = 1; j < jmaxLocal + 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, double* m, char* filename) { - int imax = solver->imax; - int jmax = solver->jmax; + int imax = solver->imax; + int jmax = solver->jmax; double* p = solver->p; - FILE *fp; - fp= fopen(filename, "w"); + FILE* fp; + fp = fopen(filename, "w"); - if (fp== NULL) { + if (fp == NULL) { printf("Error!\n"); exit(EXIT_FAILURE); } - for( int j=0; jomega); } void solveRB(Solver* solver) @@ -168,7 +169,8 @@ void solveRB(Solver* solver) it++; } - printf("Solver took %d iterations to reach %f\n", it, sqrt(res)); + // printf("Solver took %d iterations to reach %f\n", it, sqrt(res)); + printf("%d, %f\n", it, solver->omega); } void solveRBA(Solver* solver) @@ -232,8 +234,9 @@ void solveRBA(Solver* solver) it++; } - printf("Final omega: %f\n", omega); - printf("Solver took %d iterations to reach %f\n", it, sqrt(res)); + // printf("Final omega: %f\n", omega); + // printf("Solver took %d iterations to reach %f\n", it, sqrt(res)); + printf("%d, %f\n", it, omega); } void writeResult(Solver* solver)