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

Compare commits

base: ErikFabrizzi:356c4fe4fac1178908ef4791c47b1e64cc50f438
Branches Tags
ErikFabrizzi:main
ErikFabrizzi:poisson/2D-omp-alternative
ErikFabrizzi:poisson/2D-mpi-sq
ErikFabrizzi:poisson/2D-omp
moebiusband:main
..
compare: ErikFabrizzi:92e9ed764fc9fa30545986349e3d085b3667d911
Branches Tags
ErikFabrizzi:poisson/2D-omp-alternative
ErikFabrizzi:main
ErikFabrizzi:poisson/2D-mpi-sq
ErikFabrizzi:poisson/2D-omp
moebiusband:main

No commits in common. "356c4fe4fac1178908ef4791c47b1e64cc50f438" and "92e9ed764fc9fa30545986349e3d085b3667d911" have entirely different histories.
356c4fe4fa ... 92e9ed764f

3 changed files with 355 additions and 352 deletions
Show Stats Expand all files Collapse all files

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

@ -35,7 +35,7 @@ int main(int argc, char** argv)
} }
initSolver(&solver, &params, 2); initSolver(&solver, &params, 2);
solveRB(&solver); solve(&solver);
getResult(&solver); getResult(&solver);
MPI_Finalize(); MPI_Finalize();

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

