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();

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

@ -13,427 +13,430 @@
#include "parameter.h" #include "parameter.h"
#include "solver.h" #include "solver.h"
#define PI 3.14159265358979323846 #define PI 3.14159265358979323846
#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; {
int imax = solver->imax; double* p = solver->p;
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); 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("%12.8f ", P(i, j)); }
printf("\n");
} }
printf("\n"); 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) {
int top = solver->rank + 1; int top = solver->rank + 1;
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 */
if (solver->rank > 0) { if (solver->rank > 0) {
int bottom = solver->rank - 1; int bottom = solver->rank - 1;
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; {
int *rcvCounts, *displs; double* Pall = NULL;
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 * (solver->imax + 2);
rcvCounts[0] = (solver->jmaxLocal + 1) * (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);
double *sendbuffer = solver->p + (solver->imax + 2);
if(solver->rank == 0 || solver->rank == solver->size-1) cnt= (solver->jmaxLocal+1)*(solver->imax+2);
if(solver->rank == 0) sendbuffer = solver->p;
MPI_Gatherv(sendbuffer, cnt, MPI_DOUBLE, Pall, rcvCounts, displs, MPI_DOUBLE,
0, MPI_COMM_WORLD);
if (solver->rank == 0) { int cnt = solver->jmaxLocal * (solver->imax + 2);
writeResult(solver, Pall, "p.dat"); 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) {
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)); MPI_Comm_rank(MPI_COMM_WORLD, &(solver->rank));
solver->imax = params->imax; MPI_Comm_size(MPI_COMM_WORLD, &(solver->size));
solver->jmax = params->jmax; solver->imax = params->imax;
solver->jmaxLocal = sizeOfRank(solver->rank, solver->size, solver->jmax); solver->jmax = params->jmax;
printf("RANK %d: imaxLocal : %d, jmaxLocal : %d\n", solver->rank, solver->jmaxLocal = sizeOfRank(solver->rank, solver->size, solver->jmax);
solver->imax, solver->jmaxLocal); printf("RANK %d: imaxLocal : %d, jmaxLocal : %d\n",
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;
solver->ys = solver->rank * solver->jmaxLocal * solver->dy; solver->ys = solver->rank * solver->jmaxLocal * solver->dy;
solver->eps = params->eps; solver->eps = params->eps;
solver->omega = params->omg; solver->omega = params->omg;
solver->itermax = params->itermax; solver->itermax = params->itermax;
int imax = solver->imax; int imax = solver->imax;
int jmax = solver->jmax; int jmax = solver->jmax;
int jmaxLocal = solver->jmaxLocal; int jmaxLocal = solver->jmaxLocal;
solver->p = allocate(64, (imax + 2) * (jmaxLocal + 2) * sizeof(double)); solver->p = allocate(64, (imax + 2) * (jmaxLocal + 2) * sizeof(double));
solver->rhs = allocate(64, (imax + 2) * (jmax + 2) * sizeof(double)); solver->rhs = allocate(64, (imax + 2) * (jmax + 2) * sizeof(double));
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;
for (int j = 0; j < jmaxLocal + 2; j++) { for (int j = 0; j < jmaxLocal + 2; j++) {
double y = solver->ys + j * dy; double y = solver->ys + j * dy;
for (int i = 0; i < imax + 2; i++) { for (int i = 0; i < imax + 2; i++) {
P(i, j) = sin(4.0 * PI * i * dx) + sin(4.0 * PI * y); P(i, j) = sin(4.0 * PI * i * dx) + sin(4.0 * PI * y);
}
} }
}
if (problem == 2) { if (problem == 2) {
for (int j = 0; j < jmax + 2; j++) { for (int j = 0; j < jmax + 2; j++) {
for (int i = 0; i < imax + 2; i++) { for (int i = 0; i < imax + 2; i++) {
RHS(i, j) = sin(2.0 * PI * i * dx); RHS(i, j) = sin(2.0 * PI * i * dx);
} }
}
} else {
for (int j = 0; j < jmax + 2; j++) {
for (int i = 0; i < imax + 2; i++) {
RHS(i, j) = 0.0;
}
}
} }
} else {
for (int j = 0; j < jmax + 2; j++) {
for (int i = 0; i < imax + 2; i++) {
RHS(i, j) = 0.0;
}
}
}
} }
void debug(Solver *solver) { void debug(Solver* solver)
int imax = solver->imax; {
int rank = solver->rank; int imax = solver->imax;
double *p = solver->p; int rank = solver->rank;
double* p = solver->p;
/* for( int j=0; j < solver->jmaxLocal+2; j++ ) { */ /* for( int j=0; j < solver->jmaxLocal+2; j++ ) { */
/* for( int i=0; i < solver->imax+2; i++ ) { */ /* for( int i=0; i < solver->imax+2; i++ ) { */
/* P(i, j) = (double) rank; */ /* P(i, j) = (double) rank; */
/* } */ /* } */
/* } */ /* } */
/* for ( int i=0; i < solver->size; i++) { */ /* for ( int i=0; i < solver->size; i++) { */
/* if ( i == rank ) { */ /* if ( i == rank ) { */
/* print(solver); */ /* print(solver); */
/* } */ /* } */
/* MPI_Barrier(MPI_COMM_WORLD); */ /* MPI_Barrier(MPI_COMM_WORLD); */
/* } */ /* } */
/* 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); */
for (int i = 0; i < solver->size; i++) { for (int i = 0; i < solver->size; i++) {
if (i == rank) { if (i == rank) {
print(solver); print(solver);
}
MPI_Barrier(MPI_COMM_WORLD);
} }
MPI_Barrier(MPI_COMM_WORLD);
}
} }
int solve(Solver *solver) { int solve(Solver* solver)
double r; {
int it = 0; double r;
double res, res1; int it = 0;
double res, res1;
int imax = solver->imax; int imax = solver->imax;
int jmax = solver->jmax; int jmax = solver->jmax;
int jmaxLocal = solver->jmaxLocal; int jmaxLocal = solver->jmaxLocal;
double eps = solver->eps; double eps = solver->eps;
double omega = solver->omega; double omega = solver->omega;
int itermax = solver->itermax; int itermax = solver->itermax;
double dx2 = solver->dx * solver->dx; double dx2 = solver->dx * solver->dx;
double dy2 = solver->dy * solver->dy; double dy2 = solver->dy * solver->dy;
double idx2 = 1.0 / dx2; double idx2 = 1.0 / dx2;
double idy2 = 1.0 / dy2; double idy2 = 1.0 / dy2;
double factor = omega * 0.5 * (dx2 * dy2) / (dx2 + dy2); double factor = omega * 0.5 * (dx2 * dy2) / (dx2 + dy2);
double *p = solver->p; double* p = solver->p;
double *rhs = solver->rhs; double* rhs = solver->rhs;
double epssq = eps * eps; double epssq = eps * eps;
res = eps + 1.0; res = eps + 1.0;
while ((res >= epssq) && (it < itermax)) { while ((res >= epssq) && (it < itermax)) {
res = 0.0; res = 0.0;
exchange(solver); exchange(solver);
for (int j = 1; j < jmaxLocal + 1; j++) { for (int j = 1; j < jmaxLocal + 1; j++) {
for (int i = 1; i < imax + 1; i++) { for (int i = 1; i < imax + 1; i++) {
r = RHS(i, j) - ((P(i - 1, j) - 2.0 * P(i, j) + P(i + 1, j)) * idx2 + 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 - 1) - 2.0 * P(i, j) + P(i, j + 1)) * idy2);
P(i, j) -= (factor * r); P(i, j) -= (factor * r);
res += (r * r); res += (r * r);
} }
}
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++) {
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 = sqrt(res / (imax * jmax));
#ifdef DEBUG
if (solver->rank == 0) {
printf("%d Residuum: %e\n", it, res1);
}
#endif
it++;
} }
if (solver->rank == 0) { if (solver->rank == 0) {
for (int i = 1; i < imax + 1; i++) { printf("Solver took %d iterations\n", it);
P(i, 0) = P(i, 1);
}
} }
if (res < eps) {
if (solver->rank == (solver->size - 1)) { return 1;
for (int i = 1; i < imax + 1; i++) { } else {
P(i, jmaxLocal + 1) = P(i, jmaxLocal); return 0;
}
} }
}
for (int j = 1; j < jmaxLocal + 1; j++) { int solveRB(Solver* solver)
P(0, j) = P(1, j); {
P(imax + 1, j) = P(imax, j); double r;
} int it = 0;
double res, res1;
MPI_Allreduce(&res, &res1, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); int imax = solver->imax;
res = res1; int jmax = solver->jmax;
res = sqrt(res / (imax * 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 #ifdef DEBUG
printf("%d Residuum: %e\n", it, res);
#endif
it++;
}
if (solver->rank == 0) { if (solver->rank == 0) {
printf("%d Residuum: %e\n", it, res1); printf("Solver took %d iterations\n", it);
}
if (res < eps) {
return 1;
} else {
return 0;
} }
#endif
it++;
}
if (solver->rank == 0) {
printf("Solver took %d iterations\n", it);
}
if (res < eps) {
return 1;
} else {
return 0;
}
} }
int solveRB(Solver *solver) { int solveRBA(Solver* solver)
double r; {
int it = 0; double r;
double res, res1; int it = 0;
double res;
int imax = solver->imax; int imax = solver->imax;
int jmax = solver->jmax; int jmax = solver->jmax;
int jmaxLocal = solver->jmaxLocal; int jmaxLocal = solver->jmaxLocal;
double eps = solver->eps; double eps = solver->eps;
double omega = solver->omega; double omega = solver->omega;
int itermax = solver->itermax; int itermax = solver->itermax;
double dx2 = solver->dx * solver->dx; double dx2 = solver->dx * solver->dx;
double dy2 = solver->dy * solver->dy; double dy2 = solver->dy * solver->dy;
double idx2 = 1.0 / dx2; double idx2 = 1.0 / dx2;
double idy2 = 1.0 / dy2; double idy2 = 1.0 / dy2;
double factor = omega * 0.5 * (dx2 * dy2) / (dx2 + dy2); double factor = omega * 0.5 * (dx2 * dy2) / (dx2 + dy2);
double *p = solver->p; double* p = solver->p;
double *rhs = solver->rhs; double* rhs = solver->rhs;
int pass, jsw, isw; int pass, jsw, isw;
double epssq = eps * eps; double rho = solver->rho;
double epssq = eps * eps;
res = eps + 1.0; res = eps + 1.0;
while ((res >= epssq) && (it < itermax)) { while ((res >= epssq) && (it < itermax)) {
res = 0.0; res = 0.0;
jsw = 1; jsw = 1;
for (pass = 0; pass < 2; pass++) { for (pass = 0; pass < 2; pass++) {
isw = jsw; isw = jsw;
exchange(solver); exchange(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) -= (omega * factor * r);
res += (r * 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));
} }
isw = 3 - isw;
}
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++)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++) {
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) { for (int i = 1; i < imax + 1; i++) {
printf("Solver took %d iterations\n", it); P(i, 0) = P(i, 1);
} P(i, jmaxLocal + 1) = P(i, jmaxLocal);
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++) { for (int j = 1; j < jmaxLocal + 1; j++) {
P(i, 0) = P(i, 1); P(0, j) = P(1, j);
P(i, jmaxLocal + 1) = P(i, jmaxLocal); P(imax + 1, j) = P(imax, j);
} }
for (int j = 1; j < jmaxLocal + 1; j++) { res = res / (double)(imax * jmax);
P(0, j) = P(1, j);
P(imax + 1, j) = P(imax, j);
}
res = res / (double)(imax * jmax);
#ifdef DEBUG #ifdef DEBUG
printf("%d Residuum: %e Omega: %e\n", it, res, omega); printf("%d Residuum: %e Omega: %e\n", it, res, omega);
#endif #endif
it++; 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;
double *p = solver->p;
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 ", m[j * (imax + 2) + i]);
} }
fprintf(fp, "\n");
}
fclose(fp); 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;
double* p = solver->p;
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 ", m[j * (imax + 2) + i]);
}
fprintf(fp, "\n");
}
fclose(fp);
} }

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