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

checkpoint: Cleanup

Browse Source
This commit is contained in:
indiano 2025-05-04 15:59:20 +02:00
parent 899d2f162d
commit f0eeef7450
3 changed files with 104 additions and 91 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -14,15 +14,11 @@ int main(int argc, char** argv)
{ {
MPI_Init(&argc, &argv); MPI_Init(&argc, &argv);
/******* Debug Params *********/
int rank, size;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Comm_size(MPI_COMM_WORLD, &size);
printf("Process %d of %d \n", rank, size);
double startTime, endTime; double startTime, endTime;
Parameter params; Parameter params;
Solver solver; Solver solver;
initParameter(&params); initParameter(&params);
if (argc < 2) { if (argc < 2) {
@ -36,9 +32,7 @@ int main(int argc, char** argv)
solveRB(&solver); solveRB(&solver);
endTime = getTimeStamp(); endTime = getTimeStamp();
writeResult(&solver, "p.dat"); writeResult(&solver, "p.dat");
// printf(" %.2fs\n", endTime - startTime); if(solver.rank == 0) printf(" %.2fs\n", endTime - startTime);
// writeResult(&solver, "p.dat");
MPI_Finalize(); MPI_Finalize();
return EXIT_SUCCESS; return EXIT_SUCCESS;

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

@ -11,6 +11,7 @@
#include "solver.h" #include "solver.h"
#include <alloca.h> #include <alloca.h>
#include <mpi.h> #include <mpi.h>
#include <stddef.h>
#include <stdlib.h> #include <stdlib.h>
#define PI 3.14159265358979323846 #define PI 3.14159265358979323846
@ -31,8 +32,36 @@ int get_dim_start(int dim_rank, int dim_comm_size, int dim_size)
return dim_start; return dim_start;
} }
void write_matrix_to_file(
const char* _Nonnull filename, double* matrix, size_t rows, size_t cols)
{
FILE* fp;
fp = fopen(filename, "w");
if (fp == NULL) {
printf("Error!\n");
exit(EXIT_FAILURE);
}
for (int i = 0; i < rows; i++) {
for (int j = 0; j < cols; j++) {
fprintf(fp, "%f ", matrix[i * (cols) + j]);
}
fprintf(fp, "\n");
}
fclose(fp);
}
void initSolver(Solver* solver, Parameter* params, int problem) void initSolver(Solver* solver, Parameter* params, int problem)
{ {
int x_low, x_high, y_low, y_high;
int rank, size;
MPI_Datatype contiguos_boundary, vector_boundary;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Comm_size(MPI_COMM_WORLD, &size);
solver->imax = params->imax; solver->imax = params->imax;
solver->jmax = params->jmax; solver->jmax = params->jmax;
solver->dx = params->xlength / params->imax; solver->dx = params->xlength / params->imax;
@ -41,17 +70,13 @@ void initSolver(Solver* solver, Parameter* params, int problem)
solver->omega = params->omg; solver->omega = params->omg;
solver->rho = params->rho; solver->rho = params->rho;
solver->itermax = params->itermax; solver->itermax = params->itermax;
int imax = solver->imax;
int jmax = solver->jmax;
int rank, size;
solver->dims[0] = 0; solver->dims[0] = 0;
solver->dims[1] = 0; solver->dims[1] = 0;
MPI_Comm_rank(MPI_COMM_WORLD, &rank); solver->rank = 0;
MPI_Comm_size(MPI_COMM_WORLD, &size); solver->size = 0;
solver->rank = rank; solver->rank = rank;
solver->size = size; solver->size = size;
MPI_Dims_create(size, 2, solver->dims); MPI_Dims_create(size, 2, solver->dims);
MPI_Cart_create(MPI_COMM_WORLD, MPI_Cart_create(MPI_COMM_WORLD,
2, 2,
@ -61,9 +86,9 @@ void initSolver(Solver* solver, Parameter* params, int problem)
&solver->cart_comm); &solver->cart_comm);
MPI_Cart_coords(solver->cart_comm, rank, 2, solver->coords); MPI_Cart_coords(solver->cart_comm, rank, 2, solver->coords);
solver->boundary = NB; 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, 0, 1, &x_low, &x_high);
MPI_Cart_shift(solver->cart_comm, 1, 1, &y_low, &y_high); MPI_Cart_shift(solver->cart_comm, 1, 1, &y_low, &y_high);
if (x_low == MPI_PROC_NULL) solver->boundary += BOTTOM; if (x_low == MPI_PROC_NULL) solver->boundary += BOTTOM;
if (y_low == MPI_PROC_NULL) solver->boundary += LEFT; if (y_low == MPI_PROC_NULL) solver->boundary += LEFT;
if (x_high == MPI_PROC_NULL) solver->boundary += TOP; if (x_high == MPI_PROC_NULL) solver->boundary += TOP;
@ -72,6 +97,15 @@ void initSolver(Solver* solver, Parameter* params, int problem)
solver->limax = distribute_dim(solver->coords[0], solver->dims[0], solver->imax); solver->limax = distribute_dim(solver->coords[0], solver->dims[0], solver->imax);
solver->ljmax = distribute_dim(solver->coords[1], solver->dims[1], solver->jmax); solver->ljmax = distribute_dim(solver->coords[1], solver->dims[1], solver->jmax);
MPI_Type_contiguous(solver->ljmax, MPI_DOUBLE, &contiguos_boundary);
MPI_Type_vector(solver->limax, 1, solver->ljmax + 2, MPI_DOUBLE, &vector_boundary);
MPI_Type_commit(&contiguos_boundary);
MPI_Type_commit(&vector_boundary);
solver->types[0] = contiguos_boundary;
solver->types[1] = contiguos_boundary;
solver->types[2] = vector_boundary;
solver->types[3] = vector_boundary;
solver->snd_displacements[0] = ((solver->ljmax + 2) + 1) * sizeof(double); solver->snd_displacements[0] = ((solver->ljmax + 2) + 1) * sizeof(double);
solver->snd_displacements[1] = ((solver->ljmax + 2) * (solver->limax) + 1) * solver->snd_displacements[1] = ((solver->ljmax + 2) * (solver->limax) + 1) *
sizeof(double); sizeof(double);
@ -85,53 +119,39 @@ void initSolver(Solver* solver, Parameter* params, int problem)
solver->rcv_displacements[3] = ((solver->ljmax + 2) + solver->ljmax + 1) * solver->rcv_displacements[3] = ((solver->ljmax + 2) + solver->ljmax + 1) *
sizeof(double); sizeof(double);
MPI_Datatype contiguos_boundary;
MPI_Datatype vector_boundary;
MPI_Type_contiguous(solver->ljmax, MPI_DOUBLE, &contiguos_boundary);
MPI_Type_vector(solver->limax, 1, solver->ljmax + 2, MPI_DOUBLE, &vector_boundary);
MPI_Type_commit(&contiguos_boundary);
MPI_Type_commit(&vector_boundary);
solver->types[0] = contiguos_boundary;
solver->types[1] = contiguos_boundary;
solver->types[2] = vector_boundary;
solver->types[3] = vector_boundary;
size_t bytesize = (solver->limax + 2) * (solver->ljmax + 2) * sizeof(double); size_t bytesize = (solver->limax + 2) * (solver->ljmax + 2) * sizeof(double);
solver->p = malloc(bytesize); solver->p = allocate(64, bytesize);
solver->rhs = malloc(bytesize); solver->rhs = allocate(64, bytesize);
double dx = solver->dx; double dx = solver->dx;
double dy = solver->dy; double dy = solver->dy;
double* p = solver->p; double* p = solver->p;
double* rhs = solver->rhs; double* rhs = solver->rhs;
int i_start = get_dim_start(solver->coords[0], solver->dims[0], solver->imax); solver->i_start = get_dim_start(solver->coords[0], solver->dims[0], solver->imax);
int j_start = get_dim_start(solver->coords[1], solver->dims[1], solver->jmax); solver->j_start = get_dim_start(solver->coords[1], solver->dims[1], solver->jmax);
#ifdef DEBUG
printf("Prc %d; {%d,%d}, {%d,%d}, {%d,%d}\n", printf("Prc %d; {%d,%d}, {%d,%d}, {%d,%d}\n",
rank, rank,
solver->coords[0], solver->coords[0],
solver->coords[1], solver->coords[1],
solver->limax, solver->limax,
solver->ljmax, solver->ljmax,
i_start, solver->i_start,
j_start); solver->j_start);
fflush(stdout); fflush(stdout);
// for (int i = 0; i < solver->limax + 2; i++) { #endif /* ifdef DEBUG */
// 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 i = 0; i < solver->limax + 2; i++) {
for (int j = 0; j < solver->ljmax + 2; j++) { for (int j = 0; j < solver->ljmax + 2; j++) {
P(i, j) = sin(2.0 * PI * (i_start + i) * dx * 2.0) + P(i, j) = sin(2.0 * PI * (solver->i_start + i) * dx * 2.0) +
sin(2.0 * PI * (j_start + j) * dy * 2.0); sin(2.0 * PI * (solver->j_start + j) * dy * 2.0);
} }
} }
if (problem == 2) { if (problem == 2) {
for (int i = 0; i < solver->limax + 2; i++) { for (int i = 0; i < solver->limax + 2; i++) {
for (int j = 0; j < solver->ljmax + 2; j++) { for (int j = 0; j < solver->ljmax + 2; j++) {
RHS(i, j) = sin(2.0 * PI * (i_start + i) * dx); RHS(i, j) = sin(2.0 * PI * (solver->i_start + i) * dx);
} }
} }
} else { } else {
@ -160,25 +180,25 @@ void solveRB(Solver* solver)
int it = 0; int it = 0;
double res = 1.0; double res = 1.0;
int pass, jsw, isw; int pass, jsw, isw;
while ((res >= epssq) && (it < itermax)) { while ((res >= epssq) && (it < itermax)) {
res = 0.0; res = 0.0;
jsw = 1; jsw = solver->j_start % 2 ? 1 : 2;
int count[] = { 1, 1, 1, 1 };
MPI_Neighbor_alltoallw(solver->p, MPI_Neighbor_alltoallw(solver->p,
count, (int[4]) { 1, 1, 1, 1 },
solver->snd_displacements, solver->snd_displacements,
solver->types, solver->types,
solver->p, solver->p,
count, (int[4]) { 1, 1, 1, 1 },
solver->rcv_displacements, solver->rcv_displacements,
solver->types, solver->types,
solver->cart_comm); solver->cart_comm);
for (pass = 0; pass < 2; pass++) {
for (pass = 0; pass < 2; pass++) { // why ??
isw = jsw; isw = jsw;
for (int i = 1; i < solver->limax + 1; i++) {
for (int j = 1; j < solver->ljmax + 1; j++) { for (int j = isw; j < solver->ljmax + 1; j += 2) {
for (int i = isw; i < solver->limax + 1; i += 2) {
double r = RHS(i, j) - double r = RHS(i, j) -
((P(i + 1, j) - 2.0 * P(i, j) + P(i - 1, j)) * idx2 + ((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 + 1) - 2.0 * P(i, j) + P(i, j - 1)) * idy2);
@ -190,7 +210,9 @@ void solveRB(Solver* solver)
} }
jsw = 3 - jsw; jsw = 3 - jsw;
} }
MPI_Allreduce(MPI_IN_PLACE, &res, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); MPI_Allreduce(MPI_IN_PLACE, &res, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
if (solver->boundary & BOTTOM) if (solver->boundary & BOTTOM)
for (int j = 1; j < solver->ljmax + 1; j++) for (int j = 1; j < solver->ljmax + 1; j++)
P(0, j) = P(1, j); P(0, j) = P(1, j);
@ -205,31 +227,40 @@ void solveRB(Solver* solver)
P(i, solver->ljmax + 1) = P(i, solver->ljmax); P(i, solver->ljmax + 1) = P(i, solver->ljmax);
res = res / (double)(imax * jmax); res = res / (double)(imax * jmax);
#ifdef DEBUG #ifdef DEBUG
printf("%d Residuum: %e\n", it, res); printf("%d Residuum: %e\n", it, res);
#endif #endif
it++; it++;
} }
if(solver->rank==0) printf("Solver took %d iterations to reach %f\n", it, sqrt(res)); if (solver->rank == 0)
printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
} }
void writeResult(Solver* solver, char* filename) void writeResult(Solver* solver, char* filename)
{ {
int imax = solver->imax; int imax = solver->imax;
int jmax = solver->jmax; int jmax = solver->jmax;
double* p = solver->p; double* p = solver->p;
int* rcv_count = NULL; int* rcv_count = NULL;
int* displacements = NULL; int* displacements = NULL;
double* p_total = NULL; double* p_total = NULL;
double* send_buffer = NULL;
int local_size = solver->limax * solver->ljmax;
if (solver->rank == 0) { if (solver->rank == 0) {
rcv_count = malloc(solver->size * sizeof(int));
displacements = calloc(solver->size, sizeof(int)); rcv_count = allocate(64, solver->size * sizeof(int));
p_total = malloc(sizeof(double) * solver->imax * solver->jmax); displacements = allocate(64, solver->size * sizeof(int));
double* p_total_ordered = malloc(sizeof(double) * solver->imax * solver->jmax); p_total = allocate(64, sizeof(double) * solver->imax * solver->jmax);
int local_size = solver->limax * solver->ljmax; double* p_total_ordered = allocate(64,
sizeof(double) * solver->imax * solver->jmax);
MPI_Gather(&local_size, 1, MPI_INT, rcv_count, 1, MPI_INT, 0, MPI_COMM_WORLD); MPI_Gather(&local_size, 1, MPI_INT, rcv_count, 1, MPI_INT, 0, MPI_COMM_WORLD);
displacements[0] = 0; displacements[0] = 0;
for (int i = 0; i < solver->limax; i++) for (int i = 0; i < solver->limax; i++)
for (int j = 0; j < solver->ljmax; j++) for (int j = 0; j < solver->ljmax; j++)
p_total_ordered[i * (solver->ljmax) + j] = p_total_ordered[i * (solver->ljmax) + j] =
@ -237,14 +268,7 @@ void writeResult(Solver* solver, char* filename)
for (int i = 1; i < solver->size; i++) for (int i = 1; i < solver->size; i++)
displacements[i] = displacements[i - 1] + rcv_count[i - 1]; displacements[i] = displacements[i - 1] + rcv_count[i - 1];
printf("{ ");
for (int i = 0; i < solver->size; i++)
printf("%d ", rcv_count[i]);
printf("}\n");
printf("{ ");
for (int i = 0; i < solver->size; i++)
printf("%d ", displacements[i]);
printf("}\n");
MPI_Gatherv(p_total_ordered, MPI_Gatherv(p_total_ordered,
local_size, local_size,
MPI_DOUBLE, MPI_DOUBLE,
@ -254,43 +278,37 @@ void writeResult(Solver* solver, char* filename)
MPI_DOUBLE, MPI_DOUBLE,
0, 0,
MPI_COMM_WORLD); MPI_COMM_WORLD);
for (int process = 0; process < solver->size; process++) { for (int process = 0; process < solver->size; process++) {
int coords[2] = { 0 }; int coords[2] = { 0 };
MPI_Cart_coords(solver->cart_comm, process, 2, coords); MPI_Cart_coords(solver->cart_comm, process, 2, coords);
int i_max_local = distribute_dim(coords[0], solver->dims[0], solver->imax); int i_max_local = distribute_dim(coords[0], solver->dims[0], solver->imax);
int j_max_local = distribute_dim(coords[1], solver->dims[1], solver->jmax); int j_max_local = distribute_dim(coords[1], solver->dims[1], solver->jmax);
int i_start_local = get_dim_start(coords[0], solver->dims[0], solver->imax); int i_start_local = get_dim_start(coords[0], solver->dims[0], solver->imax);
int j_start_local = get_dim_start(coords[1], solver->dims[1], solver->jmax); int j_start_local = get_dim_start(coords[1], solver->dims[1], solver->jmax);
for (int i = 0; i < i_max_local; i++) for (int i = 0; i < i_max_local; i++)
for (int j = 0; j < j_max_local; j++) for (int j = 0; j < j_max_local; j++)
p_total_ordered[(i + i_start_local) * (solver->jmax) + p_total_ordered[(i + i_start_local) * (solver->jmax) +
(j + j_start_local)] = (j + j_start_local)] =
p_total[displacements[process] + i * (j_max_local) + j]; p_total[displacements[process] + i * (j_max_local) + j];
} }
FILE* fp;
fp = fopen(filename, "w");
if (fp == NULL) { write_matrix_to_file(filename, p_total_ordered, solver->imax, solver->jmax);
printf("Error!\n");
exit(EXIT_FAILURE);
}
for (int j = 0; j < jmax; j++) {
for (int i = 0; i < imax; i++) {
fprintf(fp, "%f ", p_total_ordered[i * (solver->jmax) + j]);
}
fprintf(fp, "\n");
}
fclose(fp);
} else { } else {
int local_size = solver->limax * solver->ljmax;
double* send_buffer = malloc(local_size * sizeof(double)); send_buffer = allocate(64, local_size * sizeof(double));
for (int i = 0; i < solver->limax; i++) for (int i = 0; i < solver->limax; i++)
for (int j = 0; j < solver->ljmax; j++) for (int j = 0; j < solver->ljmax; j++)
send_buffer[i * (solver->ljmax) + j] = send_buffer[i * (solver->ljmax) + j] =
solver->p[(i + 1) * (solver->ljmax + 2) + (j + 1)]; solver->p[(i + 1) * (solver->ljmax + 2) + (j + 1)];
MPI_Gather(&local_size, 1, MPI_INT, rcv_count, 1, MPI_INT, 0, MPI_COMM_WORLD); MPI_Gather(&local_size, 1, MPI_INT, rcv_count, 1, MPI_INT, 0, MPI_COMM_WORLD);
MPI_Gatherv(send_buffer, MPI_Gatherv(send_buffer,
local_size, local_size,
MPI_DOUBLE, MPI_DOUBLE,

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

@ -23,6 +23,7 @@ typedef struct {
double dx, dy; double dx, dy;
int imax, jmax; int imax, jmax;
int limax, ljmax; int limax, ljmax;
int i_start, j_start;
double *p, *rhs; double *p, *rhs;
double eps, omega, rho; double eps, omega, rho;
int itermax; int itermax;