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

checkpoint: Working Neigh com and write result

Browse Source
This commit is contained in:
indiano 2025-05-03 23:52:52 +02:00
parent 5148069114
commit 8a218ad681
3 changed files with 139 additions and 68 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -29,17 +29,16 @@ int main(int argc, char** argv)
printf("Usage: %s <configFile>\n", argv[0]);
exit(EXIT_SUCCESS);
}
readParameter(&params, argv[1]);
initSolver(&solver, &params, 2);
// writeResult(&solver, "p-0.dat");
// startTime = getTimeStamp();
// solveRB(&solver);
// endTime = getTimeStamp();
writeResult(&solver, "p.dat");
// printf(" %.2fs\n", endTime - startTime);
// writeResult(&solver, "p.dat");
if (rank == 0) {
}
MPI_Finalize();
return EXIT_SUCCESS;

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

@ -9,11 +9,13 @@
#include "allocate.h"
#include "parameter.h"
#include "solver.h"
#include <alloca.h>
#include <mpi.h>
#include <stdlib.h>
#define PI 3.14159265358979323846
#define P(i, j) p[(j) * (solver->limax + 2) + (i)]
#define RHS(i, j) rhs[(j) * (solver->limax + 2) + (i)]
#define P(i, j) p[(i) * (solver->ljmax + 2) + (j)]
#define RHS(i, j) rhs[(i) * (solver->ljmax + 2) + (j)]
int distribute_dim(int dim_rank, int dim_comm_size, int dim_size)
{
@ -48,6 +50,8 @@ void initSolver(Solver* solver, Parameter* params, int problem)
solver->dims[1] = 0;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Comm_size(MPI_COMM_WORLD, &size);
solver->rank = rank;
solver->size = size;
MPI_Dims_create(size, 2, solver->dims);
MPI_Cart_create(MPI_COMM_WORLD,
2,
@ -59,51 +63,41 @@ void initSolver(Solver* solver, Parameter* params, int problem)
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->snd_displacements[0] = ((solver->limax + 2) + 1) * sizeof(double);
solver->snd_displacements[1] = ((solver->limax + 2) * (solver->ljmax) + 1) *
solver->snd_displacements[0] = ((solver->ljmax + 2) + 1) * sizeof(double);
solver->snd_displacements[1] = ((solver->ljmax + 2) * (solver->limax) + 1) *
sizeof(double);
solver->snd_displacements[2] = ((solver->limax + 2) + 1) * sizeof(double);
solver->snd_displacements[3] = ((solver->limax + 2) + solver->limax) * sizeof(double);
solver->snd_displacements[2] = ((solver->ljmax + 2) + 1) * sizeof(double);
solver->snd_displacements[3] = ((solver->ljmax + 2) + solver->ljmax) * sizeof(double);
solver->rcv_displacements[0] = (1) * sizeof(double);
solver->rcv_displacements[1] = ((solver->limax + 2) * (solver->ljmax + 1) + 1) *
solver->rcv_displacements[1] = ((solver->ljmax + 2) * (solver->limax + 1) + 1) *
sizeof(double);
solver->rcv_displacements[2] = ((solver->limax + 2)) * sizeof(double);
solver->rcv_displacements[3] = ((solver->limax + 2) + solver->limax + 1) *
solver->rcv_displacements[2] = ((solver->ljmax + 2)) * sizeof(double);
solver->rcv_displacements[3] = ((solver->ljmax + 2) + solver->ljmax + 1) *
sizeof(double);
MPI_Datatype contiguos_boundary;
MPI_Datatype vector_boundary;
MPI_Datatype internal;
MPI_Type_contiguous(solver->limax, MPI_DOUBLE, &contiguos_boundary);
MPI_Type_vector(solver->ljmax, 1, solver->limax + 2, MPI_DOUBLE, &vector_boundary);
MPI_Type_create_subarray(2,
(int[2]){solver->limax+2,solver->ljmax+2},
(int[2]){solver->limax,solver->ljmax},
(int[2]){1,1},
MPI_ORDER_C,
MPI_DOUBLE,
&internal);
MPI_Type_commit(&internal);
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->types[4] = internal;
size_t bytesize = (solver->limax + 2) * (solver->ljmax + 2) * sizeof(double);
solver->p = allocate(64, bytesize);
solver->rhs = allocate(64, bytesize);
solver->p = malloc(bytesize);
solver->rhs = malloc(bytesize);
double dx = solver->dx;
double dy = solver->dy;
double* p = solver->p;
double* rhs = solver->rhs;
int 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);
printf("Prc %d; {%d,%d}, {%d,%d}, {%d,%d}\n",
rank,
solver->coords[0],
@ -113,23 +107,27 @@ void initSolver(Solver* solver, Parameter* params, int problem)
i_start,
j_start);
fflush(stdout);
for (int j = 0; j < solver->ljmax + 2; j++) {
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++) {
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) = sin(2.0 * PI * (i_start + i) * dx * 2.0) +
sin(2.0 * PI * (j_start + j) * dy * 2.0);
}
}
if (problem == 2) {
for (int j = 0; j < solver->ljmax + 2; j++) {
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++) {
RHS(i, j) = sin(2.0 * PI * (i_start + i) * dx);
}
}
} else {
for (int j = 0; j < solver->ljmax + 2; j++) {
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++) {
RHS(i, j) = 0.0;
}
}
@ -153,16 +151,24 @@ void solveRB(Solver* solver)
int it = 0;
double res = 1.0;
int pass, jsw, isw;
while ((res >= epssq) && (it < itermax)) {
res = 0.0;
jsw = 1;
res = 0.0;
jsw = 1;
int count[] = { 1, 1, 1, 1 };
MPI_Neighbor_alltoallw(solver->p,
count,
solver->snd_displacements,
solver->types,
solver->p,
count,
solver->rcv_displacements,
solver->types,
solver->cart_comm);
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) {
for (int j = 1; j < solver->ljmax + 1; j++) {
for (int i = isw; i < solver->limax + 1; i += 2) {
double r = RHS(i, j) -
((P(i + 1, j) - 2.0 * P(i, j) + P(i - 1, j)) * idx2 +
@ -175,16 +181,16 @@ 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);
}
//
// 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
@ -199,24 +205,89 @@ void solveRB(Solver* solver)
void writeResult(Solver* solver, char* filename)
{
int imax = solver->imax;
int jmax = solver->jmax;
double* p = solver->p;
int imax = solver->imax;
int jmax = solver->jmax;
double* p = solver->p;
int* rcv_count = NULL;
int* displacements = NULL;
double* p_total = NULL;
if (solver->rank == 0) {
rcv_count = malloc(solver->size * sizeof(int));
displacements = calloc(solver->size, sizeof(int));
p_total = malloc(sizeof(double) * solver->imax * solver->jmax);
double* p_total_ordered = malloc(sizeof(double) * solver->imax * solver->jmax);
int local_size = solver->limax * solver->ljmax;
MPI_Gather(&local_size, 1, MPI_INT, rcv_count, 1, MPI_INT, 0, MPI_COMM_WORLD);
displacements[0] = 0;
for (int i = 0; i < solver->limax; i++)
for (int j = 0; j < solver->ljmax; j++)
p_total_ordered[i * (solver->ljmax) + j] =
solver->p[(i + 1) * (solver->ljmax + 2) + (j + 1)];
FILE* fp;
fp = fopen(filename, "w");
if (fp == NULL) {
printf("Error!\n");
exit(EXIT_FAILURE);
}
for (int j = 0; j < jmax + 2; j++) {
for (int i = 0; i < imax + 2; i++) {
fprintf(fp, "%f ", P(i, j));
for (int i = 1; i < solver->size; i++)
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,
local_size,
MPI_DOUBLE,
p_total,
rcv_count,
displacements,
MPI_DOUBLE,
0,
MPI_COMM_WORLD);
for (int process = 0; process < solver->size; process++) {
int coords[2] = { 0 };
MPI_Cart_coords(solver->cart_comm, process, 2, coords);
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 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);
for (int i = 0; i < i_max_local; i++)
for (int j = 0; j < j_max_local; j++)
p_total_ordered[(i + i_start_local) * (solver->jmax) +
(j + j_start_local)] =
p_total[displacements[process] + i * (j_max_local) + j];
}
fprintf(fp, "\n");
}
FILE* fp;
fp = fopen(filename, "w");
fclose(fp);
if (fp == NULL) {
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 {
int local_size = solver->limax * solver->ljmax;
double* send_buffer = malloc(local_size * sizeof(double));
for (int i = 0; i < solver->limax; i++)
for (int j = 0; j < solver->ljmax; j++)
send_buffer[i * (solver->ljmax) + j] =
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_Gatherv(send_buffer,
local_size,
MPI_DOUBLE,
p_total,
rcv_count,
displacements,
MPI_DOUBLE,
0,
MPI_COMM_WORLD);
}
}

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

@ -15,11 +15,12 @@ typedef struct {
double eps, omega, rho;
int itermax;
MPI_Comm cart_comm;
MPI_Datatype types[5];
MPI_Datatype types[4];
MPI_Aint rcv_displacements[4];
MPI_Aint snd_displacements[4];
int coords[2];
int dims[2];
int rank, size;
} Solver;
extern void initSolver(Solver*, Parameter*, int problem);