@ -17,16 +17,17 @@
#define P(i, j) p[(j) * (imax + 2) + (i)] #define P(i, j) p[(j) * (imax + 2) + (i)]
#define RHS(i, j) rhs[(j) * (imax + 2) + (i)] #define RHS(i, j) rhs[(j) * (imax + 2) + (i)]
static int sizeOfRank(int rank, int size, int N) { 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) { static void print(Solver* solver)
{
double* p = solver->p; double* p = solver->p;
int imax = solver->imax; int imax = solver->imax;
printf( printf("### RANK %d #######################################################\n",
"### RANK %d #######################################################\n",
solver->rank); solver->rank);
for (int j = 0; j < solver->jmaxLocal + 2; j++) { for (int j = 0; j < solver->jmaxLocal + 2; j++) {
printf("%02d: ", j); printf("%02d: ", j);
@ -38,9 +39,12 @@ static void print(Solver *solver) {
fflush(stdout); fflush(stdout);
} }
static void exchange(Solver *solver) { 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 */ /* exchange ghost cells with top neighbor */
if (solver->rank + 1 < solver->size) { if (solver->rank + 1 < solver->size) {
@ -48,10 +52,8 @@ static void exchange(Solver *solver) {
double* src = solver->p + (solver->jmaxLocal) * (solver->imax + 2) + 1; double* src = solver->p + (solver->jmaxLocal) * (solver->imax + 2) + 1;
double* dst = solver->p + (solver->jmaxLocal + 1) * (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, MPI_Isend(src, solver->imax, MPI_DOUBLE, top, 1, MPI_COMM_WORLD, &requests[0]);
&requests[0]); MPI_Irecv(dst, solver->imax, MPI_DOUBLE, top, 2, MPI_COMM_WORLD, &requests[1]);
MPI_Irecv(dst, solver->imax, MPI_DOUBLE, top, 2, MPI_COMM_WORLD,
&requests[1]);
} }
/* exchange ghost cells with bottom neighbor */ /* exchange ghost cells with bottom neighbor */
@ -60,58 +62,61 @@ static void exchange(Solver *solver) {
double* src = solver->p + (solver->imax + 2) + 1; double* src = solver->p + (solver->imax + 2) + 1;
double* dst = solver->p + 1; double* dst = solver->p + 1;
MPI_Isend(src, solver->imax, MPI_DOUBLE, bottom, 2, MPI_COMM_WORLD, MPI_Isend(src, solver->imax, MPI_DOUBLE, bottom, 2, MPI_COMM_WORLD, &requests[2]);
&requests[2]); MPI_Irecv(dst, solver->imax, MPI_DOUBLE, bottom, 1, MPI_COMM_WORLD, &requests[3]);
MPI_Irecv(dst, solver->imax, MPI_DOUBLE, bottom, 1, MPI_COMM_WORLD,
&requests[3]);
} }
MPI_Waitall(4, requests, MPI_STATUSES_IGNORE); MPI_Waitall(4, requests, MPI_STATUSES_IGNORE);
} }
void getResult(Solver *solver) { void getResult(Solver* solver)
{
double* Pall = NULL; double* Pall = NULL;
int *rcvCounts, *displs; int *rcvCounts, *displs;
if (solver->rank == 0) { if (solver->rank == 0) {
Pall = Pall = allocate(64, (solver->imax + 2) * (solver->jmax + 2) * sizeof(double));
allocate(64, (solver->imax + 2) * (solver->jmax + 2) * sizeof(double));
rcvCounts = (int*)malloc(solver->size * sizeof(int)); rcvCounts = (int*)malloc(solver->size * sizeof(int));
displs = (int*)malloc(solver->size * sizeof(int)); displs = (int*)malloc(solver->size * sizeof(int));
rcvCounts[0] = (solver->jmaxLocal + 1) * (solver->imax + 2); rcvCounts[0] = solver->jmaxLocal * (solver->imax + 2);
displs[0] = 0; displs[0] = 0;
int cursor = rcvCounts[0]; int cursor = rcvCounts[0];
for (int i = 1; i < solver->size; i++) { for (int i = 1; i < solver->size; i++) {
rcvCounts[i] = rcvCounts[i] = sizeOfRank(i, solver->size, solver->jmax) * (solver->imax + 2);
sizeOfRank(i, solver->size, solver->jmax) * (solver->imax + 2);
displs[i] = cursor; displs[i] = cursor;
cursor += rcvCounts[i]; cursor += rcvCounts[i];
} }
rcvCounts[solver->size - 1] =
(sizeOfRank(solver->size - 1, solver->size, solver->jmax)+1) *
(solver->imax + 2);
} }
int cnt = solver->jmaxLocal * (solver->imax + 2); int cnt = solver->jmaxLocal * (solver->imax + 2);
double* sendbuffer = solver->p + (solver->imax + 2); double* sendbuffer = solver->p + (solver->imax + 2);
if(solver->rank == 0 || solver->rank == solver->size-1) cnt= (solver->jmaxLocal+1)*(solver->imax+2); MPI_Gatherv(sendbuffer,
if(solver->rank == 0) sendbuffer = solver->p; cnt,
MPI_Gatherv(sendbuffer, cnt, MPI_DOUBLE, Pall, rcvCounts, displs, MPI_DOUBLE, MPI_DOUBLE,
0, MPI_COMM_WORLD); Pall,
rcvCounts,
displs,
MPI_DOUBLE,
0,
MPI_COMM_WORLD);
if (solver->rank == 0) { if (solver->rank == 0) {
writeResult(solver, Pall, "p.dat"); 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_rank(MPI_COMM_WORLD, &(solver->rank));
MPI_Comm_size(MPI_COMM_WORLD, &(solver->size)); MPI_Comm_size(MPI_COMM_WORLD, &(solver->size));
solver->imax = params->imax; solver->imax = params->imax;
solver->jmax = params->jmax; solver->jmax = params->jmax;
solver->jmaxLocal = sizeOfRank(solver->rank, solver->size, solver->jmax); solver->jmaxLocal = sizeOfRank(solver->rank, solver->size, solver->jmax);
printf("RANK %d: imaxLocal : %d, jmaxLocal : %d\n", solver->rank, printf("RANK %d: imaxLocal : %d, jmaxLocal : %d\n",
solver->imax, solver->jmaxLocal); solver->rank,
solver->imax,
solver->jmaxLocal);
solver->dx = params->xlength / params->imax; solver->dx = params->xlength / params->imax;
solver->dy = params->ylength / params->jmax; solver->dy = params->ylength / params->jmax;
@ -153,7 +158,8 @@ void initSolver(Solver *solver, Parameter *params, int problem) {
} }
} }
void debug(Solver *solver) { void debug(Solver* solver)
{
int imax = solver->imax; int imax = solver->imax;
int rank = solver->rank; int rank = solver->rank;
double* p = solver->p; double* p = solver->p;
@ -172,11 +178,9 @@ void debug(Solver *solver) {
/* } */ /* } */
/* if ( rank == 0 ) { */ /* if ( rank == 0 ) { */
/* printf("##########################################################\n"); /* printf("##########################################################\n"); */
*/
/* printf("## Exchange ghost layers\n"); */ /* printf("## Exchange ghost layers\n"); */
/* printf("##########################################################\n"); /* printf("##########################################################\n"); */
*/
/* } */ /* } */
/* exchange(solver); */ /* exchange(solver); */
@ -188,7 +192,8 @@ void debug(Solver *solver) {
} }
} }
int solve(Solver *solver) { int solve(Solver* solver)
{
double r; double r;
int it = 0; int it = 0;
double res, res1; double res, res1;
@ -264,7 +269,8 @@ int solve(Solver *solver) {
} }
} }
int solveRB(Solver *solver) { int solveRB(Solver* solver)
{
double r; double r;
int it = 0; int it = 0;
double res, res1; double res, res1;
@ -299,8 +305,8 @@ int solveRB(Solver *solver) {
for (int j = 1; j < jmaxLocal + 1; j++) { for (int j = 1; j < jmaxLocal + 1; j++) {
for (int i = isw; i < imax + 1; i += 2) { for (int i = isw; i < imax + 1; i += 2) {
double r = double r = RHS(i, j) -
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);
P(i, j) -= (factor * r); P(i, j) -= (factor * r);
@ -310,16 +316,11 @@ int solveRB(Solver *solver) {
} }
jsw = 3 - jsw; jsw = 3 - jsw;
} }
if(solver->rank==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; i < imax + 1; i++) {
for(int i = 1; i < imax +1; i++)P(i, jmaxLocal + 1) = P(i, jmaxLocal); P(i, 0) = P(i, 1);
P(i, jmaxLocal + 1) = P(i, jmaxLocal);
// 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++) { for (int j = 1; j < jmaxLocal + 1; j++) {
P(0, j) = P(1, j); P(0, j) = P(1, j);
@ -344,7 +345,8 @@ int solveRB(Solver *solver) {
} }
} }
int solveRBA(Solver *solver) { int solveRBA(Solver* solver)
{
double r; double r;
int it = 0; int it = 0;
double res; double res;
@ -380,8 +382,8 @@ int solveRBA(Solver *solver) {
for (int j = 1; j < jmaxLocal + 1; j++) { for (int j = 1; j < jmaxLocal + 1; j++) {
for (int i = isw; i < imax + 1; i += 2) { for (int i = isw; i < imax + 1; i += 2) {
double r = double r = RHS(i, j) -
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);
P(i, j) -= (omega * factor * r); P(i, j) -= (omega * factor * r);
@ -415,7 +417,8 @@ int solveRBA(Solver *solver) {
printf("Solver took %d iterations to reach %f\n", it, sqrt(res)); printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
} }
void writeResult(Solver *solver, double *m, char *filename) { void writeResult(Solver* solver, double* m, 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;

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

@ -25,5 +25,5 @@ extern void debug(Solver*);
extern void initSolver(Solver*, Parameter*, int problem); extern void initSolver(Solver*, Parameter*, int problem);
extern void getResult(Solver*); extern void getResult(Solver*);
extern void writeResult(Solver*, double*, char*); extern void writeResult(Solver*, double*, char*);
extern int solveRB(Solver*); extern int solve(Solver*);
#endif #endif