diff --git a/BasicSolver/2D-mpi/src/discretization.c b/BasicSolver/2D-mpi/src/discretization.c new file mode 100644 index 0000000..e9fc421 --- /dev/null +++ b/BasicSolver/2D-mpi/src/discretization.c @@ -0,0 +1,453 @@ +/* + * Copyright (C) NHR@FAU, University Erlangen-Nuremberg. + * All rights reserved. This file is part of nusif-solver. + * Use of this source code is governed by a MIT style + * license that can be found in the LICENSE file. + */ +#include +#include +#include +#include +#include + +#include "allocate.h" +#include "comm.h" +#include "discretization.h" +#include "parameter.h" +#include "util.h" + +static void printConfig(Discretization* s) +{ + if (commIsMaster(&s->comm)) { + printf("Parameters for #%s#\n", s->problem); + printf("BC Left:%d Right:%d Bottom:%d Top:%d\n", + s->bcLeft, + s->bcRight, + s->bcBottom, + s->bcTop); + printf("\tReynolds number: %.2f\n", s->re); + printf("\tGx Gy: %.2f %.2f\n", s->gx, s->gy); + printf("Geometry data:\n"); + printf("\tDomain box size (x, y): %.2f, %.2f\n", + s->grid.xlength, + s->grid.ylength); + printf("\tCells (x, y): %d, %d\n", s->grid.imax, s->grid.jmax); + printf("\tCell size (dx, dy): %f, %f\n", s->grid.dx, s->grid.dy); + printf("Timestep parameters:\n"); + printf("\tDefault stepsize: %.2f, Final time %.2f\n", s->dt, s->te); + printf("\tdt bound: %.6f\n", s->dtBound); + printf("\tTau factor: %.2f\n", s->tau); + printf("Iterative s parameters:\n"); + printf("\tgamma factor: %f\n", s->gamma); + } + commPrintConfig(&s->comm); +} + +void initDiscretiztion(Discretization* s, Parameter* params) +{ + s->problem = params->name; + s->bcLeft = params->bcLeft; + s->bcRight = params->bcRight; + s->bcBottom = params->bcBottom; + s->bcTop = params->bcTop; + s->grid.imax = params->imax; + s->grid.jmax = params->jmax; + s->grid.xlength = params->xlength; + s->grid.ylength = params->ylength; + s->grid.dx = params->xlength / params->imax; + s->grid.dy = params->ylength / params->jmax; + s->re = params->re; + s->gx = params->gx; + s->gy = params->gy; + s->dt = params->dt; + s->te = params->te; + s->tau = params->tau; + s->gamma = params->gamma; + + /* allocate arrays */ + int imaxLocal = s->comm.imaxLocal; + int jmaxLocal = s->comm.jmaxLocal; + size_t size = (imaxLocal + 2) * (jmaxLocal + 2); + + s->u = allocate(64, size * sizeof(double)); + s->v = allocate(64, size * sizeof(double)); + s->p = allocate(64, size * sizeof(double)); + s->rhs = allocate(64, size * sizeof(double)); + s->f = allocate(64, size * sizeof(double)); + s->g = allocate(64, size * sizeof(double)); + + for (int i = 0; i < size; i++) { + s->u[i] = params->u_init; + s->v[i] = params->v_init; + s->p[i] = params->p_init; + s->rhs[i] = 0.0; + s->f[i] = 0.0; + s->g[i] = 0.0; + } + + double dx = s->grid.dx; + double dy = s->grid.dy; + + double invSqrSum = 1.0 / (dx * dx) + 1.0 / (dy * dy); + s->dtBound = 0.5 * s->re * 1.0 / invSqrSum; +#ifdef VERBOSE + printConfig(s); +#endif +} + +void computeRHS(Discretization* s) +{ + int imaxLocal = s->comm.imaxLocal; + int jmaxLocal = s->comm.jmaxLocal; + double idx = 1.0 / s->grid.dx; + double idy = 1.0 / s->grid.dy; + double idt = 1.0 / s->dt; + double* rhs = s->rhs; + double* f = s->f; + double* g = s->g; + + commShift(&s->comm, f, g); + + for (int j = 1; j < jmaxLocal + 1; j++) { + for (int i = 1; i < imaxLocal + 1; i++) { + RHS(i, j) = ((F(i, j) - F(i - 1, j)) * idx + (G(i, j) - G(i, j - 1)) * idy) * + idt; + } + } +} + +static double maxElement(Discretization* s, double* m) +{ + int imaxLocal = s->comm.imaxLocal; + int jmaxLocal = s->comm.jmaxLocal; + int size = (imaxLocal + 2) * (jmaxLocal + 2); + double maxval = DBL_MIN; + + for (int i = 0; i < size; i++) { + maxval = MAX(maxval, fabs(m[i])); + } + + commReduction(&maxval, MAX); + return maxval; +} + +void computeTimestep(Discretization* s) +{ + double dt = s->dtBound; + double dx = s->grid.dx; + double dy = s->grid.dy; + double umax = maxElement(s, s->u); + double vmax = maxElement(s, s->v); + + if (umax > 0) { + dt = (dt > dx / umax) ? dx / umax : dt; + } + if (vmax > 0) { + dt = (dt > dy / vmax) ? dy / vmax : dt; + } + + s->dt = dt * s->tau; +} + +void setBoundaryConditions(Discretization* s) +{ + int imaxLocal = s->comm.imaxLocal; + int jmaxLocal = s->comm.jmaxLocal; + double* u = s->u; + double* v = s->v; + + if (commIsBoundary(&s->comm, TOP)) { + switch (s->bcTop) { + case NOSLIP: + for (int i = 1; i < imaxLocal + 1; i++) { + V(i, jmaxLocal) = 0.0; + U(i, jmaxLocal + 1) = -U(i, jmaxLocal); + } + break; + case SLIP: + for (int i = 1; i < imaxLocal + 1; i++) { + V(i, jmaxLocal) = 0.0; + U(i, jmaxLocal + 1) = U(i, jmaxLocal); + } + break; + case OUTFLOW: + for (int i = 1; i < imaxLocal + 1; i++) { + U(i, jmaxLocal + 1) = U(i, jmaxLocal); + V(i, jmaxLocal) = V(i, jmaxLocal - 1); + } + break; + case PERIODIC: + break; + } + } + + if (commIsBoundary(&s->comm, BOTTOM)) { + switch (s->bcBottom) { + case NOSLIP: + for (int i = 1; i < imaxLocal + 1; i++) { + V(i, 0) = 0.0; + U(i, 0) = -U(i, 1); + } + break; + case SLIP: + for (int i = 1; i < imaxLocal + 1; i++) { + V(i, 0) = 0.0; + U(i, 0) = U(i, 1); + } + break; + case OUTFLOW: + for (int i = 1; i < imaxLocal + 1; i++) { + U(i, 0) = U(i, 1); + V(i, 0) = V(i, 1); + } + break; + case PERIODIC: + break; + } + } + + if (commIsBoundary(&s->comm, RIGHT)) { + switch (s->bcRight) { + case NOSLIP: + for (int j = 1; j < jmaxLocal + 1; j++) { + U(imaxLocal, j) = 0.0; + V(imaxLocal + 1, j) = -V(imaxLocal, j); + } + break; + case SLIP: + for (int j = 1; j < jmaxLocal + 1; j++) { + U(imaxLocal, j) = 0.0; + V(imaxLocal + 1, j) = V(imaxLocal, j); + } + break; + case OUTFLOW: + for (int j = 1; j < jmaxLocal + 1; j++) { + U(imaxLocal, j) = U(imaxLocal - 1, j); + V(imaxLocal + 1, j) = V(imaxLocal, j); + } + break; + case PERIODIC: + break; + } + } + + if (commIsBoundary(&s->comm, LEFT)) { + switch (s->bcLeft) { + case NOSLIP: + for (int j = 1; j < jmaxLocal + 1; j++) { + U(0, j) = 0.0; + V(0, j) = -V(1, j); + } + break; + case SLIP: + for (int j = 1; j < jmaxLocal + 1; j++) { + U(0, j) = 0.0; + V(0, j) = V(1, j); + } + break; + case OUTFLOW: + for (int j = 1; j < jmaxLocal + 1; j++) { + U(0, j) = U(1, j); + V(0, j) = V(1, j); + } + break; + case PERIODIC: + break; + } + } +} + +void setSpecialBoundaryCondition(Discretization* s) +{ + int imaxLocal = s->comm.imaxLocal; + int jmaxLocal = s->comm.jmaxLocal; + double* u = s->u; + + if (strcmp(s->problem, "dcavity") == 0) { + if (commIsBoundary(&s->comm, TOP)) { + for (int i = 1; i < imaxLocal + 1; i++) { + U(i, jmaxLocal + 1) = 2.0 - U(i, jmaxLocal); + } + } + } else if (strcmp(s->problem, "canal") == 0) { + if (commIsBoundary(&s->comm, LEFT)) { + double ylength = s->grid.ylength; + double dy = s->grid.dy; + int rest = s->grid.jmax % s->comm.size; + int yc = s->comm.rank * (s->grid.jmax / s->comm.size) + + MIN(rest, s->comm.rank); + double ys = dy * (yc + 0.5); + double y; + + /* printf("RANK %d yc: %d ys: %f\n", solver->rank, yc, ys); */ + + for (int j = 1; j < jmaxLocal + 1; j++) { + y = ys + dy * (j - 0.5); + U(0, j) = y * (ylength - y) * 4.0 / (ylength * ylength); + } + } + } + /* print(solver, solver->u); */ +} + +void computeFG(Discretization* s) +{ + double* u = s->u; + double* v = s->v; + double* f = s->f; + double* g = s->g; + + int imaxLocal = s->comm.imaxLocal; + int jmaxLocal = s->comm.jmaxLocal; + + double gx = s->gx; + double gy = s->gy; + double gamma = s->gamma; + double dt = s->dt; + double inverseRe = 1.0 / s->re; + double inverseDx = 1.0 / s->grid.dx; + double inverseDy = 1.0 / s->grid.dy; + double du2dx, dv2dy, duvdx, duvdy; + double du2dx2, du2dy2, dv2dx2, dv2dy2; + + commExchange(&s->comm, u); + commExchange(&s->comm, v); + + for (int j = 1; j < jmaxLocal + 1; j++) { + for (int i = 1; i < imaxLocal + 1; i++) { + du2dx = inverseDx * 0.25 * + ((U(i, j) + U(i + 1, j)) * (U(i, j) + U(i + 1, j)) - + (U(i, j) + U(i - 1, j)) * (U(i, j) + U(i - 1, j))) + + gamma * inverseDx * 0.25 * + (fabs(U(i, j) + U(i + 1, j)) * (U(i, j) - U(i + 1, j)) + + fabs(U(i, j) + U(i - 1, j)) * (U(i, j) - U(i - 1, j))); + + duvdy = inverseDy * 0.25 * + ((V(i, j) + V(i + 1, j)) * (U(i, j) + U(i, j + 1)) - + (V(i, j - 1) + V(i + 1, j - 1)) * (U(i, j) + U(i, j - 1))) + + gamma * inverseDy * 0.25 * + (fabs(V(i, j) + V(i + 1, j)) * (U(i, j) - U(i, j + 1)) + + fabs(V(i, j - 1) + V(i + 1, j - 1)) * + (U(i, j) - U(i, j - 1))); + + du2dx2 = inverseDx * inverseDx * (U(i + 1, j) - 2.0 * U(i, j) + U(i - 1, j)); + du2dy2 = inverseDy * inverseDy * (U(i, j + 1) - 2.0 * U(i, j) + U(i, j - 1)); + F(i, j) = U(i, j) + dt * (inverseRe * (du2dx2 + du2dy2) - du2dx - duvdy + gx); + + duvdx = inverseDx * 0.25 * + ((U(i, j) + U(i, j + 1)) * (V(i, j) + V(i + 1, j)) - + (U(i - 1, j) + U(i - 1, j + 1)) * (V(i, j) + V(i - 1, j))) + + gamma * inverseDx * 0.25 * + (fabs(U(i, j) + U(i, j + 1)) * (V(i, j) - V(i + 1, j)) + + fabs(U(i - 1, j) + U(i - 1, j + 1)) * + (V(i, j) - V(i - 1, j))); + + dv2dy = inverseDy * 0.25 * + ((V(i, j) + V(i, j + 1)) * (V(i, j) + V(i, j + 1)) - + (V(i, j) + V(i, j - 1)) * (V(i, j) + V(i, j - 1))) + + gamma * inverseDy * 0.25 * + (fabs(V(i, j) + V(i, j + 1)) * (V(i, j) - V(i, j + 1)) + + fabs(V(i, j) + V(i, j - 1)) * (V(i, j) - V(i, j - 1))); + + dv2dx2 = inverseDx * inverseDx * (V(i + 1, j) - 2.0 * V(i, j) + V(i - 1, j)); + dv2dy2 = inverseDy * inverseDy * (V(i, j + 1) - 2.0 * V(i, j) + V(i, j - 1)); + G(i, j) = V(i, j) + dt * (inverseRe * (dv2dx2 + dv2dy2) - duvdx - dv2dy + gy); + } + } + + /* ----------------------------- boundary of F --------------------------- */ + if (commIsBoundary(&s->comm, LEFT)) { + for (int j = 1; j < jmaxLocal + 1; j++) { + F(0, j) = U(0, j); + } + } + + if (commIsBoundary(&s->comm, RIGHT)) { + for (int j = 1; j < jmaxLocal + 1; j++) { + F(imaxLocal, j) = U(imaxLocal, j); + } + } + + /* ----------------------------- boundary of G --------------------------- */ + if (commIsBoundary(&s->comm, BOTTOM)) { + for (int i = 1; i < imaxLocal + 1; i++) { + G(i, 0) = V(i, 0); + } + } + + if (commIsBoundary(&s->comm, TOP)) { + for (int i = 1; i < imaxLocal + 1; i++) { + G(i, jmaxLocal) = V(i, jmaxLocal); + } + } +} + +void adaptUV(Discretization* s) +{ + int imaxLocal = s->comm.imaxLocal; + int jmaxLocal = s->comm.jmaxLocal; + + double* p = s->p; + double* u = s->u; + double* v = s->v; + double* f = s->f; + double* g = s->g; + + double factorX = s->dt / s->grid.dx; + double factorY = s->dt / s->grid.dy; + + for (int j = 1; j < jmaxLocal + 1; j++) { + for (int i = 1; i < imaxLocal + 1; i++) { + U(i, j) = F(i, j) - (P(i + 1, j) - P(i, j)) * factorX; + V(i, j) = G(i, j) - (P(i, j + 1) - P(i, j)) * factorY; + } + } +} + +void writeResult(Discretization* s, double* u, double* v, double* p) +{ + int imax = s->grid.imax; + int jmax = s->grid.jmax; + double dx = s->grid.dx; + double dy = s->grid.dy; + double x = 0.0, y = 0.0; + + FILE* fp; + fp = fopen("pressure.dat", "w"); + + if (fp == NULL) { + printf("Error!\n"); + exit(EXIT_FAILURE); + } + + for (int j = 1; j <= jmax; j++) { + y = (double)(j - 0.5) * dy; + for (int i = 1; i <= imax; i++) { + x = (double)(i - 0.5) * dx; + fprintf(fp, "%.2f %.2f %f\n", x, y, p[j * (imax + 2) + i]); + } + fprintf(fp, "\n"); + } + + fclose(fp); + + fp = fopen("velocity.dat", "w"); + + if (fp == NULL) { + printf("Error!\n"); + exit(EXIT_FAILURE); + } + + for (int j = 1; j <= jmax; j++) { + y = dy * (j - 0.5); + for (int i = 1; i <= imax; i++) { + x = dx * (i - 0.5); + double velU = (u[j * (imax + 2) + i] + u[j * (imax + 2) + (i - 1)]) / 2.0; + double velV = (v[j * (imax + 2) + i] + v[(j - 1) * (imax + 2) + i]) / 2.0; + double len = sqrt((velU * velU) + (velV * velV)); + fprintf(fp, "%.2f %.2f %f %f %f\n", x, y, velU, velV, len); + } + } + + fclose(fp); +} diff --git a/BasicSolver/2D-mpi/src/discretization.h b/BasicSolver/2D-mpi/src/discretization.h new file mode 100644 index 0000000..a20e295 --- /dev/null +++ b/BasicSolver/2D-mpi/src/discretization.h @@ -0,0 +1,43 @@ +/* + * Copyright (C) NHR@FAU, University Erlangen-Nuremberg. + * All rights reserved. This file is part of nusif-solver. + * Use of this source code is governed by a MIT style + * license that can be found in the LICENSE file. + */ +#ifndef __DISCRETIZATION_H_ +#define __DISCRETIZATION_H_ +#include "comm.h" +#include "grid.h" +#include "parameter.h" + +enum BC { NOSLIP = 1, SLIP, OUTFLOW, PERIODIC }; + +typedef struct { + /* geometry and grid information */ + Grid grid; + /* arrays */ + double *p, *rhs; + double *f, *g; + double *u, *v; + /* parameters */ + double re, tau, gamma; + double gx, gy; + /* time stepping */ + double dt, te; + double dtBound; + char* problem; + int bcLeft, bcRight, bcBottom, bcTop; + /* communication */ + Comm comm; +} Discretization; + +void initDiscretiztion(Discretization*, Parameter*); +void computeRHS(Discretization*); +void normalizePressure(Discretization*); +void computeTimestep(Discretization*); +void setBoundaryConditions(Discretization*); +void setSpecialBoundaryCondition(Discretization*); +void computeFG(Discretization*); +void adaptUV(Discretization*); +void writeResult(Discretization* s, double* u, double* v, double* p); +#endif diff --git a/BasicSolver/2D-mpi/src/grid.h b/BasicSolver/2D-mpi/src/grid.h new file mode 100644 index 0000000..7e58d1a --- /dev/null +++ b/BasicSolver/2D-mpi/src/grid.h @@ -0,0 +1,16 @@ +/* + * Copyright (C) NHR@FAU, University Erlangen-Nuremberg. + * All rights reserved. This file is part of nusif-solver. + * Use of this source code is governed by a MIT style + * license that can be found in the LICENSE file. + */ +#ifndef __GRID_H_ +#define __GRID_H_ + +typedef struct { + double dx, dy; + int imax, jmax; + double xlength, ylength; +} Grid; + +#endif // __GRID_H_ diff --git a/BasicSolver/2D-mpi/src/main.c b/BasicSolver/2D-mpi/src/main.c index 616f0a1..c7f2d2f 100644 --- a/BasicSolver/2D-mpi/src/main.c +++ b/BasicSolver/2D-mpi/src/main.c @@ -10,21 +10,22 @@ #include "allocate.h" #include "comm.h" +#include "discretization.h" #include "parameter.h" #include "progress.h" #include "solver.h" #include "timing.h" -static void writeResults(Solver* s) +static void writeResults(Discretization* s) { #ifdef _MPI - size_t bytesize = (s->imax + 2) * (s->jmax + 2) * sizeof(double); + size_t bytesize = (s->grid.imax + 2) * (s->grid.jmax + 2) * sizeof(double); double* ug = allocate(64, bytesize); double* vg = allocate(64, bytesize); double* pg = allocate(64, bytesize); - commCollectResult(&s->comm, ug, vg, pg, s->u, s->v, s->p, s->imax, s->jmax); + commCollectResult(&s->comm, ug, vg, pg, s->u, s->v, s->p, s->grid.imax, s->grid.jmax); writeResult(s, ug, vg, pg); free(ug); @@ -40,9 +41,10 @@ int main(int argc, char** argv) int rank; double timeStart, timeStop; Parameter p; + Discretization d; Solver s; - commInit(&s.comm, argc, argv); + commInit(&d.comm, argc, argv); initParameter(&p); if (argc != 2) { @@ -51,44 +53,45 @@ int main(int argc, char** argv) } readParameter(&p, argv[1]); - commPartition(&s.comm, p.jmax, p.imax); - if (commIsMaster(&s.comm)) { + commPartition(&d.comm, p.jmax, p.imax); + if (commIsMaster(&d.comm)) { printParameter(&p); } - initSolver(&s, &p); + initDiscretiztion(&d, &p); + initSolver(&s, &d, &p); #ifdef TEST - commPrintConfig(&s.comm); - commTestInit(&s.comm, s.p, s.f, s.g); - commExchange(&s.comm, s.p); - commShift(&s.comm, s.f, s.g); - commTestWrite(&s.comm, s.p, s.f, s.g); - writeResults(&s); - commFinalize(&s.comm); + commPrintConfig(&d.comm); + commTestInit(&d.comm, d.p, d.f, d.g); + commExchange(&d.comm, d.p); + commShift(&d.comm, d.f, d.g); + commTestWrite(&d.comm, d.p, d.f, d.g); + writeResults(&d); + commFinalize(&d.comm); exit(EXIT_SUCCESS); #endif #ifndef VERBOSE - initProgress(s.te); + initProgress(d.te); #endif - double tau = s.tau; - double te = s.te; + double tau = d.tau; + double te = d.te; double t = 0.0; timeStart = getTimeStamp(); while (t <= te) { - if (tau > 0.0) computeTimestep(&s); - setBoundaryConditions(&s); - setSpecialBoundaryCondition(&s); - computeFG(&s); - computeRHS(&s); - solve(&s); - adaptUV(&s); - t += s.dt; + if (tau > 0.0) computeTimestep(&d); + setBoundaryConditions(&d); + setSpecialBoundaryCondition(&d); + computeFG(&d); + computeRHS(&d); + solve(&s, d.p, d.rhs); + adaptUV(&d); + t += d.dt; #ifdef VERBOSE - if (commIsMaster(&s.comm)) { - printf("TIME %f , TIMESTEP %f\n", t, s.dt); + if (commIsMaster(s.comm)) { + printf("TIME %f , TIMESTEP %f\n", t, d.dt); } #else printProgress(t); @@ -98,11 +101,11 @@ int main(int argc, char** argv) #ifndef VERBOSE stopProgress(); #endif - if (commIsMaster(&s.comm)) { + if (commIsMaster(s.comm)) { printf("Solution took %.2fs\n", timeStop - timeStart); } - writeResults(&s); - commFinalize(&s.comm); + writeResults(&d); + commFinalize(s.comm); return EXIT_SUCCESS; } diff --git a/BasicSolver/2D-mpi/src/solver.c b/BasicSolver/2D-mpi/src/solver.c index dcd6906..654f59d 100644 --- a/BasicSolver/2D-mpi/src/solver.c +++ b/BasicSolver/2D-mpi/src/solver.c @@ -4,144 +4,38 @@ * Use of this source code is governed by a MIT style * license that can be found in the LICENSE file. */ -#include #include #include -#include -#include #include "allocate.h" #include "comm.h" +#include "discretization.h" #include "parameter.h" #include "solver.h" #include "util.h" -#define P(i, j) p[(j) * (imaxLocal + 2) + (i)] -#define F(i, j) f[(j) * (imaxLocal + 2) + (i)] -#define G(i, j) g[(j) * (imaxLocal + 2) + (i)] -#define U(i, j) u[(j) * (imaxLocal + 2) + (i)] -#define V(i, j) v[(j) * (imaxLocal + 2) + (i)] -#define RHS(i, j) rhs[(j) * (imaxLocal + 2) + (i)] - -static void printConfig(Solver* s) +void initSolver(Solver* s, Discretization* d, Parameter* p) { - if (commIsMaster(&s->comm)) { - printf("Parameters for #%s#\n", s->problem); - printf("BC Left:%d Right:%d Bottom:%d Top:%d\n", - s->bcLeft, - s->bcRight, - s->bcBottom, - s->bcTop); - printf("\tReynolds number: %.2f\n", s->re); - printf("\tGx Gy: %.2f %.2f\n", s->gx, s->gy); - printf("Geometry data:\n"); - printf("\tDomain box size (x, y): %.2f, %.2f\n", s->xlength, s->ylength); - printf("\tCells (x, y): %d, %d\n", s->imax, s->jmax); - printf("\tCell size (dx, dy): %f, %f\n", s->dx, s->dy); - printf("Timestep parameters:\n"); - printf("\tDefault stepsize: %.2f, Final time %.2f\n", s->dt, s->te); - printf("\tdt bound: %.6f\n", s->dtBound); - printf("\tTau factor: %.2f\n", s->tau); - printf("Iterative s parameters:\n"); - printf("\tMax iterations: %d\n", s->itermax); - printf("\tepsilon (stopping tolerance) : %f\n", s->eps); - printf("\tgamma factor: %f\n", s->gamma); - printf("\tomega (SOR relaxation): %f\n", s->omega); - } - commPrintConfig(&s->comm); + s->grid = &d->grid; + s->eps = p->eps; + s->omega = p->omg; + s->itermax = p->itermax; + s->comm = &d->comm; } -void initSolver(Solver* s, Parameter* params) +void solve(Solver* s, double* p, double* rhs) { - s->problem = params->name; - s->bcLeft = params->bcLeft; - s->bcRight = params->bcRight; - s->bcBottom = params->bcBottom; - s->bcTop = params->bcTop; - s->imax = params->imax; - s->jmax = params->jmax; - s->xlength = params->xlength; - s->ylength = params->ylength; - s->dx = params->xlength / params->imax; - s->dy = params->ylength / params->jmax; - s->eps = params->eps; - s->omega = params->omg; - s->itermax = params->itermax; - s->re = params->re; - s->gx = params->gx; - s->gy = params->gy; - s->dt = params->dt; - s->te = params->te; - s->tau = params->tau; - s->gamma = params->gamma; - - /* allocate arrays */ - int imaxLocal = s->comm.imaxLocal; - int jmaxLocal = s->comm.jmaxLocal; - size_t size = (imaxLocal + 2) * (jmaxLocal + 2); - - s->u = allocate(64, size * sizeof(double)); - s->v = allocate(64, size * sizeof(double)); - s->p = allocate(64, size * sizeof(double)); - s->rhs = allocate(64, size * sizeof(double)); - s->f = allocate(64, size * sizeof(double)); - s->g = allocate(64, size * sizeof(double)); - - for (int i = 0; i < size; i++) { - s->u[i] = params->u_init; - s->v[i] = params->v_init; - s->p[i] = params->p_init; - s->rhs[i] = 0.0; - s->f[i] = 0.0; - s->g[i] = 0.0; - } - - double dx = s->dx; - double dy = s->dy; - - double invSqrSum = 1.0 / (dx * dx) + 1.0 / (dy * dy); - s->dtBound = 0.5 * s->re * 1.0 / invSqrSum; -#ifdef VERBOSE - printConfig(s); -#endif -} - -void computeRHS(Solver* s) -{ - int imaxLocal = s->comm.imaxLocal; - int jmaxLocal = s->comm.jmaxLocal; - double idx = 1.0 / s->dx; - double idy = 1.0 / s->dy; - double idt = 1.0 / s->dt; - double* rhs = s->rhs; - double* f = s->f; - double* g = s->g; - - commShift(&s->comm, f, g); - - for (int j = 1; j < jmaxLocal + 1; j++) { - for (int i = 1; i < imaxLocal + 1; i++) { - RHS(i, j) = ((F(i, j) - F(i - 1, j)) * idx + (G(i, j) - G(i, j - 1)) * idy) * - idt; - } - } -} - -void solve(Solver* s) -{ - int imax = s->imax; - int jmax = s->jmax; - int imaxLocal = s->comm.imaxLocal; - int jmaxLocal = s->comm.jmaxLocal; + int imax = s->grid->imax; + int jmax = s->grid->jmax; + int imaxLocal = s->comm->imaxLocal; + int jmaxLocal = s->comm->jmaxLocal; double eps = s->eps; int itermax = s->itermax; - double dx2 = s->dx * s->dx; - double dy2 = s->dy * s->dy; + double dx2 = s->grid->dx * s->grid->dx; + double dy2 = s->grid->dy * s->grid->dy; double idx2 = 1.0 / dx2; double idy2 = 1.0 / dy2; double factor = s->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2); - double* p = s->p; - double* rhs = s->rhs; double epssq = eps * eps; int pass, jsw, isw; int it = 0; @@ -151,7 +45,7 @@ void solve(Solver* s) jsw = 1; for (pass = 0; pass < 2; pass++) { isw = jsw; - commExchange(&s->comm, p); + commExchange(s->comm, p); for (int j = 1; j < jmaxLocal + 1; j++) { for (int i = isw; i < imaxLocal + 1; i += 2) { @@ -168,25 +62,25 @@ void solve(Solver* s) jsw = 3 - jsw; } - if (commIsBoundary(&s->comm, BOTTOM)) { // set bottom bc + if (commIsBoundary(s->comm, BOTTOM)) { // set bottom bc for (int i = 1; i < imaxLocal + 1; i++) { P(i, 0) = P(i, 1); } } - if (commIsBoundary(&s->comm, TOP)) { // set top bc + if (commIsBoundary(s->comm, TOP)) { // set top bc for (int i = 1; i < imaxLocal + 1; i++) { P(i, jmaxLocal + 1) = P(i, jmaxLocal); } } - if (commIsBoundary(&s->comm, LEFT)) { // set left bc + if (commIsBoundary(s->comm, LEFT)) { // set left bc for (int j = 1; j < jmaxLocal + 1; j++) { P(0, j) = P(1, j); } } - if (commIsBoundary(&s->comm, RIGHT)) { // set right bc + if (commIsBoundary(s->comm, RIGHT)) { // set right bc for (int j = 1; j < jmaxLocal + 1; j++) { P(imaxLocal + 1, j) = P(imaxLocal, j); } @@ -195,7 +89,7 @@ void solve(Solver* s) commReduction(&res, SUM); res = res / (double)(imax * jmax); #ifdef DEBUG - if (commIsMaster(&s->comm)) { + if (commIsMaster(s->comm)) { printf("%d Residuum: %e\n", it, res); } #endif @@ -203,343 +97,8 @@ void solve(Solver* s) } #ifdef VERBOSE - if (commIsMaster(&s->comm)) { + if (commIsMaster(s->comm)) { printf("Solver took %d iterations to reach %f\n", it, sqrt(res)); } #endif } - -static double maxElement(Solver* s, double* m) -{ - int imaxLocal = s->comm.imaxLocal; - int jmaxLocal = s->comm.jmaxLocal; - int size = (imaxLocal + 2) * (jmaxLocal + 2); - double maxval = DBL_MIN; - - for (int i = 0; i < size; i++) { - maxval = MAX(maxval, fabs(m[i])); - } - - commReduction(&maxval, MAX); - return maxval; -} - -void computeTimestep(Solver* s) -{ - double dt = s->dtBound; - double dx = s->dx; - double dy = s->dy; - double umax = maxElement(s, s->u); - double vmax = maxElement(s, s->v); - - if (umax > 0) { - dt = (dt > dx / umax) ? dx / umax : dt; - } - if (vmax > 0) { - dt = (dt > dy / vmax) ? dy / vmax : dt; - } - - s->dt = dt * s->tau; -} - -void setBoundaryConditions(Solver* s) -{ - int imaxLocal = s->comm.imaxLocal; - int jmaxLocal = s->comm.jmaxLocal; - double* u = s->u; - double* v = s->v; - - if (commIsBoundary(&s->comm, TOP)) { - switch (s->bcTop) { - case NOSLIP: - for (int i = 1; i < imaxLocal + 1; i++) { - V(i, jmaxLocal) = 0.0; - U(i, jmaxLocal + 1) = -U(i, jmaxLocal); - } - break; - case SLIP: - for (int i = 1; i < imaxLocal + 1; i++) { - V(i, jmaxLocal) = 0.0; - U(i, jmaxLocal + 1) = U(i, jmaxLocal); - } - break; - case OUTFLOW: - for (int i = 1; i < imaxLocal + 1; i++) { - U(i, jmaxLocal + 1) = U(i, jmaxLocal); - V(i, jmaxLocal) = V(i, jmaxLocal - 1); - } - break; - case PERIODIC: - break; - } - } - - if (commIsBoundary(&s->comm, BOTTOM)) { - switch (s->bcBottom) { - case NOSLIP: - for (int i = 1; i < imaxLocal + 1; i++) { - V(i, 0) = 0.0; - U(i, 0) = -U(i, 1); - } - break; - case SLIP: - for (int i = 1; i < imaxLocal + 1; i++) { - V(i, 0) = 0.0; - U(i, 0) = U(i, 1); - } - break; - case OUTFLOW: - for (int i = 1; i < imaxLocal + 1; i++) { - U(i, 0) = U(i, 1); - V(i, 0) = V(i, 1); - } - break; - case PERIODIC: - break; - } - } - - if (commIsBoundary(&s->comm, RIGHT)) { - switch (s->bcRight) { - case NOSLIP: - for (int j = 1; j < jmaxLocal + 1; j++) { - U(imaxLocal, j) = 0.0; - V(imaxLocal + 1, j) = -V(imaxLocal, j); - } - break; - case SLIP: - for (int j = 1; j < jmaxLocal + 1; j++) { - U(imaxLocal, j) = 0.0; - V(imaxLocal + 1, j) = V(imaxLocal, j); - } - break; - case OUTFLOW: - for (int j = 1; j < jmaxLocal + 1; j++) { - U(imaxLocal, j) = U(imaxLocal - 1, j); - V(imaxLocal + 1, j) = V(imaxLocal, j); - } - break; - case PERIODIC: - break; - } - } - - if (commIsBoundary(&s->comm, LEFT)) { - switch (s->bcLeft) { - case NOSLIP: - for (int j = 1; j < jmaxLocal + 1; j++) { - U(0, j) = 0.0; - V(0, j) = -V(1, j); - } - break; - case SLIP: - for (int j = 1; j < jmaxLocal + 1; j++) { - U(0, j) = 0.0; - V(0, j) = V(1, j); - } - break; - case OUTFLOW: - for (int j = 1; j < jmaxLocal + 1; j++) { - U(0, j) = U(1, j); - V(0, j) = V(1, j); - } - break; - case PERIODIC: - break; - } - } -} - -void setSpecialBoundaryCondition(Solver* s) -{ - int imaxLocal = s->comm.imaxLocal; - int jmaxLocal = s->comm.jmaxLocal; - double* u = s->u; - - if (strcmp(s->problem, "dcavity") == 0) { - if (commIsBoundary(&s->comm, TOP)) { - for (int i = 1; i < imaxLocal + 1; i++) { - U(i, jmaxLocal + 1) = 2.0 - U(i, jmaxLocal); - } - } - } else if (strcmp(s->problem, "canal") == 0) { - if (commIsBoundary(&s->comm, LEFT)) { - double ylength = s->ylength; - double dy = s->dy; - int rest = s->jmax % s->comm.size; - int yc = s->comm.rank * (s->jmax / s->comm.size) + MIN(rest, s->comm.rank); - double ys = dy * (yc + 0.5); - double y; - - /* printf("RANK %d yc: %d ys: %f\n", solver->rank, yc, ys); */ - - for (int j = 1; j < jmaxLocal + 1; j++) { - y = ys + dy * (j - 0.5); - U(0, j) = y * (ylength - y) * 4.0 / (ylength * ylength); - } - } - } - /* print(solver, solver->u); */ -} - -void computeFG(Solver* s) -{ - double* u = s->u; - double* v = s->v; - double* f = s->f; - double* g = s->g; - - int imaxLocal = s->comm.imaxLocal; - int jmaxLocal = s->comm.jmaxLocal; - - double gx = s->gx; - double gy = s->gy; - double gamma = s->gamma; - double dt = s->dt; - double inverseRe = 1.0 / s->re; - double inverseDx = 1.0 / s->dx; - double inverseDy = 1.0 / s->dy; - double du2dx, dv2dy, duvdx, duvdy; - double du2dx2, du2dy2, dv2dx2, dv2dy2; - - commExchange(&s->comm, u); - commExchange(&s->comm, v); - - for (int j = 1; j < jmaxLocal + 1; j++) { - for (int i = 1; i < imaxLocal + 1; i++) { - du2dx = inverseDx * 0.25 * - ((U(i, j) + U(i + 1, j)) * (U(i, j) + U(i + 1, j)) - - (U(i, j) + U(i - 1, j)) * (U(i, j) + U(i - 1, j))) + - gamma * inverseDx * 0.25 * - (fabs(U(i, j) + U(i + 1, j)) * (U(i, j) - U(i + 1, j)) + - fabs(U(i, j) + U(i - 1, j)) * (U(i, j) - U(i - 1, j))); - - duvdy = inverseDy * 0.25 * - ((V(i, j) + V(i + 1, j)) * (U(i, j) + U(i, j + 1)) - - (V(i, j - 1) + V(i + 1, j - 1)) * (U(i, j) + U(i, j - 1))) + - gamma * inverseDy * 0.25 * - (fabs(V(i, j) + V(i + 1, j)) * (U(i, j) - U(i, j + 1)) + - fabs(V(i, j - 1) + V(i + 1, j - 1)) * - (U(i, j) - U(i, j - 1))); - - du2dx2 = inverseDx * inverseDx * (U(i + 1, j) - 2.0 * U(i, j) + U(i - 1, j)); - du2dy2 = inverseDy * inverseDy * (U(i, j + 1) - 2.0 * U(i, j) + U(i, j - 1)); - F(i, j) = U(i, j) + dt * (inverseRe * (du2dx2 + du2dy2) - du2dx - duvdy + gx); - - duvdx = inverseDx * 0.25 * - ((U(i, j) + U(i, j + 1)) * (V(i, j) + V(i + 1, j)) - - (U(i - 1, j) + U(i - 1, j + 1)) * (V(i, j) + V(i - 1, j))) + - gamma * inverseDx * 0.25 * - (fabs(U(i, j) + U(i, j + 1)) * (V(i, j) - V(i + 1, j)) + - fabs(U(i - 1, j) + U(i - 1, j + 1)) * - (V(i, j) - V(i - 1, j))); - - dv2dy = inverseDy * 0.25 * - ((V(i, j) + V(i, j + 1)) * (V(i, j) + V(i, j + 1)) - - (V(i, j) + V(i, j - 1)) * (V(i, j) + V(i, j - 1))) + - gamma * inverseDy * 0.25 * - (fabs(V(i, j) + V(i, j + 1)) * (V(i, j) - V(i, j + 1)) + - fabs(V(i, j) + V(i, j - 1)) * (V(i, j) - V(i, j - 1))); - - dv2dx2 = inverseDx * inverseDx * (V(i + 1, j) - 2.0 * V(i, j) + V(i - 1, j)); - dv2dy2 = inverseDy * inverseDy * (V(i, j + 1) - 2.0 * V(i, j) + V(i, j - 1)); - G(i, j) = V(i, j) + dt * (inverseRe * (dv2dx2 + dv2dy2) - duvdx - dv2dy + gy); - } - } - - /* ----------------------------- boundary of F --------------------------- */ - if (commIsBoundary(&s->comm, LEFT)) { - for (int j = 1; j < jmaxLocal + 1; j++) { - F(0, j) = U(0, j); - } - } - - if (commIsBoundary(&s->comm, RIGHT)) { - for (int j = 1; j < jmaxLocal + 1; j++) { - F(imaxLocal, j) = U(imaxLocal, j); - } - } - - /* ----------------------------- boundary of G --------------------------- */ - if (commIsBoundary(&s->comm, BOTTOM)) { - for (int i = 1; i < imaxLocal + 1; i++) { - G(i, 0) = V(i, 0); - } - } - - if (commIsBoundary(&s->comm, TOP)) { - for (int i = 1; i < imaxLocal + 1; i++) { - G(i, jmaxLocal) = V(i, jmaxLocal); - } - } -} - -void adaptUV(Solver* s) -{ - int imaxLocal = s->comm.imaxLocal; - int jmaxLocal = s->comm.jmaxLocal; - - double* p = s->p; - double* u = s->u; - double* v = s->v; - double* f = s->f; - double* g = s->g; - - double factorX = s->dt / s->dx; - double factorY = s->dt / s->dy; - - for (int j = 1; j < jmaxLocal + 1; j++) { - for (int i = 1; i < imaxLocal + 1; i++) { - U(i, j) = F(i, j) - (P(i + 1, j) - P(i, j)) * factorX; - V(i, j) = G(i, j) - (P(i, j + 1) - P(i, j)) * factorY; - } - } -} - -void writeResult(Solver* s, double* u, double* v, double* p) -{ - int imax = s->imax; - int jmax = s->jmax; - double dx = s->dx; - double dy = s->dy; - double x = 0.0, y = 0.0; - - FILE* fp; - fp = fopen("pressure.dat", "w"); - - if (fp == NULL) { - printf("Error!\n"); - exit(EXIT_FAILURE); - } - - for (int j = 1; j <= jmax; j++) { - y = (double)(j - 0.5) * dy; - for (int i = 1; i <= imax; i++) { - x = (double)(i - 0.5) * dx; - fprintf(fp, "%.2f %.2f %f\n", x, y, p[j * (imax + 2) + i]); - } - fprintf(fp, "\n"); - } - - fclose(fp); - - fp = fopen("velocity.dat", "w"); - - if (fp == NULL) { - printf("Error!\n"); - exit(EXIT_FAILURE); - } - - for (int j = 1; j <= jmax; j++) { - y = dy * (j - 0.5); - for (int i = 1; i <= imax; i++) { - x = dx * (i - 0.5); - double velU = (u[j * (imax + 2) + i] + u[j * (imax + 2) + (i - 1)]) / 2.0; - double velV = (v[j * (imax + 2) + i] + v[(j - 1) * (imax + 2) + i]) / 2.0; - double len = sqrt((velU * velU) + (velV * velV)); - fprintf(fp, "%.2f %.2f %f %f %f\n", x, y, velU, velV, len); - } - } - - fclose(fp); -} diff --git a/BasicSolver/2D-mpi/src/solver.h b/BasicSolver/2D-mpi/src/solver.h index 2f38522..4ceb40b 100644 --- a/BasicSolver/2D-mpi/src/solver.h +++ b/BasicSolver/2D-mpi/src/solver.h @@ -7,41 +7,20 @@ #ifndef __SOLVER_H_ #define __SOLVER_H_ #include "comm.h" +#include "discretization.h" +#include "grid.h" #include "parameter.h" -enum BC { NOSLIP = 1, SLIP, OUTFLOW, PERIODIC }; - typedef struct { /* geometry and grid information */ - double dx, dy; - int imax, jmax; - double xlength, ylength; - /* arrays */ - double *p, *rhs; - double *f, *g; - double *u, *v; + Grid* grid; /* parameters */ double eps, omega; - double re, tau, gamma; - double gx, gy; - /* time stepping */ int itermax; - double dt, te; - double dtBound; - char* problem; - int bcLeft, bcRight, bcBottom, bcTop; /* communication */ - Comm comm; + Comm* comm; } Solver; -void initSolver(Solver*, Parameter*); -void computeRHS(Solver*); -void solve(Solver*); -void normalizePressure(Solver*); -void computeTimestep(Solver*); -void setBoundaryConditions(Solver*); -void setSpecialBoundaryCondition(Solver*); -void computeFG(Solver*); -void adaptUV(Solver*); -void writeResult(Solver* s, double* u, double* v, double* p); +void initSolver(Solver*, Discretization*, Parameter*); +void solve(Solver*, double*, double*); #endif diff --git a/BasicSolver/2D-mpi/src/util.h b/BasicSolver/2D-mpi/src/util.h index c27b2ba..e36948a 100644 --- a/BasicSolver/2D-mpi/src/util.h +++ b/BasicSolver/2D-mpi/src/util.h @@ -19,4 +19,11 @@ #define ABS(a) ((a) >= 0 ? (a) : -(a)) #endif +#define P(i, j) p[(j) * (imaxLocal + 2) + (i)] +#define F(i, j) f[(j) * (imaxLocal + 2) + (i)] +#define G(i, j) g[(j) * (imaxLocal + 2) + (i)] +#define U(i, j) u[(j) * (imaxLocal + 2) + (i)] +#define V(i, j) v[(j) * (imaxLocal + 2) + (i)] +#define RHS(i, j) rhs[(j) * (imaxLocal + 2) + (i)] + #endif // __UTIL_H_ diff --git a/BasicSolver/2D-seq/src/grid.h b/BasicSolver/2D-seq/src/grid.h new file mode 100644 index 0000000..7e58d1a --- /dev/null +++ b/BasicSolver/2D-seq/src/grid.h @@ -0,0 +1,16 @@ +/* + * Copyright (C) NHR@FAU, University Erlangen-Nuremberg. + * All rights reserved. This file is part of nusif-solver. + * Use of this source code is governed by a MIT style + * license that can be found in the LICENSE file. + */ +#ifndef __GRID_H_ +#define __GRID_H_ + +typedef struct { + double dx, dy; + int imax, jmax; + double xlength, ylength; +} Grid; + +#endif // __GRID_H_