From f48927ce6f3bf799d4ffbd351a7625d91707f826 Mon Sep 17 00:00:00 2001 From: Aditya ujeniya Date: Wed, 10 Jul 2024 18:31:51 +0200 Subject: [PATCH] Completed 3D-seq and mpi refactoring in BasicSolver --- BasicSolver/3D-mpi/Makefile | 1 + BasicSolver/3D-mpi/canal.par | 11 +- BasicSolver/3D-mpi/config.mk | 4 +- BasicSolver/3D-mpi/dcavity.par | 7 + BasicSolver/3D-mpi/src/comm.c | 48 ++- BasicSolver/3D-mpi/src/comm.h | 4 + .../3D-mpi/src/{solver.c => discretization.c} | 336 +++++---------- BasicSolver/3D-mpi/src/discretization.h | 44 ++ BasicSolver/3D-mpi/src/main.c | 77 ++-- BasicSolver/3D-mpi/src/parameter.c | 28 +- BasicSolver/3D-mpi/src/parameter.h | 1 + BasicSolver/3D-mpi/src/solver-mg.c | 401 ++++++++++++++++++ BasicSolver/3D-mpi/src/solver-rb.c | 173 ++++++++ BasicSolver/3D-mpi/src/solver.h | 20 +- BasicSolver/3D-mpi/src/test.c | 22 +- BasicSolver/3D-mpi/src/util.h | 9 + BasicSolver/3D-mpi/src/vtkWriter-mpi.c | 25 +- BasicSolver/3D-mpi/src/vtkWriter-seq.c | 24 +- BasicSolver/3D-mpi/src/vtkWriter.h | 2 +- BasicSolver/3D-seq/canal.par | 7 + BasicSolver/3D-seq/config.mk | 6 +- BasicSolver/3D-seq/dcavity.par | 10 +- BasicSolver/3D-seq/src/discretization.c | 2 +- BasicSolver/3D-seq/src/main.c | 2 +- BasicSolver/3D-seq/src/parameter.c | 28 +- BasicSolver/3D-seq/src/parameter.h | 1 + BasicSolver/3D-seq/src/solver-mg.c | 86 +++- .../3D-seq/src/{solver-sor.c => solver-rb.c} | 0 BasicSolver/3D-seq/src/solver.h | 1 + BasicSolver/3D-seq/src/vtkWriter.c | 1 + 30 files changed, 1011 insertions(+), 370 deletions(-) rename BasicSolver/3D-mpi/src/{solver.c => discretization.c} (86%) create mode 100644 BasicSolver/3D-mpi/src/discretization.h create mode 100644 BasicSolver/3D-mpi/src/solver-mg.c create mode 100644 BasicSolver/3D-mpi/src/solver-rb.c rename BasicSolver/3D-seq/src/{solver-sor.c => solver-rb.c} (100%) diff --git a/BasicSolver/3D-mpi/Makefile b/BasicSolver/3D-mpi/Makefile index f61416c..d870e05 100644 --- a/BasicSolver/3D-mpi/Makefile +++ b/BasicSolver/3D-mpi/Makefile @@ -22,6 +22,7 @@ SRC = $(filter-out $(wildcard $(SRC_DIR)/*-*.c),$(wildcard $(SRC_DIR)/*.c) ASM = $(patsubst $(SRC_DIR)/%.c, $(BUILD_DIR)/%.s, $(SRC)) OBJ = $(patsubst $(SRC_DIR)/%.c, $(BUILD_DIR)/%.o, $(SRC)) OBJ += $(BUILD_DIR)/vtkWriter-$(VTK_OUTPUT_FMT).o +OBJ += $(BUILD_DIR)/solver-$(SOLVER).o SOURCES = $(SRC) $(wildcard $(SRC_DIR)/*.h) ifeq ($(VTK_OUTPUT_FMT),mpi) DEFINES += -D_VTK_WRITER_MPI diff --git a/BasicSolver/3D-mpi/canal.par b/BasicSolver/3D-mpi/canal.par index 3ff4a5f..597c5b5 100644 --- a/BasicSolver/3D-mpi/canal.par +++ b/BasicSolver/3D-mpi/canal.par @@ -38,15 +38,22 @@ kmax 50 # number of interior cells in z-direction # Time Data: # --------- -te 100.0 # final time +te 60.0 # final time dt 0.02 # time stepsize tau 0.5 # safety factor for time stepsize control (<0 constant delt) +# Multigrid data: +# --------- + +levels 3 # Multigrid levels +presmooth 5 # Pre-smoothning iterations +postsmooth 5 # Post-smoothning iterations + # Pressure Iteration Data: # ----------------------- itermax 500 # maximal number of pressure iteration in one time step eps 0.0001 # stopping tolerance for pressure iteration -omg 1.3 # relaxation parameter for SOR iteration +omg 1.7 # relaxation parameter for SOR iteration gamma 0.9 # upwind differencing factor gamma #=============================================================================== diff --git a/BasicSolver/3D-mpi/config.mk b/BasicSolver/3D-mpi/config.mk index 52154f4..0adc565 100644 --- a/BasicSolver/3D-mpi/config.mk +++ b/BasicSolver/3D-mpi/config.mk @@ -1,7 +1,9 @@ # Supported: GCC, CLANG, ICC -TAG ?= CLANG +TAG ?= ICC ENABLE_MPI ?= true ENABLE_OPENMP ?= false +# Supported: rb, mg +SOLVER ?= mg # Supported: seq, mpi VTK_OUTPUT_FMT ?= seq diff --git a/BasicSolver/3D-mpi/dcavity.par b/BasicSolver/3D-mpi/dcavity.par index c5cb201..303c292 100644 --- a/BasicSolver/3D-mpi/dcavity.par +++ b/BasicSolver/3D-mpi/dcavity.par @@ -42,6 +42,13 @@ te 10.0 # final time dt 0.02 # time stepsize tau 0.5 # safety factor for time stepsize control (<0 constant delt) +# Multigrid data: +# --------- + +levels 3 # Multigrid levels +presmooth 10 # Pre-smoothning iterations +postsmooth 5 # Post-smoothning iterations + # Pressure Iteration Data: # ----------------------- diff --git a/BasicSolver/3D-mpi/src/comm.c b/BasicSolver/3D-mpi/src/comm.c index 327bbbe..b11e9fd 100644 --- a/BasicSolver/3D-mpi/src/comm.c +++ b/BasicSolver/3D-mpi/src/comm.c @@ -336,7 +336,6 @@ void commCollectResult(Comm* c, int imaxLocal = c->imaxLocal; int jmaxLocal = c->jmaxLocal; int kmaxLocal = c->kmaxLocal; - #if defined(_MPI) int offset[c->size * NDIMS]; int imaxLocalAll[c->size]; @@ -562,9 +561,9 @@ void commPartition(Comm* c, int kmax, int jmax, int imax) MPI_Cart_shift(c->comm, KCORD, 1, &c->neighbours[FRONT], &c->neighbours[BACK]); MPI_Cart_get(c->comm, NCORDS, c->dims, periods, c->coords); - c->imaxLocal = sizeOfRank(c->rank, dims[ICORD], imax); - c->jmaxLocal = sizeOfRank(c->rank, dims[JCORD], jmax); - c->kmaxLocal = sizeOfRank(c->rank, dims[KCORD], kmax); + c->imaxLocal = sizeOfRank(c->coords[IDIM], dims[ICORD], imax); + c->jmaxLocal = sizeOfRank(c->coords[JDIM], dims[JCORD], jmax); + c->kmaxLocal = sizeOfRank(c->coords[KDIM], dims[KCORD], kmax); // setup buffer types for communication setupCommunication(c, LEFT, BULK); @@ -597,3 +596,44 @@ void commFinalize(Comm* c) MPI_Finalize(); #endif } + +void commUpdateDatatypes( + Comm* oldcomm, Comm* newcomm, int imaxLocal, int jmaxLocal, int kmaxLocal) +{ +#if defined _MPI + + int result = MPI_Comm_dup(oldcomm->comm, &newcomm->comm); + + if (result == MPI_ERR_COMM) { + printf("\nNull communicator. Duplication failed !!\n"); + } + + newcomm->imaxLocal = imaxLocal; + newcomm->jmaxLocal = jmaxLocal; + newcomm->kmaxLocal = kmaxLocal; + + setupCommunication(newcomm, LEFT, BULK); + setupCommunication(newcomm, LEFT, HALO); + setupCommunication(newcomm, RIGHT, BULK); + setupCommunication(newcomm, RIGHT, HALO); + setupCommunication(newcomm, BOTTOM, BULK); + setupCommunication(newcomm, BOTTOM, HALO); + setupCommunication(newcomm, TOP, BULK); + setupCommunication(newcomm, TOP, HALO); + setupCommunication(newcomm, FRONT, BULK); + setupCommunication(newcomm, FRONT, HALO); + setupCommunication(newcomm, BACK, BULK); + setupCommunication(newcomm, BACK, HALO); +#else + newcomm->imaxLocal = imaxLocal; + newcomm->jmaxLocal = jmaxLocal; + newcomm->kmaxLocal = kmaxLocal; +#endif +} + +void commFreeCommunicator(Comm* comm) +{ +#ifdef _MPI + MPI_Comm_free(&comm->comm); +#endif +} \ No newline at end of file diff --git a/BasicSolver/3D-mpi/src/comm.h b/BasicSolver/3D-mpi/src/comm.h index 48102fb..58ed4d7 100644 --- a/BasicSolver/3D-mpi/src/comm.h +++ b/BasicSolver/3D-mpi/src/comm.h @@ -23,6 +23,7 @@ typedef enum dimension { KDIM = 0, JDIM, IDIM, NDIMS } Dimension; enum layer { HALO = 0, BULK }; enum op { MAX = 0, SUM }; + typedef struct { int rank; int size; @@ -45,6 +46,9 @@ extern void commShift(Comm* c, double* f, double* g, double* h); extern void commReduction(double* v, int op); extern int commIsBoundary(Comm* c, Direction direction); extern void commGetOffsets(Comm* c, int offsets[], int kmax, int jmax, int imax); +extern void commFreeCommunicator(Comm* comm); +extern void commUpdateDatatypes( + Comm* oldcomm, Comm* newcomm, int imaxLocal, int jmaxLocal, int kmaxLocal); extern void commCollectResult(Comm* c, double* ug, double* vg, diff --git a/BasicSolver/3D-mpi/src/solver.c b/BasicSolver/3D-mpi/src/discretization.c similarity index 86% rename from BasicSolver/3D-mpi/src/solver.c rename to BasicSolver/3D-mpi/src/discretization.c index 121b4dc..138a037 100644 --- a/BasicSolver/3D-mpi/src/solver.c +++ b/BasicSolver/3D-mpi/src/discretization.c @@ -1,5 +1,5 @@ /* - * Copyright (C) NHR@FAU, University Erlangen-Nuremberg. + * 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. @@ -7,13 +7,11 @@ #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, k) \ @@ -33,7 +31,7 @@ #define RHS(i, j, k) \ rhs[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)] -static void printConfig(Solver* s) +static void printConfig(Discretization* s) { if (commIsMaster(&s->comm)) { printf("Parameters for #%s#\n", s->problem); @@ -72,7 +70,7 @@ static void printConfig(Solver* s) commPrintConfig(&s->comm); } -void initSolver(Solver* s, Parameter* params) +void initDiscretization(Discretization* s, Parameter* params) { s->problem = params->name; s->bcLeft = params->bcLeft; @@ -142,226 +140,7 @@ void initSolver(Solver* s, Parameter* params) #endif /* VERBOSE */ } -void computeRHS(Solver* s) -{ - int imaxLocal = s->comm.imaxLocal; - int jmaxLocal = s->comm.jmaxLocal; - int kmaxLocal = s->comm.kmaxLocal; - - double idx = 1.0 / s->grid.dx; - double idy = 1.0 / s->grid.dy; - double idz = 1.0 / s->grid.dz; - double idt = 1.0 / s->dt; - - double* rhs = s->rhs; - double* f = s->f; - double* g = s->g; - double* h = s->h; - - commShift(&s->comm, f, g, h); - - for (int k = 1; k < kmaxLocal + 1; k++) { - for (int j = 1; j < jmaxLocal + 1; j++) { - for (int i = 1; i < imaxLocal + 1; i++) { - RHS(i, j, k) = ((F(i, j, k) - F(i - 1, j, k)) * idx + - (G(i, j, k) - G(i, j - 1, k)) * idy + - (H(i, j, k) - H(i, j, k - 1)) * idz) * - idt; - } - } - } -} - -void solve(Solver* s) -{ - int imaxLocal = s->comm.imaxLocal; - int jmaxLocal = s->comm.jmaxLocal; - int kmaxLocal = s->comm.kmaxLocal; - - int imax = s->grid.imax; - int jmax = s->grid.jmax; - int kmax = s->grid.kmax; - - double eps = s->eps; - int itermax = s->itermax; - double dx2 = s->grid.dx * s->grid.dx; - double dy2 = s->grid.dy * s->grid.dy; - double dz2 = s->grid.dz * s->grid.dz; - double idx2 = 1.0 / dx2; - double idy2 = 1.0 / dy2; - double idz2 = 1.0 / dz2; - - double factor = s->omega * 0.5 * (dx2 * dy2 * dz2) / - (dy2 * dz2 + dx2 * dz2 + dx2 * dy2); - double* p = s->p; - double* rhs = s->rhs; - double epssq = eps * eps; - int it = 0; - double res = 1.0; - int pass, ksw, jsw, isw; - - while ((res >= epssq) && (it < itermax)) { - ksw = 1; - - for (pass = 0; pass < 2; pass++) { - jsw = ksw; - commExchange(&s->comm, p); - - for (int k = 1; k < kmaxLocal + 1; k++) { - isw = jsw; - for (int j = 1; j < jmaxLocal + 1; j++) { - for (int i = isw; i < imaxLocal + 1; i += 2) { - - double r = - RHS(i, j, k) - - ((P(i + 1, j, k) - 2.0 * P(i, j, k) + P(i - 1, j, k)) * idx2 + - (P(i, j + 1, k) - 2.0 * P(i, j, k) + P(i, j - 1, k)) * - idy2 + - (P(i, j, k + 1) - 2.0 * P(i, j, k) + P(i, j, k - 1)) * - idz2); - - P(i, j, k) -= (factor * r); - res += (r * r); - } - isw = 3 - isw; - } - jsw = 3 - jsw; - } - ksw = 3 - ksw; - } - - if (commIsBoundary(&s->comm, FRONT)) { - for (int j = 1; j < jmaxLocal + 1; j++) { - for (int i = 1; i < imaxLocal + 1; i++) { - P(i, j, 0) = P(i, j, 1); - } - } - } - - if (commIsBoundary(&s->comm, BACK)) { - for (int j = 1; j < jmaxLocal + 1; j++) { - for (int i = 1; i < imaxLocal + 1; i++) { - P(i, j, kmaxLocal + 1) = P(i, j, kmaxLocal); - } - } - } - - if (commIsBoundary(&s->comm, BOTTOM)) { - for (int k = 1; k < kmaxLocal + 1; k++) { - for (int i = 1; i < imaxLocal + 1; i++) { - P(i, 0, k) = P(i, 1, k); - } - } - } - - if (commIsBoundary(&s->comm, TOP)) { - for (int k = 1; k < kmaxLocal + 1; k++) { - for (int i = 1; i < imaxLocal + 1; i++) { - P(i, jmaxLocal + 1, k) = P(i, jmaxLocal, k); - } - } - } - - if (commIsBoundary(&s->comm, LEFT)) { - for (int k = 1; k < kmaxLocal + 1; k++) { - for (int j = 1; j < jmaxLocal + 1; j++) { - P(0, j, k) = P(1, j, k); - } - } - } - - if (commIsBoundary(&s->comm, RIGHT)) { - for (int k = 1; k < kmaxLocal + 1; k++) { - for (int j = 1; j < jmaxLocal + 1; j++) { - P(imaxLocal + 1, j, k) = P(imaxLocal, j, k); - } - } - } - - commReduction(&res, SUM); - res = res / (double)(imax * jmax * kmax); -#ifdef DEBUG - if (commIsMaster(&s->comm)) { - printf("%d Residuum: %e\n", it, res); - } -#endif - commExchange(&s->comm, p); - it++; - } - -#ifdef VERBOSE - 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 size = (s->comm.imaxLocal + 2) * (s->comm.jmaxLocal + 2) * - (s->comm.kmaxLocal + 2); - double maxval = DBL_MIN; - - for (int i = 0; i < size; i++) { - maxval = MAX(maxval, fabs(m[i])); - } - commReduction(&maxval, MAX); - return maxval; -} - -void normalizePressure(Solver* s) -{ - int imaxLocal = s->comm.imaxLocal; - int jmaxLocal = s->comm.jmaxLocal; - int kmaxLocal = s->comm.kmaxLocal; - - double* p = s->p; - double avgP = 0.0; - - for (int k = 1; k < kmaxLocal + 1; k++) { - for (int j = 1; j < jmaxLocal + 1; j++) { - for (int i = 1; i < imaxLocal + 1; i++) { - avgP += P(i, j, k); - } - } - } - commReduction(&avgP, SUM); - avgP /= (s->grid.imax * s->grid.jmax * s->grid.kmax); - - for (int k = 1; k < kmaxLocal + 1; k++) { - for (int j = 1; j < jmaxLocal + 1; j++) { - for (int i = 1; i < imaxLocal + 1; i++) { - P(i, j, k) = P(i, j, k) - avgP; - } - } - } -} - -void computeTimestep(Solver* s) -{ - double dt = s->dtBound; - double dx = s->grid.dx; - double dy = s->grid.dy; - double dz = s->grid.dz; - - double umax = maxElement(s, s->u); - double vmax = maxElement(s, s->v); - double wmax = maxElement(s, s->w); - - if (umax > 0) { - dt = (dt > dx / umax) ? dx / umax : dt; - } - if (vmax > 0) { - dt = (dt > dy / vmax) ? dy / vmax : dt; - } - if (wmax > 0) { - dt = (dt > dz / wmax) ? dz / wmax : dt; - } - - s->dt = dt * s->tau; -} - -void setBoundaryConditions(Solver* s) +void setBoundaryConditions(Discretization* s) { int imaxLocal = s->comm.imaxLocal; int jmaxLocal = s->comm.jmaxLocal; @@ -576,7 +355,37 @@ void setBoundaryConditions(Solver* s) } } -void setSpecialBoundaryCondition(Solver* s) +void computeRHS(Discretization* s) +{ + int imaxLocal = s->comm.imaxLocal; + int jmaxLocal = s->comm.jmaxLocal; + int kmaxLocal = s->comm.kmaxLocal; + + double idx = 1.0 / s->grid.dx; + double idy = 1.0 / s->grid.dy; + double idz = 1.0 / s->grid.dz; + double idt = 1.0 / s->dt; + + double* rhs = s->rhs; + double* f = s->f; + double* g = s->g; + double* h = s->h; + + commShift(&s->comm, f, g, h); + + for (int k = 1; k < kmaxLocal + 1; k++) { + for (int j = 1; j < jmaxLocal + 1; j++) { + for (int i = 1; i < imaxLocal + 1; i++) { + RHS(i, j, k) = ((F(i, j, k) - F(i - 1, j, k)) * idx + + (G(i, j, k) - G(i, j - 1, k)) * idy + + (H(i, j, k) - H(i, j, k - 1)) * idz) * + idt; + } + } + } +} + +void setSpecialBoundaryCondition(Discretization* s) { int imaxLocal = s->comm.imaxLocal; int jmaxLocal = s->comm.jmaxLocal; @@ -603,7 +412,74 @@ void setSpecialBoundaryCondition(Solver* s) } } -void computeFG(Solver* s) +static double maxElement(Discretization* s, double* m) +{ + int size = (s->comm.imaxLocal + 2) * (s->comm.jmaxLocal + 2) * + (s->comm.kmaxLocal + 2); + double maxval = DBL_MIN; + + for (int i = 0; i < size; i++) { + maxval = MAX(maxval, fabs(m[i])); + } + commReduction(&maxval, MAX); + return maxval; +} + +void normalizePressure(Discretization* s) +{ + int imaxLocal = s->comm.imaxLocal; + int jmaxLocal = s->comm.jmaxLocal; + int kmaxLocal = s->comm.kmaxLocal; + + double* p = s->p; + double avgP = 0.0; + + for (int k = 1; k < kmaxLocal + 1; k++) { + for (int j = 1; j < jmaxLocal + 1; j++) { + for (int i = 1; i < imaxLocal + 1; i++) { + avgP += P(i, j, k); + } + } + } + commReduction(&avgP, SUM); + avgP /= (s->grid.imax * s->grid.jmax * s->grid.kmax); + + for (int k = 1; k < kmaxLocal + 1; k++) { + for (int j = 1; j < jmaxLocal + 1; j++) { + for (int i = 1; i < imaxLocal + 1; i++) { + P(i, j, k) = P(i, j, k) - avgP; + } + } + } +} + +void computeTimestep(Discretization* s) +{ + double dt = s->dtBound; + double dx = s->grid.dx; + double dy = s->grid.dy; + double dz = s->grid.dz; + + double umax = maxElement(s, s->u); + double vmax = maxElement(s, s->v); + double wmax = maxElement(s, s->w); + + if (umax > 0) { + dt = (dt > dx / umax) ? dx / umax : dt; + } + if (vmax > 0) { + dt = (dt > dy / vmax) ? dy / vmax : dt; + } + if (wmax > 0) { + dt = (dt > dz / wmax) ? dz / wmax : dt; + } + + s->dt = dt * s->tau; +} + + + +void computeFG(Discretization* s) { int imaxLocal = s->comm.imaxLocal; int jmaxLocal = s->comm.jmaxLocal; @@ -823,7 +699,7 @@ void computeFG(Solver* s) } } -void adaptUV(Solver* s) +void adaptUV(Discretization* s) { int imaxLocal = s->comm.imaxLocal; int jmaxLocal = s->comm.jmaxLocal; @@ -850,4 +726,4 @@ void adaptUV(Solver* s) } } } -} +} \ No newline at end of file diff --git a/BasicSolver/3D-mpi/src/discretization.h b/BasicSolver/3D-mpi/src/discretization.h new file mode 100644 index 0000000..2c04800 --- /dev/null +++ b/BasicSolver/3D-mpi/src/discretization.h @@ -0,0 +1,44 @@ +/* + * 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 "grid.h" +#include "parameter.h" +#include "comm.h" + +enum BC { NOSLIP = 1, SLIP, OUTFLOW, PERIODIC }; + +typedef struct { + /* geometry and grid information */ + Grid grid; + /* arrays */ + double *p, *rhs; + double *f, *g, *h; + double *u, *v, *w; + /* parameters */ + double eps, omega; + double re, tau, gamma; + double gx, gy, gz; + /* time stepping */ + int itermax; + double dt, te; + double dtBound; + char* problem; + int bcLeft, bcRight, bcBottom, bcTop, bcFront, bcBack; + Comm comm; +} Discretization; + +extern void initDiscretization(Discretization*, Parameter*); +extern void computeRHS(Discretization*); +extern void normalizePressure(Discretization*); +extern void computeTimestep(Discretization*); +extern void setBoundaryConditions(Discretization*); +extern void setSpecialBoundaryCondition(Discretization*); +extern void computeFG(Discretization*); +extern void adaptUV(Discretization*); +#endif diff --git a/BasicSolver/3D-mpi/src/main.c b/BasicSolver/3D-mpi/src/main.c index 153abbf..944af32 100644 --- a/BasicSolver/3D-mpi/src/main.c +++ b/BasicSolver/3D-mpi/src/main.c @@ -9,7 +9,7 @@ #include #include "allocate.h" -#include "comm.h" +#include "discretization.h" #include "parameter.h" #include "progress.h" #include "solver.h" @@ -21,8 +21,9 @@ int main(int argc, char** argv) double timeStart, timeStop; Parameter p; Solver s; + Discretization d; - commInit(&s.comm, argc, argv); + commInit(&d.comm, argc, argv); initParameter(&p); if (argc != 2) { @@ -31,33 +32,40 @@ int main(int argc, char** argv) } readParameter(&p, argv[1]); - commPartition(&s.comm, p.kmax, p.jmax, p.imax); - if (commIsMaster(&s.comm)) { + commPartition(&d.comm, p.kmax, p.jmax, p.imax); + + if (commIsMaster(&d.comm)) { printParameter(&p); } - initSolver(&s, &p); + initDiscretization(&d, &p); + + initSolver(&s, &d, &p); + #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; + int nt = 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); + if (nt % 100 == 0) normalizePressure(&d); + solve(&s, d.p, d.rhs); + adaptUV(&d); + t += d.dt; + nt++; #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); @@ -67,7 +75,7 @@ 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); } @@ -75,14 +83,14 @@ int main(int argc, char** argv) #ifdef _VTK_WRITER_MPI VtkOptions opts = { .grid = s.grid, .comm = s.comm }; vtkOpen(&opts, s.problem); - vtkScalar(&opts, "pressure", s.p); - vtkVector(&opts, "velocity", (VtkVector) { s.u, s.v, s.w }); + vtkScalar(&opts, "pressure", d.p); + vtkVector(&opts, "velocity", (VtkVector) { d.u, d.v, d.w }); vtkClose(&opts); #else double *pg, *ug, *vg, *wg; - if (commIsMaster(&s.comm)) { - size_t bytesize = s.grid.imax * s.grid.jmax * s.grid.kmax * sizeof(double); + if (commIsMaster(s.comm)) { + size_t bytesize = s.grid->imax * s.grid->jmax * s.grid->kmax * sizeof(double); pg = allocate(64, bytesize); ug = allocate(64, bytesize); @@ -90,34 +98,35 @@ int main(int argc, char** argv) wg = allocate(64, bytesize); } - commCollectResult(&s.comm, + commCollectResult(s.comm, ug, vg, wg, pg, - s.u, - s.v, - s.w, - s.p, - s.grid.kmax, - s.grid.jmax, - s.grid.imax); + d.u, + d.v, + d.w, + d.p, + s.grid->kmax, + s.grid->jmax, + s.grid->imax); - if (commIsMaster(&s.comm)) { + if (commIsMaster(s.comm)) { VtkOptions opts = { .grid = s.grid }; vtkOpen(&opts, s.problem); vtkScalar(&opts, "pressure", pg); vtkVector(&opts, "velocity", (VtkVector) { ug, vg, wg }); vtkClose(&opts); } + #endif timeStop = getTimeStamp(); - if (commIsMaster(&s.comm)) { + if (commIsMaster(s.comm)) { printf("Result output took %.2fs\n", timeStop - timeStart); } - commFinalize(&s.comm); + commFinalize(s.comm); return EXIT_SUCCESS; } diff --git a/BasicSolver/3D-mpi/src/parameter.c b/BasicSolver/3D-mpi/src/parameter.c index 2128c2a..c6cde26 100644 --- a/BasicSolver/3D-mpi/src/parameter.c +++ b/BasicSolver/3D-mpi/src/parameter.c @@ -14,18 +14,22 @@ void initParameter(Parameter* param) { - param->xlength = 1.0; - param->ylength = 1.0; - param->zlength = 1.0; - param->imax = 100; - param->jmax = 100; - param->kmax = 100; - param->itermax = 1000; - param->eps = 0.0001; - param->omg = 1.7; - param->re = 100.0; - param->gamma = 0.9; - param->tau = 0.5; + param->xlength = 1.0; + param->ylength = 1.0; + param->zlength = 1.0; + param->imax = 100; + param->jmax = 100; + param->kmax = 100; + param->itermax = 1000; + param->eps = 0.0001; + param->omg = 1.7; + param->re = 100.0; + param->gamma = 0.9; + param->tau = 0.5; + param->levels = 5; + param->presmooth = 5; + param->postsmooth = 5; + } void readParameter(Parameter* param, const char* filename) diff --git a/BasicSolver/3D-mpi/src/parameter.h b/BasicSolver/3D-mpi/src/parameter.h index d946d5c..e40b9a4 100644 --- a/BasicSolver/3D-mpi/src/parameter.h +++ b/BasicSolver/3D-mpi/src/parameter.h @@ -18,6 +18,7 @@ typedef struct { char* name; int bcLeft, bcRight, bcBottom, bcTop, bcFront, bcBack; double u_init, v_init, w_init, p_init; + int levels, presmooth, postsmooth; } Parameter; void initParameter(Parameter*); diff --git a/BasicSolver/3D-mpi/src/solver-mg.c b/BasicSolver/3D-mpi/src/solver-mg.c new file mode 100644 index 0000000..d3f532a --- /dev/null +++ b/BasicSolver/3D-mpi/src/solver-mg.c @@ -0,0 +1,401 @@ +/* + * 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 "allocate.h" +#include "solver.h" +#include "util.h" + +#define FINEST_LEVEL 0 +#define COARSEST_LEVEL (s->levels - 1) +#define S(i, j, k) \ + s[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)] +#define E(i, j, k) \ + e[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)] +#define R(i, j, k) \ + r[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)] +#define OLD(i, j, k) \ + old[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)] + +static void restrictMG(Solver* s, int level, Comm* comm) +{ + double* r = s->r[level + 1]; + double* old = s->r[level]; + + int imaxLocal = comm->imaxLocal; + int jmaxLocal = comm->jmaxLocal; + int kmaxLocal = comm->kmaxLocal; + + commExchange(comm, old); + + for (int k = 1; k < comm->kmaxLocal + 1; k++) { + for (int j = 1; j < comm->jmaxLocal + 1; j++) { + for (int i = 1; i < comm->imaxLocal + 1; ++i) { + R(i, j, k) = (OLD(2 * i - 1, 2 * j - 1, 2 * k) + + OLD(2 * i, 2 * j - 1, 2 * k) * 2 + + OLD(2 * i + 1, 2 * j - 1, 2 * k) + + OLD(2 * i - 1, 2 * j, 2 * k) * 2 + + OLD(2 * i, 2 * j, 2 * k) * 8 + + OLD(2 * i + 1, 2 * j, 2 * k) * 2 + + OLD(2 * i - 1, 2 * j + 1, 2 * k) + + OLD(2 * i, 2 * j + 1, 2 * k) * 2 + + OLD(2 * i + 1, 2 * j + 1, 2 * k) + + + OLD(2 * i - 1, 2 * j - 1, 2 * k - 1) + + OLD(2 * i, 2 * j - 1, 2 * k - 1) * 2 + + OLD(2 * i + 1, 2 * j - 1, 2 * k - 1) + + OLD(2 * i - 1, 2 * j, 2 * k - 1) * 2 + + OLD(2 * i, 2 * j, 2 * k - 1) * 4 + + OLD(2 * i + 1, 2 * j, 2 * k - 1) * 2 + + OLD(2 * i - 1, 2 * j + 1, 2 * k - 1) + + OLD(2 * i, 2 * j + 1, 2 * k - 1) * 2 + + OLD(2 * i + 1, 2 * j + 1, 2 * k - 1) + + + OLD(2 * i - 1, 2 * j - 1, 2 * k + 1) + + OLD(2 * i, 2 * j - 1, 2 * k + 1) * 2 + + OLD(2 * i + 1, 2 * j - 1, 2 * k + 1) + + OLD(2 * i - 1, 2 * j, 2 * k + 1) * 2 + + OLD(2 * i, 2 * j, 2 * k + 1) * 4 + + OLD(2 * i + 1, 2 * j, 2 * k + 1) * 2 + + OLD(2 * i - 1, 2 * j + 1, 2 * k + 1) + + OLD(2 * i, 2 * j + 1, 2 * k + 1) * 2 + + OLD(2 * i + 1, 2 * j + 1, 2 * k + 1)) / + 64.0; + } + } + } +} + +static void prolongate(Solver* s, int level, Comm* comm) +{ + double* old = s->r[level + 1]; + double* e = s->r[level]; + + int imaxLocal = comm->imaxLocal; + int jmaxLocal = comm->jmaxLocal; + int kmaxLocal = comm->kmaxLocal; + + for (int k = 2; k < kmaxLocal + 1; k += 2) { + for (int j = 2; j < jmaxLocal + 1; j += 2) { + for (int i = 2; i < imaxLocal + 1; i += 2) { + E(i, j, k) = OLD(i / 2, j / 2, k / 2); + } + } + } +} + +static void correct(Solver* s, double* p, int level, Comm* comm) +{ + double* e = s->e[level]; + + int imaxLocal = comm->imaxLocal; + int jmaxLocal = comm->jmaxLocal; + int kmaxLocal = comm->kmaxLocal; + + for (int k = 1; k < kmaxLocal + 1; ++k) { + for (int j = 1; j < jmaxLocal + 1; ++j) { + for (int i = 1; i < imaxLocal + 1; ++i) { + P(i, j, k) += E(i, j, k); + } + } + } +} + +static void setBoundaryCondition( + Solver* s, double* p, int imaxLocal, int jmaxLocal, int kmaxLocal) +{ +#ifdef _MPI + if (commIsBoundary(s->comm, FRONT)) { + for (int j = 1; j < jmaxLocal + 1; j++) { + for (int i = 1; i < imaxLocal + 1; i++) { + P(i, j, 0) = P(i, j, 1); + } + } + } + + if (commIsBoundary(s->comm, BACK)) { + for (int j = 1; j < jmaxLocal + 1; j++) { + for (int i = 1; i < imaxLocal + 1; i++) { + P(i, j, kmaxLocal + 1) = P(i, j, kmaxLocal); + } + } + } + + if (commIsBoundary(s->comm, BOTTOM)) { + for (int k = 1; k < kmaxLocal + 1; k++) { + for (int i = 1; i < imaxLocal + 1; i++) { + P(i, 0, k) = P(i, 1, k); + } + } + } + + if (commIsBoundary(s->comm, TOP)) { + for (int k = 1; k < kmaxLocal + 1; k++) { + for (int i = 1; i < imaxLocal + 1; i++) { + P(i, jmaxLocal + 1, k) = P(i, jmaxLocal, k); + } + } + } + + if (commIsBoundary(s->comm, LEFT)) { + for (int k = 1; k < kmaxLocal + 1; k++) { + for (int j = 1; j < jmaxLocal + 1; j++) { + P(0, j, k) = P(1, j, k); + } + } + } + + if (commIsBoundary(s->comm, RIGHT)) { + for (int k = 1; k < kmaxLocal + 1; k++) { + for (int j = 1; j < jmaxLocal + 1; j++) { + P(imaxLocal + 1, j, k) = P(imaxLocal, j, k); + } + } + } +#else + for (int j = 1; j < jmaxLocal + 1; j++) { + for (int i = 1; i < imaxLocal + 1; i++) { + P(i, j, 0) = P(i, j, 1); + P(i, j, kmaxLocal + 1) = P(i, j, kmaxLocal); + } + } + + for (int k = 1; k < kmaxLocal + 1; k++) { + for (int i = 1; i < imaxLocal + 1; i++) { + P(i, 0, k) = P(i, 1, k); + P(i, jmaxLocal + 1, k) = P(i, jmaxLocal, k); + } + } + + for (int k = 1; k < kmaxLocal + 1; k++) { + for (int j = 1; j < jmaxLocal + 1; j++) { + P(0, j, k) = P(1, j, k); + P(imaxLocal + 1, j, k) = P(imaxLocal, j, k); + } + } +#endif +} + +static double smooth(Solver* s, double* p, double* rhs, int level, Comm* comm) +{ + int imaxLocal = comm->imaxLocal; + int jmaxLocal = comm->jmaxLocal; + int kmaxLocal = comm->kmaxLocal; + + double eps = s->eps; + int itermax = s->itermax; + double dx2 = s->grid->dx * s->grid->dx; + double dy2 = s->grid->dy * s->grid->dy; + double dz2 = s->grid->dz * s->grid->dz; + double idx2 = 1.0 / dx2; + double idy2 = 1.0 / dy2; + double idz2 = 1.0 / dz2; + double factor = s->omega * 0.5 * (dx2 * dy2 * dz2) / + (dy2 * dz2 + dx2 * dz2 + dx2 * dy2); + double* r = s->r[level]; + double epssq = eps * eps; + int it = 0; + int pass, ksw, jsw, isw; + double res = 1.0; + + ksw = 1; + + for (pass = 0; pass < 2; pass++) { + jsw = ksw; + + commExchange(comm, p); + + for (int k = 1; k < kmaxLocal + 1; k++) { + isw = jsw; + for (int j = 1; j < jmaxLocal + 1; j++) { + for (int i = isw; i < imaxLocal + 1; i += 2) { + + P(i, j, k) -= + factor * + (RHS(i, j, k) - + ((P(i + 1, j, k) - 2.0 * P(i, j, k) + P(i - 1, j, k)) * idx2 + + (P(i, j + 1, k) - 2.0 * P(i, j, k) + P(i, j - 1, k)) * + idy2 + + (P(i, j, k + 1) - 2.0 * P(i, j, k) + P(i, j, k - 1)) * + idz2)); + } + isw = 3 - isw; + } + jsw = 3 - jsw; + } + ksw = 3 - ksw; + } +} + +static double calculateResidual(Solver* s, double* p, double* rhs, int level, Comm* comm) +{ + + int imaxLocal = comm->imaxLocal; + int jmaxLocal = comm->jmaxLocal; + int kmaxLocal = comm->kmaxLocal; + + double eps = s->eps; + int itermax = s->itermax; + double dx2 = s->grid->dx * s->grid->dx; + double dy2 = s->grid->dy * s->grid->dy; + double dz2 = s->grid->dz * s->grid->dz; + double idx2 = 1.0 / dx2; + double idy2 = 1.0 / dy2; + double idz2 = 1.0 / dz2; + double factor = s->omega * 0.5 * (dx2 * dy2 * dz2) / + (dy2 * dz2 + dx2 * dz2 + dx2 * dy2); + double* r = s->r[level]; + double epssq = eps * eps; + int it = 0; + int pass, ksw, jsw, isw; + double res = 1.0; + + ksw = 1; + + for (pass = 0; pass < 2; pass++) { + jsw = ksw; + + commExchange(comm, p); + + for (int k = 1; k < kmaxLocal + 1; k++) { + isw = jsw; + for (int j = 1; j < jmaxLocal + 1; j++) { + for (int i = isw; i < imaxLocal + 1; i += 2) { + + R(i, + j, + k) = (RHS(i, j, k) - + ((P(i + 1, j, k) - 2.0 * P(i, j, k) + P(i - 1, j, k)) * + idx2 + + (P(i, j + 1, k) - 2.0 * P(i, j, k) + P(i, j - 1, k)) * + idy2 + + (P(i, j, k + 1) - 2.0 * P(i, j, k) + P(i, j, k - 1)) * + idz2)); + + res += (R(i, j, k) * R(i, j, k)); + } + isw = 3 - isw; + } + jsw = 3 - jsw; + } + ksw = 3 - ksw; + } + + commReduction(&res, SUM); + + res = res / (double)(imaxLocal * jmaxLocal * kmaxLocal); +#ifdef DEBUG + if (commIsMaster(s->comm)) { + printf("%d Residuum: %e\n", it, res); + } +#endif + return res; +} + +static double multiGrid(Solver* s, double* p, double* rhs, int level, Comm* comm) +{ + int imaxLocal = comm->imaxLocal; + int jmaxLocal = comm->jmaxLocal; + int kmaxLocal = comm->kmaxLocal; + + double res = 0.0; + + // coarsest level + if (level == COARSEST_LEVEL) { + for (int i = 0; i < 5; i++) { + smooth(s, p, rhs, level, comm); + } + return res; + } + + // pre-smoothing + for (int i = 0; i < s->presmooth; i++) { + smooth(s, p, rhs, level, comm); + if (level == FINEST_LEVEL) + setBoundaryCondition(s, p, imaxLocal, jmaxLocal, kmaxLocal); + } + + res = calculateResidual(s, p, rhs, level, comm); + + // restrict + restrictMG(s, level, comm); + + Comm newcomm; + commUpdateDatatypes(s->comm, + &newcomm, + comm->imaxLocal / 2, + comm->jmaxLocal / 2, + comm->kmaxLocal / 2); + + // MGSolver on residual and error. + multiGrid(s, s->e[level + 1], s->r[level + 1], level + 1, &newcomm); + + commFreeCommunicator(&newcomm); + + // prolongate + prolongate(s, level, comm); + + // correct p on finer level using residual + correct(s, p, level, comm); + if (level == FINEST_LEVEL) + setBoundaryCondition(s, p, imaxLocal, jmaxLocal, kmaxLocal); + + // post-smoothing + for (int i = 0; i < s->postsmooth; i++) { + smooth(s, p, rhs, level, comm); + if (level == FINEST_LEVEL) + setBoundaryCondition(s, p, imaxLocal, jmaxLocal, kmaxLocal); + } + + return res; +} + +void initSolver(Solver* s, Discretization* d, Parameter* p) +{ + s->eps = p->eps; + s->omega = p->omg; + s->itermax = p->itermax; + s->levels = p->levels; + s->grid = &d->grid; + s->presmooth = p->presmooth; + s->postsmooth = p->postsmooth; + s->comm = &d->comm; + s->problem = p->name; + + int imax = s->grid->imax; + int jmax = s->grid->jmax; + int kmax = s->grid->kmax; + int levels = s->levels; + printf("Using Multigrid solver with %d levels\n", levels); + + s->r = malloc(levels * sizeof(double*)); + s->e = malloc(levels * sizeof(double*)); + + size_t size = (imax + 2) * (jmax + 2) * (kmax + 2); + + for (int j = 0; j < levels; j++) { + s->r[j] = allocate(64, size * sizeof(double)); + s->e[j] = allocate(64, size * sizeof(double)); + + for (size_t i = 0; i < size; i++) { + s->r[j][i] = 0.0; + s->e[j][i] = 0.0; + } + } +} + +void solve(Solver* s, double* p, double* rhs) +{ + double res = multiGrid(s, p, rhs, 0, s->comm); + +#ifdef VERBOSE + if (commIsMaster(s->comm)) { + printf("Residuum: %.6f\n", res); + } +#endif +} diff --git a/BasicSolver/3D-mpi/src/solver-rb.c b/BasicSolver/3D-mpi/src/solver-rb.c new file mode 100644 index 0000000..54148e9 --- /dev/null +++ b/BasicSolver/3D-mpi/src/solver-rb.c @@ -0,0 +1,173 @@ +/* + * 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 "parameter.h" +#include "solver.h" +#include "util.h" + +void solve(Solver* s, double* p, double* rhs) +{ + int imaxLocal = s->comm->imaxLocal; + int jmaxLocal = s->comm->jmaxLocal; + int kmaxLocal = s->comm->kmaxLocal; + + int imax = s->grid->imax; + int jmax = s->grid->jmax; + int kmax = s->grid->kmax; + + double eps = s->eps; + int itermax = s->itermax; + double dx2 = s->grid->dx * s->grid->dx; + double dy2 = s->grid->dy * s->grid->dy; + double dz2 = s->grid->dz * s->grid->dz; + double idx2 = 1.0 / dx2; + double idy2 = 1.0 / dy2; + double idz2 = 1.0 / dz2; + + double factor = s->omega * 0.5 * (dx2 * dy2 * dz2) / + (dy2 * dz2 + dx2 * dz2 + dx2 * dy2); + double epssq = eps * eps; + int it = 0; + double res = 1.0; + int pass, ksw, jsw, isw; + + while ((res >= epssq) && (it < itermax)) { + ksw = 1; + + for (pass = 0; pass < 2; pass++) { + jsw = ksw; + commExchange(s->comm, p); + + for (int k = 1; k < kmaxLocal + 1; k++) { + isw = jsw; + for (int j = 1; j < jmaxLocal + 1; j++) { + for (int i = isw; i < imaxLocal + 1; i += 2) { + + double r = + RHS(i, j, k) - + ((P(i + 1, j, k) - 2.0 * P(i, j, k) + P(i - 1, j, k)) * idx2 + + (P(i, j + 1, k) - 2.0 * P(i, j, k) + P(i, j - 1, k)) * + idy2 + + (P(i, j, k + 1) - 2.0 * P(i, j, k) + P(i, j, k - 1)) * + idz2); + + P(i, j, k) -= (factor * r); + res += (r * r); + } + isw = 3 - isw; + } + jsw = 3 - jsw; + } + ksw = 3 - ksw; + } +#ifdef _MPI + if (commIsBoundary(s->comm, FRONT)) { + for (int j = 1; j < jmaxLocal + 1; j++) { + for (int i = 1; i < imaxLocal + 1; i++) { + P(i, j, 0) = P(i, j, 1); + } + } + } + + if (commIsBoundary(s->comm, BACK)) { + for (int j = 1; j < jmaxLocal + 1; j++) { + for (int i = 1; i < imaxLocal + 1; i++) { + P(i, j, kmaxLocal + 1) = P(i, j, kmaxLocal); + } + } + } + + if (commIsBoundary(s->comm, BOTTOM)) { + for (int k = 1; k < kmaxLocal + 1; k++) { + for (int i = 1; i < imaxLocal + 1; i++) { + P(i, 0, k) = P(i, 1, k); + } + } + } + + if (commIsBoundary(s->comm, TOP)) { + for (int k = 1; k < kmaxLocal + 1; k++) { + for (int i = 1; i < imaxLocal + 1; i++) { + P(i, jmaxLocal + 1, k) = P(i, jmaxLocal, k); + } + } + } + + if (commIsBoundary(s->comm, LEFT)) { + for (int k = 1; k < kmaxLocal + 1; k++) { + for (int j = 1; j < jmaxLocal + 1; j++) { + P(0, j, k) = P(1, j, k); + } + } + } + + if (commIsBoundary(s->comm, RIGHT)) { + for (int k = 1; k < kmaxLocal + 1; k++) { + for (int j = 1; j < jmaxLocal + 1; j++) { + P(imaxLocal + 1, j, k) = P(imaxLocal, j, k); + } + } + } +#else + for (int j = 1; j < jmax + 1; j++) { + for (int i = 1; i < imax + 1; i++) { + P(i, j, 0) = P(i, j, 1); + P(i, j, kmax + 1) = P(i, j, kmax); + } + } + + for (int k = 1; k < kmax + 1; k++) { + for (int i = 1; i < imax + 1; i++) { + P(i, 0, k) = P(i, 1, k); + P(i, jmax + 1, k) = P(i, jmax, k); + } + } + + for (int k = 1; k < kmax + 1; k++) { + for (int j = 1; j < jmax + 1; j++) { + P(0, j, k) = P(1, j, k); + P(imax + 1, j, k) = P(imax, j, k); + } + } +#endif + commReduction(&res, SUM); + res = res / (double)(imax * jmax * kmax); +#ifdef DEBUG + if (commIsMaster(&s->comm)) { + printf("%d Residuum: %e\n", it, res); + } +#endif + commExchange(s->comm, p); + it++; + } + +#ifdef VERBOSE + if (commIsMaster(s->comm)) { + printf("Solver took %d iterations to reach %f\n", it, sqrt(res)); + } +#endif +} + +void initSolver(Solver* s, Discretization* d, Parameter* p) +{ + s->eps = p->eps; + s->omega = p->omg; + s->itermax = p->itermax; + s->levels = p->levels; + s->grid = &d->grid; + s->presmooth = p->presmooth; + s->postsmooth = p->postsmooth; + s->comm = &d->comm; + s->problem = p->name; +} \ No newline at end of file diff --git a/BasicSolver/3D-mpi/src/solver.h b/BasicSolver/3D-mpi/src/solver.h index ae734f7..09b1374 100644 --- a/BasicSolver/3D-mpi/src/solver.h +++ b/BasicSolver/3D-mpi/src/solver.h @@ -9,12 +9,11 @@ #include "comm.h" #include "grid.h" #include "parameter.h" - -enum BC { NOSLIP = 1, SLIP, OUTFLOW, PERIODIC }; +#include "discretization.h" typedef struct { /* geometry and grid information */ - Grid grid; + Grid* grid; /* arrays */ double *p, *rhs; double *f, *g, *h; @@ -30,16 +29,11 @@ typedef struct { char* problem; int bcLeft, bcRight, bcBottom, bcTop, bcFront, bcBack; /* communication */ - Comm comm; + double **r, **e; + int levels, presmooth, postsmooth; + Comm* comm; } Solver; -extern void initSolver(Solver*, Parameter*); -extern void computeRHS(Solver*); -extern void solve(Solver*); -extern void normalizePressure(Solver*); -extern void computeTimestep(Solver*); -extern void setBoundaryConditions(Solver*); -extern void setSpecialBoundaryCondition(Solver*); -extern void computeFG(Solver*); -extern void adaptUV(Solver*); +extern void solve(Solver* , double* , double* ); +extern void initSolver(Solver*, Discretization*, Parameter*); #endif diff --git a/BasicSolver/3D-mpi/src/test.c b/BasicSolver/3D-mpi/src/test.c index 4f27776..ba6bc78 100644 --- a/BasicSolver/3D-mpi/src/test.c +++ b/BasicSolver/3D-mpi/src/test.c @@ -14,10 +14,10 @@ void testInit(Solver* s) { - int imaxLocal = s->comm.imaxLocal; - int jmaxLocal = s->comm.jmaxLocal; - int kmaxLocal = s->comm.kmaxLocal; - int myrank = s->comm.rank; + int imaxLocal = s->comm->imaxLocal; + int jmaxLocal = s->comm->jmaxLocal; + int kmaxLocal = s->comm->kmaxLocal; + int myrank = s->comm->rank; double* p = s->p; double* f = s->f; double* g = s->g; @@ -76,11 +76,11 @@ static char* direction2String(Direction dir) static void printPlane(Solver* s, double* a, int ymax, int xmax, Direction dir) { - int imaxLocal = s->comm.imaxLocal; - int jmaxLocal = s->comm.jmaxLocal; - int kmaxLocal = s->comm.kmaxLocal; + int imaxLocal = s->comm->imaxLocal; + int jmaxLocal = s->comm->jmaxLocal; + int kmaxLocal = s->comm->kmaxLocal; char filename[50]; - snprintf(filename, 50, "halo-%s-r%d.txt", direction2String(dir), s->comm.rank); + snprintf(filename, 50, "halo-%s-r%d.txt", direction2String(dir), s->comm->rank); FILE* fh = fopen(filename, "w"); for (int y = 0; y < ymax; y++) { @@ -116,9 +116,9 @@ static void printPlane(Solver* s, double* a, int ymax, int xmax, Direction dir) void testPrintHalo(Solver* s, double* a) { - int imaxLocal = s->comm.imaxLocal; - int jmaxLocal = s->comm.jmaxLocal; - int kmaxLocal = s->comm.kmaxLocal; + int imaxLocal = s->comm->imaxLocal; + int jmaxLocal = s->comm->jmaxLocal; + int kmaxLocal = s->comm->kmaxLocal; printPlane(s, a, kmaxLocal + 2, imaxLocal + 2, BOTTOM); printPlane(s, a, kmaxLocal + 2, imaxLocal + 2, TOP); diff --git a/BasicSolver/3D-mpi/src/util.h b/BasicSolver/3D-mpi/src/util.h index c27b2ba..68b87f0 100644 --- a/BasicSolver/3D-mpi/src/util.h +++ b/BasicSolver/3D-mpi/src/util.h @@ -19,4 +19,13 @@ #define ABS(a) ((a) >= 0 ? (a) : -(a)) #endif +#define P(i, j, k) p[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)] +#define F(i, j, k) f[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)] +#define G(i, j, k) g[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)] +#define H(i, j, k) h[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)] +#define U(i, j, k) u[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)] +#define V(i, j, k) v[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)] +#define W(i, j, k) w[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)] +#define RHS(i, j, k) rhs[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)] + #endif // __UTIL_H_ diff --git a/BasicSolver/3D-mpi/src/vtkWriter-mpi.c b/BasicSolver/3D-mpi/src/vtkWriter-mpi.c index 8b29c56..d9670d5 100644 --- a/BasicSolver/3D-mpi/src/vtkWriter-mpi.c +++ b/BasicSolver/3D-mpi/src/vtkWriter-mpi.c @@ -46,18 +46,18 @@ static void writeHeader(VtkOptions* o) cursor += sprintf(cursor, "DATASET STRUCTURED_POINTS\n"); cursor += sprintf(cursor, "DIMENSIONS %d %d %d\n", - o->grid.imax, - o->grid.jmax, - o->grid.kmax); + o->grid->imax, + o->grid->jmax, + o->grid->kmax); cursor += sprintf(cursor, "ORIGIN %f %f %f\n", - o->grid.dx * 0.5, - o->grid.dy * 0.5, - o->grid.dz * 0.5); - cursor += sprintf(cursor, "SPACING %f %f %f\n", o->grid.dx, o->grid.dy, o->grid.dz); + o->grid->dx * 0.5, + o->grid->dy * 0.5, + o->grid->dz * 0.5); + cursor += sprintf(cursor, "SPACING %f %f %f\n", o->grid->dx, o->grid->dy, o->grid->dz); cursor += sprintf(cursor, "POINT_DATA %d\n", - o->grid.imax * o->grid.jmax * o->grid.kmax); + o->grid->imax * o->grid->jmax * o->grid->kmax); if (commIsMaster(&o->comm)) { MPI_File_write(o->fh, header, (int)strlen(header), MPI_CHAR, MPI_STATUS_IGNORE); @@ -67,7 +67,6 @@ static void writeHeader(VtkOptions* o) void vtkOpen(VtkOptions* o, char* problem) { char filename[50]; - snprintf(filename, 50, "%s-p%d.vtk", problem, o->comm.size); MPI_File_open(o->comm.comm, filename, @@ -98,7 +97,7 @@ void vtkScalar(VtkOptions* o, char* name, double* s) } int offsets[NDIMS]; - commGetOffsets(&o->comm, offsets, o->grid.kmax, o->grid.jmax, o->grid.imax); + commGetOffsets(&o->comm, offsets, o->grid->kmax, o->grid->jmax, o->grid->imax); // set global view in file MPI_Offset disp; @@ -108,7 +107,7 @@ void vtkScalar(VtkOptions* o, char* name, double* s) MPI_File_get_size(o->fh, &disp); MPI_Type_create_subarray(NDIMS, - (int[NDIMS]) { o->grid.kmax, o->grid.jmax, o->grid.imax }, + (int[NDIMS]) { o->grid->kmax, o->grid->jmax, o->grid->imax }, (int[NDIMS]) { o->comm.kmaxLocal, o->comm.jmaxLocal, o->comm.imaxLocal }, offsets, MPI_ORDER_C, @@ -177,7 +176,7 @@ void vtkVector(VtkOptions* o, char* name, VtkVector vec) } int offsets[NDIMS]; - commGetOffsets(&o->comm, offsets, o->grid.kmax, o->grid.jmax, o->grid.imax); + commGetOffsets(&o->comm, offsets, o->grid->kmax, o->grid->jmax, o->grid->imax); // set global view in file MPI_Offset disp; @@ -190,7 +189,7 @@ void vtkVector(VtkOptions* o, char* name, VtkVector vec) MPI_Type_commit(&vectorType); MPI_Type_create_subarray(NDIMS, - (int[NDIMS]) { o->grid.kmax, o->grid.jmax, o->grid.imax }, + (int[NDIMS]) { o->grid->kmax, o->grid->jmax, o->grid->imax }, (int[NDIMS]) { kmaxLocal, jmaxLocal, imaxLocal }, offsets, MPI_ORDER_C, diff --git a/BasicSolver/3D-mpi/src/vtkWriter-seq.c b/BasicSolver/3D-mpi/src/vtkWriter-seq.c index 7bd8a14..4bb4fbb 100644 --- a/BasicSolver/3D-mpi/src/vtkWriter-seq.c +++ b/BasicSolver/3D-mpi/src/vtkWriter-seq.c @@ -41,14 +41,14 @@ static void writeHeader(VtkOptions* o) } fprintf(o->fh, "DATASET STRUCTURED_POINTS\n"); - fprintf(o->fh, "DIMENSIONS %d %d %d\n", o->grid.imax, o->grid.jmax, o->grid.kmax); + fprintf(o->fh, "DIMENSIONS %d %d %d\n", o->grid->imax, o->grid->jmax, o->grid->kmax); fprintf(o->fh, "ORIGIN %f %f %f\n", - o->grid.dx * 0.5, - o->grid.dy * 0.5, - o->grid.dz * 0.5); - fprintf(o->fh, "SPACING %f %f %f\n", o->grid.dx, o->grid.dy, o->grid.dz); - fprintf(o->fh, "POINT_DATA %d\n", o->grid.imax * o->grid.jmax * o->grid.kmax); + o->grid->dx * 0.5, + o->grid->dy * 0.5, + o->grid->dz * 0.5); + fprintf(o->fh, "SPACING %f %f %f\n", o->grid->dx, o->grid->dy, o->grid->dz); + fprintf(o->fh, "POINT_DATA %d\n", o->grid->imax * o->grid->jmax * o->grid->kmax); } void vtkOpen(VtkOptions* o, char* problem) @@ -64,9 +64,9 @@ void vtkOpen(VtkOptions* o, char* problem) static void writeScalar(VtkOptions* o, double* s) { - int imax = o->grid.imax; - int jmax = o->grid.jmax; - int kmax = o->grid.kmax; + int imax = o->grid->imax; + int jmax = o->grid->jmax; + int kmax = o->grid->kmax; for (int k = 0; k < kmax; k++) { for (int j = 0; j < jmax; j++) { @@ -105,9 +105,9 @@ void vtkScalar(VtkOptions* o, char* name, double* s) static void writeVector(VtkOptions* o, VtkVector vec) { - int imax = o->grid.imax; - int jmax = o->grid.jmax; - int kmax = o->grid.kmax; + int imax = o->grid->imax; + int jmax = o->grid->jmax; + int kmax = o->grid->kmax; for (int k = 0; k < kmax; k++) { for (int j = 0; j < jmax; j++) { diff --git a/BasicSolver/3D-mpi/src/vtkWriter.h b/BasicSolver/3D-mpi/src/vtkWriter.h index 4b4c2fe..e93835c 100644 --- a/BasicSolver/3D-mpi/src/vtkWriter.h +++ b/BasicSolver/3D-mpi/src/vtkWriter.h @@ -14,7 +14,7 @@ typedef enum VtkFormat { ASCII = 0, BINARY } VtkFormat; typedef struct VtkOptions { - Grid grid; + Grid* grid; #ifdef _VTK_WRITER_MPI MPI_File fh; #else diff --git a/BasicSolver/3D-seq/canal.par b/BasicSolver/3D-seq/canal.par index 30e1f41..118802b 100644 --- a/BasicSolver/3D-seq/canal.par +++ b/BasicSolver/3D-seq/canal.par @@ -42,6 +42,13 @@ te 100.0 # final time dt 0.02 # time stepsize tau 0.5 # safety factor for time stepsize control (<0 constant delt) +# Multigrid data: +# --------- + +levels 3 # Multigrid levels +presmooth 5 # Pre-smoothning iterations +postsmooth 5 # Post-smoothning iterations + # Pressure Iteration Data: # ----------------------- diff --git a/BasicSolver/3D-seq/config.mk b/BasicSolver/3D-seq/config.mk index ccdfa2c..2cf33f9 100644 --- a/BasicSolver/3D-seq/config.mk +++ b/BasicSolver/3D-seq/config.mk @@ -1,12 +1,12 @@ # Supported: GCC, CLANG, ICC -TAG ?= CLANG +TAG ?= ICC ENABLE_OPENMP ?= false -# Supported: sor, mg +# Supported: rb, mg SOLVER ?= mg # Run in debug settings DEBUG ?= false #Feature options OPTIONS += -DARRAY_ALIGNMENT=64 -#OPTIONS += -DVERBOSE +OPTIONS += -DVERBOSE #OPTIONS += -DDEBUG diff --git a/BasicSolver/3D-seq/dcavity.par b/BasicSolver/3D-seq/dcavity.par index 8224887..c3263dd 100644 --- a/BasicSolver/3D-seq/dcavity.par +++ b/BasicSolver/3D-seq/dcavity.par @@ -38,10 +38,17 @@ kmax 128 # number of interior cells in z-direction # Time Data: # --------- -te 2.0 # final time +te 2.0 # final time dt 0.02 # time stepsize tau 0.5 # safety factor for time stepsize control (<0 constant delt) +# Multigrid data: +# --------- + +levels 3 # Multigrid levels +presmooth 10 # Pre-smoothning iterations +postsmooth 5 # Post-smoothning iterations + # Pressure Iteration Data: # ----------------------- @@ -50,5 +57,4 @@ eps 0.001 # stopping tolerance for pressure iteration rho 0.5 omg 1.7 # relaxation parameter for SOR iteration gamma 0.9 # upwind differencing factor gamma -levels 5 # Multigrid levels #=============================================================================== diff --git a/BasicSolver/3D-seq/src/discretization.c b/BasicSolver/3D-seq/src/discretization.c index eacac17..8821e7e 100644 --- a/BasicSolver/3D-seq/src/discretization.c +++ b/BasicSolver/3D-seq/src/discretization.c @@ -106,7 +106,7 @@ void initDiscretization(Discretization* d, Parameter* p) d->dtBound = 0.5 * d->re * 1.0 / invSqrSum; #ifdef VERBOSE - printConfig(s); + printConfig(d); #endif /* VERBOSE */ } diff --git a/BasicSolver/3D-seq/src/main.c b/BasicSolver/3D-seq/src/main.c index c289631..129993a 100644 --- a/BasicSolver/3D-seq/src/main.c +++ b/BasicSolver/3D-seq/src/main.c @@ -105,7 +105,7 @@ int main(int argc, char** argv) nt++; #ifdef VERBOSE - printf("TIME %f , TIMESTEP %f\n", t, solver.dt); + printf("TIME %f , TIMESTEP %f\n", t, d.dt); #else printProgress(t); #endif diff --git a/BasicSolver/3D-seq/src/parameter.c b/BasicSolver/3D-seq/src/parameter.c index 005b69a..13a2d8a 100644 --- a/BasicSolver/3D-seq/src/parameter.c +++ b/BasicSolver/3D-seq/src/parameter.c @@ -14,19 +14,21 @@ void initParameter(Parameter* param) { - param->xlength = 1.0; - param->ylength = 1.0; - param->zlength = 1.0; - param->imax = 100; - param->jmax = 100; - param->kmax = 100; - param->itermax = 1000; - param->eps = 0.0001; - param->omg = 1.7; - param->re = 100.0; - param->gamma = 0.9; - param->tau = 0.5; - param->levels = 5; + param->xlength = 1.0; + param->ylength = 1.0; + param->zlength = 1.0; + param->imax = 100; + param->jmax = 100; + param->kmax = 100; + param->itermax = 1000; + param->eps = 0.0001; + param->omg = 1.7; + param->re = 100.0; + param->gamma = 0.9; + param->tau = 0.5; + param->levels = 5; + param->presmooth = 5; + param->postsmooth = 5; } void readParameter(Parameter* param, const char* filename) diff --git a/BasicSolver/3D-seq/src/parameter.h b/BasicSolver/3D-seq/src/parameter.h index 83917c0..eb2b753 100644 --- a/BasicSolver/3D-seq/src/parameter.h +++ b/BasicSolver/3D-seq/src/parameter.h @@ -18,6 +18,7 @@ typedef struct { char* name; int bcLeft, bcRight, bcBottom, bcTop, bcFront, bcBack; double u_init, v_init, w_init, p_init; + int presmooth, postsmooth; } Parameter; void initParameter(Parameter*); diff --git a/BasicSolver/3D-seq/src/solver-mg.c b/BasicSolver/3D-seq/src/solver-mg.c index 9644f59..668312e 100644 --- a/BasicSolver/3D-seq/src/solver-mg.c +++ b/BasicSolver/3D-seq/src/solver-mg.c @@ -141,13 +141,62 @@ static double smooth( for (int j = 1; j < jmax + 1; j++) { for (int i = isw; i < imax + 1; i += 2) { - R(i, j, k) = - RHS(i, j, k) - - ((P(i + 1, j, k) - 2.0 * P(i, j, k) + P(i - 1, j, k)) * idx2 + - (P(i, j + 1, k) - 2.0 * P(i, j, k) + P(i, j - 1, k)) * idy2 + - (P(i, j, k + 1) - 2.0 * P(i, j, k) + P(i, j, k - 1)) * idz2); + P(i, j, k) -= + factor * + (RHS(i, j, k) - + ((P(i + 1, j, k) - 2.0 * P(i, j, k) + P(i - 1, j, k)) * idx2 + + (P(i, j + 1, k) - 2.0 * P(i, j, k) + P(i, j - 1, k)) * + idy2 + + (P(i, j, k + 1) - 2.0 * P(i, j, k) + P(i, j, k - 1)) * + idz2)); + } + isw = 3 - isw; + } + jsw = 3 - jsw; + } + ksw = 3 - ksw; + } +} + +static double calculateResidual( + Solver* s, double* p, double* rhs, int level, int imax, int jmax, int kmax) +{ + double eps = s->eps; + int itermax = s->itermax; + double dx2 = s->grid->dx * s->grid->dx; + double dy2 = s->grid->dy * s->grid->dy; + double dz2 = s->grid->dz * s->grid->dz; + double idx2 = 1.0 / dx2; + double idy2 = 1.0 / dy2; + double idz2 = 1.0 / dz2; + double factor = s->omega * 0.5 * (dx2 * dy2 * dz2) / + (dy2 * dz2 + dx2 * dz2 + dx2 * dy2); + double* r = s->r[level]; + double epssq = eps * eps; + int it = 0; + int pass, ksw, jsw, isw; + double res = 1.0; + + ksw = 1; + + for (pass = 0; pass < 2; pass++) { + jsw = ksw; + + for (int k = 1; k < kmax + 1; k++) { + isw = jsw; + for (int j = 1; j < jmax + 1; j++) { + for (int i = isw; i < imax + 1; i += 2) { + + R(i, + j, + k) = (RHS(i, j, k) - + ((P(i + 1, j, k) - 2.0 * P(i, j, k) + P(i - 1, j, k)) * + idx2 + + (P(i, j + 1, k) - 2.0 * P(i, j, k) + P(i, j - 1, k)) * + idy2 + + (P(i, j, k + 1) - 2.0 * P(i, j, k) + P(i, j, k - 1)) * + idz2)); - P(i, j, k) -= (factor * R(i, j, k)); res += (R(i, j, k) * R(i, j, k)); } isw = 3 - isw; @@ -167,7 +216,7 @@ static double multiGrid( { double res = 0.0; - // coarsest level TODO: Use direct solver? + // coarsest level if (level == COARSEST_LEVEL) { for (int i = 0; i < 5; i++) { smooth(s, p, rhs, level, imax, jmax, kmax); @@ -175,17 +224,18 @@ static double multiGrid( return res; } - // pre-smoothing TODO: Make smoothing steps configurable? - for (int i = 0; i < 5; i++) { + // pre-smoothing + for (int i = 0; i < s->presmooth; i++) { smooth(s, p, rhs, level, imax, jmax, kmax); if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax, kmax); } + res = calculateResidual(s, p, rhs, level, imax, jmax, kmax); + // restrict restrictMG(s, level, imax, jmax, kmax); // MGSolver on residual and error. - // TODO: What if there is a rest? multiGrid(s, s->e[level + 1], s->r[level + 1], @@ -202,8 +252,8 @@ static double multiGrid( if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax, kmax); // post-smoothing - for (int i = 0; i < 5; i++) { - res = smooth(s, p, rhs, level, imax, jmax, kmax); + for (int i = 0; i < s->postsmooth; i++) { + smooth(s, p, rhs, level, imax, jmax, kmax); if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax, kmax); } @@ -212,11 +262,13 @@ static double multiGrid( void initSolver(Solver* s, Discretization* d, Parameter* p) { - s->eps = p->eps; - s->omega = p->omg; - s->itermax = p->itermax; - s->levels = p->levels; - s->grid = &d->grid; + s->eps = p->eps; + s->omega = p->omg; + s->itermax = p->itermax; + s->levels = p->levels; + s->grid = &d->grid; + s->presmooth = p->presmooth; + s->postsmooth = p->postsmooth; int imax = s->grid->imax; int jmax = s->grid->jmax; diff --git a/BasicSolver/3D-seq/src/solver-sor.c b/BasicSolver/3D-seq/src/solver-rb.c similarity index 100% rename from BasicSolver/3D-seq/src/solver-sor.c rename to BasicSolver/3D-seq/src/solver-rb.c diff --git a/BasicSolver/3D-seq/src/solver.h b/BasicSolver/3D-seq/src/solver.h index d7ce513..9a6a2f6 100644 --- a/BasicSolver/3D-seq/src/solver.h +++ b/BasicSolver/3D-seq/src/solver.h @@ -18,6 +18,7 @@ typedef struct { int itermax; int levels; double **r, **e; + int presmooth, postsmooth; } Solver; extern void initSolver(Solver*, Discretization*, Parameter*); diff --git a/BasicSolver/3D-seq/src/vtkWriter.c b/BasicSolver/3D-seq/src/vtkWriter.c index 8668b44..2b8f92f 100644 --- a/BasicSolver/3D-seq/src/vtkWriter.c +++ b/BasicSolver/3D-seq/src/vtkWriter.c @@ -70,6 +70,7 @@ void vtkScalar(VtkOptions* o, char* name, double* s) exit(EXIT_FAILURE); } fprintf(o->fh, "SCALARS %s float\n", name); + fprintf(o->fh, "LOOKUP_TABLE default\n"); for (int k = 0; k < kmax; k++) { for (int j = 0; j < jmax; j++) {