From d9361a8086a6606960235825659cf48716ca8b4e Mon Sep 17 00:00:00 2001 From: Aditya ujeniya Date: Wed, 24 Jul 2024 22:56:05 +0200 Subject: [PATCH] Fully working BasicSolver --- BasicSolver/2D-mpi/Makefile | 15 ++++++++++-- BasicSolver/2D-mpi/config.mk | 3 ++- BasicSolver/2D-mpi/dcavity.par | 8 +++---- BasicSolver/2D-mpi/residual.plot | 9 ++++++++ BasicSolver/2D-mpi/src/comm-v1.c | 36 ++++++++++++++++++++++++++--- BasicSolver/2D-mpi/src/comm-v2.c | 33 +++++++++++++++----------- BasicSolver/2D-mpi/src/comm-v3.c | 24 +++++++++++-------- BasicSolver/2D-mpi/src/main.c | 12 ++++++++-- BasicSolver/2D-mpi/src/progress.c | 21 ++++++++++++++++- BasicSolver/2D-mpi/src/progress.h | 3 ++- BasicSolver/2D-mpi/src/solver-mg.c | 18 +++++++-------- BasicSolver/2D-mpi/src/solver-rb.c | 4 +++- BasicSolver/2D-mpi/src/solver.h | 2 +- BasicSolver/2D-seq/Makefile | 17 ++++++++++++-- BasicSolver/2D-seq/residual.plot | 9 ++++++++ BasicSolver/2D-seq/src/main.c | 15 ++++++++---- BasicSolver/2D-seq/src/progress.c | 19 +++++++++++++++ BasicSolver/2D-seq/src/progress.h | 2 ++ BasicSolver/2D-seq/src/solver-mg.c | 4 +++- BasicSolver/2D-seq/src/solver-rb.c | 4 +++- BasicSolver/2D-seq/src/solver-sor.c | 4 +++- BasicSolver/2D-seq/src/solver.h | 2 +- BasicSolver/3D-mpi/Makefile | 14 +++++++++-- BasicSolver/3D-mpi/canal.par | 8 +++---- BasicSolver/3D-mpi/config.mk | 1 + BasicSolver/3D-mpi/dcavity.par | 4 ++-- BasicSolver/3D-mpi/residual.plot | 9 ++++++++ BasicSolver/3D-mpi/src/comm.c | 28 ++++++++++++++-------- BasicSolver/3D-mpi/src/main.c | 12 ++++++++-- BasicSolver/3D-mpi/src/parameter.c | 3 +++ BasicSolver/3D-mpi/src/progress.c | 19 +++++++++++++-- BasicSolver/3D-mpi/src/progress.h | 7 ++++-- BasicSolver/3D-mpi/src/solver-mg.c | 12 ++++++---- BasicSolver/3D-mpi/src/solver-rb.c | 28 +++++++++++----------- BasicSolver/3D-mpi/src/solver.h | 2 +- BasicSolver/3D-seq/Makefile | 12 ++++++++++ BasicSolver/3D-seq/canal.par | 2 +- BasicSolver/3D-seq/dcavity.par | 4 ++-- BasicSolver/3D-seq/residual.plot | 9 ++++++++ BasicSolver/3D-seq/src/main.c | 12 ++++++++-- BasicSolver/3D-seq/src/progress.c | 19 +++++++++++++++ BasicSolver/3D-seq/src/progress.h | 3 ++- BasicSolver/3D-seq/src/solver-mg.c | 4 +++- BasicSolver/3D-seq/src/solver-rb.c | 4 +++- BasicSolver/3D-seq/src/solver.h | 2 +- 45 files changed, 372 insertions(+), 110 deletions(-) create mode 100644 BasicSolver/2D-mpi/residual.plot create mode 100644 BasicSolver/2D-seq/residual.plot create mode 100644 BasicSolver/3D-mpi/residual.plot create mode 100644 BasicSolver/3D-seq/residual.plot diff --git a/BasicSolver/2D-mpi/Makefile b/BasicSolver/2D-mpi/Makefile index 44863f9..eb5349c 100644 --- a/BasicSolver/2D-mpi/Makefile +++ b/BasicSolver/2D-mpi/Makefile @@ -39,9 +39,20 @@ $(BUILD_DIR)/%.s: %.c $(info ===> GENERATE ASM $@) $(CC) -S $(CPPFLAGS) $(CFLAGS) $< -o $@ -.PHONY: clean distclean tags info asm format +.PHONY: clean distclean vis vis_clean tags info asm format -clean: +vis: + $(info ===> GENERATE VISUALIZATION) + @gnuplot -e "filename='pressure.dat'" ./surface.plot + @gnuplot -e "filename='velocity.dat'" ./vector.plot + @gnuplot -e "filename='residual.dat'" ./residual.plot + +vis_clean: + $(info ===> CLEAN VISUALIZATION) + @rm -f *.dat + @rm -f *.png + +clean: vis_clean $(info ===> CLEAN) @rm -rf $(BUILD_DIR) @rm -f tags diff --git a/BasicSolver/2D-mpi/config.mk b/BasicSolver/2D-mpi/config.mk index d635a3c..eb5546a 100644 --- a/BasicSolver/2D-mpi/config.mk +++ b/BasicSolver/2D-mpi/config.mk @@ -1,10 +1,11 @@ # Supported: GCC, CLANG, ICC TAG ?= ICC +# Supported: true, false ENABLE_MPI ?= true ENABLE_OPENMP ?= false # Supported: rb, mg SOLVER ?= mg -# Run in debug settings ?= mg +# Supported: v1, v2, v3 COMM_TYPE ?= v3 #Feature options diff --git a/BasicSolver/2D-mpi/dcavity.par b/BasicSolver/2D-mpi/dcavity.par index 37f195e..8f229ed 100644 --- a/BasicSolver/2D-mpi/dcavity.par +++ b/BasicSolver/2D-mpi/dcavity.par @@ -33,14 +33,14 @@ jmax 128 # number of interior cells in y-direction # --------- te 10.0 # final time -dt 0.02 # time stepsize -tau 0.5 # safety factor for time stepsize control (<0 constant delt) +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 +presmooth 20 # Pre-smoothning iterations postsmooth 5 # Post-smoothning iterations # Pressure Iteration Data: @@ -48,6 +48,6 @@ postsmooth 5 # Post-smoothning iterations itermax 1000 # maximal number of pressure iteration in one time step eps 0.001 # stopping tolerance for pressure iteration -omg 1.9 # relaxation parameter for SOR iteration +omg 1.7 # relaxation parameter for SOR iteration gamma 0.9 # upwind differencing factor gamma #=============================================================================== diff --git a/BasicSolver/2D-mpi/residual.plot b/BasicSolver/2D-mpi/residual.plot new file mode 100644 index 0000000..36fb011 --- /dev/null +++ b/BasicSolver/2D-mpi/residual.plot @@ -0,0 +1,9 @@ +set terminal png size 1800,768 enhanced font ,12 +set output 'residual.png' +set datafile separator whitespace +set xlabel "Timestep" +set ylabel "Residual" + +set logscale y 2 + +plot 'residual.dat' using 1:2 title "Residual" \ No newline at end of file diff --git a/BasicSolver/2D-mpi/src/comm-v1.c b/BasicSolver/2D-mpi/src/comm-v1.c index 69d6600..16a0422 100644 --- a/BasicSolver/2D-mpi/src/comm-v1.c +++ b/BasicSolver/2D-mpi/src/comm-v1.c @@ -14,7 +14,7 @@ static int sum(int* sizes, int position) { int sum = 0; - for (int i = 0; i < position; i++) { + for (int i = 0; i < position; i += position) { sum += sizes[i]; } @@ -67,6 +67,7 @@ int commIsBoundary(Comm* c, int direction) void commExchange(Comm* c, double* grid) { + // printf("Rank : %d In exchange \n", c->rank); #ifdef _MPI MPI_Request requests[4] = { MPI_REQUEST_NULL, MPI_REQUEST_NULL, @@ -174,6 +175,17 @@ void commPartition(Comm* c, int jmax, int imax) #ifdef _MPI c->imaxLocal = imax; c->jmaxLocal = sizeOfRank(c->coords[JDIM], c->size, jmax); + + c->neighbours[BOTTOM] = c->rank == 0 ? -1 : c->rank - 1; + c->neighbours[TOP] = c->rank == (c->size - 1) ? -1 : c->rank + 1; + c->neighbours[LEFT] = -1; + c->neighbours[RIGHT] = -1; + + c->coords[IDIM] = 0; + c->coords[JDIM] = c->rank; + + c->dims[IDIM] = 1; + c->dims[JDIM] = c->size; #else c->imaxLocal = imax; c->jmaxLocal = jmax; @@ -182,15 +194,33 @@ void commPartition(Comm* c, int jmax, int imax) 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); + int result = MPI_Comm_dup(MPI_COMM_WORLD, &newcomm->comm); if (result == MPI_ERR_COMM) { printf("\nNull communicator. Duplication failed !!\n"); } + + newcomm->rank = oldcomm->rank; + newcomm->size = oldcomm->size; + + newcomm->imaxLocal = imaxLocal / 2; + newcomm->jmaxLocal = jmaxLocal / 2; + + newcomm->neighbours[BOTTOM] = newcomm->rank == 0 ? -1 : newcomm->rank - 1; + newcomm->neighbours[TOP] = newcomm->rank == (newcomm->size - 1) ? -1 : newcomm->rank + 1; + newcomm->neighbours[LEFT] = -1; + newcomm->neighbours[RIGHT] = -1; + + newcomm->coords[IDIM] = 0; + newcomm->coords[JDIM] = newcomm->rank; + + newcomm->dims[IDIM] = 1; + newcomm->dims[JDIM] = newcomm->size; + + #endif newcomm->imaxLocal = imaxLocal; newcomm->jmaxLocal = jmaxLocal; diff --git a/BasicSolver/2D-mpi/src/comm-v2.c b/BasicSolver/2D-mpi/src/comm-v2.c index 90d8e43..5ace203 100644 --- a/BasicSolver/2D-mpi/src/comm-v2.c +++ b/BasicSolver/2D-mpi/src/comm-v2.c @@ -4,18 +4,18 @@ * Use of this source code is governed by a MIT style * license that can be found in the LICENSE file. */ +#include "comm.h" #include #include - -#include "comm.h" +#include #ifdef _MPI // subroutines local to this module -static int sum(int* sizes, int position) +static int sum(int* sizes, int init, int offset, int coord) { int sum = 0; - for (int i = 0; i < position; i++) { + for (int i = init - offset; coord > 0; i -= offset, --coord) { sum += sizes[i]; } @@ -79,8 +79,8 @@ static void assembleResult(Comm* c, double* src, double* dst, int imax, int jmax int newSizes[NDIMS] = { newSizesJ[i], newSizesI[i] }; int coords[NDIMS]; MPI_Cart_coords(c->comm, i, NDIMS, coords); - int starts[NDIMS] = { sum(newSizesJ, coords[JDIM]), - sum(newSizesI, coords[IDIM]) }; + int starts[NDIMS] = { sum(newSizesJ, i, 1, coords[JDIM]), + sum(newSizesI, i, c->dims[JDIM], coords[IDIM]) }; printf( "Rank: %d, Coords(i,j): %d %d, Size(i,j): %d %d, Target Size(i,j): %d %d " "Starts(i,j): %d %d\n", @@ -290,14 +290,21 @@ void commUpdateDatatypes(Comm* oldcomm, Comm* newcomm, int imaxLocal, int jmaxLo { #if defined _MPI newcomm->comm = MPI_COMM_NULL; - int result = MPI_Comm_dup(oldcomm->comm, &newcomm->comm); + 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->rank = oldcomm->rank; + newcomm->size = oldcomm->size; + + memcpy(&newcomm->neighbours, &oldcomm->neighbours, sizeof(oldcomm->neighbours)); + memcpy(&newcomm->coords, &oldcomm->coords, sizeof(oldcomm->coords)); + memcpy(&newcomm->dims, &oldcomm->dims, sizeof(oldcomm->dims)); + + newcomm->imaxLocal = imaxLocal/2; + newcomm->jmaxLocal = jmaxLocal/2; MPI_Datatype jBufferType; MPI_Type_contiguous(imaxLocal, MPI_DOUBLE, &jBufferType); @@ -329,7 +336,7 @@ void commUpdateDatatypes(Comm* oldcomm, Comm* newcomm, int imaxLocal, int jmaxLo void commFreeCommunicator(Comm* comm) { - #ifdef _MPI - MPI_Comm_free(&comm->comm); - #endif +#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 1fd8315..a1e6fb2 100644 --- a/BasicSolver/2D-mpi/src/comm-v3.c +++ b/BasicSolver/2D-mpi/src/comm-v3.c @@ -11,11 +11,11 @@ #ifdef _MPI // subroutines local to this module -static int sum(int* sizes, int position) +static int sum(int* sizes, int init, int offset, int coord) { int sum = 0; - for (int i = 0; i < position; i++) { + for (int i = init - offset; coord > 0; i -= offset, --coord) { sum += sizes[i]; } @@ -79,8 +79,8 @@ static void assembleResult(Comm* c, double* src, double* dst, int imax, int jmax int newSizes[NDIMS] = { newSizesJ[i], newSizesI[i] }; int coords[NDIMS]; MPI_Cart_coords(c->comm, i, NDIMS, coords); - int starts[NDIMS] = { sum(newSizesJ, coords[JDIM]), - sum(newSizesI, coords[IDIM]) }; + int starts[NDIMS] = { sum(newSizesJ, i, 1, coords[JDIM]), + sum(newSizesI, i, c->dims[JDIM], coords[IDIM]) }; printf( "Rank: %d, Coords(i,j): %d %d, Size(i,j): %d %d, Target Size(i,j): %d %d " "Starts(i,j): %d %d\n", @@ -271,14 +271,18 @@ void commUpdateDatatypes(Comm* oldcomm, Comm* newcomm, int imaxLocal, int jmaxLo { #if defined _MPI - int result = MPI_Comm_dup(oldcomm->comm, &newcomm->comm); + 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->rank = oldcomm->rank; + newcomm->size = oldcomm->size; + + + newcomm->imaxLocal = imaxLocal / 2; + newcomm->jmaxLocal = jmaxLocal / 2; MPI_Datatype jBufferType; MPI_Type_contiguous(imaxLocal, MPI_DOUBLE, &jBufferType); @@ -310,7 +314,7 @@ void commUpdateDatatypes(Comm* oldcomm, Comm* newcomm, int imaxLocal, int jmaxLo void commFreeCommunicator(Comm* comm) { - #ifdef _MPI - MPI_Comm_free(&comm->comm); - #endif +#ifdef _MPI + MPI_Comm_free(&comm->comm); +#endif } \ No newline at end of file diff --git a/BasicSolver/2D-mpi/src/main.c b/BasicSolver/2D-mpi/src/main.c index b08f0ac..8f93767 100644 --- a/BasicSolver/2D-mpi/src/main.c +++ b/BasicSolver/2D-mpi/src/main.c @@ -49,6 +49,9 @@ int main(int argc, char** argv) commInit(&d.comm, argc, argv); initParameter(&p); + FILE* fp; + if (commIsMaster(&d.comm)) fp = initResidualWriter(); + if (argc != 2) { printf("Usage: %s \n", argv[0]); exit(EXIT_SUCCESS); @@ -79,16 +82,21 @@ int main(int argc, char** argv) double tau = d.tau; double te = d.te; double t = 0.0; + double res = 0.0; timeStart = getTimeStamp(); while (t <= te) { + if (tau > 0.0) computeTimestep(&d); setBoundaryConditions(&d); setSpecialBoundaryCondition(&d); computeFG(&d); computeRHS(&d); - solve(&s, d.p, d.rhs); + res = solve(&s, d.p, d.rhs); adaptUV(&d); + + if (commIsMaster(&d.comm)) writeResidual(fp, t, res); + t += d.dt; #ifdef VERBOSE @@ -106,7 +114,7 @@ int main(int argc, char** argv) if (commIsMaster(s.comm)) { printf("Solution took %.2fs\n", timeStop - timeStart); } - + if (commIsMaster(&d.comm)) fclose(fp); writeResults(&d); commFinalize(s.comm); return EXIT_SUCCESS; diff --git a/BasicSolver/2D-mpi/src/progress.c b/BasicSolver/2D-mpi/src/progress.c index be98fce..4f3b22e 100644 --- a/BasicSolver/2D-mpi/src/progress.c +++ b/BasicSolver/2D-mpi/src/progress.c @@ -8,7 +8,7 @@ #include #include #include - +#include #include "progress.h" static double _end; @@ -49,3 +49,22 @@ void stopProgress() printf("\n"); fflush(stdout); } + +FILE* initResidualWriter() +{ + FILE* fp; + fp = fopen("residual.dat", "w"); + + if (fp == NULL) { + printf("Error!\n"); + exit(EXIT_FAILURE); + } + + return fp; + +} + +void writeResidual(FILE* fp, double ts, double res) +{ + fprintf(fp, "%f, %f\n", ts, res); +} diff --git a/BasicSolver/2D-mpi/src/progress.h b/BasicSolver/2D-mpi/src/progress.h index 1cdb9a3..6803315 100644 --- a/BasicSolver/2D-mpi/src/progress.h +++ b/BasicSolver/2D-mpi/src/progress.h @@ -10,5 +10,6 @@ extern void initProgress(double); extern void printProgress(double); extern void stopProgress(); - +extern FILE* initResidualWriter(void); +extern void writeResidual(FILE*, double, double); #endif diff --git a/BasicSolver/2D-mpi/src/solver-mg.c b/BasicSolver/2D-mpi/src/solver-mg.c index 44cd039..23d64f3 100644 --- a/BasicSolver/2D-mpi/src/solver-mg.c +++ b/BasicSolver/2D-mpi/src/solver-mg.c @@ -14,9 +14,9 @@ #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)] +#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) { @@ -109,7 +109,6 @@ static void setBoundaryCondition(Solver* s, double* p, int imaxLocal, int jmaxLo #endif } - static double smooth(Solver* s, double* p, double* rhs, int level, Comm* comm) { int imaxLocal = comm->imaxLocal; @@ -144,8 +143,7 @@ static double smooth(Solver* s, double* p, double* rhs, int level, Comm* comm) 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)); + (P(i, j + 1) - 2.0 * P(i, j) + P(i, j - 1)) * idy2)); } isw = 3 - isw; } @@ -231,10 +229,10 @@ static double multiGrid(Solver* s, double* p, double* rhs, int level, Comm* comm restrictMG(s, level, comm); Comm newcomm; - commUpdateDatatypes(s->comm, &newcomm, comm->imaxLocal / 2, comm->jmaxLocal / 2); + commUpdateDatatypes(s->comm, &newcomm, comm->imaxLocal, comm->jmaxLocal); + // 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); @@ -290,7 +288,7 @@ void initSolver(Solver* s, Discretization* d, Parameter* p) } } -void solve(Solver* s, double* p, double* rhs) +double solve(Solver* s, double* p, double* rhs) { double res = multiGrid(s, p, rhs, 0, s->comm); @@ -299,4 +297,6 @@ void solve(Solver* s, double* p, double* rhs) printf("Residuum: %.6f\n", res); } #endif + + return res; } diff --git a/BasicSolver/2D-mpi/src/solver-rb.c b/BasicSolver/2D-mpi/src/solver-rb.c index 654f59d..622a11e 100644 --- a/BasicSolver/2D-mpi/src/solver-rb.c +++ b/BasicSolver/2D-mpi/src/solver-rb.c @@ -23,7 +23,7 @@ void initSolver(Solver* s, Discretization* d, Parameter* p) s->comm = &d->comm; } -void solve(Solver* s, double* p, double* rhs) +double solve(Solver* s, double* p, double* rhs) { int imax = s->grid->imax; int jmax = s->grid->jmax; @@ -101,4 +101,6 @@ void solve(Solver* s, double* p, double* rhs) printf("Solver took %d iterations to reach %f\n", it, sqrt(res)); } #endif + + return res; } diff --git a/BasicSolver/2D-mpi/src/solver.h b/BasicSolver/2D-mpi/src/solver.h index 890e664..dcca9bd 100644 --- a/BasicSolver/2D-mpi/src/solver.h +++ b/BasicSolver/2D-mpi/src/solver.h @@ -25,5 +25,5 @@ typedef struct { } Solver; void initSolver(Solver*, Discretization*, Parameter*); -void solve(Solver*, double*, double*); +double solve(Solver*, double*, double*); #endif diff --git a/BasicSolver/2D-seq/Makefile b/BasicSolver/2D-seq/Makefile index 979def5..d33dd17 100644 --- a/BasicSolver/2D-seq/Makefile +++ b/BasicSolver/2D-seq/Makefile @@ -38,9 +38,22 @@ $(BUILD_DIR)/%.s: %.c $(info ===> GENERATE ASM $@) $(CC) -S $(CPPFLAGS) $(CFLAGS) $< -o $@ -.PHONY: clean distclean tags info asm format +.PHONY: clean distclean vis vis_clean tags info asm format -clean: +vis: + $(info ===> GENERATE VISUALIZATION) + @gnuplot -e "filename='pressure.dat'" ./surface.plot + @gnuplot -e "filename='velocity.dat'" ./vector.plot + @gnuplot -e "filename='residual.dat'" ./residual.plot + +vis_clean: + $(info ===> CLEAN VISUALIZATION) + @rm -f *.dat + @rm -f *.png + @rm -f ./vis_files/*.dat + @rm -f ./vis_files/*.gif + +clean: vis_clean $(info ===> CLEAN) @rm -rf $(BUILD_DIR) @rm -f tags diff --git a/BasicSolver/2D-seq/residual.plot b/BasicSolver/2D-seq/residual.plot new file mode 100644 index 0000000..36fb011 --- /dev/null +++ b/BasicSolver/2D-seq/residual.plot @@ -0,0 +1,9 @@ +set terminal png size 1800,768 enhanced font ,12 +set output 'residual.png' +set datafile separator whitespace +set xlabel "Timestep" +set ylabel "Residual" + +set logscale y 2 + +plot 'residual.dat' using 1:2 title "Residual" \ No newline at end of file diff --git a/BasicSolver/2D-seq/src/main.c b/BasicSolver/2D-seq/src/main.c index 51b9877..5f617fa 100644 --- a/BasicSolver/2D-seq/src/main.c +++ b/BasicSolver/2D-seq/src/main.c @@ -20,8 +20,10 @@ int main(int argc, char** argv) Parameter p; Discretization d; Solver s; - + initParameter(&p); + FILE* fp; + fp = initResidualWriter(); if (argc != 2) { printf("Usage: %s \n", argv[0]); @@ -32,7 +34,7 @@ int main(int argc, char** argv) printParameter(&p); initDiscretization(&d, &p); initSolver(&s, &d, &p); - + #ifndef VERBOSE initProgress(d.te); #endif @@ -41,8 +43,9 @@ int main(int argc, char** argv) double te = d.te; double t = 0.0; int nt = 0; + double res = 0.0; - timeStart = getTimeStamp(); + timeStart = getTimeStamp(); while (t <= te) { if (tau > 0.0) computeTimestep(&d); setBoundaryConditions(&d); @@ -50,8 +53,11 @@ int main(int argc, char** argv) computeFG(&d); computeRHS(&d); if (nt % 100 == 0) normalizePressure(&d); - solve(&s, d.p, d.rhs); + res = solve(&s, d.p, d.rhs); adaptUV(&d); + + writeResidual(fp, t, res); + t += d.dt; nt++; @@ -61,6 +67,7 @@ int main(int argc, char** argv) printProgress(t); #endif } + fclose(fp); timeStop = getTimeStamp(); stopProgress(); printf("Solution took %.2fs\n", timeStop - timeStart); diff --git a/BasicSolver/2D-seq/src/progress.c b/BasicSolver/2D-seq/src/progress.c index d5393c4..bcf6e5f 100644 --- a/BasicSolver/2D-seq/src/progress.c +++ b/BasicSolver/2D-seq/src/progress.c @@ -49,3 +49,22 @@ void stopProgress() printf("\n"); fflush(stdout); } + +FILE* initResidualWriter() +{ + FILE* fp; + fp = fopen("residual.dat", "w"); + + if (fp == NULL) { + printf("Error!\n"); + exit(EXIT_FAILURE); + } + + return fp; + +} + +void writeResidual(FILE* fp, double ts, double res) +{ + fprintf(fp, "%f, %f\n", ts, res); +} diff --git a/BasicSolver/2D-seq/src/progress.h b/BasicSolver/2D-seq/src/progress.h index 314e921..a79ecfe 100644 --- a/BasicSolver/2D-seq/src/progress.h +++ b/BasicSolver/2D-seq/src/progress.h @@ -10,5 +10,7 @@ extern void initProgress(double); extern void printProgress(double); extern void stopProgress(void); +extern FILE* initResidualWriter(void); +extern void writeResidual(FILE*, double, double); #endif diff --git a/BasicSolver/2D-seq/src/solver-mg.c b/BasicSolver/2D-seq/src/solver-mg.c index 92537b8..acd3d3b 100644 --- a/BasicSolver/2D-seq/src/solver-mg.c +++ b/BasicSolver/2D-seq/src/solver-mg.c @@ -209,11 +209,13 @@ void initSolver(Solver* s, Discretization* d, Parameter* p) } } -void solve(Solver* s, double* p, double* rhs) +double solve(Solver* s, double* p, double* rhs) { double res = multiGrid(s, p, rhs, 0, s->grid->imax, s->grid->jmax); #ifdef VERBOSE printf("Residuum: %.6f\n", res); #endif + +return res; } diff --git a/BasicSolver/2D-seq/src/solver-rb.c b/BasicSolver/2D-seq/src/solver-rb.c index 265456a..a9facf4 100644 --- a/BasicSolver/2D-seq/src/solver-rb.c +++ b/BasicSolver/2D-seq/src/solver-rb.c @@ -15,7 +15,7 @@ void initSolver(Solver* s, Discretization* d, Parameter* p) s->omega = p->omg; } -void solve(Solver* solver, double* p, double* rhs) +double solve(Solver* solver, double* p, double* rhs) { int imax = solver->grid->imax; int jmax = solver->grid->jmax; @@ -73,4 +73,6 @@ void solve(Solver* solver, double* p, double* rhs) #ifdef VERBOSE printf("Solver took %d iterations to reach %f\n", it, sqrt(res)); #endif + + return res; } diff --git a/BasicSolver/2D-seq/src/solver-sor.c b/BasicSolver/2D-seq/src/solver-sor.c index 3d21ea5..21af4eb 100644 --- a/BasicSolver/2D-seq/src/solver-sor.c +++ b/BasicSolver/2D-seq/src/solver-sor.c @@ -15,7 +15,7 @@ void initSolver(Solver* s, Discretization* d, Parameter* p) s->omega = p->omg; } -void solve(Solver* solver, double* p, double* rhs) +double solve(Solver* solver, double* p, double* rhs) { int imax = solver->grid->imax; int jmax = solver->grid->jmax; @@ -65,4 +65,6 @@ void solve(Solver* solver, double* p, double* rhs) #ifdef VERBOSE printf("Solver took %d iterations to reach %f\n", it, sqrt(res)); #endif + + return res; } diff --git a/BasicSolver/2D-seq/src/solver.h b/BasicSolver/2D-seq/src/solver.h index 955b88d..0343f6e 100644 --- a/BasicSolver/2D-seq/src/solver.h +++ b/BasicSolver/2D-seq/src/solver.h @@ -22,6 +22,6 @@ typedef struct { } Solver; extern void initSolver(Solver*, Discretization*, Parameter*); -extern void solve(Solver*, double*, double*); +extern double solve(Solver*, double*, double*); #endif diff --git a/BasicSolver/3D-mpi/Makefile b/BasicSolver/3D-mpi/Makefile index d870e05..883bb88 100644 --- a/BasicSolver/3D-mpi/Makefile +++ b/BasicSolver/3D-mpi/Makefile @@ -43,9 +43,19 @@ $(BUILD_DIR)/%.s: %.c $(info ===> GENERATE ASM $@) $(CC) -S $(CPPFLAGS) $(CFLAGS) $< -o $@ -.PHONY: clean distclean tags info asm format +.PHONY: clean distclean vis vis_clean tags info asm format -clean: +vis: + $(info ===> GENERATE VISUALIZATION) + @gnuplot -e "filename='residual.dat'" ./residual.plot + +vis_clean: + $(info ===> CLEAN VISUALIZATION) + @rm -f *.dat + @rm -f *.vtk + @rm -f *.png + +clean: vis_clean $(info ===> CLEAN) @rm -rf $(BUILD_DIR) @rm -f tags diff --git a/BasicSolver/3D-mpi/canal.par b/BasicSolver/3D-mpi/canal.par index 597c5b5..b8d35c0 100644 --- a/BasicSolver/3D-mpi/canal.par +++ b/BasicSolver/3D-mpi/canal.par @@ -32,8 +32,8 @@ xlength 30.0 # domain size in x-direction ylength 4.0 # domain size in y-direction zlength 4.0 # domain size in z-direction imax 200 # number of interior cells in x-direction -jmax 50 # number of interior cells in y-direction -kmax 50 # number of interior cells in z-direction +jmax 40 # number of interior cells in y-direction +kmax 40 # number of interior cells in z-direction # Time Data: # --------- @@ -45,8 +45,8 @@ tau 0.5 # safety factor for time stepsize control (<0 constant delt) # Multigrid data: # --------- -levels 3 # Multigrid levels -presmooth 5 # Pre-smoothning iterations +levels 2 # Multigrid levels +presmooth 20 # Pre-smoothning iterations postsmooth 5 # Post-smoothning iterations # Pressure Iteration Data: diff --git a/BasicSolver/3D-mpi/config.mk b/BasicSolver/3D-mpi/config.mk index 0adc565..0dbe440 100644 --- a/BasicSolver/3D-mpi/config.mk +++ b/BasicSolver/3D-mpi/config.mk @@ -1,5 +1,6 @@ # Supported: GCC, CLANG, ICC TAG ?= ICC +# Supported: true, false ENABLE_MPI ?= true ENABLE_OPENMP ?= false # Supported: rb, mg diff --git a/BasicSolver/3D-mpi/dcavity.par b/BasicSolver/3D-mpi/dcavity.par index 303c292..77f1991 100644 --- a/BasicSolver/3D-mpi/dcavity.par +++ b/BasicSolver/3D-mpi/dcavity.par @@ -18,7 +18,7 @@ gx 0.0 # Body forces (e.g. gravity) gy 0.0 # gz 0.0 # -re 1000.0 # Reynolds number +re 1000.0 # Reynolds number u_init 0.0 # initial value for velocity in x-direction v_init 0.0 # initial value for velocity in y-direction @@ -46,7 +46,7 @@ tau 0.5 # safety factor for time stepsize control (<0 constant delt) # --------- levels 3 # Multigrid levels -presmooth 10 # Pre-smoothning iterations +presmooth 20 # Pre-smoothning iterations postsmooth 5 # Post-smoothning iterations # Pressure Iteration Data: diff --git a/BasicSolver/3D-mpi/residual.plot b/BasicSolver/3D-mpi/residual.plot new file mode 100644 index 0000000..36fb011 --- /dev/null +++ b/BasicSolver/3D-mpi/residual.plot @@ -0,0 +1,9 @@ +set terminal png size 1800,768 enhanced font ,12 +set output 'residual.png' +set datafile separator whitespace +set xlabel "Timestep" +set ylabel "Residual" + +set logscale y 2 + +plot 'residual.dat' using 1:2 title "Residual" \ No newline at end of file diff --git a/BasicSolver/3D-mpi/src/comm.c b/BasicSolver/3D-mpi/src/comm.c index b11e9fd..df9f795 100644 --- a/BasicSolver/3D-mpi/src/comm.c +++ b/BasicSolver/3D-mpi/src/comm.c @@ -167,11 +167,12 @@ static void assembleResult(Comm* c, MPI_Waitall(numRequests, requests, MPI_STATUSES_IGNORE); } -static int sum(int* sizes, int position) +// subroutines local to this module +static int sum(int* sizes, int init, int offset, int coord) { int sum = 0; - for (int i = 0; i < position; i++) { + for (int i = init - offset; coord > 0; i -= offset, --coord) { sum += sizes[i]; } @@ -350,9 +351,13 @@ void commCollectResult(Comm* c, for (int i = 0; i < c->size; i++) { int coords[NCORDS]; MPI_Cart_coords(c->comm, i, NDIMS, coords); - offset[i * NDIMS + IDIM] = sum(imaxLocalAll, coords[ICORD]); - offset[i * NDIMS + JDIM] = sum(jmaxLocalAll, coords[JCORD]); - offset[i * NDIMS + KDIM] = sum(kmaxLocalAll, coords[KCORD]); + offset[i * NDIMS + IDIM] = sum(imaxLocalAll, + i, + c->dims[IDIM] * c->dims[JDIM], + coords[ICORD]); + offset[i * NDIMS + JDIM] = sum(jmaxLocalAll, i, c->dims[IDIM], coords[JCORD]); + offset[i * NDIMS + KDIM] = sum(kmaxLocalAll, i, 1, coords[KCORD]); + printf("Rank: %d, Coords(k,j,i): %d %d %d, Size(k,j,i): %d %d %d, " "Offset(k,j,i): %d %d %d\n", i, @@ -561,9 +566,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->coords[IDIM], dims[ICORD], imax); + c->imaxLocal = sizeOfRank(c->coords[KDIM], dims[ICORD], imax); c->jmaxLocal = sizeOfRank(c->coords[JDIM], dims[JCORD], jmax); - c->kmaxLocal = sizeOfRank(c->coords[KDIM], dims[KCORD], kmax); + c->kmaxLocal = sizeOfRank(c->coords[IDIM], dims[KCORD], kmax); // setup buffer types for communication setupCommunication(c, LEFT, BULK); @@ -608,9 +613,12 @@ void commUpdateDatatypes( printf("\nNull communicator. Duplication failed !!\n"); } - newcomm->imaxLocal = imaxLocal; - newcomm->jmaxLocal = jmaxLocal; - newcomm->kmaxLocal = kmaxLocal; + newcomm->rank = oldcomm->rank; + newcomm->size = oldcomm->size; + + newcomm->imaxLocal = imaxLocal / 2; + newcomm->jmaxLocal = jmaxLocal / 2; + newcomm->kmaxLocal = kmaxLocal / 2; setupCommunication(newcomm, LEFT, BULK); setupCommunication(newcomm, LEFT, HALO); diff --git a/BasicSolver/3D-mpi/src/main.c b/BasicSolver/3D-mpi/src/main.c index 944af32..3b68f70 100644 --- a/BasicSolver/3D-mpi/src/main.c +++ b/BasicSolver/3D-mpi/src/main.c @@ -25,6 +25,8 @@ int main(int argc, char** argv) commInit(&d.comm, argc, argv); initParameter(&p); + FILE* fp; + if (commIsMaster(&d.comm)) fp = initResidualWriter(); if (argc != 2) { printf("Usage: %s \n", argv[0]); @@ -37,8 +39,8 @@ int main(int argc, char** argv) if (commIsMaster(&d.comm)) { printParameter(&p); } + initDiscretization(&d, &p); - initSolver(&s, &d, &p); #ifndef VERBOSE @@ -49,6 +51,7 @@ int main(int argc, char** argv) double te = d.te; double t = 0.0; int nt = 0; + double res = 0.0; timeStart = getTimeStamp(); while (t <= te) { @@ -58,8 +61,11 @@ int main(int argc, char** argv) computeFG(&d); computeRHS(&d); if (nt % 100 == 0) normalizePressure(&d); - solve(&s, d.p, d.rhs); + res = solve(&s, d.p, d.rhs); adaptUV(&d); + + if (commIsMaster(&d.comm)) writeResidual(fp, t, res); + t += d.dt; nt++; @@ -87,6 +93,8 @@ int main(int argc, char** argv) vtkVector(&opts, "velocity", (VtkVector) { d.u, d.v, d.w }); vtkClose(&opts); #else + if (commIsMaster(&d.comm)) fclose(fp); + double *pg, *ug, *vg, *wg; if (commIsMaster(s.comm)) { diff --git a/BasicSolver/3D-mpi/src/parameter.c b/BasicSolver/3D-mpi/src/parameter.c index c6cde26..3f98561 100644 --- a/BasicSolver/3D-mpi/src/parameter.c +++ b/BasicSolver/3D-mpi/src/parameter.c @@ -69,6 +69,9 @@ void readParameter(Parameter* param, const char* filename) PARSE_INT(jmax); PARSE_INT(kmax); PARSE_INT(itermax); + PARSE_INT(levels); + PARSE_INT(presmooth); + PARSE_INT(postsmooth); PARSE_REAL(eps); PARSE_REAL(omg); PARSE_REAL(re); diff --git a/BasicSolver/3D-mpi/src/progress.c b/BasicSolver/3D-mpi/src/progress.c index f03dc02..8249ed1 100644 --- a/BasicSolver/3D-mpi/src/progress.c +++ b/BasicSolver/3D-mpi/src/progress.c @@ -4,12 +4,12 @@ * Use of this source code is governed by a MIT style * license that can be found in the LICENSE file. */ +#include "progress.h" #include #include +#include #include -#include "progress.h" - static double _end; static int _current; @@ -48,3 +48,18 @@ void stopProgress() printf("\n"); fflush(stdout); } + +FILE* initResidualWriter() +{ + FILE* fp; + fp = fopen("residual.dat", "w"); + + if (fp == NULL) { + printf("Error!\n"); + exit(EXIT_FAILURE); + } + + return fp; +} + +void writeResidual(FILE* fp, double ts, double res) { fprintf(fp, "%f, %f\n", ts, res); } \ No newline at end of file diff --git a/BasicSolver/3D-mpi/src/progress.h b/BasicSolver/3D-mpi/src/progress.h index 314e921..d612ebe 100644 --- a/BasicSolver/3D-mpi/src/progress.h +++ b/BasicSolver/3D-mpi/src/progress.h @@ -4,11 +4,14 @@ * Use of this source code is governed by a MIT-style * license that can be found in the LICENSE file. */ + +#include + #ifndef __PROGRESS_H_ #define __PROGRESS_H_ - extern void initProgress(double); extern void printProgress(double); extern void stopProgress(void); - +extern FILE* initResidualWriter(void); +extern void writeResidual(FILE*, double, double); #endif diff --git a/BasicSolver/3D-mpi/src/solver-mg.c b/BasicSolver/3D-mpi/src/solver-mg.c index d3f532a..27b1b18 100644 --- a/BasicSolver/3D-mpi/src/solver-mg.c +++ b/BasicSolver/3D-mpi/src/solver-mg.c @@ -328,9 +328,9 @@ static double multiGrid(Solver* s, double* p, double* rhs, int level, Comm* comm Comm newcomm; commUpdateDatatypes(s->comm, &newcomm, - comm->imaxLocal / 2, - comm->jmaxLocal / 2, - comm->kmaxLocal / 2); + imaxLocal, + jmaxLocal, + kmaxLocal); // MGSolver on residual and error. multiGrid(s, s->e[level + 1], s->r[level + 1], level + 1, &newcomm); @@ -371,7 +371,7 @@ void initSolver(Solver* s, Discretization* d, Parameter* p) int jmax = s->grid->jmax; int kmax = s->grid->kmax; int levels = s->levels; - printf("Using Multigrid solver with %d levels\n", levels); + if (commIsMaster(s->comm)) printf("Using Multigrid solver with %d levels\n", levels); s->r = malloc(levels * sizeof(double*)); s->e = malloc(levels * sizeof(double*)); @@ -389,7 +389,7 @@ void initSolver(Solver* s, Discretization* d, Parameter* p) } } -void solve(Solver* s, double* p, double* rhs) +double solve(Solver* s, double* p, double* rhs) { double res = multiGrid(s, p, rhs, 0, s->comm); @@ -398,4 +398,6 @@ void solve(Solver* s, double* p, double* rhs) printf("Residuum: %.6f\n", res); } #endif + + return res; } diff --git a/BasicSolver/3D-mpi/src/solver-rb.c b/BasicSolver/3D-mpi/src/solver-rb.c index 54148e9..d64fecd 100644 --- a/BasicSolver/3D-mpi/src/solver-rb.c +++ b/BasicSolver/3D-mpi/src/solver-rb.c @@ -16,7 +16,20 @@ #include "solver.h" #include "util.h" -void solve(Solver* s, double* p, double* rhs) +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; +} + +double solve(Solver* s, double* p, double* rhs) { int imaxLocal = s->comm->imaxLocal; int jmaxLocal = s->comm->jmaxLocal; @@ -157,17 +170,6 @@ void solve(Solver* s, double* p, double* rhs) 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; +return res; } \ No newline at end of file diff --git a/BasicSolver/3D-mpi/src/solver.h b/BasicSolver/3D-mpi/src/solver.h index 09b1374..d9a0c39 100644 --- a/BasicSolver/3D-mpi/src/solver.h +++ b/BasicSolver/3D-mpi/src/solver.h @@ -34,6 +34,6 @@ typedef struct { Comm* comm; } Solver; -extern void solve(Solver* , double* , double* ); +extern double solve(Solver* , double* , double* ); extern void initSolver(Solver*, Discretization*, Parameter*); #endif diff --git a/BasicSolver/3D-seq/Makefile b/BasicSolver/3D-seq/Makefile index 64735f6..65e6d8d 100644 --- a/BasicSolver/3D-seq/Makefile +++ b/BasicSolver/3D-seq/Makefile @@ -40,6 +40,18 @@ $(BUILD_DIR)/%.s: %.c .PHONY: clean distclean tags info asm format +vis: + $(info ===> GENERATE VISUALIZATION) + @gnuplot -e "filename='residual.dat'" ./residual.plot + +vis_clean: + $(info ===> CLEAN VISUALIZATION) + @rm -f *.dat + @rm -f *.vtk + @rm -f *.png + +clean: vis_clean + clean: $(info ===> CLEAN) @rm -rf $(BUILD_DIR) diff --git a/BasicSolver/3D-seq/canal.par b/BasicSolver/3D-seq/canal.par index 118802b..597c5b5 100644 --- a/BasicSolver/3D-seq/canal.par +++ b/BasicSolver/3D-seq/canal.par @@ -38,7 +38,7 @@ 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) diff --git a/BasicSolver/3D-seq/dcavity.par b/BasicSolver/3D-seq/dcavity.par index c3263dd..b668f49 100644 --- a/BasicSolver/3D-seq/dcavity.par +++ b/BasicSolver/3D-seq/dcavity.par @@ -38,7 +38,7 @@ kmax 128 # number of interior cells in z-direction # Time Data: # --------- -te 2.0 # final time +te 10.0 # final time dt 0.02 # time stepsize tau 0.5 # safety factor for time stepsize control (<0 constant delt) @@ -46,7 +46,7 @@ tau 0.5 # safety factor for time stepsize control (<0 constant delt) # --------- levels 3 # Multigrid levels -presmooth 10 # Pre-smoothning iterations +presmooth 20 # Pre-smoothning iterations postsmooth 5 # Post-smoothning iterations # Pressure Iteration Data: diff --git a/BasicSolver/3D-seq/residual.plot b/BasicSolver/3D-seq/residual.plot new file mode 100644 index 0000000..36fb011 --- /dev/null +++ b/BasicSolver/3D-seq/residual.plot @@ -0,0 +1,9 @@ +set terminal png size 1800,768 enhanced font ,12 +set output 'residual.png' +set datafile separator whitespace +set xlabel "Timestep" +set ylabel "Residual" + +set logscale y 2 + +plot 'residual.dat' using 1:2 title "Residual" \ No newline at end of file diff --git a/BasicSolver/3D-seq/src/main.c b/BasicSolver/3D-seq/src/main.c index 129993a..f453fdd 100644 --- a/BasicSolver/3D-seq/src/main.c +++ b/BasicSolver/3D-seq/src/main.c @@ -73,6 +73,9 @@ int main(int argc, char** argv) Solver s; initParameter(&p); + FILE* fp; + fp = initResidualWriter(); + if (argc != 2) { printf("Usage: %s \n", argv[0]); exit(EXIT_SUCCESS); @@ -90,6 +93,7 @@ int main(int argc, char** argv) double te = d.te; double t = 0.0; int nt = 0; + double res = 0.0; timeStart = getTimeStamp(); while (t <= te) { @@ -99,8 +103,11 @@ int main(int argc, char** argv) computeFG(&d); computeRHS(&d); if (nt % 100 == 0) normalizePressure(&d); - solve(&s, d.p, d.rhs); + res = solve(&s, d.p, d.rhs); adaptUV(&d); + + writeResidual(fp, t, res); + t += d.dt; nt++; @@ -123,7 +130,8 @@ int main(int argc, char** argv) ug = allocate(64, bytesize); vg = allocate(64, bytesize); wg = allocate(64, bytesize); - + + fclose(fp); createBulkArrays(&d, pg, ug, vg, wg); VtkOptions opts = { .grid = d.grid }; vtkOpen(&opts, d.problem); diff --git a/BasicSolver/3D-seq/src/progress.c b/BasicSolver/3D-seq/src/progress.c index d5393c4..3523426 100644 --- a/BasicSolver/3D-seq/src/progress.c +++ b/BasicSolver/3D-seq/src/progress.c @@ -49,3 +49,22 @@ void stopProgress() printf("\n"); fflush(stdout); } + +FILE* initResidualWriter() +{ + FILE* fp; + fp = fopen("residual.dat", "w"); + + if (fp == NULL) { + printf("Error!\n"); + exit(EXIT_FAILURE); + } + + return fp; + +} + +void writeResidual(FILE* fp, double ts, double res) +{ + fprintf(fp, "%f, %f\n", ts, res); +} \ No newline at end of file diff --git a/BasicSolver/3D-seq/src/progress.h b/BasicSolver/3D-seq/src/progress.h index 314e921..240c279 100644 --- a/BasicSolver/3D-seq/src/progress.h +++ b/BasicSolver/3D-seq/src/progress.h @@ -10,5 +10,6 @@ extern void initProgress(double); extern void printProgress(double); extern void stopProgress(void); - +extern FILE* initResidualWriter(void); +extern void writeResidual(FILE*, double, double); #endif diff --git a/BasicSolver/3D-seq/src/solver-mg.c b/BasicSolver/3D-seq/src/solver-mg.c index 668312e..743292c 100644 --- a/BasicSolver/3D-seq/src/solver-mg.c +++ b/BasicSolver/3D-seq/src/solver-mg.c @@ -292,11 +292,13 @@ void initSolver(Solver* s, Discretization* d, Parameter* p) } } -void solve(Solver* s, double* p, double* rhs) +double solve(Solver* s, double* p, double* rhs) { double res = multiGrid(s, p, rhs, 0, s->grid->imax, s->grid->jmax, s->grid->kmax); #ifdef VERBOSE printf("Residuum: %.6f\n", res); #endif + +return res; } diff --git a/BasicSolver/3D-seq/src/solver-rb.c b/BasicSolver/3D-seq/src/solver-rb.c index 6f954f5..6332cce 100644 --- a/BasicSolver/3D-seq/src/solver-rb.c +++ b/BasicSolver/3D-seq/src/solver-rb.c @@ -15,7 +15,7 @@ void initSolver(Solver* s, Discretization* d, Parameter* p) s->omega = p->omg; } -void solve(Solver* s, double* p, double* rhs) +double solve(Solver* s, double* p, double* rhs) { int imax = s->grid->imax; int jmax = s->grid->jmax; @@ -96,4 +96,6 @@ void solve(Solver* s, double* p, double* rhs) #ifdef VERBOSE printf("Solver took %d iterations to reach %f\n", it, sqrt(res)); #endif + +return res; } diff --git a/BasicSolver/3D-seq/src/solver.h b/BasicSolver/3D-seq/src/solver.h index 9a6a2f6..e5f1b38 100644 --- a/BasicSolver/3D-seq/src/solver.h +++ b/BasicSolver/3D-seq/src/solver.h @@ -22,6 +22,6 @@ typedef struct { } Solver; extern void initSolver(Solver*, Discretization*, Parameter*); -extern void solve(Solver*, double*, double*); +extern double solve(Solver*, double*, double*); #endif