From 5a872d05331041c19020e4c7ffda45f333ca9dd0 Mon Sep 17 00:00:00 2001 From: Jan Eitzinger Date: Mon, 4 Mar 2024 14:29:49 +0100 Subject: [PATCH] Separate discretization and solver. Port Multigrid solver. --- BasicSolver/2D-mpi/Makefile | 2 +- BasicSolver/2D-seq/Makefile | 3 +- BasicSolver/2D-seq/config.mk | 10 +- BasicSolver/2D-seq/dcavity.par | 4 +- BasicSolver/2D-seq/include_CLANG.mk | 12 +- BasicSolver/2D-seq/src/discretization.c | 435 ++++++++++++++++++ BasicSolver/2D-seq/src/discretization.h | 40 ++ BasicSolver/2D-seq/src/main.c | 51 +-- BasicSolver/2D-seq/src/parameter.c | 6 + BasicSolver/2D-seq/src/parameter.h | 4 +- BasicSolver/2D-seq/src/solver-mg.c | 192 ++++++++ BasicSolver/2D-seq/src/solver-sor.c | 129 ++++++ BasicSolver/2D-seq/src/solver.c | 564 ------------------------ BasicSolver/2D-seq/src/solver.h | 38 +- BasicSolver/2D-seq/src/util.h | 12 +- 15 files changed, 863 insertions(+), 639 deletions(-) create mode 100644 BasicSolver/2D-seq/src/discretization.c create mode 100644 BasicSolver/2D-seq/src/discretization.h create mode 100644 BasicSolver/2D-seq/src/solver-mg.c create mode 100644 BasicSolver/2D-seq/src/solver-sor.c delete mode 100644 BasicSolver/2D-seq/src/solver.c diff --git a/BasicSolver/2D-mpi/Makefile b/BasicSolver/2D-mpi/Makefile index 3d972dc..2b23422 100644 --- a/BasicSolver/2D-mpi/Makefile +++ b/BasicSolver/2D-mpi/Makefile @@ -1,5 +1,5 @@ #======================================================================================= -# Copyright (C) NHR@FAU, University Erlangen-Nuremberg. +# Copyright (C) NHR@FAU, University Erlangen-Nuremberg. # All rights reserved. # Use of this source code is governed by a MIT-style # license that can be found in the LICENSE file. diff --git a/BasicSolver/2D-seq/Makefile b/BasicSolver/2D-seq/Makefile index 60f4df1..979def5 100644 --- a/BasicSolver/2D-seq/Makefile +++ b/BasicSolver/2D-seq/Makefile @@ -1,5 +1,5 @@ #======================================================================================= -# Copyright (C) NHR@FAU, University Erlangen-Nuremberg. +# Copyright (C) NHR@FAU, University Erlangen-Nuremberg. # All rights reserved. # Use of this source code is governed by a MIT-style # license that can be found in the LICENSE file. @@ -21,6 +21,7 @@ VPATH = $(SRC_DIR) 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)/solver-$(SOLVER).o SOURCES = $(SRC) $(wildcard $(SRC_DIR)/*.h) CPPFLAGS := $(CPPFLAGS) $(DEFINES) $(OPTIONS) $(INCLUDES) diff --git a/BasicSolver/2D-seq/config.mk b/BasicSolver/2D-seq/config.mk index d228f66..a782baf 100644 --- a/BasicSolver/2D-seq/config.mk +++ b/BasicSolver/2D-seq/config.mk @@ -1,12 +1,12 @@ # Supported: GCC, CLANG, ICC TAG ?= CLANG ENABLE_OPENMP ?= false +# Supported: sor, mg +SOLVER ?= sor +# Run in debug settings +DEBUG ?= false #Feature options OPTIONS += -DARRAY_ALIGNMENT=64 -# OPTIONS += -DVERBOSE +#OPTIONS += -DVERBOSE #OPTIONS += -DDEBUG -#OPTIONS += -DBOUNDCHECK -#OPTIONS += -DVERBOSE_AFFINITY -#OPTIONS += -DVERBOSE_DATASIZE -#OPTIONS += -DVERBOSE_TIMER diff --git a/BasicSolver/2D-seq/dcavity.par b/BasicSolver/2D-seq/dcavity.par index 67287ef..b1c238d 100644 --- a/BasicSolver/2D-seq/dcavity.par +++ b/BasicSolver/2D-seq/dcavity.par @@ -36,11 +36,13 @@ te 10.0 # final time dt 0.02 # time stepsize tau 0.5 # safety factor for time stepsize control (<0 constant delt) -# Pressure Iteration Data: +# Solver Data: # ----------------------- itermax 1000 # maximal number of pressure iteration in one time step 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/2D-seq/include_CLANG.mk b/BasicSolver/2D-seq/include_CLANG.mk index 466f441..a124053 100644 --- a/BasicSolver/2D-seq/include_CLANG.mk +++ b/BasicSolver/2D-seq/include_CLANG.mk @@ -2,16 +2,18 @@ CC = clang GCC = cc LINKER = $(CC) -ifeq ($(ENABLE_OPENMP),true) +ifeq ($(strip $(ENABLE_OPENMP)),true) OPENMP = -fopenmp #OPENMP = -Xpreprocessor -fopenmp #required on Macos with homebrew libomp LIBS = # -lomp endif +ifeq ($(strip $(DEBUG)),true) +CFLAGS = -O0 -g -std=c17 +else +CFLAGS = -O3 -std=c17 $(OPENMP) +endif VERSION = --version -# CFLAGS = -O3 -std=c17 $(OPENMP) -CFLAGS = -Ofast -std=c17 -#CFLAGS = -Ofast -fnt-store=aggressive -std=c99 $(OPENMP) #AMD CLANG LFLAGS = $(OPENMP) -lm -DEFINES = -D_GNU_SOURCE# -DDEBUG +DEFINES = -D_GNU_SOURCE INCLUDES = diff --git a/BasicSolver/2D-seq/src/discretization.c b/BasicSolver/2D-seq/src/discretization.c new file mode 100644 index 0000000..184d9ac --- /dev/null +++ b/BasicSolver/2D-seq/src/discretization.c @@ -0,0 +1,435 @@ +/* + * 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 "discretization.h" +#include "parameter.h" +#include "util.h" + +static void print(Discretization* d, double* grid) +{ + int imax = d->grid.imax; + + for (int j = 0; j < d->grid.jmax + 2; j++) { + printf("%02d: ", j); + for (int i = 0; i < d->grid.imax + 2; i++) { + printf("%12.8f ", grid[j * (imax + 2) + i]); + } + printf("\n"); + } + fflush(stdout); +} + +static void printConfig(Discretization* d) +{ + printf("Parameters for #%s#\n", d->problem); + printf("Boundary conditions Left:%d Right:%d Bottom:%d Top:%d\n", + d->bcLeft, + d->bcRight, + d->bcBottom, + d->bcTop); + printf("\tReynolds number: %.2f\n", d->re); + printf("\tGx Gy: %.2f %.2f\n", d->gx, d->gy); + printf("Geometry data:\n"); + printf("\tDomain box size (x, y): %.2f, %.2f\n", d->grid.xlength, d->grid.ylength); + printf("\tCells (x, y): %d, %d\n", d->grid.imax, d->grid.jmax); + printf("Timestep parameters:\n"); + printf("\tDefault stepsize: %.2f, Final time %.2f\n", d->dt, d->te); + printf("\tdt bound: %.6f\n", d->dtBound); + printf("\tTau factor: %.2f\n", d->tau); + printf("Iterative d parameters:\n"); + printf("\tgamma factor: %f\n", d->gamma); +} + +void initDiscretization(Discretization* d, Parameter* p) +{ + d->problem = p->name; + d->bcLeft = p->bcLeft; + d->bcRight = p->bcRight; + d->bcBottom = p->bcBottom; + d->bcTop = p->bcTop; + d->grid.imax = p->imax; + d->grid.jmax = p->jmax; + d->grid.xlength = p->xlength; + d->grid.ylength = p->ylength; + d->grid.dx = p->xlength / p->imax; + d->grid.dy = p->ylength / p->jmax; + d->re = p->re; + d->gx = p->gx; + d->gy = p->gy; + d->dt = p->dt; + d->te = p->te; + d->tau = p->tau; + d->gamma = p->gamma; + + int imax = d->grid.imax; + int jmax = d->grid.jmax; + size_t size = (imax + 2) * (jmax + 2) * sizeof(double); + d->u = allocate(64, size); + d->v = allocate(64, size); + d->p = allocate(64, size); + d->rhs = allocate(64, size); + d->f = allocate(64, size); + d->g = allocate(64, size); + + for (int i = 0; i < (imax + 2) * (jmax + 2); i++) { + d->u[i] = p->u_init; + d->v[i] = p->v_init; + d->p[i] = p->p_init; + d->rhs[i] = 0.0; + d->f[i] = 0.0; + d->g[i] = 0.0; + } + + double dx = d->grid.dx; + double dy = d->grid.dy; + double invSqrSum = 1.0 / (dx * dx) + 1.0 / (dy * dy); + d->dtBound = 0.5 * d->re * 1.0 / invSqrSum; +#ifdef VERBOSE + printConfig(d); +#endif +} + +void computeRHS(Discretization* d) +{ + int imax = d->grid.imax; + int jmax = d->grid.jmax; + double idx = 1.0 / d->grid.dx; + double idy = 1.0 / d->grid.dy; + double idt = 1.0 / d->dt; + double* rhs = d->rhs; + double* f = d->f; + double* g = d->g; + + for (int j = 1; j < jmax + 1; j++) { + for (int i = 1; i < imax + 1; i++) { + RHS(i, j) = idt * + ((F(i, j) - F(i - 1, j)) * idx + (G(i, j) - G(i, j - 1)) * idy); + } + } +} + +static double maxElement(Discretization* d, double* m) +{ + int size = (d->grid.imax + 2) * (d->grid.jmax + 2); + double maxval = DBL_MIN; + + for (int i = 0; i < size; i++) { + maxval = MAX(maxval, fabs(m[i])); + } + + return maxval; +} + +void normalizePressure(Discretization* d) +{ + int size = (d->grid.imax + 2) * (d->grid.jmax + 2); + double* p = d->p; + double avgP = 0.0; + + for (int i = 0; i < size; i++) { + avgP += p[i]; + } + avgP /= size; + + for (int i = 0; i < size; i++) { + p[i] = p[i] - avgP; + } +} + +void computeTimestep(Discretization* d) +{ + double dt = d->dtBound; + double dx = d->grid.dx; + double dy = d->grid.dy; + double umax = maxElement(d, d->u); + double vmax = maxElement(d, d->v); + + if (umax > 0) { + dt = (dt > dx / umax) ? dx / umax : dt; + } + if (vmax > 0) { + dt = (dt > dy / vmax) ? dy / vmax : dt; + } + + d->dt = dt * d->tau; +} + +void setBoundaryConditions(Discretization* d) +{ + int imax = d->grid.imax; + int jmax = d->grid.jmax; + double* u = d->u; + double* v = d->v; + + // Left boundary + switch (d->bcLeft) { + case NOSLIP: + for (int j = 1; j < jmax + 1; j++) { + U(0, j) = 0.0; + V(0, j) = -V(1, j); + } + break; + case SLIP: + for (int j = 1; j < jmax + 1; j++) { + U(0, j) = 0.0; + V(0, j) = V(1, j); + } + break; + case OUTFLOW: + for (int j = 1; j < jmax + 1; j++) { + U(0, j) = U(1, j); + V(0, j) = V(1, j); + } + break; + case PERIODIC: + break; + } + + // Right boundary + switch (d->bcRight) { + case NOSLIP: + for (int j = 1; j < jmax + 1; j++) { + U(imax, j) = 0.0; + V(imax + 1, j) = -V(imax, j); + } + break; + case SLIP: + for (int j = 1; j < jmax + 1; j++) { + U(imax, j) = 0.0; + V(imax + 1, j) = V(imax, j); + } + break; + case OUTFLOW: + for (int j = 1; j < jmax + 1; j++) { + U(imax, j) = U(imax - 1, j); + V(imax + 1, j) = V(imax, j); + } + break; + case PERIODIC: + break; + } + + // Bottom boundary + switch (d->bcBottom) { + case NOSLIP: + for (int i = 1; i < imax + 1; i++) { + V(i, 0) = 0.0; + U(i, 0) = -U(i, 1); + } + break; + case SLIP: + for (int i = 1; i < imax + 1; i++) { + V(i, 0) = 0.0; + U(i, 0) = U(i, 1); + } + break; + case OUTFLOW: + for (int i = 1; i < imax + 1; i++) { + U(i, 0) = U(i, 1); + V(i, 0) = V(i, 1); + } + break; + case PERIODIC: + break; + } + + // Top boundary + switch (d->bcTop) { + case NOSLIP: + for (int i = 1; i < imax + 1; i++) { + V(i, jmax) = 0.0; + U(i, jmax + 1) = -U(i, jmax); + } + break; + case SLIP: + for (int i = 1; i < imax + 1; i++) { + V(i, jmax) = 0.0; + U(i, jmax + 1) = U(i, jmax); + } + break; + case OUTFLOW: + for (int i = 1; i < imax + 1; i++) { + U(i, jmax + 1) = U(i, jmax); + V(i, jmax) = V(i, jmax - 1); + } + break; + case PERIODIC: + break; + } +} + +void setSpecialBoundaryCondition(Discretization* d) +{ + int imax = d->grid.imax; + int jmax = d->grid.jmax; + double mDy = d->grid.dy; + double* u = d->u; + + if (strcmp(d->problem, "dcavity") == 0) { + for (int i = 1; i < imax; i++) { + U(i, jmax + 1) = 2.0 - U(i, jmax); + } + } else if (strcmp(d->problem, "canal") == 0) { + double ylength = d->grid.ylength; + double y; + + for (int j = 1; j < jmax + 1; j++) { + y = mDy * (j - 0.5); + U(0, j) = y * (ylength - y) * 4.0 / (ylength * ylength); + } + } +} + +void computeFG(Discretization* d) +{ + double* u = d->u; + double* v = d->v; + double* f = d->f; + double* g = d->g; + int imax = d->grid.imax; + int jmax = d->grid.jmax; + double gx = d->gx; + double gy = d->gy; + double gamma = d->gamma; + double dt = d->dt; + double inverseRe = 1.0 / d->re; + double inverseDx = 1.0 / d->grid.dx; + double inverseDy = 1.0 / d->grid.dy; + double du2dx, dv2dy, duvdx, duvdy; + double du2dx2, du2dy2, dv2dx2, dv2dy2; + + for (int j = 1; j < jmax + 1; j++) { + for (int i = 1; i < imax + 1; i++) { + du2dx = inverseDx * 0.25 * + ((U(i, j) + U(i + 1, j)) * (U(i, j) + U(i + 1, j)) - + (U(i, j) + U(i - 1, j)) * (U(i, j) + U(i - 1, j))) + + gamma * inverseDx * 0.25 * + (fabs(U(i, j) + U(i + 1, j)) * (U(i, j) - U(i + 1, j)) + + fabs(U(i, j) + U(i - 1, j)) * (U(i, j) - U(i - 1, j))); + + duvdy = inverseDy * 0.25 * + ((V(i, j) + V(i + 1, j)) * (U(i, j) + U(i, j + 1)) - + (V(i, j - 1) + V(i + 1, j - 1)) * (U(i, j) + U(i, j - 1))) + + gamma * inverseDy * 0.25 * + (fabs(V(i, j) + V(i + 1, j)) * (U(i, j) - U(i, j + 1)) + + fabs(V(i, j - 1) + V(i + 1, j - 1)) * + (U(i, j) - U(i, j - 1))); + + du2dx2 = inverseDx * inverseDx * (U(i + 1, j) - 2.0 * U(i, j) + U(i - 1, j)); + du2dy2 = inverseDy * inverseDy * (U(i, j + 1) - 2.0 * U(i, j) + U(i, j - 1)); + F(i, j) = U(i, j) + dt * (inverseRe * (du2dx2 + du2dy2) - du2dx - duvdy + gx); + + duvdx = inverseDx * 0.25 * + ((U(i, j) + U(i, j + 1)) * (V(i, j) + V(i + 1, j)) - + (U(i - 1, j) + U(i - 1, j + 1)) * (V(i, j) + V(i - 1, j))) + + gamma * inverseDx * 0.25 * + (fabs(U(i, j) + U(i, j + 1)) * (V(i, j) - V(i + 1, j)) + + fabs(U(i - 1, j) + U(i - 1, j + 1)) * + (V(i, j) - V(i - 1, j))); + + dv2dy = inverseDy * 0.25 * + ((V(i, j) + V(i, j + 1)) * (V(i, j) + V(i, j + 1)) - + (V(i, j) + V(i, j - 1)) * (V(i, j) + V(i, j - 1))) + + gamma * inverseDy * 0.25 * + (fabs(V(i, j) + V(i, j + 1)) * (V(i, j) - V(i, j + 1)) + + fabs(V(i, j) + V(i, j - 1)) * (V(i, j) - V(i, j - 1))); + + dv2dx2 = inverseDx * inverseDx * (V(i + 1, j) - 2.0 * V(i, j) + V(i - 1, j)); + dv2dy2 = inverseDy * inverseDy * (V(i, j + 1) - 2.0 * V(i, j) + V(i, j - 1)); + G(i, j) = V(i, j) + dt * (inverseRe * (dv2dx2 + dv2dy2) - duvdx - dv2dy + gy); + } + } + + /* ---------------------- boundary of F --------------------------- */ + for (int j = 1; j < jmax + 1; j++) { + F(0, j) = U(0, j); + F(imax, j) = U(imax, j); + } + + /* ---------------------- boundary of G --------------------------- */ + for (int i = 1; i < imax + 1; i++) { + G(i, 0) = V(i, 0); + G(i, jmax) = V(i, jmax); + } +} + +void adaptUV(Discretization* d) +{ + int imax = d->grid.imax; + int jmax = d->grid.jmax; + double* p = d->p; + double* u = d->u; + double* v = d->v; + double* f = d->f; + double* g = d->g; + double factorX = d->dt / d->grid.dx; + double factorY = d->dt / d->grid.dy; + + for (int j = 1; j < jmax + 1; j++) { + for (int i = 1; i < imax + 1; i++) { + U(i, j) = F(i, j) - (P(i + 1, j) - P(i, j)) * factorX; + V(i, j) = G(i, j) - (P(i, j + 1) - P(i, j)) * factorY; + } + } +} + +void writeResult(Discretization* d) +{ + int imax = d->grid.imax; + int jmax = d->grid.jmax; + double dx = d->grid.dx; + double dy = d->grid.dy; + double* p = d->p; + double* u = d->u; + double* v = d->v; + double x = 0.0, y = 0.0; + + FILE* fp; + fp = fopen("pressure.dat", "w"); + + if (fp == NULL) { + printf("Error!\n"); + exit(EXIT_FAILURE); + } + + for (int j = 1; j < jmax + 1; j++) { + y = (double)(j - 0.5) * dy; + for (int i = 1; i < imax + 1; i++) { + x = (double)(i - 0.5) * dx; + fprintf(fp, "%.2f %.2f %f\n", x, y, P(i, j)); + } + fprintf(fp, "\n"); + } + + fclose(fp); + + fp = fopen("velocity.dat", "w"); + + if (fp == NULL) { + printf("Error!\n"); + exit(EXIT_FAILURE); + } + + for (int j = 1; j < jmax + 1; j++) { + y = dy * (j - 0.5); + for (int i = 1; i < imax + 1; i++) { + x = dx * (i - 0.5); + double velU = (U(i, j) + U(i - 1, j)) / 2.0; + double velV = (V(i, j) + V(i, j - 1)) / 2.0; + double len = sqrt((velU * velU) + (velV * velV)); + fprintf(fp, "%.2f %.2f %f %f %f\n", x, y, velU, velV, len); + } + } + + fclose(fp); +} diff --git a/BasicSolver/2D-seq/src/discretization.h b/BasicSolver/2D-seq/src/discretization.h new file mode 100644 index 0000000..3007093 --- /dev/null +++ b/BasicSolver/2D-seq/src/discretization.h @@ -0,0 +1,40 @@ +/* + * 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" + +enum BC { NOSLIP = 1, SLIP, OUTFLOW, PERIODIC }; + +typedef struct { + /* geometry and grid information */ + Grid grid; + /* arrays */ + double *p, *rhs; + double *f, *g; + double *u, *v; + /* parameters */ + double re, tau, gamma; + double gx, gy; + /* time stepping */ + double dt, te; + double dtBound; + char* problem; + int bcLeft, bcRight, bcBottom, bcTop; +} 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*); +extern void writeResult(Discretization*); +#endif diff --git a/BasicSolver/2D-seq/src/main.c b/BasicSolver/2D-seq/src/main.c index fe9fee2..a284745 100644 --- a/BasicSolver/2D-seq/src/main.c +++ b/BasicSolver/2D-seq/src/main.c @@ -4,12 +4,11 @@ * 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 "discretization.h" #include "parameter.h" #include "progress.h" #include "solver.h" @@ -17,39 +16,41 @@ int main(int argc, char** argv) { - double S, E; - Parameter params; - Solver solver; - initParameter(¶ms); + double timeStart, timeStop; + Parameter p; + Discretization d; + Solver s; + initParameter(&p); if (argc != 2) { printf("Usage: %s \n", argv[0]); exit(EXIT_SUCCESS); } - readParameter(¶ms, argv[1]); - printParameter(¶ms); - initSolver(&solver, ¶ms); + readParameter(&p, argv[1]); + printParameter(&p); + initDiscretization(&d, &p); + initSolver(&s, &d, &p); #ifndef VERBOSE - initProgress(solver.te); + initProgress(d.te); #endif - double tau = solver.tau; - double te = solver.te; + double tau = d.tau; + double te = d.te; double t = 0.0; int nt = 0; - S = getTimeStamp(); + timeStart = getTimeStamp(); while (t <= te) { - if (tau > 0.0) computeTimestep(&solver); - setBoundaryConditions(&solver); - setSpecialBoundaryCondition(&solver); - computeFG(&solver); - computeRHS(&solver); - if (nt % 100 == 0) normalizePressure(&solver); - solveRB(&solver); - adaptUV(&solver); - t += solver.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 @@ -58,9 +59,9 @@ int main(int argc, char** argv) printProgress(t); #endif } - E = getTimeStamp(); + timeStop = getTimeStamp(); stopProgress(); - printf("Solution took %.2fs\n", E - S); - writeResult(&solver); + printf("Solution took %.2fs\n", timeStop - timeStart); + writeResult(&d); return EXIT_SUCCESS; } diff --git a/BasicSolver/2D-seq/src/parameter.c b/BasicSolver/2D-seq/src/parameter.c index 8385199..777a274 100644 --- a/BasicSolver/2D-seq/src/parameter.c +++ b/BasicSolver/2D-seq/src/parameter.c @@ -24,6 +24,8 @@ void initParameter(Parameter* param) param->re = 100.0; param->gamma = 0.9; param->tau = 0.5; + param->rho = 0.99; + param->levels = 5; } void readParameter(Parameter* param, const char* filename) @@ -61,6 +63,7 @@ void readParameter(Parameter* param, const char* filename) PARSE_INT(imax); PARSE_INT(jmax); PARSE_INT(itermax); + PARSE_INT(levels); PARSE_REAL(eps); PARSE_REAL(omg); PARSE_REAL(re); @@ -78,6 +81,7 @@ void readParameter(Parameter* param, const char* filename) PARSE_REAL(u_init); PARSE_REAL(v_init); PARSE_REAL(p_init); + PARSE_REAL(rho); } } @@ -108,4 +112,6 @@ void printParameter(Parameter* param) printf("\tepsilon (stopping tolerance) : %f\n", param->eps); printf("\tgamma (stopping tolerance) : %f\n", param->gamma); printf("\tomega (SOR relaxation): %f\n", param->omg); + printf("\trho (SOR relaxation): %f\n", param->rho); + printf("\tMultiGrid levels : %d\n", param->levels); } diff --git a/BasicSolver/2D-seq/src/parameter.h b/BasicSolver/2D-seq/src/parameter.h index be28e82..5098794 100644 --- a/BasicSolver/2D-seq/src/parameter.h +++ b/BasicSolver/2D-seq/src/parameter.h @@ -10,8 +10,8 @@ typedef struct { double xlength, ylength; int imax, jmax; - int itermax; - double eps, omg; + int itermax, levels; + double eps, omg, rho; double re, tau, gamma; double te, dt; double gx, gy; diff --git a/BasicSolver/2D-seq/src/solver-mg.c b/BasicSolver/2D-seq/src/solver-mg.c new file mode 100644 index 0000000..721be56 --- /dev/null +++ b/BasicSolver/2D-seq/src/solver-mg.c @@ -0,0 +1,192 @@ +/* + * 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 "discretization.h" +#include "parameter.h" +#include "solver.h" +#include "util.h" + +#define FINEST_LEVEL 0 +#define S(i, j) s[(j) * (imax + 2) + (i)] +#define E(i, j) e[(j) * (imax + 2) + (i)] +#define R(i, j) r[(j) * (imax + 2) + (i)] +#define OLD(i, j) old[(j) * (imax + 2) + (i)] + +static void restrictMG(Solver* s, int level, int imax, int jmax) +{ + double* r = s->r[level + 1]; + double* old = s->r[level]; + + for (int j = 1; j < jmax + 1; j++) { + for (int i = 1; i < imax + 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, int imax, int jmax) +{ + double* old = s->r[level + 1]; + double* e = s->r[level]; + + for (int j = 2; j < jmax + 1; j += 2) { + for (int i = 2; i < imax + 1; i += 2) { + E(i, j) = OLD(i / 2, j / 2); + } + } +} + +static void correct(Solver* s, double* p, int level, int imax, int jmax) +{ + double* e = s->e[level]; + for (int j = 1; j < jmax + 1; ++j) { + for (int i = 1; i < imax + 1; ++i) { + P(i, j) += E(i, j); + } + } +} + +static void setBoundaryCondition(double* p, int imax, int jmax) +{ + for (int i = 1; i < imax + 1; i++) { + P(i, 0) = P(i, 1); + P(i, jmax + 1) = P(i, jmax); + } + + for (int j = 1; j < jmax + 1; j++) { + P(0, j) = P(1, j); + P(imax + 1, j) = P(imax, j); + } +} + +static double smooth(Solver* s, double* p, double* rhs, int level, int imax, int 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; + + for (int j = 1; j < jmax + 1; j++) { + for (int i = isw; i < imax + 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); + + P(i, j) -= (factor * R(i, j)); + res += (R(i, j) * R(i, j)); + } + isw = 3 - isw; + } + jsw = 3 - jsw; + } + + res = res / (double)(imax * jmax); + return res; +} + +void initSolver(Solver* s, Discretization* d, Parameter* p) +{ + s->eps = p->eps; + s->omega = p->omg; + s->itermax = p->itermax; + s->rho = p->rho; + s->levels = p->levels; + s->grid = &d->grid; + + 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; + } + } +} + +double multiGrid(Solver* solver, double* p, double* rhs, int level, int imax, int jmax) +{ + double res = 0.0; + + // coarsest level TODO: Use direct solver? + if (level == (solver->levels - 1)) { + for (int i = 0; i < 5; i++) { + smooth(solver, p, rhs, level, imax, jmax); + } + return res; + } + + // pre-smoothing TODO: Make smoothing steps configurable? + for (int i = 0; i < 5; i++) { + smooth(solver, p, rhs, level, imax, jmax); + if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax); + } + + // restrict + restrictMG(solver, level, imax, jmax); + + // MGSolver on residual and error. + // TODO: What if there is a rest? + multiGrid(solver, + solver->e[level + 1], + solver->r[level], + level + 1, + imax / 2, + jmax / 2); + + // prolongate + prolongate(solver, level, imax, jmax); + + // correct p on finer level using residual + correct(solver, p, level, imax, jmax); + if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax); + + // post-smoothing + for (int i = 0; i < 5; i++) { + res = smooth(solver, p, rhs, level, imax, jmax); + if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax); + } + + return res; +} + +void 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 +} diff --git a/BasicSolver/2D-seq/src/solver-sor.c b/BasicSolver/2D-seq/src/solver-sor.c new file mode 100644 index 0000000..b516e5b --- /dev/null +++ b/BasicSolver/2D-seq/src/solver-sor.c @@ -0,0 +1,129 @@ +/* + * 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 "discretization.h" +#include "solver.h" +#include "util.h" + +void initSolver(Solver* s, Discretization* d, Parameter* p) +{ + s->grid = &d->grid; + s->itermax = p->itermax; + s->eps = p->eps; + s->omega = p->omg; +} + +void solveSOR(Solver* solver, double* p, double* rhs) +{ + int imax = solver->grid->imax; + int jmax = solver->grid->jmax; + double eps = solver->eps; + int itermax = solver->itermax; + double dx2 = solver->grid->dx * solver->grid->dx; + double dy2 = solver->grid->dy * solver->grid->dy; + double idx2 = 1.0 / dx2; + double idy2 = 1.0 / dy2; + double factor = solver->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2); + double epssq = eps * eps; + int it = 0; + double res = 1.0; + + while ((res >= epssq) && (it < itermax)) { + res = 0.0; + + for (int j = 1; j < jmax + 1; j++) { + for (int i = 1; i < imax + 1; i++) { + + double r = RHS(i, j) - + ((P(i + 1, j) - 2.0 * P(i, j) + P(i - 1, j)) * idx2 + + (P(i, j + 1) - 2.0 * P(i, j) + P(i, j - 1)) * idy2); + + P(i, j) -= (factor * r); + res += (r * r); + } + } + + for (int i = 1; i < imax + 1; i++) { + P(i, 0) = P(i, 1); + P(i, jmax + 1) = P(i, jmax); + } + + for (int j = 1; j < jmax + 1; j++) { + P(0, j) = P(1, j); + P(imax + 1, j) = P(imax, j); + } + + res = res / (double)(imax * jmax); +#ifdef DEBUG + printf("%d Residuum: %e\n", it, res); +#endif + it++; + } + +#ifdef VERBOSE + printf("Solver took %d iterations to reach %f\n", it, sqrt(res)); +#endif +} + +void solve(Solver* solver, double* p, double* rhs) +{ + int imax = solver->grid->imax; + int jmax = solver->grid->jmax; + double eps = solver->eps; + int itermax = solver->itermax; + double dx2 = solver->grid->dx * solver->grid->dx; + double dy2 = solver->grid->dy * solver->grid->dy; + double idx2 = 1.0 / dx2; + double idy2 = 1.0 / dy2; + double factor = solver->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2); + double epssq = eps * eps; + int it = 0; + double res = 1.0; + int pass, jsw, isw; + + while ((res >= epssq) && (it < itermax)) { + res = 0.0; + jsw = 1; + + for (pass = 0; pass < 2; pass++) { + isw = jsw; + + for (int j = 1; j < jmax + 1; j++) { + for (int i = isw; i < imax + 1; i += 2) { + + double r = RHS(i, j) - + ((P(i + 1, j) - 2.0 * P(i, j) + P(i - 1, j)) * idx2 + + (P(i, j + 1) - 2.0 * P(i, j) + P(i, j - 1)) * idy2); + + P(i, j) -= (factor * r); + res += (r * r); + } + isw = 3 - isw; + } + jsw = 3 - jsw; + } + + for (int i = 1; i < imax + 1; i++) { + P(i, 0) = P(i, 1); + P(i, jmax + 1) = P(i, jmax); + } + + for (int j = 1; j < jmax + 1; j++) { + P(0, j) = P(1, j); + P(imax + 1, j) = P(imax, j); + } + + res = res / (double)(imax * jmax); +#ifdef DEBUG + printf("%d Residuum: %e\n", it, res); +#endif + it++; + } + +#ifdef VERBOSE + printf("Solver took %d iterations to reach %f\n", it, sqrt(res)); +#endif +} diff --git a/BasicSolver/2D-seq/src/solver.c b/BasicSolver/2D-seq/src/solver.c deleted file mode 100644 index e42c950..0000000 --- a/BasicSolver/2D-seq/src/solver.c +++ /dev/null @@ -1,564 +0,0 @@ -/* - * 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 "parameter.h" -#include "solver.h" -#include "util.h" - -#define P(i, j) p[(j) * (imax + 2) + (i)] -#define F(i, j) f[(j) * (imax + 2) + (i)] -#define G(i, j) g[(j) * (imax + 2) + (i)] -#define U(i, j) u[(j) * (imax + 2) + (i)] -#define V(i, j) v[(j) * (imax + 2) + (i)] -#define RHS(i, j) rhs[(j) * (imax + 2) + (i)] - -static void print(Solver* solver, double* grid) -{ - int imax = solver->imax; - - for (int j = 0; j < solver->jmax + 2; j++) { - printf("%02d: ", j); - for (int i = 0; i < solver->imax + 2; i++) { - printf("%12.8f ", grid[j * (imax + 2) + i]); - } - printf("\n"); - } - fflush(stdout); -} - -static void printConfig(Solver* solver) -{ - printf("Parameters for #%s#\n", solver->problem); - printf("Boundary conditions Left:%d Right:%d Bottom:%d Top:%d\n", - solver->bcLeft, - solver->bcRight, - solver->bcBottom, - solver->bcTop); - printf("\tReynolds number: %.2f\n", solver->re); - printf("\tGx Gy: %.2f %.2f\n", solver->gx, solver->gy); - printf("Geometry data:\n"); - printf("\tDomain box size (x, y): %.2f, %.2f\n", solver->xlength, solver->ylength); - printf("\tCells (x, y): %d, %d\n", solver->imax, solver->jmax); - printf("Timestep parameters:\n"); - printf("\tDefault stepsize: %.2f, Final time %.2f\n", solver->dt, solver->te); - printf("\tdt bound: %.6f\n", solver->dtBound); - printf("\tTau factor: %.2f\n", solver->tau); - printf("Iterative solver parameters:\n"); - printf("\tMax iterations: %d\n", solver->itermax); - printf("\tepsilon (stopping tolerance) : %f\n", solver->eps); - printf("\tgamma factor: %f\n", solver->gamma); - printf("\tomega (SOR relaxation): %f\n", solver->omega); -} - -void initSolver(Solver* solver, Parameter* params) -{ - solver->problem = params->name; - solver->bcLeft = params->bcLeft; - solver->bcRight = params->bcRight; - solver->bcBottom = params->bcBottom; - solver->bcTop = params->bcTop; - solver->imax = params->imax; - solver->jmax = params->jmax; - solver->xlength = params->xlength; - solver->ylength = params->ylength; - solver->dx = params->xlength / params->imax; - solver->dy = params->ylength / params->jmax; - solver->eps = params->eps; - solver->omega = params->omg; - solver->itermax = params->itermax; - solver->re = params->re; - solver->gx = params->gx; - solver->gy = params->gy; - solver->dt = params->dt; - solver->te = params->te; - solver->tau = params->tau; - solver->gamma = params->gamma; - - int imax = solver->imax; - int jmax = solver->jmax; - size_t size = (imax + 2) * (jmax + 2) * sizeof(double); - solver->u = allocate(64, size); - solver->v = allocate(64, size); - solver->p = allocate(64, size); - solver->rhs = allocate(64, size); - solver->f = allocate(64, size); - solver->g = allocate(64, size); - - for (int i = 0; i < (imax + 2) * (jmax + 2); i++) { - solver->u[i] = params->u_init; - solver->v[i] = params->v_init; - solver->p[i] = params->p_init; - solver->rhs[i] = 0.0; - solver->f[i] = 0.0; - solver->g[i] = 0.0; - } - - double dx = solver->dx; - double dy = solver->dy; - double invSqrSum = 1.0 / (dx * dx) + 1.0 / (dy * dy); - solver->dtBound = 0.5 * solver->re * 1.0 / invSqrSum; -#ifdef VERBOSE - printConfig(solver); -#endif -} - -void computeRHS(Solver* solver) -{ - int imax = solver->imax; - int jmax = solver->jmax; - double idx = 1.0 / solver->dx; - double idy = 1.0 / solver->dy; - double idt = 1.0 / solver->dt; - double* rhs = solver->rhs; - double* f = solver->f; - double* g = solver->g; - - for (int j = 1; j < jmax + 1; j++) { - for (int i = 1; i < imax + 1; i++) { - RHS(i, j) = idt * - ((F(i, j) - F(i - 1, j)) * idx + (G(i, j) - G(i, j - 1)) * idy); - } - } -} - -void solve(Solver* solver) -{ - int imax = solver->imax; - int jmax = solver->jmax; - double eps = solver->eps; - int itermax = solver->itermax; - double dx2 = solver->dx * solver->dx; - double dy2 = solver->dy * solver->dy; - double idx2 = 1.0 / dx2; - double idy2 = 1.0 / dy2; - double factor = solver->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2); - double* p = solver->p; - double* rhs = solver->rhs; - double epssq = eps * eps; - int it = 0; - double res = 1.0; - - while ((res >= epssq) && (it < itermax)) { - res = 0.0; - - for (int j = 1; j < jmax + 1; j++) { - for (int i = 1; i < imax + 1; i++) { - - double r = RHS(i, j) - - ((P(i + 1, j) - 2.0 * P(i, j) + P(i - 1, j)) * idx2 + - (P(i, j + 1) - 2.0 * P(i, j) + P(i, j - 1)) * idy2); - - P(i, j) -= (factor * r); - res += (r * r); - } - } - - for (int i = 1; i < imax + 1; i++) { - P(i, 0) = P(i, 1); - P(i, jmax + 1) = P(i, jmax); - } - - for (int j = 1; j < jmax + 1; j++) { - P(0, j) = P(1, j); - P(imax + 1, j) = P(imax, j); - } - - res = res / (double)(imax * jmax); -#ifdef DEBUG - printf("%d Residuum: %e\n", it, res); -#endif - it++; - } - -#ifdef VERBOSE - printf("Solver took %d iterations to reach %f\n", it, sqrt(res)); -#endif -} - -void solveRB(Solver* solver) -{ - int imax = solver->imax; - int jmax = solver->jmax; - double eps = solver->eps; - int itermax = solver->itermax; - double dx2 = solver->dx * solver->dx; - double dy2 = solver->dy * solver->dy; - double idx2 = 1.0 / dx2; - double idy2 = 1.0 / dy2; - double factor = solver->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2); - double* p = solver->p; - double* rhs = solver->rhs; - double epssq = eps * eps; - int it = 0; - double res = 1.0; - int pass, jsw, isw; - - while ((res >= epssq) && (it < itermax)) { - res = 0.0; - jsw = 1; - - for (pass = 0; pass < 2; pass++) { - isw = jsw; - - for (int j = 1; j < jmax + 1; j++) { - for (int i = isw; i < imax + 1; i += 2) { - - double r = RHS(i, j) - - ((P(i + 1, j) - 2.0 * P(i, j) + P(i - 1, j)) * idx2 + - (P(i, j + 1) - 2.0 * P(i, j) + P(i, j - 1)) * idy2); - - P(i, j) -= (factor * r); - res += (r * r); - } - isw = 3 - isw; - } - jsw = 3 - jsw; - } - - for (int i = 1; i < imax + 1; i++) { - P(i, 0) = P(i, 1); - P(i, jmax + 1) = P(i, jmax); - } - - for (int j = 1; j < jmax + 1; j++) { - P(0, j) = P(1, j); - P(imax + 1, j) = P(imax, j); - } - - res = res / (double)(imax * jmax); -#ifdef DEBUG - printf("%d Residuum: %e\n", it, res); -#endif - it++; - } - -#ifdef VERBOSE - printf("Solver took %d iterations to reach %f\n", it, sqrt(res)); -#endif -} - -static double maxElement(Solver* solver, double* m) -{ - int size = (solver->imax + 2) * (solver->jmax + 2); - double maxval = DBL_MIN; - - for (int i = 0; i < size; i++) { - maxval = MAX(maxval, fabs(m[i])); - } - - return maxval; -} - -void normalizePressure(Solver* solver) -{ - int size = (solver->imax + 2) * (solver->jmax + 2); - double* p = solver->p; - double avgP = 0.0; - - for (int i = 0; i < size; i++) { - avgP += p[i]; - } - avgP /= size; - - for (int i = 0; i < size; i++) { - p[i] = p[i] - avgP; - } -} - -void computeTimestep(Solver* solver) -{ - double dt = solver->dtBound; - double dx = solver->dx; - double dy = solver->dy; - double umax = maxElement(solver, solver->u); - double vmax = maxElement(solver, solver->v); - - if (umax > 0) { - dt = (dt > dx / umax) ? dx / umax : dt; - } - if (vmax > 0) { - dt = (dt > dy / vmax) ? dy / vmax : dt; - } - - solver->dt = dt * solver->tau; -} - -void setBoundaryConditions(Solver* solver) -{ - int imax = solver->imax; - int jmax = solver->jmax; - double* u = solver->u; - double* v = solver->v; - - // Left boundary - switch (solver->bcLeft) { - case NOSLIP: - for (int j = 1; j < jmax + 1; j++) { - U(0, j) = 0.0; - V(0, j) = -V(1, j); - } - break; - case SLIP: - for (int j = 1; j < jmax + 1; j++) { - U(0, j) = 0.0; - V(0, j) = V(1, j); - } - break; - case OUTFLOW: - for (int j = 1; j < jmax + 1; j++) { - U(0, j) = U(1, j); - V(0, j) = V(1, j); - } - break; - case PERIODIC: - break; - } - - // Right boundary - switch (solver->bcRight) { - case NOSLIP: - for (int j = 1; j < jmax + 1; j++) { - U(imax, j) = 0.0; - V(imax + 1, j) = -V(imax, j); - } - break; - case SLIP: - for (int j = 1; j < jmax + 1; j++) { - U(imax, j) = 0.0; - V(imax + 1, j) = V(imax, j); - } - break; - case OUTFLOW: - for (int j = 1; j < jmax + 1; j++) { - U(imax, j) = U(imax - 1, j); - V(imax + 1, j) = V(imax, j); - } - break; - case PERIODIC: - break; - } - - // Bottom boundary - switch (solver->bcBottom) { - case NOSLIP: - for (int i = 1; i < imax + 1; i++) { - V(i, 0) = 0.0; - U(i, 0) = -U(i, 1); - } - break; - case SLIP: - for (int i = 1; i < imax + 1; i++) { - V(i, 0) = 0.0; - U(i, 0) = U(i, 1); - } - break; - case OUTFLOW: - for (int i = 1; i < imax + 1; i++) { - U(i, 0) = U(i, 1); - V(i, 0) = V(i, 1); - } - break; - case PERIODIC: - break; - } - - // Top boundary - switch (solver->bcTop) { - case NOSLIP: - for (int i = 1; i < imax + 1; i++) { - V(i, jmax) = 0.0; - U(i, jmax + 1) = -U(i, jmax); - } - break; - case SLIP: - for (int i = 1; i < imax + 1; i++) { - V(i, jmax) = 0.0; - U(i, jmax + 1) = U(i, jmax); - } - break; - case OUTFLOW: - for (int i = 1; i < imax + 1; i++) { - U(i, jmax + 1) = U(i, jmax); - V(i, jmax) = V(i, jmax - 1); - } - break; - case PERIODIC: - break; - } -} - -void setSpecialBoundaryCondition(Solver* solver) -{ - int imax = solver->imax; - int jmax = solver->jmax; - double mDy = solver->dy; - double* u = solver->u; - - if (strcmp(solver->problem, "dcavity") == 0) { - for (int i = 1; i < imax; i++) { - U(i, jmax + 1) = 2.0 - U(i, jmax); - } - } else if (strcmp(solver->problem, "canal") == 0) { - double ylength = solver->ylength; - double y; - - for (int j = 1; j < jmax + 1; j++) { - y = mDy * (j - 0.5); - U(0, j) = y * (ylength - y) * 4.0 / (ylength * ylength); - } - } -} - -void computeFG(Solver* solver) -{ - double* u = solver->u; - double* v = solver->v; - double* f = solver->f; - double* g = solver->g; - int imax = solver->imax; - int jmax = solver->jmax; - double gx = solver->gx; - double gy = solver->gy; - double gamma = solver->gamma; - double dt = solver->dt; - double inverseRe = 1.0 / solver->re; - double inverseDx = 1.0 / solver->dx; - double inverseDy = 1.0 / solver->dy; - double du2dx, dv2dy, duvdx, duvdy; - double du2dx2, du2dy2, dv2dx2, dv2dy2; - - for (int j = 1; j < jmax + 1; j++) { - for (int i = 1; i < imax + 1; i++) { - du2dx = inverseDx * 0.25 * - ((U(i, j) + U(i + 1, j)) * (U(i, j) + U(i + 1, j)) - - (U(i, j) + U(i - 1, j)) * (U(i, j) + U(i - 1, j))) + - gamma * inverseDx * 0.25 * - (fabs(U(i, j) + U(i + 1, j)) * (U(i, j) - U(i + 1, j)) + - fabs(U(i, j) + U(i - 1, j)) * (U(i, j) - U(i - 1, j))); - - duvdy = inverseDy * 0.25 * - ((V(i, j) + V(i + 1, j)) * (U(i, j) + U(i, j + 1)) - - (V(i, j - 1) + V(i + 1, j - 1)) * (U(i, j) + U(i, j - 1))) + - gamma * inverseDy * 0.25 * - (fabs(V(i, j) + V(i + 1, j)) * (U(i, j) - U(i, j + 1)) + - fabs(V(i, j - 1) + V(i + 1, j - 1)) * - (U(i, j) - U(i, j - 1))); - - du2dx2 = inverseDx * inverseDx * (U(i + 1, j) - 2.0 * U(i, j) + U(i - 1, j)); - du2dy2 = inverseDy * inverseDy * (U(i, j + 1) - 2.0 * U(i, j) + U(i, j - 1)); - F(i, j) = U(i, j) + dt * (inverseRe * (du2dx2 + du2dy2) - du2dx - duvdy + gx); - - duvdx = inverseDx * 0.25 * - ((U(i, j) + U(i, j + 1)) * (V(i, j) + V(i + 1, j)) - - (U(i - 1, j) + U(i - 1, j + 1)) * (V(i, j) + V(i - 1, j))) + - gamma * inverseDx * 0.25 * - (fabs(U(i, j) + U(i, j + 1)) * (V(i, j) - V(i + 1, j)) + - fabs(U(i - 1, j) + U(i - 1, j + 1)) * - (V(i, j) - V(i - 1, j))); - - dv2dy = inverseDy * 0.25 * - ((V(i, j) + V(i, j + 1)) * (V(i, j) + V(i, j + 1)) - - (V(i, j) + V(i, j - 1)) * (V(i, j) + V(i, j - 1))) + - gamma * inverseDy * 0.25 * - (fabs(V(i, j) + V(i, j + 1)) * (V(i, j) - V(i, j + 1)) + - fabs(V(i, j) + V(i, j - 1)) * (V(i, j) - V(i, j - 1))); - - dv2dx2 = inverseDx * inverseDx * (V(i + 1, j) - 2.0 * V(i, j) + V(i - 1, j)); - dv2dy2 = inverseDy * inverseDy * (V(i, j + 1) - 2.0 * V(i, j) + V(i, j - 1)); - G(i, j) = V(i, j) + dt * (inverseRe * (dv2dx2 + dv2dy2) - duvdx - dv2dy + gy); - } - } - - /* ---------------------- boundary of F --------------------------- */ - for (int j = 1; j < jmax + 1; j++) { - F(0, j) = U(0, j); - F(imax, j) = U(imax, j); - } - - /* ---------------------- boundary of G --------------------------- */ - for (int i = 1; i < imax + 1; i++) { - G(i, 0) = V(i, 0); - G(i, jmax) = V(i, jmax); - } -} - -void adaptUV(Solver* solver) -{ - int imax = solver->imax; - int jmax = solver->jmax; - double* p = solver->p; - double* u = solver->u; - double* v = solver->v; - double* f = solver->f; - double* g = solver->g; - double factorX = solver->dt / solver->dx; - double factorY = solver->dt / solver->dy; - - for (int j = 1; j < jmax + 1; j++) { - for (int i = 1; i < imax + 1; i++) { - U(i, j) = F(i, j) - (P(i + 1, j) - P(i, j)) * factorX; - V(i, j) = G(i, j) - (P(i, j + 1) - P(i, j)) * factorY; - } - } -} - -void writeResult(Solver* solver) -{ - int imax = solver->imax; - int jmax = solver->jmax; - double dx = solver->dx; - double dy = solver->dy; - double* p = solver->p; - double* u = solver->u; - double* v = solver->v; - double x = 0.0, y = 0.0; - - FILE* fp; - fp = fopen("pressure.dat", "w"); - - if (fp == NULL) { - printf("Error!\n"); - exit(EXIT_FAILURE); - } - - for (int j = 1; j < jmax + 1; j++) { - y = (double)(j - 0.5) * dy; - for (int i = 1; i < imax + 1; i++) { - x = (double)(i - 0.5) * dx; - fprintf(fp, "%.2f %.2f %f\n", x, y, P(i, j)); - } - fprintf(fp, "\n"); - } - - fclose(fp); - - fp = fopen("velocity.dat", "w"); - - if (fp == NULL) { - printf("Error!\n"); - exit(EXIT_FAILURE); - } - - for (int j = 1; j < jmax + 1; j++) { - y = dy * (j - 0.5); - for (int i = 1; i < imax + 1; i++) { - x = dx * (i - 0.5); - double vel_u = (U(i, j) + U(i - 1, j)) / 2.0; - double vel_v = (V(i, j) + V(i, j - 1)) / 2.0; - double len = sqrt((vel_u * vel_u) + (vel_v * vel_v)); - fprintf(fp, "%.2f %.2f %f %f %f\n", x, y, vel_u, vel_v, len); - } - } - - fclose(fp); -} diff --git a/BasicSolver/2D-seq/src/solver.h b/BasicSolver/2D-seq/src/solver.h index 4525799..d7ce513 100644 --- a/BasicSolver/2D-seq/src/solver.h +++ b/BasicSolver/2D-seq/src/solver.h @@ -6,41 +6,21 @@ */ #ifndef __SOLVER_H_ #define __SOLVER_H_ +#include "discretization.h" +#include "grid.h" #include "parameter.h" -enum BC { NOSLIP = 1, SLIP, OUTFLOW, PERIODIC }; - typedef struct { /* geometry and grid information */ - double dx, dy; - int imax, jmax; - double xlength, ylength; - /* arrays */ - double *p, *rhs; - double *f, *g; - double *u, *v; + Grid* grid; /* parameters */ - double eps, omega; - double re, tau, gamma; - double gx, gy; - /* time stepping */ + double eps, omega, rho; int itermax; - double dt, te; - double dtBound; - char* problem; - int bcLeft, bcRight, bcBottom, bcTop; + int levels; + double **r, **e; } Solver; -extern void initSolver(Solver*, Parameter*); -extern void computeRHS(Solver*); -extern void solve(Solver*); -extern void solveRB(Solver*); -extern void solveRBA(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 writeResult(Solver*); +extern void initSolver(Solver*, Discretization*, Parameter*); +extern void solve(Solver*, double*, double*); + #endif diff --git a/BasicSolver/2D-seq/src/util.h b/BasicSolver/2D-seq/src/util.h index e36948a..e4a93c2 100644 --- a/BasicSolver/2D-seq/src/util.h +++ b/BasicSolver/2D-seq/src/util.h @@ -19,11 +19,11 @@ #define ABS(a) ((a) >= 0 ? (a) : -(a)) #endif -#define P(i, j) p[(j) * (imaxLocal + 2) + (i)] -#define F(i, j) f[(j) * (imaxLocal + 2) + (i)] -#define G(i, j) g[(j) * (imaxLocal + 2) + (i)] -#define U(i, j) u[(j) * (imaxLocal + 2) + (i)] -#define V(i, j) v[(j) * (imaxLocal + 2) + (i)] -#define RHS(i, j) rhs[(j) * (imaxLocal + 2) + (i)] +#define P(i, j) p[(j) * (imax + 2) + (i)] +#define F(i, j) f[(j) * (imax + 2) + (i)] +#define G(i, j) g[(j) * (imax + 2) + (i)] +#define U(i, j) u[(j) * (imax + 2) + (i)] +#define V(i, j) v[(j) * (imax + 2) + (i)] +#define RHS(i, j) rhs[(j) * (imax + 2) + (i)] #endif // __UTIL_H_