From ae6303b519762f8f8339d43196ea50f2aba7d73d Mon Sep 17 00:00:00 2001 From: Aditya ujeniya Date: Tue, 9 Jul 2024 17:34:23 +0200 Subject: [PATCH] Porting working 2D-seq and 2D-mpi with MG --- BasicSolver/2D-mpi/Makefile | 1 + BasicSolver/2D-mpi/canal.par | 19 +- BasicSolver/2D-mpi/config.mk | 5 +- BasicSolver/2D-mpi/dcavity.par | 11 +- BasicSolver/2D-mpi/src/comm-v1.c | 25 +- BasicSolver/2D-mpi/src/comm-v2.c | 54 +++- BasicSolver/2D-mpi/src/comm-v3.c | 53 ++- BasicSolver/2D-mpi/src/comm.h | 2 + BasicSolver/2D-mpi/src/discretization.c | 7 +- BasicSolver/2D-mpi/src/main.c | 4 +- BasicSolver/2D-mpi/src/parameter.c | 20 +- BasicSolver/2D-mpi/src/parameter.h | 1 + BasicSolver/2D-mpi/src/solver-mg.c | 302 ++++++++++++++++++ .../2D-mpi/src/{solver.c => solver-rb.c} | 0 BasicSolver/2D-mpi/src/solver.h | 3 + BasicSolver/2D-seq/canal.par | 1 + BasicSolver/2D-seq/config.mk | 6 +- BasicSolver/2D-seq/dcavity.par | 4 +- BasicSolver/2D-seq/src/main.c | 2 +- 19 files changed, 488 insertions(+), 32 deletions(-) create mode 100644 BasicSolver/2D-mpi/src/solver-mg.c rename BasicSolver/2D-mpi/src/{solver.c => solver-rb.c} (100%) diff --git a/BasicSolver/2D-mpi/Makefile b/BasicSolver/2D-mpi/Makefile index 2b23422..44863f9 100644 --- a/BasicSolver/2D-mpi/Makefile +++ b/BasicSolver/2D-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)/comm-$(COMM_TYPE).o +OBJ += $(BUILD_DIR)/solver-$(SOLVER).o SOURCES = $(SRC) $(wildcard $(SRC_DIR)/*.h) CPPFLAGS := $(CPPFLAGS) $(DEFINES) $(OPTIONS) $(INCLUDES) diff --git a/BasicSolver/2D-mpi/canal.par b/BasicSolver/2D-mpi/canal.par index 80dfa19..37b26fa 100644 --- a/BasicSolver/2D-mpi/canal.par +++ b/BasicSolver/2D-mpi/canal.par @@ -7,10 +7,10 @@ name canal # name of flow setup -bcN 1 # flags for boundary conditions -bcE 3 # 1 = no-slip 3 = outflow -bcS 1 # 2 = free-slip 4 = periodic -bcW 3 # +bcTop 1 # flags for boundary conditions +bcBottom 1 # 1 = no-slip 3 = outflow +bcLeft 3 # 2 = free-slip 4 = periodic +bcRight 3 # gx 0.0 # Body forces (e.g. gravity) gy 0.0 # @@ -27,15 +27,22 @@ p_init 0.0 # initial value for pressure xlength 30.0 # domain size in x-direction ylength 4.0 # domain size in y-direction imax 200 # number of interior cells in x-direction -jmax 50 # number of interior cells in y-direction +jmax 40 # number of interior cells in y-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 2 # Multigrid levels +presmooth 5 # Pre-smoothning iterations +postsmooth 5 # Post-smoothning iterations + # Pressure Iteration Data: # ----------------------- diff --git a/BasicSolver/2D-mpi/config.mk b/BasicSolver/2D-mpi/config.mk index a7b7925..d635a3c 100644 --- a/BasicSolver/2D-mpi/config.mk +++ b/BasicSolver/2D-mpi/config.mk @@ -1,7 +1,10 @@ # Supported: GCC, CLANG, ICC -TAG ?= CLANG +TAG ?= ICC ENABLE_MPI ?= true ENABLE_OPENMP ?= false +# Supported: rb, mg +SOLVER ?= mg +# Run in debug settings ?= mg COMM_TYPE ?= v3 #Feature options diff --git a/BasicSolver/2D-mpi/dcavity.par b/BasicSolver/2D-mpi/dcavity.par index 2c81044..37f195e 100644 --- a/BasicSolver/2D-mpi/dcavity.par +++ b/BasicSolver/2D-mpi/dcavity.par @@ -26,8 +26,8 @@ p_init 0.0 # initial value for pressure xlength 1.0 # domain size in x-direction ylength 1.0 # domain size in y-direction -imax 80 # number of interior cells in x-direction -jmax 80 # number of interior cells in y-direction +imax 128 # number of interior cells in x-direction +jmax 128 # number of interior cells in y-direction # Time Data: # --------- @@ -36,6 +36,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 2 # Multigrid levels +presmooth 20 # Pre-smoothning iterations +postsmooth 5 # Post-smoothning iterations + # Pressure Iteration Data: # ----------------------- diff --git a/BasicSolver/2D-mpi/src/comm-v1.c b/BasicSolver/2D-mpi/src/comm-v1.c index 8f0ea3e..69d6600 100644 --- a/BasicSolver/2D-mpi/src/comm-v1.c +++ b/BasicSolver/2D-mpi/src/comm-v1.c @@ -173,9 +173,32 @@ void commPartition(Comm* c, int jmax, int imax) { #ifdef _MPI c->imaxLocal = imax; - c->jmaxLocal = sizeOfRank(c->rank, c->size, jmax); + c->jmaxLocal = sizeOfRank(c->coords[JDIM], c->size, jmax); #else c->imaxLocal = imax; c->jmaxLocal = jmax; #endif } + +void commUpdateDatatypes(Comm* oldcomm, Comm* newcomm, int imaxLocal, int jmaxLocal) +{ + Comm* newcomm; + +#if defined _MPI + newcomm->comm = MPI_COMM_NULL; + int result = MPI_Comm_dup(oldcomm->comm, &newcomm->comm); + + if (result == MPI_ERR_COMM) { + printf("\nNull communicator. Duplication failed !!\n"); + } +#endif + newcomm->imaxLocal = imaxLocal; + newcomm->jmaxLocal = jmaxLocal; +} + +void commFreeCommunicator(Comm* comm) +{ +#ifdef _MPI + MPI_Comm_free(&comm->comm); +#endif +} \ No newline at end of file diff --git a/BasicSolver/2D-mpi/src/comm-v2.c b/BasicSolver/2D-mpi/src/comm-v2.c index c50a343..c0285ca 100644 --- a/BasicSolver/2D-mpi/src/comm-v2.c +++ b/BasicSolver/2D-mpi/src/comm-v2.c @@ -252,8 +252,8 @@ void commPartition(Comm* c, int jmax, int imax) MPI_Cart_shift(c->comm, JDIM, 1, &c->neighbours[BOTTOM], &c->neighbours[TOP]); MPI_Cart_get(c->comm, NDIMS, c->dims, periods, c->coords); - int imaxLocal = sizeOfRank(c->rank, dims[IDIM], imax); - int jmaxLocal = sizeOfRank(c->rank, dims[JDIM], jmax); + int imaxLocal = sizeOfRank(c->coords[IDIM], dims[IDIM], imax); + int jmaxLocal = sizeOfRank(c->coords[JDIM], dims[JDIM], jmax); c->imaxLocal = imaxLocal; c->jmaxLocal = jmaxLocal; @@ -285,3 +285,53 @@ void commPartition(Comm* c, int jmax, int imax) c->jmaxLocal = jmax; #endif } + +void commUpdateDatatypes(Comm* oldcomm, Comm* newcomm, int imaxLocal, int jmaxLocal) +{ + Comm* newcomm; + +#if defined _MPI + newcomm->comm = MPI_COMM_NULL; + 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; + + MPI_Datatype jBufferType; + MPI_Type_contiguous(imaxLocal, MPI_DOUBLE, &jBufferType); + MPI_Type_commit(&jBufferType); + + MPI_Datatype iBufferType; + MPI_Type_vector(jmaxLocal, 1, imaxLocal + 2, MPI_DOUBLE, &iBufferType); + MPI_Type_commit(&iBufferType); + + newcomm->bufferTypes[LEFT] = iBufferType; + newcomm->bufferTypes[RIGHT] = iBufferType; + newcomm->bufferTypes[BOTTOM] = jBufferType; + newcomm->bufferTypes[TOP] = jBufferType; + + newcomm->sdispls[LEFT] = (imaxLocal + 2) + 1; + newcomm->sdispls[RIGHT] = (imaxLocal + 2) + imaxLocal; + newcomm->sdispls[BOTTOM] = (imaxLocal + 2) + 1; + newcomm->sdispls[TOP] = jmaxLocal * (imaxLocal + 2) + 1; + + newcomm->rdispls[LEFT] = (imaxLocal + 2); + newcomm->rdispls[RIGHT] = (imaxLocal + 2) + (imaxLocal + 1); + newcomm->rdispls[BOTTOM] = 1; + newcomm->rdispls[TOP] = (jmaxLocal + 1) * (imaxLocal + 2) + 1; +#else + newcomm->imaxLocal = imaxLocal; + newcomm->jmaxLocal = jmaxLocal; +#endif +} + +void commFreeCommunicator(Comm* comm) +{ + #ifdef _MPI + MPI_Comm_free(&comm->comm); + #endif +} \ No newline at end of file diff --git a/BasicSolver/2D-mpi/src/comm-v3.c b/BasicSolver/2D-mpi/src/comm-v3.c index 6f0d3b7..1fd8315 100644 --- a/BasicSolver/2D-mpi/src/comm-v3.c +++ b/BasicSolver/2D-mpi/src/comm-v3.c @@ -139,7 +139,6 @@ void commExchange(Comm* c, double* grid) { #ifdef _MPI int counts[NDIRS] = { 1, 1, 1, 1 }; - MPI_Neighbor_alltoallw(grid, counts, c->sdispls, @@ -233,8 +232,8 @@ void commPartition(Comm* c, int jmax, int imax) MPI_Cart_shift(c->comm, JDIM, 1, &c->neighbours[BOTTOM], &c->neighbours[TOP]); MPI_Cart_get(c->comm, NDIMS, c->dims, periods, c->coords); - int imaxLocal = sizeOfRank(c->rank, dims[IDIM], imax); - int jmaxLocal = sizeOfRank(c->rank, dims[JDIM], jmax); + int imaxLocal = sizeOfRank(c->coords[IDIM], dims[IDIM], imax); + int jmaxLocal = sizeOfRank(c->coords[JDIM], dims[JDIM], jmax); c->imaxLocal = imaxLocal; c->jmaxLocal = jmaxLocal; @@ -267,3 +266,51 @@ void commPartition(Comm* c, int jmax, int imax) c->jmaxLocal = jmax; #endif } + +void commUpdateDatatypes(Comm* oldcomm, Comm* newcomm, int imaxLocal, int jmaxLocal) +{ +#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; + + MPI_Datatype jBufferType; + MPI_Type_contiguous(imaxLocal, MPI_DOUBLE, &jBufferType); + MPI_Type_commit(&jBufferType); + + MPI_Datatype iBufferType; + MPI_Type_vector(jmaxLocal, 1, imaxLocal + 2, MPI_DOUBLE, &iBufferType); + MPI_Type_commit(&iBufferType); + + newcomm->bufferTypes[LEFT] = iBufferType; + newcomm->bufferTypes[RIGHT] = iBufferType; + newcomm->bufferTypes[BOTTOM] = jBufferType; + newcomm->bufferTypes[TOP] = jBufferType; + + newcomm->sdispls[LEFT] = (imaxLocal + 2) + 1; + newcomm->sdispls[RIGHT] = (imaxLocal + 2) + imaxLocal; + newcomm->sdispls[BOTTOM] = (imaxLocal + 2) + 1; + newcomm->sdispls[TOP] = jmaxLocal * (imaxLocal + 2) + 1; + + newcomm->rdispls[LEFT] = (imaxLocal + 2); + newcomm->rdispls[RIGHT] = (imaxLocal + 2) + (imaxLocal + 1); + newcomm->rdispls[BOTTOM] = 1; + newcomm->rdispls[TOP] = (jmaxLocal + 1) * (imaxLocal + 2) + 1; +#else + newcomm->imaxLocal = imaxLocal; + newcomm->jmaxLocal = jmaxLocal; +#endif +} + +void commFreeCommunicator(Comm* comm) +{ + #ifdef _MPI + MPI_Comm_free(&comm->comm); + #endif +} \ No newline at end of file diff --git a/BasicSolver/2D-mpi/src/comm.h b/BasicSolver/2D-mpi/src/comm.h index 85b4f70..30eaf7b 100644 --- a/BasicSolver/2D-mpi/src/comm.h +++ b/BasicSolver/2D-mpi/src/comm.h @@ -41,6 +41,8 @@ extern void commExchange(Comm*, double*); extern void commShift(Comm* c, double* f, double* g); extern void commReduction(double* v, int op); extern int commIsBoundary(Comm* c, int direction); +extern void commUpdateDatatypes(Comm*, Comm*, int, int); +extern void commFreeCommunicator(Comm*); extern void commCollectResult(Comm* c, double* ug, double* vg, diff --git a/BasicSolver/2D-mpi/src/discretization.c b/BasicSolver/2D-mpi/src/discretization.c index e9fc421..69d3582 100644 --- a/BasicSolver/2D-mpi/src/discretization.c +++ b/BasicSolver/2D-mpi/src/discretization.c @@ -273,17 +273,18 @@ void setSpecialBoundaryCondition(Discretization* s) 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) + + int rest = s->grid.jmax % s->comm.dims[JDIM]; + int yc = s->comm.rank * (s->grid.jmax / s->comm.dims[JDIM]) + 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); */ + // printf("RANK %d yc: %d ys: %f\n", s->comm.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); + } } } diff --git a/BasicSolver/2D-mpi/src/main.c b/BasicSolver/2D-mpi/src/main.c index c7f2d2f..b08f0ac 100644 --- a/BasicSolver/2D-mpi/src/main.c +++ b/BasicSolver/2D-mpi/src/main.c @@ -26,7 +26,9 @@ static void writeResults(Discretization* s) double* pg = allocate(64, bytesize); commCollectResult(&s->comm, ug, vg, pg, s->u, s->v, s->p, s->grid.imax, s->grid.jmax); - writeResult(s, ug, vg, pg); + if (commIsMaster(&s->comm)) { + writeResult(s, ug, vg, pg); + } free(ug); free(vg); diff --git a/BasicSolver/2D-mpi/src/parameter.c b/BasicSolver/2D-mpi/src/parameter.c index 7bac1e3..e5e7029 100644 --- a/BasicSolver/2D-mpi/src/parameter.c +++ b/BasicSolver/2D-mpi/src/parameter.c @@ -14,13 +14,16 @@ void initParameter(Parameter* param) { - param->xlength = 1.0; - param->ylength = 1.0; - param->imax = 100; - param->jmax = 100; - param->itermax = 1000; - param->eps = 0.0001; - param->omg = 1.8; + param->xlength = 1.0; + param->ylength = 1.0; + param->imax = 100; + param->jmax = 100; + param->itermax = 1000; + param->eps = 0.0001; + param->omg = 1.8; + param->levels = 5; + param->presmooth = 5; + param->postsmooth = 5; } void readParameter(Parameter* param, const char* filename) @@ -72,6 +75,9 @@ void readParameter(Parameter* param, const char* filename) PARSE_INT(bcRight); PARSE_INT(bcBottom); PARSE_INT(bcTop); + PARSE_INT(levels); + PARSE_INT(presmooth); + PARSE_INT(postsmooth); PARSE_REAL(u_init); PARSE_REAL(v_init); PARSE_REAL(p_init); diff --git a/BasicSolver/2D-mpi/src/parameter.h b/BasicSolver/2D-mpi/src/parameter.h index be28e82..ed765eb 100644 --- a/BasicSolver/2D-mpi/src/parameter.h +++ b/BasicSolver/2D-mpi/src/parameter.h @@ -18,6 +18,7 @@ typedef struct { char* name; int bcLeft, bcRight, bcBottom, bcTop; double u_init, v_init, p_init; + int levels, presmooth, postsmooth; } Parameter; void initParameter(Parameter*); diff --git a/BasicSolver/2D-mpi/src/solver-mg.c b/BasicSolver/2D-mpi/src/solver-mg.c new file mode 100644 index 0000000..b165214 --- /dev/null +++ b/BasicSolver/2D-mpi/src/solver-mg.c @@ -0,0 +1,302 @@ +/* + * 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) s[(j) * (imaxLocal + 2) + (i)] +#define E(i, j) e[(j) * (imaxLocal + 2) + (i)] +#define R(i, j) r[(j) * (imaxLocal + 2) + (i)] +#define OLD(i, j) old[(j) * (imaxLocal + 2) + (i)] + +static void restrictMG(Solver* s, int level, Comm* comm) +{ + int imaxLocal = comm->imaxLocal; + int jmaxLocal = comm->jmaxLocal; + + double* r = s->r[level + 1]; + double* old = s->r[level]; + +#ifdef _MPI + commExchange(comm, old); +#endif + + for (int j = 1; j < (jmaxLocal / 2) + 1; j++) { + for (int i = 1; i < (imaxLocal / 2) + 1; i++) { + R(i, j) = (OLD(2 * i - 1, 2 * j - 1) + OLD(2 * i, 2 * j - 1) * 2 + + OLD(2 * i + 1, 2 * j - 1) + OLD(2 * i - 1, 2 * j) * 2 + + OLD(2 * i, 2 * j) * 4 + OLD(2 * i + 1, 2 * j) * 2 + + OLD(2 * i - 1, 2 * j + 1) + OLD(2 * i, 2 * j + 1) * 2 + + OLD(2 * i + 1, 2 * j + 1)) / + 16.0; + } + } +} + +static void prolongate(Solver* s, int level, Comm* comm) +{ + int imaxLocal = comm->imaxLocal; + int jmaxLocal = comm->jmaxLocal; + + double* old = s->r[level + 1]; + double* e = s->r[level]; + + for (int j = 2; j < jmaxLocal + 1; j += 2) { + for (int i = 2; i < imaxLocal + 1; i += 2) { + E(i, j) = OLD(i / 2, j / 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; + + for (int j = 1; j < jmaxLocal + 1; ++j) { + for (int i = 1; i < imaxLocal + 1; ++i) { + P(i, j) += E(i, j); + } + } +} + +static void setBoundaryCondition(Solver* s, double* p, int imaxLocal, int jmaxLocal) +{ +#ifdef _MPI + 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 + for (int i = 1; i < imaxLocal + 1; i++) { + P(i, jmaxLocal + 1) = P(i, jmaxLocal); + } + } + + 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 + for (int j = 1; j < jmaxLocal + 1; j++) { + P(imaxLocal + 1, j) = P(imaxLocal, j); + } + } +#else + for (int i = 1; i < imaxLocal + 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(imaxLocal + 1, j) = P(imaxLocal, j); + } +#endif +} + + +static double smooth(Solver* s, double* p, double* rhs, int level, Comm* comm) +{ + int imaxLocal = comm->imaxLocal; + int jmaxLocal = comm->jmaxLocal; + + int imax = s->grid->imax; + int jmax = s->grid->jmax; + + 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* r = s->r[level]; + + double res = 1.0; + int pass, jsw, isw; + + jsw = 1; + + for (pass = 0; pass < 2; pass++) { + isw = jsw; + +#ifdef _MPI + commExchange(comm, p); +#endif + + for (int j = 1; j < jmaxLocal + 1; j++) { + for (int i = isw; i < imaxLocal + 1; i += 2) { + + P(i, j) -= factor * + (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)); + } + isw = 3 - isw; + } + jsw = 3 - jsw; + } +} + +static double calculateResidual(Solver* s, double* p, double* rhs, int level, Comm* comm) +{ + int imax = s->grid->imax; + int jmax = s->grid->jmax; + int imaxLocal = comm->imaxLocal; + int jmaxLocal = comm->jmaxLocal; + + 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* r = s->r[level]; + double res = 1.0; + int pass, jsw, isw; + + jsw = 1; + + for (pass = 0; pass < 2; pass++) { + isw = jsw; + +#ifdef _MPI + commExchange(comm, p); +#endif + + for (int j = 1; j < jmaxLocal + 1; j++) { + for (int i = isw; i < imaxLocal + 1; i += 2) { + + R(i, j) = 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); + + res += (R(i, j) * R(i, j)); + } + isw = 3 - isw; + } + jsw = 3 - jsw; + } + +#ifdef _MPI + commReduction(&res, SUM); +#endif + + res = res / (double)(imax * jmax); +#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) +{ + 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, comm->imaxLocal, comm->jmaxLocal); + } + + // calculate residuals + res = calculateResidual(s, p, rhs, level, comm); + + // restrict + restrictMG(s, level, comm); + + Comm newcomm; + commUpdateDatatypes(s->comm, &newcomm, comm->imaxLocal / 2, comm->jmaxLocal / 2); + + // MGSolver on residual and error. + // TODO: What if there is a rest? + 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, comm->imaxLocal, comm->jmaxLocal); + + // post-smoothing + for (int i = 0; i < s->postsmooth; i++) { + smooth(s, p, rhs, level, comm); + if (level == FINEST_LEVEL) + setBoundaryCondition(s, p, comm->imaxLocal, comm->jmaxLocal); + } + + 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->comm = &d->comm; + s->presmooth = p->presmooth; + s->postsmooth = p->postsmooth; + + int imax = s->grid->imax; + int jmax = s->grid->jmax; + 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) * sizeof(double); + + for (int j = 0; j < levels; j++) { + s->r[j] = allocate(64, size); + s->e[j] = allocate(64, size); + + for (int i = 0; i < (imax + 2) * (jmax + 2); 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/2D-mpi/src/solver.c b/BasicSolver/2D-mpi/src/solver-rb.c similarity index 100% rename from BasicSolver/2D-mpi/src/solver.c rename to BasicSolver/2D-mpi/src/solver-rb.c diff --git a/BasicSolver/2D-mpi/src/solver.h b/BasicSolver/2D-mpi/src/solver.h index 4ceb40b..890e664 100644 --- a/BasicSolver/2D-mpi/src/solver.h +++ b/BasicSolver/2D-mpi/src/solver.h @@ -9,6 +9,7 @@ #include "comm.h" #include "discretization.h" #include "grid.h" +#include "mpi.h" #include "parameter.h" typedef struct { @@ -17,6 +18,8 @@ typedef struct { /* parameters */ double eps, omega; int itermax; + int levels, presmooth, postsmooth; + double **r, **e; /* communication */ Comm* comm; } Solver; diff --git a/BasicSolver/2D-seq/canal.par b/BasicSolver/2D-seq/canal.par index 8b9fd9c..f113e4a 100644 --- a/BasicSolver/2D-seq/canal.par +++ b/BasicSolver/2D-seq/canal.par @@ -43,4 +43,5 @@ itermax 500 # maximal number of pressure iteration in one time step eps 0.00001 # stopping tolerance for pressure iteration omg 1.8 # relaxation parameter for SOR iteration gamma 0.9 # upwind differencing factor gamma +levels 5 # Multigrid levels #=============================================================================== diff --git a/BasicSolver/2D-seq/config.mk b/BasicSolver/2D-seq/config.mk index dc874ee..0b2c509 100644 --- a/BasicSolver/2D-seq/config.mk +++ b/BasicSolver/2D-seq/config.mk @@ -1,12 +1,12 @@ # Supported: GCC, CLANG, ICC -TAG ?= CLANG +TAG ?= ICC ENABLE_OPENMP ?= false # Supported: sor, rb, mg -SOLVER ?= rb +SOLVER ?= mg # Run in debug settings DEBUG ?= false #Feature options OPTIONS += -DARRAY_ALIGNMENT=64 -#OPTIONS += -DVERBOSE +OPTIONS += -DVERBOSE #OPTIONS += -DDEBUG diff --git a/BasicSolver/2D-seq/dcavity.par b/BasicSolver/2D-seq/dcavity.par index 798819f..d8e311e 100644 --- a/BasicSolver/2D-seq/dcavity.par +++ b/BasicSolver/2D-seq/dcavity.par @@ -26,8 +26,8 @@ p_init 0.0 # initial value for pressure xlength 1.0 # domain size in x-direction ylength 1.0 # domain size in y-direction -imax 128 # number of interior cells in x-direction -jmax 128 # number of interior cells in y-direction +imax 800 # number of interior cells in x-direction +jmax 800 # number of interior cells in y-direction # Time Data: # --------- diff --git a/BasicSolver/2D-seq/src/main.c b/BasicSolver/2D-seq/src/main.c index 822e32a..51b9877 100644 --- a/BasicSolver/2D-seq/src/main.c +++ b/BasicSolver/2D-seq/src/main.c @@ -56,7 +56,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