From fe57042556780d1027a41bc62c669200443463c3 Mon Sep 17 00:00:00 2001 From: Jan Eitzinger Date: Thu, 15 Feb 2024 09:44:06 +0100 Subject: [PATCH] Integrate sequential enhanced solver --- EnhancedSolver/2D-seq/Makefile | 73 ++ EnhancedSolver/2D-seq/README.md | 78 ++ EnhancedSolver/2D-seq/backstep.par | 72 ++ EnhancedSolver/2D-seq/canal.par | 73 ++ EnhancedSolver/2D-seq/config.mk | 12 + EnhancedSolver/2D-seq/dcavity.par | 73 ++ EnhancedSolver/2D-seq/include_CLANG.mk | 17 + EnhancedSolver/2D-seq/include_GCC.mk | 14 + EnhancedSolver/2D-seq/include_ICC.mk | 14 + EnhancedSolver/2D-seq/karman.par | 73 ++ EnhancedSolver/2D-seq/src/affinity.c | 61 ++ EnhancedSolver/2D-seq/src/affinity.h | 14 + EnhancedSolver/2D-seq/src/allocate.c | 35 + EnhancedSolver/2D-seq/src/allocate.h | 13 + EnhancedSolver/2D-seq/src/grid.h | 16 + EnhancedSolver/2D-seq/src/likwid-marker.h | 54 ++ EnhancedSolver/2D-seq/src/main.c | 83 ++ EnhancedSolver/2D-seq/src/parameter.c | 147 +++ EnhancedSolver/2D-seq/src/parameter.h | 35 + EnhancedSolver/2D-seq/src/particletracing.c | 251 +++++ EnhancedSolver/2D-seq/src/particletracing.h | 46 + EnhancedSolver/2D-seq/src/progress.c | 51 + EnhancedSolver/2D-seq/src/progress.h | 14 + EnhancedSolver/2D-seq/src/solver.c | 900 ++++++++++++++++++ EnhancedSolver/2D-seq/src/solver.h | 67 ++ EnhancedSolver/2D-seq/src/timing.c | 24 + EnhancedSolver/2D-seq/src/timing.h | 14 + EnhancedSolver/2D-seq/src/util.h | 23 + EnhancedSolver/2D-seq/src/vtkWriter.c | 125 +++ EnhancedSolver/2D-seq/src/vtkWriter.h | 25 + EnhancedSolver/2D-seq/surface.plot | 7 + EnhancedSolver/2D-seq/vector.plot | 7 + EnhancedSolver/2D-seq/vis_files/animate.plot | 10 + .../2D-seq/vis_files/backstep_animate.plot | 14 + .../2D-seq/vis_files/karman_animate.plot | 13 + 35 files changed, 2548 insertions(+) create mode 100644 EnhancedSolver/2D-seq/Makefile create mode 100644 EnhancedSolver/2D-seq/README.md create mode 100644 EnhancedSolver/2D-seq/backstep.par create mode 100644 EnhancedSolver/2D-seq/canal.par create mode 100644 EnhancedSolver/2D-seq/config.mk create mode 100644 EnhancedSolver/2D-seq/dcavity.par create mode 100644 EnhancedSolver/2D-seq/include_CLANG.mk create mode 100644 EnhancedSolver/2D-seq/include_GCC.mk create mode 100644 EnhancedSolver/2D-seq/include_ICC.mk create mode 100644 EnhancedSolver/2D-seq/karman.par create mode 100644 EnhancedSolver/2D-seq/src/affinity.c create mode 100644 EnhancedSolver/2D-seq/src/affinity.h create mode 100644 EnhancedSolver/2D-seq/src/allocate.c create mode 100644 EnhancedSolver/2D-seq/src/allocate.h create mode 100644 EnhancedSolver/2D-seq/src/grid.h create mode 100644 EnhancedSolver/2D-seq/src/likwid-marker.h create mode 100644 EnhancedSolver/2D-seq/src/main.c create mode 100644 EnhancedSolver/2D-seq/src/parameter.c create mode 100644 EnhancedSolver/2D-seq/src/parameter.h create mode 100644 EnhancedSolver/2D-seq/src/particletracing.c create mode 100644 EnhancedSolver/2D-seq/src/particletracing.h create mode 100644 EnhancedSolver/2D-seq/src/progress.c create mode 100644 EnhancedSolver/2D-seq/src/progress.h create mode 100644 EnhancedSolver/2D-seq/src/solver.c create mode 100644 EnhancedSolver/2D-seq/src/solver.h create mode 100644 EnhancedSolver/2D-seq/src/timing.c create mode 100644 EnhancedSolver/2D-seq/src/timing.h create mode 100644 EnhancedSolver/2D-seq/src/util.h create mode 100644 EnhancedSolver/2D-seq/src/vtkWriter.c create mode 100644 EnhancedSolver/2D-seq/src/vtkWriter.h create mode 100644 EnhancedSolver/2D-seq/surface.plot create mode 100644 EnhancedSolver/2D-seq/vector.plot create mode 100644 EnhancedSolver/2D-seq/vis_files/animate.plot create mode 100644 EnhancedSolver/2D-seq/vis_files/backstep_animate.plot create mode 100644 EnhancedSolver/2D-seq/vis_files/karman_animate.plot diff --git a/EnhancedSolver/2D-seq/Makefile b/EnhancedSolver/2D-seq/Makefile new file mode 100644 index 0000000..aab8407 --- /dev/null +++ b/EnhancedSolver/2D-seq/Makefile @@ -0,0 +1,73 @@ +#======================================================================================= +# 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. +#======================================================================================= + +#CONFIGURE BUILD SYSTEM +TARGET = exe-$(TAG) +BUILD_DIR = ./$(TAG) +SRC_DIR = ./src +MAKE_DIR = ./ +Q ?= @ + +#DO NOT EDIT BELOW +include $(MAKE_DIR)/config.mk +include $(MAKE_DIR)/include_$(TAG).mk +INCLUDES += -I$(SRC_DIR) -I$(BUILD_DIR) + +VPATH = $(SRC_DIR) +SRC = $(wildcard $(SRC_DIR)/*.c) +ASM = $(patsubst $(SRC_DIR)/%.c, $(BUILD_DIR)/%.s, $(SRC)) +OBJ = $(patsubst $(SRC_DIR)/%.c, $(BUILD_DIR)/%.o, $(SRC)) +SOURCES = $(SRC) $(wildcard $(SRC_DIR)/*.h) +CPPFLAGS := $(CPPFLAGS) $(DEFINES) $(OPTIONS) $(INCLUDES) + +${TARGET}: $(BUILD_DIR) $(OBJ) + $(info ===> LINKING $(TARGET)) + $(Q)${LINKER} ${LFLAGS} -o $(TARGET) $(OBJ) $(LIBS) + +$(BUILD_DIR)/%.o: %.c $(MAKE_DIR)/include_$(TAG).mk $(MAKE_DIR)/config.mk + $(info ===> COMPILE $@) + $(CC) -c $(CPPFLAGS) $(CFLAGS) $< -o $@ + $(Q)$(GCC) $(CPPFLAGS) -MT $(@:.d=.o) -MM $< > $(BUILD_DIR)/$*.d + +$(BUILD_DIR)/%.s: %.c + $(info ===> GENERATE ASM $@) + $(CC) -S $(CPPFLAGS) $(CFLAGS) $< -o $@ + +.PHONY: clean distclean tags info asm format + +clean: + $(info ===> CLEAN) + @rm -rf $(BUILD_DIR) + @rm -f tags + +distclean: clean + $(info ===> DIST CLEAN) + @rm -f $(TARGET) + @rm -f *.dat + @rm -f *.png + +info: + $(info $(CFLAGS)) + $(Q)$(CC) $(VERSION) + +asm: $(BUILD_DIR) $(ASM) + +tags: + $(info ===> GENERATE TAGS) + $(Q)ctags -R + +format: + @for src in $(SOURCES) ; do \ + echo "Formatting $$src" ; \ + clang-format -i $$src ; \ + done + @echo "Done" + +$(BUILD_DIR): + @mkdir $(BUILD_DIR) + +-include $(OBJ:.o=.d) diff --git a/EnhancedSolver/2D-seq/README.md b/EnhancedSolver/2D-seq/README.md new file mode 100644 index 0000000..d980b54 --- /dev/null +++ b/EnhancedSolver/2D-seq/README.md @@ -0,0 +1,78 @@ +# C source skeleton + +## Build + +1. Configure the toolchain and additional options in `config.mk`: +``` +# Supported: GCC, CLANG, ICC +TAG ?= GCC +ENABLE_OPENMP ?= false + +OPTIONS += -DARRAY_ALIGNMENT=64 +#OPTIONS += -DVERBOSE +#OPTIONS += -DVERBOSE_AFFINITY +#OPTIONS += -DVERBOSE_DATASIZE +#OPTIONS += -DVERBOSE_TIMER +``` + +The verbosity options enable detailed output about solver, affinity settings, allocation sizes and timer resolution. +For debugging you may want to set the VERBOSE option: +``` +# Supported: GCC, CLANG, ICC +TAG ?= GCC +ENABLE_OPENMP ?= false + +OPTIONS += -DARRAY_ALIGNMENT=64 +OPTIONS += -DVERBOSE +#OPTIONS += -DVERBOSE_AFFINITY +#OPTIONS += -DVERBOSE_DATASIZE +#OPTIONS += -DVERBOSE_TIMER +` + +2. Build with: +``` +make +``` + +You can build multiple toolchains in the same directory, but notice that the Makefile is only acting on the one currently set. +Intermediate build results are located in the `` directory. + +To output the executed commands use: +``` +make Q= +``` + +3. Clean up with: +``` +make clean +``` +to clean intermediate build results. + +``` +make distclean +``` +to clean intermediate build results and binary. + +4. (Optional) Generate assembler: +``` +make asm +``` +The assembler files will also be located in the `` directory. + +## Usage + +You have to provide a parameter file describing the problem you want to solve: +``` +./exe-CLANG dcavity.par +``` + +Examples are given in in dcavity (a lid driven cavity test case) and canal (simulating a empty canal). + +You can plot the resulting velocity and pressure fields using gnuplot: +``` +gnuplot vector.plot +``` +and for the pressure: +``` +gnuplot surface.plot +``` diff --git a/EnhancedSolver/2D-seq/backstep.par b/EnhancedSolver/2D-seq/backstep.par new file mode 100644 index 0000000..656ac31 --- /dev/null +++ b/EnhancedSolver/2D-seq/backstep.par @@ -0,0 +1,72 @@ +#============================================================================== +# Laminar Canal Flow +#============================================================================== + +# Problem specific Data: +# --------------------- + +name backstep # name of flow setup + +bcTop 1 # flags for boundary conditions +bcBottom 1 # 1 = no-slip 3 = outflow +bcLeft 3 # 2 = free-slip 4 = periodic +bcRight 3 # + +gx 0.0 # Body forces (e.g. gravity) +gy 0.0 # + +re 36500.0 # Reynolds number + +u_init 1.0 # initial value for velocity in x-direction +v_init 0.0 # initial value for velocity in y-direction +p_init 1.0 # initial value for pressure + +# Geometry Data: +# ------------- + +xlength 7.0 # domain size in x-direction +ylength 1.5 # domain size in y-direction +imax 210 # number of interior cells in x-direction +jmax 45 # number of interior cells in y-direction + +# Time Data: +# --------- + +te 60.0 # final time +dt 0.02 # time stepsize +tau 0.5 # safety factor for time stepsize control (<0 constant delt) + +# Pressure Iteration Data: +# ----------------------- + +itermax 500 # maximal number of pressure iteration in one time step +eps 0.0001 # stopping tolerance for pressure iteration +rho 0.52 +omg 1.8 # relaxation parameter for SOR iteration +gamma 0.9 # upwind differencing factor gamma + +# Particle Tracing Data: +# ----------------------- + +numberOfParticles 200 +startTime 100 +injectTimePeriod 1.0 +writeTimePeriod 0.5 + +x1 0.0 +y1 0.5 +x2 0.0 +y2 1.5 + +# Obstacle Geometry Data: +# ----------------------- +# Shape 0 disable, 1 Rectangle/Square, 2 Circle + +shape 1 +xCenter 0.0 +yCenter 0.0 +xRectLength 2.0 +yRectLength 1.0 +circleRadius 1.0 + +#=============================================================================== diff --git a/EnhancedSolver/2D-seq/canal.par b/EnhancedSolver/2D-seq/canal.par new file mode 100644 index 0000000..ca30e51 --- /dev/null +++ b/EnhancedSolver/2D-seq/canal.par @@ -0,0 +1,73 @@ +#============================================================================== +# Laminar Canal Flow +#============================================================================== + +# Problem specific Data: +# --------------------- + +name canal # name of flow setup + +bcTop 1 # flags for boundary conditions +bcBottom 1 # 1 = no-slip 3 = outflow +bcLeft 3 # 2 = free-slip 4 = periodic +bcRight 3 # + +gx 0.0 # Body forces (e.g. gravity) +gy 0.0 # + +re 100.0 # Reynolds number + +u_init 1.0 # initial value for velocity in x-direction +v_init 0.0 # initial value for velocity in y-direction +p_init 1.0 # initial value for pressure + +# Geometry Data: +# ------------- + +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 + +# Time Data: +# --------- + +te 60.0 # final time +dt 0.02 # time stepsize +tau 0.5 # safety factor for time stepsize control (<0 constant delt) + +# Pressure Iteration Data: +# ----------------------- + +itermax 500 # maximal number of pressure iteration in one time step +eps 0.0001 # stopping tolerance for pressure iteration +rho 0.52 +omg 1.8 # relaxation parameter for SOR iteration +gamma 0.9 # upwind differencing factor gamma +levels 5 # Multigrid levels + +# Particle Tracing Data: +# ----------------------- + +numberOfParticles 60 +startTime 100 +injectTimePeriod 2.0 +writeTimePeriod 0.5 + +x1 1.0 +y1 0.0 +x2 1.0 +y2 4.0 + +# Obstacle Geometry Data: +# ----------------------- +# Shape 0 disable, 1 Rectangle/Square, 2 Circle + +shape 0 +xCenter 10.0 +yCenter 2 +xRectLength 6.0 +yRectLength 1.0 +circleRadius 1.0 + +#=============================================================================== diff --git a/EnhancedSolver/2D-seq/config.mk b/EnhancedSolver/2D-seq/config.mk new file mode 100644 index 0000000..d228f66 --- /dev/null +++ b/EnhancedSolver/2D-seq/config.mk @@ -0,0 +1,12 @@ +# Supported: GCC, CLANG, ICC +TAG ?= CLANG +ENABLE_OPENMP ?= false + +#Feature options +OPTIONS += -DARRAY_ALIGNMENT=64 +# 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 new file mode 100644 index 0000000..11dab71 --- /dev/null +++ b/EnhancedSolver/2D-seq/dcavity.par @@ -0,0 +1,73 @@ +#============================================================================== +# Driven Cavity +#============================================================================== + +# Problem specific Data: +# --------------------- + +name dcavity # name of flow setup + +bcTop 1 # flags for boundary conditions +bcBottom 1 # 1 = no-slip 3 = outflow +bcLeft 1 # 2 = free-slip 4 = periodic +bcRight 1 # + +gx 0.0 # Body forces (e.g. gravity) +gy 0.0 # + +re 10.0 # Reynolds number + +u_init 1.0 # initial value for velocity in x-direction +v_init 0.0 # initial value for velocity in y-direction +p_init 0.0 # initial value for pressure + +# Geometry Data: +# ------------- + +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 + +# Time Data: +# --------- + +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: +# ----------------------- + +itermax 1000 # maximal number of pressure iteration in one time step +eps 0.001 # stopping tolerance for pressure iteration +rho 0.5 +omg 1.8 # relaxation parameter for SOR iteration +gamma 0.9 # upwind differencing factor gamma +levels 5 # Multigrid levels + +# Particle Tracing Data: +# ----------------------- + +numberOfParticles 200 +startTime 100 +injectTimePeriod 0.5 +writeTimePeriod 0.2 + +x1 0.1 +y1 0.9 +x2 0.9 +y2 0.9 + + +# Obstacle Geometry Data: +# ----------------------- +# Shape 0 disable, 1 Rectangle/Square, 2 Circle + +shape 0 +xCenter 0.5 +yCenter 0.5 +xRectLength 0.5 +yRectLength 0.5 +circleRadius 0.5 +#=============================================================================== diff --git a/EnhancedSolver/2D-seq/include_CLANG.mk b/EnhancedSolver/2D-seq/include_CLANG.mk new file mode 100644 index 0000000..466f441 --- /dev/null +++ b/EnhancedSolver/2D-seq/include_CLANG.mk @@ -0,0 +1,17 @@ +CC = clang +GCC = cc +LINKER = $(CC) + +ifeq ($(ENABLE_OPENMP),true) +OPENMP = -fopenmp +#OPENMP = -Xpreprocessor -fopenmp #required on Macos with homebrew libomp +LIBS = # -lomp +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 +INCLUDES = diff --git a/EnhancedSolver/2D-seq/include_GCC.mk b/EnhancedSolver/2D-seq/include_GCC.mk new file mode 100644 index 0000000..427e798 --- /dev/null +++ b/EnhancedSolver/2D-seq/include_GCC.mk @@ -0,0 +1,14 @@ +CC = gcc +GCC = gcc +LINKER = $(CC) + +ifeq ($(ENABLE_OPENMP),true) +OPENMP = -fopenmp +endif + +VERSION = --version +CFLAGS = -Ofast -ffreestanding -std=c99 $(OPENMP) +LFLAGS = $(OPENMP) +DEFINES = -D_GNU_SOURCE +INCLUDES = +LIBS = diff --git a/EnhancedSolver/2D-seq/include_ICC.mk b/EnhancedSolver/2D-seq/include_ICC.mk new file mode 100644 index 0000000..d0c2a1b --- /dev/null +++ b/EnhancedSolver/2D-seq/include_ICC.mk @@ -0,0 +1,14 @@ +CC = icx +GCC = gcc +LINKER = $(CC) + +ifeq ($(ENABLE_OPENMP),true) +OPENMP = -qopenmp +endif + +VERSION = --version +CFLAGS = -O3 -xHost -qopt-zmm-usage=high -std=c99 $(OPENMP) +LFLAGS = $(OPENMP) +DEFINES = -D_GNU_SOURCE +INCLUDES = +LIBS = diff --git a/EnhancedSolver/2D-seq/karman.par b/EnhancedSolver/2D-seq/karman.par new file mode 100644 index 0000000..60664b6 --- /dev/null +++ b/EnhancedSolver/2D-seq/karman.par @@ -0,0 +1,73 @@ +#============================================================================== +# Laminar Canal Flow +#============================================================================== + +# Problem specific Data: +# --------------------- + +name karman # name of flow setup + +bcTop 1 # flags for boundary conditions +bcBottom 1 # 1 = no-slip 3 = outflow +bcLeft 3 # 2 = free-slip 4 = periodic +bcRight 3 # + +gx 0.0 # Body forces (e.g. gravity) +gy 0.0 # + +re 5050.0 # Reynolds number + +u_init 1.0 # initial value for velocity in x-direction +v_init 0.0 # initial value for velocity in y-direction +p_init 0.0 # initial value for pressure + +# Geometry Data: +# ------------- + +xlength 30.0 # domain size in x-direction +ylength 8.0 # domain size in y-direction +imax 400 # number of interior cells in x-direction +jmax 200 # number of interior cells in y-direction + +# Time Data: +# --------- + +te 150.0 # final time +dt 0.02 # time stepsize +tau 0.5 # safety factor for time stepsize control (<0 constant delt) + +# Pressure Iteration Data: +# ----------------------- + +itermax 200 # maximal number of pressure iteration in one time step +eps 0.001 # stopping tolerance for pressure iteration +rho 0.52 +omg 1.75 # relaxation parameter for SOR iteration +gamma 0.9 # upwind differencing factor gamma +levels 5 # Multigrid levels + +# Particle Tracing Data: +# ----------------------- + +numberOfParticles 200 +startTime 201 +injectTimePeriod 1.0 +writeTimePeriod 0.5 + +x1 0.0 +y1 3.8 +x2 0.0 +y2 4.1 + +# Obstacle Geometry Data: +# ----------------------- +# Shape 0 disable, 1 Rectangle/Square, 2 Circle + +shape 2 +xCenter 5.0 +yCenter 4.0 +xRectLength 2.0 +yRectLength 1.0 +circleRadius 1.0 + +#=============================================================================== diff --git a/EnhancedSolver/2D-seq/src/affinity.c b/EnhancedSolver/2D-seq/src/affinity.c new file mode 100644 index 0000000..ef8355a --- /dev/null +++ b/EnhancedSolver/2D-seq/src/affinity.c @@ -0,0 +1,61 @@ +/* + * 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 new file mode 100644 index 0000000..84bf733 --- /dev/null +++ b/EnhancedSolver/2D-seq/src/affinity.h @@ -0,0 +1,14 @@ +/* + * 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 new file mode 100644 index 0000000..855cf32 --- /dev/null +++ b/EnhancedSolver/2D-seq/src/allocate.c @@ -0,0 +1,35 @@ +/* + * 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. + */ +#include +#include +#include + +void* allocate(int alignment, size_t bytesize) +{ + int errorCode; + void* ptr; + + errorCode = posix_memalign(&ptr, alignment, bytesize); + + if (errorCode) { + if (errorCode == EINVAL) { + fprintf(stderr, "Error: Alignment parameter is not a power of two\n"); + exit(EXIT_FAILURE); + } + if (errorCode == ENOMEM) { + fprintf(stderr, "Error: Insufficient memory to fulfill the request\n"); + exit(EXIT_FAILURE); + } + } + + if (ptr == NULL) { + fprintf(stderr, "Error: posix_memalign failed!\n"); + exit(EXIT_FAILURE); + } + + return ptr; +} diff --git a/EnhancedSolver/2D-seq/src/allocate.h b/EnhancedSolver/2D-seq/src/allocate.h new file mode 100644 index 0000000..1a5c8cb --- /dev/null +++ b/EnhancedSolver/2D-seq/src/allocate.h @@ -0,0 +1,13 @@ +/* + * 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 __ALLOCATE_H_ +#define __ALLOCATE_H_ +#include + +extern void* allocate(int alignment, size_t bytesize); + +#endif diff --git a/EnhancedSolver/2D-seq/src/grid.h b/EnhancedSolver/2D-seq/src/grid.h new file mode 100644 index 0000000..7e58d1a --- /dev/null +++ b/EnhancedSolver/2D-seq/src/grid.h @@ -0,0 +1,16 @@ +/* + * 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 __GRID_H_ +#define __GRID_H_ + +typedef struct { + double dx, dy; + int imax, jmax; + double xlength, ylength; +} Grid; + +#endif // __GRID_H_ diff --git a/EnhancedSolver/2D-seq/src/likwid-marker.h b/EnhancedSolver/2D-seq/src/likwid-marker.h new file mode 100644 index 0000000..eb7cc78 --- /dev/null +++ b/EnhancedSolver/2D-seq/src/likwid-marker.h @@ -0,0 +1,54 @@ +/* + * ======================================================================================= + * + * Author: Jan Eitzinger (je), jan.eitzinger@fau.de + * Copyright (c) 2020 RRZE, University Erlangen-Nuremberg + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), to + * deal in the Software without restriction, including without limitation the + * rights to use, copy, modify, merge, publish, distribute, sublicense, and/or + * sell copies of the Software, and to permit persons to whom the Software is + * furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS + * IN THE SOFTWARE. + * + * ======================================================================================= + */ +#ifndef LIKWID_MARKERS_H +#define LIKWID_MARKERS_H + +#ifdef LIKWID_PERFMON +#include +#define LIKWID_MARKER_INIT likwid_markerInit() +#define LIKWID_MARKER_THREADINIT likwid_markerThreadInit() +#define LIKWID_MARKER_SWITCH likwid_markerNextGroup() +#define LIKWID_MARKER_REGISTER(regionTag) likwid_markerRegisterRegion(regionTag) +#define LIKWID_MARKER_START(regionTag) likwid_markerStartRegion(regionTag) +#define LIKWID_MARKER_STOP(regionTag) likwid_markerStopRegion(regionTag) +#define LIKWID_MARKER_CLOSE likwid_markerClose() +#define LIKWID_MARKER_RESET(regionTag) likwid_markerResetRegion(regionTag) +#define LIKWID_MARKER_GET(regionTag, nevents, events, time, count) \ + likwid_markerGetRegion(regionTag, nevents, events, time, count) +#else /* LIKWID_PERFMON */ +#define LIKWID_MARKER_INIT +#define LIKWID_MARKER_THREADINIT +#define LIKWID_MARKER_SWITCH +#define LIKWID_MARKER_REGISTER(regionTag) +#define LIKWID_MARKER_START(regionTag) +#define LIKWID_MARKER_STOP(regionTag) +#define LIKWID_MARKER_CLOSE +#define LIKWID_MARKER_GET(regionTag, nevents, events, time, count) +#define LIKWID_MARKER_RESET(regionTag) +#endif /* LIKWID_PERFMON */ + +#endif /*LIKWID_MARKERS_H*/ diff --git a/EnhancedSolver/2D-seq/src/main.c b/EnhancedSolver/2D-seq/src/main.c new file mode 100644 index 0000000..d889c3b --- /dev/null +++ b/EnhancedSolver/2D-seq/src/main.c @@ -0,0 +1,83 @@ +/* + * 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. + */ +#include +#include +#include + +#include "parameter.h" +#include "particletracing.h" +#include "progress.h" +#include "solver.h" +#include "timing.h" + +int main(int argc, char** argv) +{ + double S, E; + Parameter params; + Solver solver; + ParticleTracer particletracer; + initParameter(¶ms); + + 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); + printParticleTracerParameters(&particletracer); + +#ifndef VERBOSE + initProgress(solver.te); +#endif + + double tau = solver.tau; + double te = solver.te; + double t = 0.0; + int nt = 0; + + S = 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); + + adaptUV(&solver); + + /* 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; + nt++; + +#ifdef VERBOSE + printf("TIME %f , TIMESTEP %f\n", t, solver.dt); +#else + printProgress(t); +#endif + } + printf("Total particles : %d\n", particletracer.totalParticles); + + E = getTimeStamp(); + stopProgress(); + + freeParticles(&particletracer); + + printf("Solution took %.2fs\n", E - S); + writeResult(&solver); + return EXIT_SUCCESS; +} diff --git a/EnhancedSolver/2D-seq/src/parameter.c b/EnhancedSolver/2D-seq/src/parameter.c new file mode 100644 index 0000000..cf40d51 --- /dev/null +++ b/EnhancedSolver/2D-seq/src/parameter.c @@ -0,0 +1,147 @@ +/* + * 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 "parameter.h" +#include "util.h" +#define MAXLINE 4096 + +void initParameter(Parameter* param) +{ + param->xlength = 1.0; + param->ylength = 1.0; + param->imax = 100; + param->jmax = 100; + param->itermax = 1000; + param->eps = 0.0001; + param->omg = 1.7; + 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) +{ + FILE* fp = fopen(filename, "r"); + char line[MAXLINE]; + int i; + + if (!fp) { + fprintf(stderr, "Could not open parameter file: %s\n", filename); + exit(EXIT_FAILURE); + } + + while (!feof(fp)) { + line[0] = '\0'; + fgets(line, MAXLINE, fp); + for (i = 0; line[i] != '\0' && line[i] != '#'; i++) + ; + line[i] = '\0'; + + char* tok = strtok(line, " "); + char* val = strtok(NULL, " "); + +#define PARSE_PARAM(p, f) \ + if (strncmp(tok, #p, sizeof(#p) / sizeof(#p[0]) - 1) == 0) { \ + param->p = f(val); \ + } +#define PARSE_STRING(p) PARSE_PARAM(p, strdup) +#define PARSE_INT(p) PARSE_PARAM(p, atoi) +#define PARSE_REAL(p) PARSE_PARAM(p, atof) + + if (tok != NULL && val != NULL) { + PARSE_REAL(xlength); + PARSE_REAL(ylength); + PARSE_INT(imax); + PARSE_INT(jmax); + PARSE_INT(itermax); + PARSE_INT(levels); + PARSE_REAL(eps); + PARSE_REAL(omg); + PARSE_REAL(re); + PARSE_REAL(tau); + PARSE_REAL(gamma); + PARSE_REAL(dt); + PARSE_REAL(te); + PARSE_REAL(gx); + PARSE_REAL(gy); + PARSE_STRING(name); + PARSE_INT(bcLeft); + PARSE_INT(bcRight); + PARSE_INT(bcBottom); + PARSE_INT(bcTop); + PARSE_REAL(u_init); + PARSE_REAL(v_init); + PARSE_REAL(p_init); + PARSE_REAL(rho); + + /* Added new particle tracing parameters */ + PARSE_INT(numberOfParticles); + PARSE_REAL(startTime); + PARSE_REAL(injectTimePeriod); + PARSE_REAL(writeTimePeriod); + PARSE_REAL(x1); + PARSE_REAL(y1); + PARSE_REAL(x2); + PARSE_REAL(y2); + + /* Added obstacle geometry parameters */ + PARSE_INT(shape); + PARSE_REAL(xCenter); + PARSE_REAL(yCenter); + PARSE_REAL(xRectLength); + PARSE_REAL(yRectLength); + PARSE_REAL(circleRadius); + } + } + + fclose(fp); +} + +void printParameter(Parameter* param) +{ + printf("Parameters for %s\n", param->name); + printf("Boundary conditions Left:%d Right:%d Bottom:%d Top:%d\n", + param->bcLeft, + param->bcRight, + param->bcBottom, + param->bcTop); + printf("\tReynolds number: %.2f\n", param->re); + printf("\tInit arrays: U:%.2f V:%.2f P:%.2f\n", + param->u_init, + param->v_init, + param->p_init); + printf("Geometry data:\n"); + printf("\tDomain box size (x, y): %.2f, %.2f\n", param->xlength, param->ylength); + printf("\tCells (x, y): %d, %d\n", param->imax, param->jmax); + printf("Timestep parameters:\n"); + printf("\tDefault stepsize: %.2f, Final time %.2f\n", param->dt, param->te); + printf("\tTau factor: %.2f\n", param->tau); + printf("Iterative solver parameters:\n"); + printf("\tMax iterations: %d\n", param->itermax); + 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); + + printf("Particle Tracing data:\n"); + printf("\tNumber of particles : %d being injected for every period of %.2f\n", + param->numberOfParticles, + param->injectTimePeriod); + printf("\tstartTime : %.2f\n", param->startTime); + printf("\t(Line along which the particles are to be injected) \n\tx1 : %.2f, y1 : " + "%.2f, x2 : %.2f, y2 : %.2f\n", + param->x1, + param->y1, + param->x2, + param->y2); +} diff --git a/EnhancedSolver/2D-seq/src/parameter.h b/EnhancedSolver/2D-seq/src/parameter.h new file mode 100644 index 0000000..9b905ca --- /dev/null +++ b/EnhancedSolver/2D-seq/src/parameter.h @@ -0,0 +1,35 @@ +/* + * 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 __PARAMETER_H_ +#define __PARAMETER_H_ + +typedef struct { + double xlength, ylength; + int imax, jmax; + int itermax, levels; + double eps, omg, rho; + double re, tau, gamma; + double te, dt; + double gx, gy; + char* name; + int bcLeft, bcRight, bcBottom, bcTop; + double u_init, v_init, p_init; + + int numberOfParticles; + double startTime, injectTimePeriod, writeTimePeriod; + + double x1, y1, x2, y2; + + int shape; + double xCenter, yCenter, xRectLength, yRectLength, circleRadius; + +} Parameter; + +void initParameter(Parameter*); +void readParameter(Parameter*, const char*); +void printParameter(Parameter*); +#endif diff --git a/EnhancedSolver/2D-seq/src/particletracing.c b/EnhancedSolver/2D-seq/src/particletracing.c new file mode 100644 index 0000000..b95a1a7 --- /dev/null +++ b/EnhancedSolver/2D-seq/src/particletracing.c @@ -0,0 +1,251 @@ +/* + * 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 "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) { + printf("Particle position X : %.2f, Y : %.2f, flag : %d\n", + p->particlePool[i].x, + p->particlePool[i].y, + p->particlePool[i].flag); + } +} +void injectParticles(ParticleTracer* p) +{ + for (int i = 0; i < p->numberOfParticles; ++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); + } +} + +void advanceParticles(ParticleTracer* p, + double* restrict u, + double* restrict v, + int* restrict s, + double time) +{ + 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; + + for (int i = 0; i < p->totalParticles; ++i) { + if (p->particlePool[i].flag == true) { + double x = p->particlePool[i].x; + double y = p->particlePool[i].y; + + int iCoord = (int)(x / dx) + 1; + int jCoord = (int)((y + 0.5 * dy) / dy) + 1; + + double x1 = (double)(iCoord - 1) * dx; + double y1 = ((double)(jCoord - 1) - 0.5) * dy; + 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 new_x = x + p->dt * u_n; + p->particlePool[i].x = new_x; + + iCoord = (int)((x + 0.5 * dx) / dx) + 1; + jCoord = (int)(y / dy) + 1; + + x1 = ((double)(iCoord - 1) - 0.5) * dx; + y1 = (double)(jCoord - 1) * dy; + 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 new_y = y + p->dt * v_n; + p->particlePool[i].y = new_y; + + // 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))) { + p->particlePool[i].flag = false; + } + int i_new = new_x / dx, j_new = new_y / dy; + if (S(i_new, j_new) != NONE) { + p->particlePool[i].flag = false; + } + } + } +} + +void freeParticles(ParticleTracer* p) +{ + if (p->particlePool != NULL) { + free(p->particlePool); + free(p->linSpaceLine); + } +} + +void writeParticles(ParticleTracer* p) +{ + VtkOptions opts = { .particletracer = p }; + + char filename[50]; + // snprintf(filename, 50, "vtk_files/particles%d.vtk", ts); + // vtkOpen(&opts, filename, ts); + // vtkParticle(&opts, "particle"); + // vtkClose(&opts); + + FILE* fp; + Particle* particlePool = p->particlePool; + + snprintf(filename, 50, "vis_files/particles_%d.dat", ts); + fp = fopen(filename, "w"); + + if (fp == NULL) { + printf("Error!\n"); + exit(EXIT_FAILURE); + } + + for (int i = 0; i < p->totalParticles; ++i) { + double x = particlePool[i].x; + double y = particlePool[i].y; + fprintf(fp, "%f %f\n", x, y); + } + fclose(fp); + + ++ts; +} + +void initParticleTracer(ParticleTracer* p, Grid* grid, Parameter* params) +{ + p->numberOfParticles = params->numberOfParticles; + p->startTime = params->startTime; + p->injectTimePeriod = params->injectTimePeriod; + p->writeTimePeriod = params->writeTimePeriod; + + p->dt = params->dt; + p->grid = grid; + + p->x1 = params->x1; + p->y1 = params->y1; + p->x2 = params->x2; + p->y2 = params->y2; + + p->lastInjectTime = params->startTime; + p->lastUpdateTime = params->startTime; + p->lastWriteTime = params->startTime; + + p->pointer = 0; + p->totalParticles = 0; + + if (params->te > params->startTime) { + p->estimatedNumParticles = ((params->te - params->startTime) + 2) * + params->numberOfParticles; + + 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; + } + } else { + p->particlePool = NULL; + p->linSpaceLine = NULL; + } +} + +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->injectTimePeriod); + printf("\tstartTime : %.2f\n", p->startTime); + printf("\t(Line along which the particles are to be injected) \n\tx1 : %.2f, y1 : " + "%.2f, x2 : %.2f, y2 : %.2f\n", + p->x1, + p->y1, + 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) +{ + 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; + + 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->pointer = totalParticles + 1; + memcpy(p->particlePool, tempPool, totalParticles * sizeof(Particle)); +} diff --git a/EnhancedSolver/2D-seq/src/particletracing.h b/EnhancedSolver/2D-seq/src/particletracing.h new file mode 100644 index 0000000..40621d4 --- /dev/null +++ b/EnhancedSolver/2D-seq/src/particletracing.h @@ -0,0 +1,46 @@ +/* + * 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 __PARTICLETRACING_H_ +#define __PARTICLETRACING_H_ +#include + +#include "grid.h" +#include "parameter.h" + +typedef enum COORD { X = 0, Y, NCOORD } COORD; + +typedef struct { + double x, y; + bool flag; +} Particle; + +typedef struct { + int numberOfParticles, totalParticles; + double startTime, injectTimePeriod, writeTimePeriod; + double lastInjectTime, lastUpdateTime, lastWriteTime; + + int estimatedNumParticles; + + double dt; + Particle* linSpaceLine; + Particle* particlePool; + + int pointer; + double x1, y1, x2, y2; + Grid* grid; +} 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*); +#endif diff --git a/EnhancedSolver/2D-seq/src/progress.c b/EnhancedSolver/2D-seq/src/progress.c new file mode 100644 index 0000000..d5393c4 --- /dev/null +++ b/EnhancedSolver/2D-seq/src/progress.c @@ -0,0 +1,51 @@ +/* + * 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 "progress.h" + +static double _end; +static int _current; + +void initProgress(double end) +{ + _end = end; + _current = 0; + + printf("[ ]"); + fflush(stdout); +} + +void printProgress(double current) +{ + int new = (int)rint((current / _end) * 10.0); + + if (new > _current) { + char progress[11]; + _current = new; + progress[0] = 0; + + for (int i = 0; i < 10; i++) { + if (i < _current) { + sprintf(progress + strlen(progress), "#"); + } else { + sprintf(progress + strlen(progress), " "); + } + } + printf("\r[%s]", progress); + } + fflush(stdout); +} + +void stopProgress() +{ + printf("\n"); + fflush(stdout); +} diff --git a/EnhancedSolver/2D-seq/src/progress.h b/EnhancedSolver/2D-seq/src/progress.h new file mode 100644 index 0000000..1cdb9a3 --- /dev/null +++ b/EnhancedSolver/2D-seq/src/progress.h @@ -0,0 +1,14 @@ +/* + * 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 __PROGRESS_H_ +#define __PROGRESS_H_ + +extern void initProgress(double); +extern void printProgress(double); +extern void stopProgress(); + +#endif diff --git a/EnhancedSolver/2D-seq/src/solver.c b/EnhancedSolver/2D-seq/src/solver.c new file mode 100644 index 0000000..bbb51f5 --- /dev/null +++ b/EnhancedSolver/2D-seq/src/solver.c @@ -0,0 +1,900 @@ +/* + * 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 new file mode 100644 index 0000000..a72e4e0 --- /dev/null +++ b/EnhancedSolver/2D-seq/src/solver.h @@ -0,0 +1,67 @@ +/* + * 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 __SOLVER_H_ +#define __SOLVER_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; + /* 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; +} 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*); + +#endif diff --git a/EnhancedSolver/2D-seq/src/timing.c b/EnhancedSolver/2D-seq/src/timing.c new file mode 100644 index 0000000..50bc72f --- /dev/null +++ b/EnhancedSolver/2D-seq/src/timing.c @@ -0,0 +1,24 @@ +/* + * 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. + */ +#include +#include + +double getTimeStamp() +{ + struct timespec ts; + clock_gettime(CLOCK_MONOTONIC, &ts); + return (double)ts.tv_sec + (double)ts.tv_nsec * 1.e-9; +} + +double getTimeResolution() +{ + 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 new file mode 100644 index 0000000..d00aea9 --- /dev/null +++ b/EnhancedSolver/2D-seq/src/timing.h @@ -0,0 +1,14 @@ +/* + * 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 __TIMING_H_ +#define __TIMING_H_ + +extern double getTimeStamp(); +extern double getTimeResolution(); +extern double getTimeStamp_(); + +#endif // __TIMING_H_ diff --git a/EnhancedSolver/2D-seq/src/util.h b/EnhancedSolver/2D-seq/src/util.h new file mode 100644 index 0000000..780c778 --- /dev/null +++ b/EnhancedSolver/2D-seq/src/util.h @@ -0,0 +1,23 @@ +/* + * 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 __UTIL_H_ +#define __UTIL_H_ +#define HLINE \ + "------------------------------------------------------------------------" \ + "----\n" + +#ifndef MIN +#define MIN(x, y) ((x) < (y) ? (x) : (y)) +#endif +#ifndef MAX +#define MAX(x, y) ((x) > (y) ? (x) : (y)) +#endif +#ifndef ABS +#define ABS(a) ((a) >= 0 ? (a) : -(a)) +#endif + +#endif // __UTIL_H_ diff --git a/EnhancedSolver/2D-seq/src/vtkWriter.c b/EnhancedSolver/2D-seq/src/vtkWriter.c new file mode 100644 index 0000000..521335d --- /dev/null +++ b/EnhancedSolver/2D-seq/src/vtkWriter.c @@ -0,0 +1,125 @@ +/* + * 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 "vtkWriter.h" + +static float floatSwap(float f) +{ + union { + float f; + char b[4]; + } dat1, dat2; + + dat1.f = f; + dat2.b[0] = dat1.b[3]; + dat2.b[1] = dat1.b[2]; + dat2.b[2] = dat1.b[1]; + dat2.b[3] = dat1.b[0]; + return dat2.f; +} + +static void writeHeader(VtkOptions* o, int ts) +{ + fprintf(o->fh, "# vtk DataFile Version 3.0\n"); + fprintf(o->fh, "PAMPI cfd solver particle tracing file\n"); + if (o->fmt == ASCII) { + fprintf(o->fh, "ASCII\n"); + } else if (o->fmt == BINARY) { + fprintf(o->fh, "BINARY\n"); + } + + fprintf(o->fh, "DATASET UNSTRUCTURED_GRID\n"); + fprintf(o->fh, "FIELD FieldData 2\n"); + fprintf(o->fh, "TIME 1 1 double\n"); + fprintf(o->fh, "%d\n", ts); + fprintf(o->fh, "CYCLE 1 1 int\n"); + fprintf(o->fh, "1\n"); +} + +void vtkOpen(VtkOptions* o, char* problem, int ts) +{ + o->fh = fopen(problem, "w"); + + if (o->fh == NULL) { + printf("vtkWriter not initialize! Call vtkOpen first!\n"); + exit(EXIT_FAILURE); + } + + writeHeader(o, ts); + + printf("Writing VTK output for %s\n", problem); +} + +void vtkParticle(VtkOptions* o, char* name) +{ + Particle* particlePool = o->particletracer->particlePool; + + int imax = o->particletracer->grid->imax; + int jmax = o->particletracer->grid->jmax; + + if (o->fh == NULL) { + printf("vtkWriter not initialize! Call vtkOpen first!\n"); + exit(EXIT_FAILURE); + } + + fprintf(o->fh, "POINTS %d float\n", o->particletracer->totalParticles); + + for (int i = 0; i < o->particletracer->totalParticles; ++i) { + double x = particlePool[i].x; + double y = particlePool[i].y; + fprintf(o->fh, "%.2f %.2f 0.0\n", x, y); + } + + fprintf(o->fh, + "CELLS %d %d\n", + o->particletracer->totalParticles, + 2 * o->particletracer->totalParticles); + + for (int i = 0; i < o->particletracer->totalParticles; ++i) { + fprintf(o->fh, "1 %d\n", i); + } + + fprintf(o->fh, "CELL_TYPES %d\n", o->particletracer->totalParticles); + + for (int i = 0; i < o->particletracer->totalParticles; ++i) { + fprintf(o->fh, "1\n"); + } + + /* + for (int k = 0; k < kmax; k++) { + for (int j = 0; j < jmax; j++) { + for (int i = 0; i < imax; i++) { + if (o->fmt == ASCII) { + fprintf(o->fh, + "%f %f %f\n", + G(vec.u, i, j, k), + G(vec.v, i, j, k), + G(vec.w, i, j, k)); + } else if (o->fmt == BINARY) { + fwrite((float[3]) { floatSwap(G(vec.u, i, j, k)), + floatSwap(G(vec.v, i, j, k)), + floatSwap(G(vec.w, i, j, k)) }, + sizeof(float), + 3, + o->fh); + } + } + } + } + if (o->fmt == BINARY) fprintf(o->fh, "\n"); + + */ +} + +void vtkClose(VtkOptions* o) +{ + fclose(o->fh); + o->fh = NULL; +} diff --git a/EnhancedSolver/2D-seq/src/vtkWriter.h b/EnhancedSolver/2D-seq/src/vtkWriter.h new file mode 100644 index 0000000..1e582fc --- /dev/null +++ b/EnhancedSolver/2D-seq/src/vtkWriter.h @@ -0,0 +1,25 @@ +/* + * 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 __VTKWRITER_H_ +#define __VTKWRITER_H_ +#include + +#include "particletracing.h" +#include "solver.h" + +typedef enum VtkFormat { ASCII = 0, BINARY } VtkFormat; + +typedef struct VtkOptions { + VtkFormat fmt; + FILE* fh; + ParticleTracer* particletracer; +} VtkOptions; + +extern void vtkOpen(VtkOptions* opts, char* filename, int ts); +extern void vtkParticle(VtkOptions* opts, char* name); +extern void vtkClose(VtkOptions* opts); +#endif // __VTKWRITER_H_ diff --git a/EnhancedSolver/2D-seq/surface.plot b/EnhancedSolver/2D-seq/surface.plot new file mode 100644 index 0000000..4f7ccd9 --- /dev/null +++ b/EnhancedSolver/2D-seq/surface.plot @@ -0,0 +1,7 @@ +set terminal png size 1024,768 enhanced font ,12 +set output 'p.png' +set datafile separator whitespace + +set grid +set hidden3d +splot 'pressure.dat' using 1:2:3 with lines diff --git a/EnhancedSolver/2D-seq/vector.plot b/EnhancedSolver/2D-seq/vector.plot new file mode 100644 index 0000000..232da08 --- /dev/null +++ b/EnhancedSolver/2D-seq/vector.plot @@ -0,0 +1,7 @@ +set terminal png size 3600,768 enhanced font ,28 +set output 'velocity.png' + +set size ratio -1 +set datafile separator whitespace + +plot 'velocity.dat' using 1:2:3:4:5 with vectors filled head size 0.01,20,60 lc palette diff --git a/EnhancedSolver/2D-seq/vis_files/animate.plot b/EnhancedSolver/2D-seq/vis_files/animate.plot new file mode 100644 index 0000000..d54ba32 --- /dev/null +++ b/EnhancedSolver/2D-seq/vis_files/animate.plot @@ -0,0 +1,10 @@ +unset border; unset tics; unset key; +set term gif animate delay 30 +set output "trace.gif" +set xrange [0:1] +set yrange [0:1] + +do for [ts=0:120] { + plot "particles_".ts.".dat" with points pointtype 7 +} +unset output diff --git a/EnhancedSolver/2D-seq/vis_files/backstep_animate.plot b/EnhancedSolver/2D-seq/vis_files/backstep_animate.plot new file mode 100644 index 0000000..2888e5d --- /dev/null +++ b/EnhancedSolver/2D-seq/vis_files/backstep_animate.plot @@ -0,0 +1,14 @@ +unset border; unset tics; unset key; +set term gif animate delay 10 +set output "trace.gif" +set xrange [0:7] +set yrange [0:1.5] +set size ratio -1 + +set object 1 rect from 0.0,0.0 to 1.0,0.5 lw 5 + + +do for [ts=0:300] { + plot "particles_".ts.".dat" with points pointtype 7 +} +unset output diff --git a/EnhancedSolver/2D-seq/vis_files/karman_animate.plot b/EnhancedSolver/2D-seq/vis_files/karman_animate.plot new file mode 100644 index 0000000..d78f0c3 --- /dev/null +++ b/EnhancedSolver/2D-seq/vis_files/karman_animate.plot @@ -0,0 +1,13 @@ +unset border; unset tics; unset key; +set term gif animate delay 10 +set output "trace.gif" +set xrange [0:30] +set yrange [0:8] +set size ratio -1 +set object 1 circle front at 5.0,4.0 size 1.0 fillcolor rgb "black" lw 2 + + +do for [ts=0:500] { + plot "particles_".ts.".dat" with points pointtype 7 pointsize 0.3 +} +unset output