ErikFabrizzi/NuSiF-Solver
1
0
Fork 0
forked from moebiusband/NuSiF-Solver
Code Pull Requests Activity

ft: 2D Poisson solver with Square domain decomposition

Browse Source
This commit is contained in:
indiano 2025-05-04 01:26:20 +02:00
parent 8a218ad681
commit 899d2f162d
3 changed files with 46 additions and 21 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
PoissonSolver/2D-mpi-sq/src/main.c
View File

@ -32,9 +32,9 @@ int main(int argc, char** argv)
readParameter(&params, argv[1]);
initSolver(&solver, &params, 2);
// startTime = getTimeStamp();
// solveRB(&solver);
// endTime = getTimeStamp();
startTime = getTimeStamp();
solveRB(&solver);
endTime = getTimeStamp();
writeResult(&solver, "p.dat");
// printf(" %.2fs\n", endTime - startTime);
// writeResult(&solver, "p.dat");

45
PoissonSolver/2D-mpi-sq/src/solver.c
View File

@ -60,6 +60,15 @@ void initSolver(Solver* solver, Parameter* params, int problem)
0,
&solver->cart_comm);
MPI_Cart_coords(solver->cart_comm, rank, 2, solver->coords);
solver->boundary = NB;
int x_low, x_high, y_low, y_high;
MPI_Cart_shift(solver->cart_comm, 0, 1, &x_low, &x_high);
MPI_Cart_shift(solver->cart_comm, 1, 1, &y_low, &y_high);
if (x_low == MPI_PROC_NULL) solver->boundary += BOTTOM;
if (y_low == MPI_PROC_NULL) solver->boundary += LEFT;
if (x_high == MPI_PROC_NULL) solver->boundary += TOP;
if (y_high == MPI_PROC_NULL) solver->boundary += RIGHT;
solver->limax = distribute_dim(solver->coords[0], solver->dims[0], solver->imax);
solver->ljmax = distribute_dim(solver->coords[1], solver->dims[1], solver->jmax);
@ -107,11 +116,11 @@ void initSolver(Solver* solver, Parameter* params, int problem)
i_start,
j_start);
fflush(stdout);
for (int i = 0; i < solver->limax + 2; i++) {
for (int j = 0; j < solver->ljmax + 2; j++) {
P(i, j) = rank;
}
}
// for (int i = 0; i < solver->limax + 2; i++) {
// for (int j = 0; j < solver->ljmax + 2; j++) {
// P(i, j) = (double)solver->boundary;
// }
// }
for (int i = 0; i < solver->limax + 2; i++) {
for (int j = 0; j < solver->ljmax + 2; j++) {
P(i, j) = sin(2.0 * PI * (i_start + i) * dx * 2.0) +
@ -181,16 +190,19 @@ void solveRB(Solver* solver)
}
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);
// }
MPI_Allreduce(MPI_IN_PLACE, &res, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
if (solver->boundary & BOTTOM)
for (int j = 1; j < solver->ljmax + 1; j++)
P(0, j) = P(1, j);
if (solver->boundary & TOP)
for (int j = 1; j < solver->ljmax + 1; j++)
P(solver->limax + 1, j) = P(solver->limax, j);
if (solver->boundary & LEFT)
for (int i = 1; i < solver->limax + 1; i++)
P(i, 0) = P(i, 1);
if (solver->boundary & RIGHT)
for (int i = 1; i < solver->limax + 1; i++)
P(i, solver->ljmax + 1) = P(i, solver->ljmax);
res = res / (double)(imax * jmax);
#ifdef DEBUG
@ -199,8 +211,7 @@ void solveRB(Solver* solver)
it++;
}
// printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
printf("%d, %f\n", it, solver->omega);
if(solver->rank==0) printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
}
void writeResult(Solver* solver, char* filename)

16
PoissonSolver/2D-mpi-sq/src/solver.h
View File

@ -6,11 +6,23 @@
#define __SOLVER_H_
#include "parameter.h"
#include <mpi.h>
typedef enum {
NB = 0b0000,
LEFT = 0b0001,
RIGHT = 0b0010,
TOP = 0b0100,
BOTTOM = 0b1000,
LFTO = LEFT + TOP,
LFBO = LEFT + BOTTOM,
RITO = RIGHT + TOP,
ROBO = RIGHT + BOTTOM
} boundary_t;
typedef struct {
double dx, dy;
int imax, jmax;
int limax,ljmax;
int limax, ljmax;
double *p, *rhs;
double eps, omega, rho;
int itermax;
@ -21,6 +33,8 @@ typedef struct {
int coords[2];
int dims[2];
int rank, size;
boundary_t boundary;
} Solver;
extern void initSolver(Solver*, Parameter*, int problem);