From 6f3d9e73ef5fae2e50a2d86c693a0a3a76cac5f4 Mon Sep 17 00:00:00 2001 From: Jan Eitzinger Date: Tue, 5 Mar 2024 21:24:45 +0100 Subject: [PATCH] Review particle tracer. Refactor. --- EnhancedSolver/2D-seq/Makefile | 5 +- EnhancedSolver/2D-seq/canal.par | 12 +- EnhancedSolver/2D-seq/config.mk | 10 +- EnhancedSolver/2D-seq/dcavity.par | 6 +- EnhancedSolver/2D-seq/include_CLANG.mk | 12 +- EnhancedSolver/2D-seq/src/affinity.c | 61 -- EnhancedSolver/2D-seq/src/affinity.h | 14 - EnhancedSolver/2D-seq/src/allocate.c | 5 +- EnhancedSolver/2D-seq/src/allocate.h | 2 +- EnhancedSolver/2D-seq/src/discretization.c | 649 ++++++++++++++ EnhancedSolver/2D-seq/src/discretization.h | 44 + EnhancedSolver/2D-seq/src/grid.h | 22 + EnhancedSolver/2D-seq/src/main.c | 66 +- EnhancedSolver/2D-seq/src/particletracing.c | 202 +++-- EnhancedSolver/2D-seq/src/particletracing.h | 11 +- EnhancedSolver/2D-seq/src/progress.h | 2 +- EnhancedSolver/2D-seq/src/solver-mg.c | 186 ++++ EnhancedSolver/2D-seq/src/solver-sor.c | 128 +++ EnhancedSolver/2D-seq/src/solver.c | 900 -------------------- EnhancedSolver/2D-seq/src/solver.h | 55 +- EnhancedSolver/2D-seq/src/timing.c | 6 +- EnhancedSolver/2D-seq/src/timing.h | 5 +- EnhancedSolver/2D-seq/src/util.h | 10 +- EnhancedSolver/2D-seq/vis_files/canal.plot | 10 + 24 files changed, 1221 insertions(+), 1202 deletions(-) delete mode 100644 EnhancedSolver/2D-seq/src/affinity.c delete mode 100644 EnhancedSolver/2D-seq/src/affinity.h create mode 100644 EnhancedSolver/2D-seq/src/discretization.c create mode 100644 EnhancedSolver/2D-seq/src/discretization.h create mode 100644 EnhancedSolver/2D-seq/src/solver-mg.c create mode 100644 EnhancedSolver/2D-seq/src/solver-sor.c delete mode 100644 EnhancedSolver/2D-seq/src/solver.c create mode 100644 EnhancedSolver/2D-seq/vis_files/canal.plot diff --git a/EnhancedSolver/2D-seq/Makefile b/EnhancedSolver/2D-seq/Makefile index aab8407..dc8a79d 100644 --- a/EnhancedSolver/2D-seq/Makefile +++ b/EnhancedSolver/2D-seq/Makefile @@ -18,9 +18,10 @@ include $(MAKE_DIR)/include_$(TAG).mk INCLUDES += -I$(SRC_DIR) -I$(BUILD_DIR) VPATH = $(SRC_DIR) -SRC = $(wildcard $(SRC_DIR)/*.c) +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) @@ -49,6 +50,8 @@ distclean: clean @rm -f $(TARGET) @rm -f *.dat @rm -f *.png + @rm -f ./vis_files/*.dat + @rm -f ./vis_files/*.gif info: $(info $(CFLAGS)) diff --git a/EnhancedSolver/2D-seq/canal.par b/EnhancedSolver/2D-seq/canal.par index ca30e51..1661db6 100644 --- a/EnhancedSolver/2D-seq/canal.par +++ b/EnhancedSolver/2D-seq/canal.par @@ -26,13 +26,13 @@ p_init 1.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 +imax 256 # number of interior cells in x-direction +jmax 64 # number of interior cells in y-direction # Time Data: # --------- -te 60.0 # final time +te 80.0 # final time dt 0.02 # time stepsize tau 0.5 # safety factor for time stepsize control (<0 constant delt) @@ -50,9 +50,9 @@ levels 5 # Multigrid levels # ----------------------- numberOfParticles 60 -startTime 100 -injectTimePeriod 2.0 -writeTimePeriod 0.5 +startTime 5.0 +injectTimePeriod 4.0 +writeTimePeriod 1.0 x1 1.0 y1 0.0 diff --git a/EnhancedSolver/2D-seq/config.mk b/EnhancedSolver/2D-seq/config.mk index d228f66..ccdfa2c 100644 --- a/EnhancedSolver/2D-seq/config.mk +++ b/EnhancedSolver/2D-seq/config.mk @@ -1,12 +1,12 @@ # Supported: GCC, CLANG, ICC TAG ?= CLANG ENABLE_OPENMP ?= false +# Supported: sor, mg +SOLVER ?= mg +# 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/EnhancedSolver/2D-seq/dcavity.par b/EnhancedSolver/2D-seq/dcavity.par index 11dab71..8ff64b6 100644 --- a/EnhancedSolver/2D-seq/dcavity.par +++ b/EnhancedSolver/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 100 # number of interior cells in x-direction -jmax 100 # 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: # --------- @@ -50,7 +50,7 @@ levels 5 # Multigrid levels # ----------------------- numberOfParticles 200 -startTime 100 +startTime 2.0 injectTimePeriod 0.5 writeTimePeriod 0.2 diff --git a/EnhancedSolver/2D-seq/include_CLANG.mk b/EnhancedSolver/2D-seq/include_CLANG.mk index 466f441..a124053 100644 --- a/EnhancedSolver/2D-seq/include_CLANG.mk +++ b/EnhancedSolver/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/EnhancedSolver/2D-seq/src/affinity.c b/EnhancedSolver/2D-seq/src/affinity.c deleted file mode 100644 index ef8355a..0000000 --- a/EnhancedSolver/2D-seq/src/affinity.c +++ /dev/null @@ -1,61 +0,0 @@ -/* - * 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. - */ -#ifdef __linux__ -#ifdef _OPENMP -#include -#include -#include -#include -#include -#include -#include - -#define MAX_NUM_THREADS 128 -#define gettid() syscall(SYS_gettid) - -static int getProcessorID(cpu_set_t* cpu_set) -{ - int processorId; - - for (processorId = 0; processorId < MAX_NUM_THREADS; processorId++) { - if (CPU_ISSET(processorId, cpu_set)) { - break; - } - } - return processorId; -} - -int affinity_getProcessorId() -{ - cpu_set_t cpu_set; - CPU_ZERO(&cpu_set); - sched_getaffinity(gettid(), sizeof(cpu_set_t), &cpu_set); - - return getProcessorID(&cpu_set); -} - -void affinity_pinThread(int processorId) -{ - cpu_set_t cpuset; - pthread_t thread; - - thread = pthread_self(); - CPU_ZERO(&cpuset); - CPU_SET(processorId, &cpuset); - pthread_setaffinity_np(thread, sizeof(cpu_set_t), &cpuset); -} - -void affinity_pinProcess(int processorId) -{ - cpu_set_t cpuset; - - CPU_ZERO(&cpuset); - CPU_SET(processorId, &cpuset); - sched_setaffinity(0, sizeof(cpu_set_t), &cpuset); -} -#endif /*_OPENMP*/ -#endif /*__linux__*/ diff --git a/EnhancedSolver/2D-seq/src/affinity.h b/EnhancedSolver/2D-seq/src/affinity.h deleted file mode 100644 index 84bf733..0000000 --- a/EnhancedSolver/2D-seq/src/affinity.h +++ /dev/null @@ -1,14 +0,0 @@ -/* - * 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. - */ -#ifndef AFFINITY_H -#define AFFINITY_H - -extern int affinity_getProcessorId(); -extern void affinity_pinProcess(int); -extern void affinity_pinThread(int); - -#endif /*AFFINITY_H*/ diff --git a/EnhancedSolver/2D-seq/src/allocate.c b/EnhancedSolver/2D-seq/src/allocate.c index 855cf32..cf2efd6 100644 --- a/EnhancedSolver/2D-seq/src/allocate.c +++ b/EnhancedSolver/2D-seq/src/allocate.c @@ -5,10 +5,13 @@ * license that can be found in the LICENSE file. */ #include +#include #include #include -void* allocate(int alignment, size_t bytesize) +#include "allocate.h" + +void* allocate(size_t alignment, size_t bytesize) { int errorCode; void* ptr; diff --git a/EnhancedSolver/2D-seq/src/allocate.h b/EnhancedSolver/2D-seq/src/allocate.h index 1a5c8cb..77f4ba0 100644 --- a/EnhancedSolver/2D-seq/src/allocate.h +++ b/EnhancedSolver/2D-seq/src/allocate.h @@ -8,6 +8,6 @@ #define __ALLOCATE_H_ #include -extern void* allocate(int alignment, size_t bytesize); +extern void* allocate(size_t alignment, size_t bytesize); #endif diff --git a/EnhancedSolver/2D-seq/src/discretization.c b/EnhancedSolver/2D-seq/src/discretization.c new file mode 100644 index 0000000..daa346d --- /dev/null +++ b/EnhancedSolver/2D-seq/src/discretization.c @@ -0,0 +1,649 @@ +/* + * 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 "grid.h" +#include "parameter.h" +#include "util.h" + +#define S(i, j) s[(j) * (imax + 2) + (i)] + +static double distance(double i, double j, double iCenter, double jCenter) +{ + return sqrt(pow(iCenter - i, 2) + pow(jCenter - j, 2) * 1.0); +} + +void print(Discretization* d, double* grid) +{ + int imax = d->grid.imax; + int jmax = d->grid.jmax; + + for (int j = 0; j < jmax + 2; j++) { + printf("%02d: ", j); + for (int i = 0; i < imax + 2; i++) { + printf("%3.2f ", grid[j * (imax + 2) + i]); + } + printf("\n"); + } + fflush(stdout); +} + +void printGrid(Discretization* d, int* grid) +{ + int imax = d->grid.imax; + int jmax = d->grid.jmax; + + for (int j = 0; j < jmax + 2; j++) { + printf("%02d: ", j); + for (int i = 0; i < imax + 2; i++) { + printf("%2d ", 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->grid.s = 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; + d->grid.s[i] = FLUID; + } + + 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; + + double xCenter = 0, yCenter = 0, radius = 0; + double x1 = 0, x2 = 0, y1 = 0, y2 = 0; + + int* s = d->grid.s; + + switch (p->shape) { + case NOSHAPE: + break; + case RECT: + x1 = p->xCenter - p->xRectLength / 2; + x2 = p->xCenter + p->xRectLength / 2; + y1 = p->yCenter - p->yRectLength / 2; + y2 = p->yCenter + p->yRectLength / 2; + + for (int j = 1; j < jmax + 1; j++) { + for (int i = 1; i < imax + 1; i++) { + if ((x1 <= (i * dx)) && ((i * dx) <= x2) && (y1 <= (j * dy)) && + ((j * dy) <= y2)) { + S(i, j) = OBSTACLE; + } + } + } + break; + case CIRCLE: + xCenter = p->xCenter; + yCenter = p->yCenter; + radius = p->circleRadius; + + for (int j = 1; j < jmax + 1; j++) { + for (int i = 1; i < imax + 1; i++) { + if (distance((i * dx), (j * dy), xCenter, yCenter) <= radius) { + S(i, j) = OBSTACLE; + } + } + } + break; + } + + if (p->shape != NOSHAPE) { + for (int j = 1; j < jmax + 1; j++) { + for (int i = 1; i < imax + 1; i++) { + + if (S(i, j - 1) == FLUID && S(i, j + 1) == OBSTACLE && + S(i, j) == OBSTACLE) + S(i, j) = BOTTOM; // TOP + if (S(i - 1, j) == FLUID && S(i + 1, j) == OBSTACLE && + S(i, j) == OBSTACLE) + S(i, j) = LEFT; // LEFT + if (S(i + 1, j) == FLUID && S(i - 1, j) == OBSTACLE && + S(i, j) == OBSTACLE) + S(i, j) = RIGHT; // RIGHT + if (S(i, j + 1) == FLUID && S(i, j - 1) == OBSTACLE && + S(i, j) == OBSTACLE) + S(i, j) = TOP; // BOTTOM + if (S(i - 1, j - 1) == FLUID && S(i, j - 1) == FLUID && + S(i - 1, j) == FLUID && S(i + 1, j + 1) == OBSTACLE && + (S(i, j) == OBSTACLE || S(i, j) == LEFT || S(i, j) == BOTTOM)) + S(i, j) = BOTTOMLEFT; // TOPLEFT + if (S(i + 1, j - 1) == FLUID && S(i, j - 1) == FLUID && + S(i + 1, j) == FLUID && S(i - 1, j + 1) == OBSTACLE && + (S(i, j) == OBSTACLE || S(i, j) == RIGHT || S(i, j) == BOTTOM)) + S(i, j) = BOTTOMRIGHT; // TOPRIGHT + if (S(i - 1, j + 1) == FLUID && S(i - 1, j) == FLUID && + S(i, j + 1) == FLUID && S(i + 1, j - 1) == OBSTACLE && + (S(i, j) == OBSTACLE || S(i, j) == LEFT || S(i, j) == TOP)) + S(i, j) = TOPLEFT; // BOTTOMLEFT + if (S(i + 1, j + 1) == FLUID && S(i + 1, j) == FLUID && + S(i, j + 1) == FLUID && S(i - 1, j - 1) == OBSTACLE && + (S(i, j) == OBSTACLE || S(i, j) == RIGHT || S(i, j) == TOP)) + S(i, j) = TOPRIGHT; // BOTTOMRIGHT + } + } + } + +#ifdef VERBOSE + printConfig(solver); +#endif +} + +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 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; + int* s = d->grid.s; + + 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 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; + int* s = d->grid.s; + + 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); + } + } else if (strcmp(d->problem, "backstep") == 0) { + for (int j = 1; j < jmax + 1; j++) { + if (S(0, j) == FLUID) U(0, j) = 1.0; + } + } else if (strcmp(d->problem, "karman") == 0) { + for (int j = 1; j < jmax + 1; j++) { + U(0, j) = 1.0; + } + } +} + +void setObjectBoundaryCondition(Discretization* d) +{ + int imax = d->grid.imax; + int jmax = d->grid.jmax; + double* u = d->u; + double* v = d->v; + int* s = d->grid.s; + + for (int j = 1; j < jmax + 1; j++) { + for (int i = 1; i < imax + 1; i++) { + switch (S(i, j)) { + case TOP: + U(i, j) = -U(i, j + 1); + U(i - 1, j) = -U(i - 1, j + 1); + V(i, j) = 0.0; + break; + case BOTTOM: + U(i, j) = -U(i, j - 1); + U(i - 1, j) = -U(i - 1, j - 1); + V(i, j) = 0.0; + break; + case LEFT: + U(i - 1, j) = 0.0; + V(i, j) = -V(i - 1, j); + V(i, j - 1) = -V(i - 1, j - 1); + break; + case RIGHT: + U(i, j) = 0.0; + V(i, j) = -V(i + 1, j); + V(i, j - 1) = -V(i + 1, j - 1); + break; + case TOPLEFT: + U(i, j) = -U(i, j + 1); + U(i - 1, j) = 0.0; + V(i, j) = 0.0; + V(i, j - 1) = -V(i - 1, j - 1); + break; + case TOPRIGHT: + U(i, j) = 0.0; + U(i - 1, j) = -U(i - 1, j + 1); + V(i, j) = 0.0; + V(i, j - 1) = -V(i + 1, j - 1); + break; + case BOTTOMLEFT: + U(i, j) = -U(i, j - 1); + U(i - 1, j) = 0.0; + V(i, j) = -V(i - 1, j); + V(i, j - 1) = 0.0; + break; + case BOTTOMRIGHT: + U(i, j) = 0.0; + U(i - 1, j) = -U(i - 1, j - 1); + V(i, j) = -V(i, j + 1); + V(i, j - 1) = 0.0; + break; + } + } + } +} + +void computeFG(Discretization* d) +{ + double* u = d->u; + double* v = d->v; + double* f = d->f; + double* g = d->g; + int* s = d->grid.s; + 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++) { + if (S(i, j) == FLUID) { + 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); + } else { + switch (S(i, j)) { + case TOP: + G(i, j) = V(i, j); + break; + case BOTTOM: + G(i, j - 1) = V(i, j - 1); + break; + case LEFT: + F(i - 1, j) = U(i - 1, j); + break; + case RIGHT: + F(i, j) = U(i, j); + break; + case TOPLEFT: + F(i - 1, j) = U(i - 1, j); + G(i, j) = V(i, j); + break; + case TOPRIGHT: + F(i, j) = U(i, j); + G(i, j) = V(i, j); + break; + case BOTTOMLEFT: + F(i - 1, j) = U(i - 1, j); + G(i, j - 1) = V(i, j - 1); + break; + case BOTTOMRIGHT: + F(i, j) = U(i, j); + G(i, j - 1) = V(i, j - 1); + break; + } + } + } + } + + /* ---------------------- 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/EnhancedSolver/2D-seq/src/discretization.h b/EnhancedSolver/2D-seq/src/discretization.h new file mode 100644 index 0000000..1b8d92d --- /dev/null +++ b/EnhancedSolver/2D-seq/src/discretization.h @@ -0,0 +1,44 @@ +/* + * Copyright (C) NHR@FAU, University Erlangen-Nuremberg. + * All rights reserved. This file is part of nusif-solver. + * Use of this source code is governed by a MIT style + * license that can be found in the LICENSE file. + */ +#ifndef __DISCRETIZATION_H_ +#define __DISCRETIZATION_H_ +#include "grid.h" +#include "parameter.h" + +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 rho; + 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 setObjectBoundaryCondition(Discretization*); +extern void computeFG(Discretization*); +extern void adaptUV(Discretization*); +extern void writeResult(Discretization*); +extern void print(Discretization*, double*); +extern void printGrid(Discretization*, int*); +#endif diff --git a/EnhancedSolver/2D-seq/src/grid.h b/EnhancedSolver/2D-seq/src/grid.h index 7e58d1a..218ffa9 100644 --- a/EnhancedSolver/2D-seq/src/grid.h +++ b/EnhancedSolver/2D-seq/src/grid.h @@ -7,10 +7,32 @@ #ifndef __GRID_H_ #define __GRID_H_ +#define S(i, j) s[(j) * (imax + 2) + (i)] + +enum OBJECTBOUNDARY { + FLUID = 0, + TOP, + BOTTOM, + LEFT, + RIGHT, + TOPLEFT, + BOTTOMLEFT, + TOPRIGHT, + BOTTOMRIGHT, + OBSTACLE +}; + +enum SHAPE { NOSHAPE = 0, RECT, CIRCLE }; + typedef struct { double dx, dy; int imax, jmax; double xlength, ylength; + int* s; } Grid; +static inline int gridIsFluid(Grid* g, int i, int j) +{ + return g->s[j * (g->imax + 2) + i]; +} #endif // __GRID_H_ diff --git a/EnhancedSolver/2D-seq/src/main.c b/EnhancedSolver/2D-seq/src/main.c index d889c3b..079d5e3 100644 --- a/EnhancedSolver/2D-seq/src/main.c +++ b/EnhancedSolver/2D-seq/src/main.c @@ -8,6 +8,7 @@ #include #include +#include "discretization.h" #include "parameter.h" #include "particletracing.h" #include "progress.h" @@ -16,52 +17,51 @@ int main(int argc, char** argv) { - double S, E; - Parameter params; - Solver solver; + double timeStart, timeStop; + Parameter p; + Discretization d; + Solver s; ParticleTracer particletracer; - initParameter(¶ms); + + initParameter(&p); if (argc != 2) { printf("Usage: %s \n", argv[0]); exit(EXIT_SUCCESS); } - readParameter(¶ms, argv[1]); - printParameter(¶ms); - initSolver(&solver, ¶ms); - printf("initsolver done\n"); - - initParticleTracer(&particletracer, &solver.grid, ¶ms); + readParameter(&p, argv[1]); + printParameter(&p); + initDiscretization(&d, &p); + initSolver(&s, &d, &p); + initParticleTracer(&particletracer, &d.grid, &p); printParticleTracerParameters(&particletracer); #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); - setObjectBoundaryCondition(&solver); - computeFG(&solver); - computeRHS(&solver); - if (nt % 100 == 0) normalizePressure(&solver); - multiGrid(&solver); + if (tau > 0.0) computeTimestep(&d); + setBoundaryConditions(&d); + setSpecialBoundaryCondition(&d); + setObjectBoundaryCondition(&d); - adaptUV(&solver); + computeFG(&d); + computeRHS(&d); + if (nt % 100 == 0) normalizePressure(&d); + solve(&s, d.p, d.rhs); + adaptUV(&d); + trace(&particletracer, d.u, d.v, d.dt, t); - /* Added function for particle tracing. Will inject and advance particles as per - * timePeriod */ - trace(&particletracer, solver.u, solver.v, solver.s, t); - - t += solver.dt; + t += d.dt; nt++; #ifdef VERBOSE @@ -70,14 +70,12 @@ int main(int argc, char** argv) printProgress(t); #endif } - printf("Total particles : %d\n", particletracer.totalParticles); - E = getTimeStamp(); + timeStop = getTimeStamp(); + stopProgress(); - freeParticles(&particletracer); - - printf("Solution took %.2fs\n", E - S); - writeResult(&solver); + printf("Solution took %.2fs\n", timeStop - timeStart); + writeResult(&d); return EXIT_SUCCESS; } diff --git a/EnhancedSolver/2D-seq/src/particletracing.c b/EnhancedSolver/2D-seq/src/particletracing.c index b95a1a7..5750b0c 100644 --- a/EnhancedSolver/2D-seq/src/particletracing.c +++ b/EnhancedSolver/2D-seq/src/particletracing.c @@ -8,14 +8,13 @@ #include #include +#include "grid.h" #include "vtkWriter.h" #define U(i, j) u[(j) * (imax + 2) + (i)] #define V(i, j) v[(j) * (imax + 2) + (i)] #define S(i, j) s[(j) * (imax + 2) + (i)] -static int ts = 0; - void printParticles(ParticleTracer* p) { for (int i = 0; i < p->totalParticles; ++i) { @@ -25,29 +24,28 @@ void printParticles(ParticleTracer* p) p->particlePool[i].flag); } } -void injectParticles(ParticleTracer* p) + +static void injectParticles(ParticleTracer* p) { - for (int i = 0; i < p->numberOfParticles; ++i) { + if (p->totalParticles + p->numParticlesInLine > p->numAllocatedParticles) { + return; + } + for (int i = 0; i < p->numParticlesInLine; ++i) { p->particlePool[p->pointer].x = p->linSpaceLine[i].x; p->particlePool[p->pointer].y = p->linSpaceLine[i].y; p->particlePool[p->pointer].flag = true; - ++(p->pointer); - ++(p->totalParticles); + p->pointer++; + p->totalParticles++; } } -void advanceParticles(ParticleTracer* p, - double* restrict u, - double* restrict v, - int* restrict s, - double time) +static void advanceParticles( + ParticleTracer* p, double* restrict u, double* restrict v, double dt) { - int imax = p->grid->imax; - int jmax = p->grid->jmax; - - double dx = p->grid->dx; - double dy = p->grid->dy; - + int imax = p->grid->imax; + int jmax = p->grid->jmax; + double dx = p->grid->dx; + double dy = p->grid->dy; double xlength = p->grid->xlength; double ylength = p->grid->ylength; @@ -64,14 +62,14 @@ void advanceParticles(ParticleTracer* p, double x2 = (double)iCoord * dx; double y2 = ((double)jCoord - 0.5) * dy; - double u_n = (1.0 / (dx * dy)) * - ((x2 - x) * (y2 - y) * U(iCoord - 1, jCoord - 1) + - (x - x1) * (y2 - y) * U(iCoord, jCoord - 1) + - (x2 - x) * (y - y1) * U(iCoord - 1, jCoord) + - (x - x1) * (y - y1) * U(iCoord, jCoord)); + double intU = (1.0 / (dx * dy)) * + ((x2 - x) * (y2 - y) * U(iCoord - 1, jCoord - 1) + + (x - x1) * (y2 - y) * U(iCoord, jCoord - 1) + + (x2 - x) * (y - y1) * U(iCoord - 1, jCoord) + + (x - x1) * (y - y1) * U(iCoord, jCoord)); - double new_x = x + p->dt * u_n; - p->particlePool[i].x = new_x; + double newX = x + dt * intU; + p->particlePool[i].x = newX; iCoord = (int)((x + 0.5 * dx) / dx) + 1; jCoord = (int)(y / dy) + 1; @@ -81,50 +79,55 @@ void advanceParticles(ParticleTracer* p, x2 = ((double)iCoord - 0.5) * dx; y2 = (double)jCoord * dy; - double v_n = (1.0 / (dx * dy)) * - ((x2 - x) * (y2 - y) * V(iCoord - 1, jCoord - 1) + - (x - x1) * (y2 - y) * V(iCoord, jCoord - 1) + - (x2 - x) * (y - y1) * V(iCoord - 1, jCoord) + - (x - x1) * (y - y1) * V(iCoord, jCoord)); + double intV = (1.0 / (dx * dy)) * + ((x2 - x) * (y2 - y) * V(iCoord - 1, jCoord - 1) + + (x - x1) * (y2 - y) * V(iCoord, jCoord - 1) + + (x2 - x) * (y - y1) * V(iCoord - 1, jCoord) + + (x - x1) * (y - y1) * V(iCoord, jCoord)); - double new_y = y + p->dt * v_n; - p->particlePool[i].y = new_y; + double newY = y + dt * intV; + p->particlePool[i].y = newY; - // printf("\tOld X : %.2f, New X : %.2f, iCoord : %d\n\tOld Y : %.2f, New Y : - // %.2f, jCoord : %d\n\n", x, new_x, iCoord, y, new_y, jCoord); - // printf("\tU(iCoord - 1, jCoord - 1) : %.2f, U(iCoord, jCoord - 1) : %.2f, - // U(iCoord - 1, jCoord) : %.2f, U(iCoord, jCoord) : %.2f\n", U(iCoord - 1, - // jCoord - 1), U(iCoord, jCoord - 1), U(iCoord - 1, jCoord), U(iCoord, - // jCoord)); printf("\tV(iCoord - 1, jCoord - 1) : %.2f, V(iCoord, jCoord - 1) - // : %.2f, V(iCoord - 1, jCoord) : %.2f, V(iCoord, jCoord) : %.2f\n\n", - // V(iCoord - 1, jCoord - 1), V(iCoord, jCoord - 1), V(iCoord - 1, jCoord), - // V(iCoord, jCoord)); printf("\t U N : %.2f, V N : %.2f\n\n", u_n, v_n); - // printf("\t j-1 * (imax + 2) + i-1 = %d with element from U : %.2f", (jCoord - // - 1) * (200 + 2) + (iCoord - 1), u[(jCoord - 1) * (imax + 2) + (iCoord - - // 1)]); printf("\nimax : %d, jmax : %d\n", imax, jmax); - - if (((new_x < 0.0) || (new_x > xlength) || (new_y < 0.0) || - (new_y > ylength))) { + if (((newX < 0.0) || (newX > xlength) || (newY < 0.0) || (newY > ylength))) { p->particlePool[i].flag = false; + p->removedParticles++; } - int i_new = new_x / dx, j_new = new_y / dy; - if (S(i_new, j_new) != NONE) { + + int newI = newX / dx, newJ = newY / dy; + + if (!gridIsFluid(p->grid, newI, newJ)) { p->particlePool[i].flag = false; + p->removedParticles++; + printf("Forbidden movement of particle into obstacle!\n"); } } } } -void freeParticles(ParticleTracer* p) +static void compress(ParticleTracer* p) { - if (p->particlePool != NULL) { - free(p->particlePool); - free(p->linSpaceLine); + Particle* memPool = p->particlePool; + Particle tempPool[p->totalParticles]; + int totalParticles = 0; + + for (int i = 0; i < p->totalParticles; i++) { + if (memPool[i].flag == 1) { + tempPool[totalParticles].x = memPool[i].x; + tempPool[totalParticles].y = memPool[i].y; + tempPool[totalParticles].flag = memPool[i].flag; + totalParticles++; + } } + + p->totalParticles = totalParticles; + p->removedParticles = 0; + p->pointer = totalParticles + 1; + memcpy(p->particlePool, tempPool, totalParticles * sizeof(Particle)); } void writeParticles(ParticleTracer* p) { + static int ts = 0; VtkOptions opts = { .particletracer = p }; char filename[50]; @@ -154,44 +157,43 @@ void writeParticles(ParticleTracer* p) ++ts; } -void initParticleTracer(ParticleTracer* p, Grid* grid, Parameter* params) +void initParticleTracer(ParticleTracer* pt, Grid* g, Parameter* p) { - p->numberOfParticles = params->numberOfParticles; - p->startTime = params->startTime; - p->injectTimePeriod = params->injectTimePeriod; - p->writeTimePeriod = params->writeTimePeriod; + pt->numParticlesInLine = p->numberOfParticles; + pt->startTime = p->startTime; + pt->injectTimePeriod = p->injectTimePeriod; + pt->writeTimePeriod = p->writeTimePeriod; + pt->grid = g; - p->dt = params->dt; - p->grid = grid; + pt->x1 = p->x1; + pt->y1 = p->y1; + pt->x2 = p->x2; + pt->y2 = p->y2; - p->x1 = params->x1; - p->y1 = params->y1; - p->x2 = params->x2; - p->y2 = params->y2; + pt->lastInjectTime = p->startTime; + pt->lastWriteTime = p->startTime; - p->lastInjectTime = params->startTime; - p->lastUpdateTime = params->startTime; - p->lastWriteTime = params->startTime; + pt->pointer = 0; + pt->removedParticles = 0; + pt->totalParticles = 0; - p->pointer = 0; - p->totalParticles = 0; + if (p->te > p->startTime) { + pt->numAllocatedParticles = ((p->te - p->startTime) / p->injectTimePeriod) * + p->numberOfParticles; + pt->numAllocatedParticles += (2 * p->numberOfParticles); - if (params->te > params->startTime) { - p->estimatedNumParticles = ((params->te - params->startTime) + 2) * - params->numberOfParticles; + pt->particlePool = malloc(sizeof(Particle) * pt->numAllocatedParticles); + pt->linSpaceLine = malloc(sizeof(Particle) * pt->numParticlesInLine); - p->particlePool = malloc(sizeof(Particle) * p->estimatedNumParticles); - p->linSpaceLine = malloc(sizeof(Particle) * p->numberOfParticles); - - for (int i = 0; i < p->numberOfParticles; ++i) { - double spacing = (double)i / (double)(p->numberOfParticles - 1); - p->linSpaceLine[i].x = spacing * p->x1 + (1.0 - spacing) * p->x2; - p->linSpaceLine[i].y = spacing * p->y1 + (1.0 - spacing) * p->y2; - p->linSpaceLine[i].flag = true; + for (int i = 0; i < pt->numParticlesInLine; ++i) { + double spacing = (double)i / (double)(pt->numParticlesInLine - 1); + pt->linSpaceLine[i].x = spacing * pt->x1 + (1.0 - spacing) * pt->x2; + pt->linSpaceLine[i].y = spacing * pt->y1 + (1.0 - spacing) * pt->y2; + pt->linSpaceLine[i].flag = true; } } else { - p->particlePool = NULL; - p->linSpaceLine = NULL; + pt->particlePool = NULL; + pt->linSpaceLine = NULL; } } @@ -199,7 +201,7 @@ void printParticleTracerParameters(ParticleTracer* p) { printf("Particle Tracing data:\n"); printf("\tNumber of particles : %d being injected for every period of %.2f\n", - p->numberOfParticles, + p->numParticlesInLine, p->injectTimePeriod); printf("\tstartTime : %.2f\n", p->startTime); printf("\t(Line along which the particles are to be injected) \n\tx1 : %.2f, y1 : " @@ -209,43 +211,33 @@ void printParticleTracerParameters(ParticleTracer* p) p->x2, p->y2); printf("\tPointer : %d, TotalParticles : %d\n", p->pointer, p->totalParticles); - printf("\tdt : %.2f, dx : %.2f, dy : %.2f\n", p->dt, p->grid->dx, p->grid->dy); } -void trace(ParticleTracer* p, double* u, double* v, int* s, double time) +void trace(ParticleTracer* p, double* u, double* v, double dt, double time) { if (time >= p->startTime) { - // printParticles(particletracer); if ((time - p->lastInjectTime) >= p->injectTimePeriod) { injectParticles(p); p->lastInjectTime = time; } + if ((time - p->lastWriteTime) >= p->writeTimePeriod) { writeParticles(p); p->lastWriteTime = time; } - advanceParticles(p, u, v, s, time); - compress(p); - p->lastUpdateTime = time; - } -} -void compress(ParticleTracer* p) -{ - Particle* memPool = p->particlePool; - Particle tempPool[p->totalParticles]; - int totalParticles = 0; + advanceParticles(p, u, v, dt); - for (int i = 0; i < p->totalParticles; ++i) { - if (memPool[i].flag == 1) { - tempPool[totalParticles].x = memPool[i].x; - tempPool[totalParticles].y = memPool[i].y; - tempPool[totalParticles].flag = memPool[i].flag; - ++totalParticles; + if (p->removedParticles > (p->totalParticles * 0.2)) { + compress(p); } } - - p->totalParticles = totalParticles; - p->pointer = totalParticles + 1; - memcpy(p->particlePool, tempPool, totalParticles * sizeof(Particle)); +} + +void freeParticles(ParticleTracer* p) +{ + if (p->particlePool != NULL) { + free(p->particlePool); + free(p->linSpaceLine); + } } diff --git a/EnhancedSolver/2D-seq/src/particletracing.h b/EnhancedSolver/2D-seq/src/particletracing.h index 40621d4..09855a5 100644 --- a/EnhancedSolver/2D-seq/src/particletracing.h +++ b/EnhancedSolver/2D-seq/src/particletracing.h @@ -19,13 +19,12 @@ typedef struct { } Particle; typedef struct { - int numberOfParticles, totalParticles; + int numParticlesInLine, removedParticles, totalParticles; double startTime, injectTimePeriod, writeTimePeriod; double lastInjectTime, lastUpdateTime, lastWriteTime; - int estimatedNumParticles; + int numAllocatedParticles; - double dt; Particle* linSpaceLine; Particle* particlePool; @@ -35,12 +34,8 @@ typedef struct { } ParticleTracer; extern void initParticleTracer(ParticleTracer*, Grid*, Parameter*); -extern void injectParticles(ParticleTracer*); -extern void advanceParticles(ParticleTracer*, double*, double*, int*, double); extern void freeParticles(ParticleTracer*); extern void writeParticles(ParticleTracer*); extern void printParticleTracerParameters(ParticleTracer*); -extern void printParticles(ParticleTracer*); -extern void trace(ParticleTracer*, double*, double*, int*, double); -extern void compress(ParticleTracer*); +extern void trace(ParticleTracer*, double*, double*, double, double); #endif diff --git a/EnhancedSolver/2D-seq/src/progress.h b/EnhancedSolver/2D-seq/src/progress.h index 1cdb9a3..314e921 100644 --- a/EnhancedSolver/2D-seq/src/progress.h +++ b/EnhancedSolver/2D-seq/src/progress.h @@ -9,6 +9,6 @@ extern void initProgress(double); extern void printProgress(double); -extern void stopProgress(); +extern void stopProgress(void); #endif diff --git a/EnhancedSolver/2D-seq/src/solver-mg.c b/EnhancedSolver/2D-seq/src/solver-mg.c new file mode 100644 index 0000000..d1e56b4 --- /dev/null +++ b/EnhancedSolver/2D-seq/src/solver-mg.c @@ -0,0 +1,186 @@ +/* + * 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) * (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; +} + +static double multiGrid(Solver* s, double* p, double* rhs, int level, int imax, int jmax) +{ + double res = 0.0; + + // coarsest level TODO: Use direct solver? + if (level == COARSEST_LEVEL) { + for (int i = 0; i < 5; i++) { + smooth(s, p, rhs, level, imax, jmax); + } + return res; + } + + // pre-smoothing TODO: Make smoothing steps configurable? + for (int i = 0; i < 5; i++) { + smooth(s, p, rhs, level, imax, jmax); + if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax); + } + + // restrict + restrictMG(s, level, imax, jmax); + + // MGSolver on residual and error. + // TODO: What if there is a rest? + multiGrid(s, s->e[level + 1], s->r[level + 1], level + 1, imax / 2, jmax / 2); + + // prolongate + prolongate(s, level, imax, jmax); + + // correct p on finer level using residual + correct(s, p, level, imax, jmax); + if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax); + + // post-smoothing + for (int i = 0; i < 5; i++) { + res = smooth(s, p, rhs, level, imax, jmax); + if (level == FINEST_LEVEL) setBoundaryCondition(p, 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->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; + } + } +} + +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/EnhancedSolver/2D-seq/src/solver-sor.c b/EnhancedSolver/2D-seq/src/solver-sor.c new file mode 100644 index 0000000..f817b5f --- /dev/null +++ b/EnhancedSolver/2D-seq/src/solver-sor.c @@ -0,0 +1,128 @@ +/* + * 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 "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/EnhancedSolver/2D-seq/src/solver.c b/EnhancedSolver/2D-seq/src/solver.c deleted file mode 100644 index bbb51f5..0000000 --- a/EnhancedSolver/2D-seq/src/solver.c +++ /dev/null @@ -1,900 +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 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 oldR(i, j) oldr[(j) * (imax + 2) + (i)] -#define oldE(i, j) olde[(j) * (imax + 2) + (i)] -#define RHS(i, j) rhs[(j) * (imax + 2) + (i)] - -static double distance(double i, double j, double iCenter, double jCenter) -{ - return sqrt(pow(iCenter - i, 2) + pow(jCenter - j, 2) * 1.0); -} - -void print(Solver* solver, double* grid) -{ - int imax = solver->grid.imax; - int jmax = solver->grid.jmax; - - for (int j = 0; j < jmax + 2; j++) { - printf("%02d: ", j); - for (int i = 0; i < imax + 2; i++) { - printf("%3.2f ", grid[j * (imax + 2) + i]); - } - printf("\n"); - } - fflush(stdout); -} - -void printGrid(Solver* solver, int* grid) -{ - int imax = solver->grid.imax; - int jmax = solver->grid.jmax; - - for (int j = 0; j < jmax + 2; j++) { - printf("%02d: ", j); - for (int i = 0; i < imax + 2; i++) { - printf("%2d ", 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->grid.xlength, - solver->grid.ylength); - printf("\tCells (x, y): %d, %d\n", solver->grid.imax, solver->grid.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->grid.imax = params->imax; - solver->grid.jmax = params->jmax; - solver->grid.xlength = params->xlength; - solver->grid.ylength = params->ylength; - solver->grid.dx = params->xlength / params->imax; - solver->grid.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; - solver->rho = params->rho; - solver->levels = params->levels; - solver->currentlevel = 0; - - int imax = solver->grid.imax; - int jmax = solver->grid.jmax; - int levels = solver->levels; - - size_t size_level = levels * (imax + 2) * (jmax + 2) * sizeof(double); - - size_t size = (imax + 2) * (jmax + 2) * sizeof(double); - solver->u = allocate(64, size); - solver->v = allocate(64, size); - solver->s = allocate(64, size); - solver->p = allocate(64, size); - solver->rhs = allocate(64, size); - solver->f = allocate(64, size); - solver->g = allocate(64, size); - - solver->r = malloc(levels * sizeof(double*)); - solver->e = malloc(levels * sizeof(double*)); - - for (int j = 0; j < levels; ++j) { - - solver->r[j] = allocate(64, size); - solver->e[j] = 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; - solver->s[i] = NONE; - for (int j = 0; j < levels; ++j) { - - solver->r[j][i] = 0.0; - solver->e[j][i] = 0.0; - } - } - - double dx = solver->grid.dx; - double dy = solver->grid.dy; - double invSqrSum = 1.0 / (dx * dx) + 1.0 / (dy * dy); - solver->dtBound = 0.5 * solver->re * 1.0 / invSqrSum; - - double xCenter = 0, yCenter = 0, radius = 0; - double x1 = 0, x2 = 0, y1 = 0, y2 = 0; - - int* s = solver->s; - - switch (params->shape) { - case NOSHAPE: - break; - case RECT: - x1 = params->xCenter - params->xRectLength / 2; - x2 = params->xCenter + params->xRectLength / 2; - y1 = params->yCenter - params->yRectLength / 2; - y2 = params->yCenter + params->yRectLength / 2; - - for (int j = 1; j < jmax + 1; j++) { - for (int i = 1; i < imax + 1; i++) { - if ((x1 <= (i * dx)) && ((i * dx) <= x2) && (y1 <= (j * dy)) && - ((j * dy) <= y2)) { - S(i, j) = LOCAL; - } - } - } - - break; - case CIRCLE: - xCenter = params->xCenter; - yCenter = params->yCenter; - radius = params->circleRadius; - - for (int j = 1; j < jmax + 1; j++) { - for (int i = 1; i < imax + 1; i++) { - if (distance((i * dx), (j * dy), xCenter, yCenter) <= radius) { - S(i, j) = LOCAL; - } - } - } - - break; - default: - break; - } - - if (params->shape != NOSHAPE) { - for (int j = 1; j < jmax + 1; j++) { - for (int i = 1; i < imax + 1; i++) { - - if (S(i, j - 1) == NONE && S(i, j + 1) == LOCAL && S(i, j) == LOCAL) - S(i, j) = BOTTOM; // TOP - if (S(i - 1, j) == NONE && S(i + 1, j) == LOCAL && S(i, j) == LOCAL) - S(i, j) = LEFT; // LEFT - if (S(i + 1, j) == NONE && S(i - 1, j) == LOCAL && S(i, j) == LOCAL) - S(i, j) = RIGHT; // RIGHT - if (S(i, j + 1) == NONE && S(i, j - 1) == LOCAL && S(i, j) == LOCAL) - S(i, j) = TOP; // BOTTOM - if (S(i - 1, j - 1) == NONE && S(i, j - 1) == NONE && - S(i - 1, j) == NONE && S(i + 1, j + 1) == LOCAL && - (S(i, j) == LOCAL || S(i, j) == LEFT || S(i, j) == BOTTOM)) - S(i, j) = BOTTOMLEFT; // TOPLEFT - if (S(i + 1, j - 1) == NONE && S(i, j - 1) == NONE && - S(i + 1, j) == NONE && S(i - 1, j + 1) == LOCAL && - (S(i, j) == LOCAL || S(i, j) == RIGHT || S(i, j) == BOTTOM)) - S(i, j) = BOTTOMRIGHT; // TOPRIGHT - if (S(i - 1, j + 1) == NONE && S(i - 1, j) == NONE && - S(i, j + 1) == NONE && S(i + 1, j - 1) == LOCAL && - (S(i, j) == LOCAL || S(i, j) == LEFT || S(i, j) == TOP)) - S(i, j) = TOPLEFT; // BOTTOMLEFT - if (S(i + 1, j + 1) == NONE && S(i + 1, j) == NONE && - S(i, j + 1) == NONE && S(i - 1, j - 1) == LOCAL && - (S(i, j) == LOCAL || S(i, j) == RIGHT || S(i, j) == TOP)) - S(i, j) = TOPRIGHT; // BOTTOMRIGHT - } - } - } - -#ifdef VERBOSE - printConfig(solver); -#endif -} - -static double maxElement(Solver* solver, double* m) -{ - int size = (solver->grid.imax + 2) * (solver->grid.jmax + 2); - double maxval = DBL_MIN; - - for (int i = 0; i < size; i++) { - maxval = MAX(maxval, fabs(m[i])); - } - - return maxval; -} - -void computeRHS(Solver* solver) -{ - int imax = solver->grid.imax; - int jmax = solver->grid.jmax; - double idx = 1.0 / solver->grid.dx; - double idy = 1.0 / solver->grid.dy; - double idt = 1.0 / solver->dt; - double* rhs = solver->rhs; - double* f = solver->f; - double* g = solver->g; - int* s = solver->s; - - 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 normalizePressure(Solver* solver) -{ - int size = (solver->grid.imax + 2) * (solver->grid.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->grid.dx; - double dy = solver->grid.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->grid.imax; - int jmax = solver->grid.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->grid.imax; - int jmax = solver->grid.jmax; - double mDy = solver->grid.dy; - double* u = solver->u; - int* s = solver->s; - - 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->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); - } - } else if (strcmp(solver->problem, "backstep") == 0) { - for (int j = 1; j < jmax + 1; j++) { - if (S(0, j) == NONE) U(0, j) = 1.0; - } - } else if (strcmp(solver->problem, "karman") == 0) { - for (int j = 1; j < jmax + 1; j++) { - U(0, j) = 1.0; - } - } -} - -void setObjectBoundaryCondition(Solver* solver) -{ - int imax = solver->grid.imax; - int jmax = solver->grid.jmax; - double* u = solver->u; - double* v = solver->v; - int* s = solver->s; - - for (int j = 1; j < jmax + 1; j++) { - for (int i = 1; i < imax + 1; i++) { - switch (S(i, j)) { - case TOP: - U(i, j) = -U(i, j + 1); - U(i - 1, j) = -U(i - 1, j + 1); - V(i, j) = 0.0; - break; - case BOTTOM: - U(i, j) = -U(i, j - 1); - U(i - 1, j) = -U(i - 1, j - 1); - V(i, j) = 0.0; - break; - case LEFT: - U(i - 1, j) = 0.0; - V(i, j) = -V(i - 1, j); - V(i, j - 1) = -V(i - 1, j - 1); - break; - case RIGHT: - U(i, j) = 0.0; - V(i, j) = -V(i + 1, j); - V(i, j - 1) = -V(i + 1, j - 1); - break; - case TOPLEFT: - U(i, j) = -U(i, j + 1); - U(i - 1, j) = 0.0; - V(i, j) = 0.0; - V(i, j - 1) = -V(i - 1, j - 1); - break; - case TOPRIGHT: - U(i, j) = 0.0; - U(i - 1, j) = -U(i - 1, j + 1); - V(i, j) = 0.0; - V(i, j - 1) = -V(i + 1, j - 1); - break; - case BOTTOMLEFT: - U(i, j) = -U(i, j - 1); - U(i - 1, j) = 0.0; - V(i, j) = -V(i - 1, j); - V(i, j - 1) = 0.0; - break; - case BOTTOMRIGHT: - U(i, j) = 0.0; - U(i - 1, j) = -U(i - 1, j - 1); - V(i, j) = -V(i, j + 1); - V(i, j - 1) = 0.0; - break; - } - } - } -} - -void computeFG(Solver* solver) -{ - double* u = solver->u; - double* v = solver->v; - double* f = solver->f; - double* g = solver->g; - int* s = solver->s; - int imax = solver->grid.imax; - int jmax = solver->grid.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->grid.dx; - double inverseDy = 1.0 / solver->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++) { - if (S(i, j) == NONE) { - 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); - } else { - switch (S(i, j)) { - case TOP: - G(i, j) = V(i, j); - break; - case BOTTOM: - G(i, j - 1) = V(i, j - 1); - break; - case LEFT: - F(i - 1, j) = U(i - 1, j); - break; - case RIGHT: - F(i, j) = U(i, j); - break; - case TOPLEFT: - F(i - 1, j) = U(i - 1, j); - G(i, j) = V(i, j); - break; - case TOPRIGHT: - F(i, j) = U(i, j); - G(i, j) = V(i, j); - break; - case BOTTOMLEFT: - F(i - 1, j) = U(i - 1, j); - G(i, j - 1) = V(i, j - 1); - break; - case BOTTOMRIGHT: - F(i, j) = U(i, j); - G(i, j - 1) = V(i, j - 1); - break; - } - } - } - } - - /* ---------------------- 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->grid.imax; - int jmax = solver->grid.jmax; - double* p = solver->p; - double* u = solver->u; - double* v = solver->v; - int* s = solver->s; - double* f = solver->f; - double* g = solver->g; - double factorX = solver->dt / solver->grid.dx; - double factorY = solver->dt / solver->grid.dy; - - double val = 0; - - 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; - } - } -} - -double smoothRB(Solver* solver) -{ - 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* p = solver->p; - double* r = solver->r[solver->currentlevel]; - double* rhs = solver->rhs; - double epssq = eps * eps; - int it = 0; - 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 multiGrid(Solver* solver) -{ - double res = 0.0; - int imax = solver->grid.imax; - int jmax = solver->grid.jmax; - if (solver->currentlevel == (solver->levels - 1)) { - for (int i = 0; i < 5; i++) { - smoothRB(solver); - } - return; - } - - for (int i = 0; i < 5; i++) { - smoothRB(solver); - if (solver->currentlevel == 0) { - - double* p = solver->p; - - 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); - } - } - } - Solver coarseSolver = copySolver(solver); - - // restrict - restrictMG(solver); - - coarseSolver.p = solver->e[coarseSolver.currentlevel]; - coarseSolver.rhs = solver->r[coarseSolver.currentlevel]; - coarseSolver.grid.imax /= 2; - coarseSolver.grid.jmax /= 2; - - // MGSolver on residual and error. - multiGrid(&coarseSolver); - - // prolongate - prolongate(solver); - - // correct p on finest level using residual - correct(solver); - - if (solver->currentlevel == 0) { - - double* p = solver->p; - - 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); - } - } - - for (int i = 0; i < 5; i++) { - res = smoothRB(solver); - if (solver->currentlevel == 0) { - - double* p = solver->p; - - 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); - } - } - } -#ifdef VERBOSE - if (solver->currentlevel == 0) { - printf("Residuum: %.6f\n", res); - } -#endif -} - -void restrictMG(Solver* solver) -{ - int imax = solver->grid.imax; - int jmax = solver->grid.jmax; - double* r = solver->r[solver->currentlevel + 1]; - double* oldr = solver->r[solver->currentlevel]; - - for (int j = 1; j < jmax + 1; j++) { - for (int i = 1; i < imax + 1; ++i) { - R(i, j) = (oldR(2 * i - 1, 2 * j - 1) + oldR(2 * i, 2 * j - 1) * 2 + - oldR(2 * i + 1, 2 * j - 1) + oldR(2 * i - 1, 2 * j) * 2 + - oldR(2 * i, 2 * j) * 4 + oldR(2 * i + 1, 2 * j) * 2 + - oldR(2 * i - 1, 2 * j + 1) + oldR(2 * i, 2 * j + 1) * 2 + - oldR(2 * i + 1, 2 * j + 1)) / - 16.0; - } - } -} - -void prolongate(Solver* solver) -{ - int imax = solver->grid.imax; - int jmax = solver->grid.jmax; - double* olde = solver->r[solver->currentlevel + 1]; - double* e = solver->r[solver->currentlevel]; - - for (int j = 2; j < jmax + 1; j += 2) { - for (int i = 2; i < imax + 1; i += 2) { - E(i, j) = oldE(i / 2, j / 2); - } - } -} - -void correct(Solver* solver) -{ - int imax = solver->grid.imax; - int jmax = solver->grid.jmax; - double* p = solver->p; - double* e = solver->e[solver->currentlevel]; - for (int j = 1; j < jmax + 1; ++j) { - for (int i = 1; i < imax + 1; ++i) { - P(i, j) += E(i, j); - } - } -} - -Solver copySolver(Solver* solver) -{ - Solver newSolver; - newSolver.problem = solver->problem; - newSolver.bcLeft = solver->bcLeft; - newSolver.bcRight = solver->bcRight; - newSolver.bcBottom = solver->bcBottom; - newSolver.bcTop = solver->bcTop; - newSolver.grid.imax = solver->grid.imax; - newSolver.grid.jmax = solver->grid.jmax; - newSolver.grid.xlength = solver->grid.xlength; - newSolver.grid.ylength = solver->grid.ylength; - newSolver.grid.dx = solver->grid.xlength / solver->grid.imax; - newSolver.grid.dy = solver->grid.ylength / solver->grid.jmax; - newSolver.eps = solver->eps; - newSolver.omega = solver->omega; - newSolver.itermax = solver->itermax; - newSolver.re = solver->re; - newSolver.gx = solver->gx; - newSolver.gy = solver->gy; - newSolver.dt = solver->dt; - newSolver.te = solver->te; - newSolver.tau = solver->tau; - newSolver.gamma = solver->gamma; - newSolver.rho = solver->rho; - newSolver.levels = solver->levels; - newSolver.currentlevel = solver->currentlevel + 1; - - newSolver.r = solver->r; - newSolver.e = solver->e; - - return newSolver; -} - -void writeResult(Solver* solver) -{ - int imax = solver->grid.imax; - int jmax = solver->grid.jmax; - double dx = solver->grid.dx; - double dy = solver->grid.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/EnhancedSolver/2D-seq/src/solver.h b/EnhancedSolver/2D-seq/src/solver.h index a72e4e0..d7ce513 100644 --- a/EnhancedSolver/2D-seq/src/solver.h +++ b/EnhancedSolver/2D-seq/src/solver.h @@ -6,62 +6,21 @@ */ #ifndef __SOLVER_H_ #define __SOLVER_H_ +#include "discretization.h" #include "grid.h" #include "parameter.h" -enum OBJECTBOUNDARY { - NONE = 0, - TOP, - BOTTOM, - LEFT, - RIGHT, - TOPLEFT, - BOTTOMLEFT, - TOPRIGHT, - BOTTOMRIGHT, - LOCAL -}; -enum BC { NOSLIP = 1, SLIP, OUTFLOW, PERIODIC }; -/// @brief -enum SHAPE { NOSHAPE = 0, RECT, CIRCLE }; - typedef struct { /* geometry and grid information */ - Grid grid; - /* arrays */ - double *p, *rhs, **r, **e; - double *f, *g; - double *u, *v; - int* s; + Grid* grid; /* parameters */ double eps, omega, rho; - double re, tau, gamma; - double gx, gy; - /* time stepping */ - int itermax, levels, currentlevel; - double dt, te; - double dtBound; - char* problem; - int bcLeft, bcRight, bcBottom, bcTop; + int itermax; + int levels; + double **r, **e; } Solver; -extern void initSolver(Solver*, Parameter*); -extern void computeRHS(Solver*); -extern double smoothRB(Solver*); -extern void restrictMG(Solver*); -extern void prolongate(Solver*); -extern void correct(Solver*); -extern Solver copySolver(Solver*); -extern void multiGrid(Solver*); -extern void normalizePressure(Solver*); -extern void computeTimestep(Solver*); -extern void setBoundaryConditions(Solver*); -extern void setSpecialBoundaryCondition(Solver*); -extern void setObjectBoundaryCondition(Solver*); -extern void computeFG(Solver*); -extern void adaptUV(Solver*); -extern void writeResult(Solver*); -extern void print(Solver*, double*); -extern void printGrid(Solver*, int*); +extern void initSolver(Solver*, Discretization*, Parameter*); +extern void solve(Solver*, double*, double*); #endif diff --git a/EnhancedSolver/2D-seq/src/timing.c b/EnhancedSolver/2D-seq/src/timing.c index 50bc72f..78b01c4 100644 --- a/EnhancedSolver/2D-seq/src/timing.c +++ b/EnhancedSolver/2D-seq/src/timing.c @@ -7,18 +7,16 @@ #include #include -double getTimeStamp() +double getTimeStamp(void) { struct timespec ts; clock_gettime(CLOCK_MONOTONIC, &ts); return (double)ts.tv_sec + (double)ts.tv_nsec * 1.e-9; } -double getTimeResolution() +double getTimeResolution(void) { struct timespec ts; clock_getres(CLOCK_MONOTONIC, &ts); return (double)ts.tv_sec + (double)ts.tv_nsec * 1.e-9; } - -double getTimeStamp_() { return getTimeStamp(); } diff --git a/EnhancedSolver/2D-seq/src/timing.h b/EnhancedSolver/2D-seq/src/timing.h index d00aea9..ed05a8c 100644 --- a/EnhancedSolver/2D-seq/src/timing.h +++ b/EnhancedSolver/2D-seq/src/timing.h @@ -7,8 +7,7 @@ #ifndef __TIMING_H_ #define __TIMING_H_ -extern double getTimeStamp(); -extern double getTimeResolution(); -extern double getTimeStamp_(); +extern double getTimeStamp(void); +extern double getTimeResolution(void); #endif // __TIMING_H_ diff --git a/EnhancedSolver/2D-seq/src/util.h b/EnhancedSolver/2D-seq/src/util.h index 780c778..e4a93c2 100644 --- a/EnhancedSolver/2D-seq/src/util.h +++ b/EnhancedSolver/2D-seq/src/util.h @@ -7,8 +7,7 @@ #ifndef __UTIL_H_ #define __UTIL_H_ #define HLINE \ - "------------------------------------------------------------------------" \ - "----\n" + "----------------------------------------------------------------------------\n" #ifndef MIN #define MIN(x, y) ((x) < (y) ? (x) : (y)) @@ -20,4 +19,11 @@ #define ABS(a) ((a) >= 0 ? (a) : -(a)) #endif +#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_ diff --git a/EnhancedSolver/2D-seq/vis_files/canal.plot b/EnhancedSolver/2D-seq/vis_files/canal.plot new file mode 100644 index 0000000..bf9ce60 --- /dev/null +++ b/EnhancedSolver/2D-seq/vis_files/canal.plot @@ -0,0 +1,10 @@ +unset border; unset tics; unset key; +set term gif animate delay 30 +set output "trace.gif" +set xrange [0:30] +set yrange [0:4] + +do for [ts=0:120] { + plot "particles_".ts.".dat" with points pointtype 7 +} +unset output