From 91cbb29f6e6c3b1bbeb8f2074ae1a709257096a3 Mon Sep 17 00:00:00 2001 From: Jan Eitzinger Date: Wed, 31 Jan 2024 17:40:36 +0100 Subject: [PATCH] Merge 3D mpi versions --- BasicSolver/3D-mpi-io/Makefile | 71 -- BasicSolver/3D-mpi-io/README.md | 78 -- BasicSolver/3D-mpi-io/canal.par | 52 -- BasicSolver/3D-mpi-io/config.mk | 13 - BasicSolver/3D-mpi-io/dcavity.par | 52 -- BasicSolver/3D-mpi-io/include_CLANG.mk | 21 - BasicSolver/3D-mpi-io/include_GCC.mk | 14 - BasicSolver/3D-mpi-io/include_ICC.mk | 14 - BasicSolver/3D-mpi-io/src/allocate.c | 38 - BasicSolver/3D-mpi-io/src/allocate.h | 13 - BasicSolver/3D-mpi-io/src/comm.c | 359 -------- BasicSolver/3D-mpi-io/src/comm.h | 51 -- BasicSolver/3D-mpi-io/src/grid.h | 16 - BasicSolver/3D-mpi-io/src/likwid-marker.h | 54 -- BasicSolver/3D-mpi-io/src/main.c | 84 -- BasicSolver/3D-mpi-io/src/parameter.c | 126 --- BasicSolver/3D-mpi-io/src/parameter.h | 26 - BasicSolver/3D-mpi-io/src/progress.c | 51 -- BasicSolver/3D-mpi-io/src/progress.h | 14 - BasicSolver/3D-mpi-io/src/solver.c | 853 ------------------ BasicSolver/3D-mpi-io/src/solver.h | 45 - BasicSolver/3D-mpi-io/src/test.c | 129 --- BasicSolver/3D-mpi-io/src/test.h | 13 - BasicSolver/3D-mpi-io/src/timing.c | 22 - BasicSolver/3D-mpi-io/src/timing.h | 13 - BasicSolver/3D-mpi-io/src/util.h | 22 - BasicSolver/3D-mpi-io/src/vtkWriter.h | 29 - BasicSolver/3D-mpi/Makefile | 19 +- BasicSolver/3D-mpi/config.mk | 2 + BasicSolver/3D-mpi/include_CLANG.mk | 4 +- BasicSolver/3D-mpi/include_GCC.mk | 8 +- BasicSolver/3D-mpi/include_ICC.mk | 8 +- BasicSolver/3D-mpi/src/allocate.c | 2 +- BasicSolver/3D-mpi/src/allocate.h | 2 +- BasicSolver/3D-mpi/src/comm.c | 38 +- BasicSolver/3D-mpi/src/comm.h | 3 +- BasicSolver/3D-mpi/src/grid.h | 2 +- BasicSolver/3D-mpi/src/main.c | 12 +- BasicSolver/3D-mpi/src/parameter.c | 2 +- BasicSolver/3D-mpi/src/parameter.h | 2 +- BasicSolver/3D-mpi/src/progress.c | 2 +- BasicSolver/3D-mpi/src/progress.h | 2 +- BasicSolver/3D-mpi/src/solver.c | 22 +- BasicSolver/3D-mpi/src/solver.h | 2 +- BasicSolver/3D-mpi/src/test.c | 21 +- BasicSolver/3D-mpi/src/test.h | 2 +- BasicSolver/3D-mpi/src/timing.c | 2 +- BasicSolver/3D-mpi/src/timing.h | 2 +- BasicSolver/3D-mpi/src/util.h | 2 +- .../src/vtkWriter-mpi.c} | 47 +- .../src/{vtkWriter.c => vtkWriter-seq.c} | 2 +- BasicSolver/3D-mpi/src/vtkWriter.h | 10 +- 52 files changed, 149 insertions(+), 2344 deletions(-) delete mode 100644 BasicSolver/3D-mpi-io/Makefile delete mode 100644 BasicSolver/3D-mpi-io/README.md delete mode 100644 BasicSolver/3D-mpi-io/canal.par delete mode 100644 BasicSolver/3D-mpi-io/config.mk delete mode 100644 BasicSolver/3D-mpi-io/dcavity.par delete mode 100644 BasicSolver/3D-mpi-io/include_CLANG.mk delete mode 100644 BasicSolver/3D-mpi-io/include_GCC.mk delete mode 100644 BasicSolver/3D-mpi-io/include_ICC.mk delete mode 100644 BasicSolver/3D-mpi-io/src/allocate.c delete mode 100644 BasicSolver/3D-mpi-io/src/allocate.h delete mode 100644 BasicSolver/3D-mpi-io/src/comm.c delete mode 100644 BasicSolver/3D-mpi-io/src/comm.h delete mode 100644 BasicSolver/3D-mpi-io/src/grid.h delete mode 100644 BasicSolver/3D-mpi-io/src/likwid-marker.h delete mode 100644 BasicSolver/3D-mpi-io/src/main.c delete mode 100644 BasicSolver/3D-mpi-io/src/parameter.c delete mode 100644 BasicSolver/3D-mpi-io/src/parameter.h delete mode 100644 BasicSolver/3D-mpi-io/src/progress.c delete mode 100644 BasicSolver/3D-mpi-io/src/progress.h delete mode 100644 BasicSolver/3D-mpi-io/src/solver.c delete mode 100644 BasicSolver/3D-mpi-io/src/solver.h delete mode 100644 BasicSolver/3D-mpi-io/src/test.c delete mode 100644 BasicSolver/3D-mpi-io/src/test.h delete mode 100644 BasicSolver/3D-mpi-io/src/timing.c delete mode 100644 BasicSolver/3D-mpi-io/src/timing.h delete mode 100644 BasicSolver/3D-mpi-io/src/util.h delete mode 100644 BasicSolver/3D-mpi-io/src/vtkWriter.h rename BasicSolver/{3D-mpi-io/src/vtkWriter.c => 3D-mpi/src/vtkWriter-mpi.c} (84%) rename BasicSolver/3D-mpi/src/{vtkWriter.c => vtkWriter-seq.c} (98%) diff --git a/BasicSolver/3D-mpi-io/Makefile b/BasicSolver/3D-mpi-io/Makefile deleted file mode 100644 index 57f99f4..0000000 --- a/BasicSolver/3D-mpi-io/Makefile +++ /dev/null @@ -1,71 +0,0 @@ -#======================================================================================= -# Copyright (C) 2022 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) - -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/BasicSolver/3D-mpi-io/README.md b/BasicSolver/3D-mpi-io/README.md deleted file mode 100644 index d980b54..0000000 --- a/BasicSolver/3D-mpi-io/README.md +++ /dev/null @@ -1,78 +0,0 @@ -# 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/BasicSolver/3D-mpi-io/canal.par b/BasicSolver/3D-mpi-io/canal.par deleted file mode 100644 index 3ff4a5f..0000000 --- a/BasicSolver/3D-mpi-io/canal.par +++ /dev/null @@ -1,52 +0,0 @@ -#============================================================================== -# Laminar Canal Flow -#============================================================================== - -# Problem specific Data: -# --------------------- - -name canal # name of flow setup - -bcLeft 3 # flags for boundary conditions -bcRight 3 # 1 = no-slip 3 = outflow -bcBottom 1 # 2 = free-slip 4 = periodic -bcTop 1 # -bcFront 1 # -bcBack 1 # - -gx 0.0 # Body forces (e.g. gravity) -gy 0.0 # -gz 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 -w_init 0.0 # initial value for velocity in z-direction -p_init 0.0 # initial value for pressure - -# Geometry Data: -# ------------- - -xlength 30.0 # domain size in x-direction -ylength 4.0 # domain size in y-direction -zlength 4.0 # domain size in z-direction -imax 200 # number of interior cells in x-direction -jmax 50 # number of interior cells in y-direction -kmax 50 # number of interior cells in z-direction - -# Time Data: -# --------- - -te 100.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 -omg 1.3 # relaxation parameter for SOR iteration -gamma 0.9 # upwind differencing factor gamma -#=============================================================================== diff --git a/BasicSolver/3D-mpi-io/config.mk b/BasicSolver/3D-mpi-io/config.mk deleted file mode 100644 index c163666..0000000 --- a/BasicSolver/3D-mpi-io/config.mk +++ /dev/null @@ -1,13 +0,0 @@ -# Supported: GCC, CLANG, ICC -TAG ?= CLANG -ENABLE_MPI ?= true -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/BasicSolver/3D-mpi-io/dcavity.par b/BasicSolver/3D-mpi-io/dcavity.par deleted file mode 100644 index c5cb201..0000000 --- a/BasicSolver/3D-mpi-io/dcavity.par +++ /dev/null @@ -1,52 +0,0 @@ -#============================================================================== -# Driven Cavity -#============================================================================== - -# Problem specific Data: -# --------------------- - -name dcavity # name of flow setup - -bcLeft 1 # flags for boundary conditions -bcRight 1 # 1 = no-slip 3 = outflow -bcBottom 1 # 2 = free-slip 4 = periodic -bcTop 1 # -bcFront 1 # -bcBack 1 # - -gx 0.0 # Body forces (e.g. gravity) -gy 0.0 # -gz 0.0 # - -re 1000.0 # Reynolds number - -u_init 0.0 # initial value for velocity in x-direction -v_init 0.0 # initial value for velocity in y-direction -w_init 0.0 # initial value for velocity in z-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 -zlength 1.0 # domain size in z-direction -imax 128 # number of interior cells in x-direction -jmax 128 # number of interior cells in y-direction -kmax 128 # number of interior cells in z-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 -omg 1.8 # relaxation parameter for SOR iteration -gamma 0.9 # upwind differencing factor gamma -#=============================================================================== diff --git a/BasicSolver/3D-mpi-io/include_CLANG.mk b/BasicSolver/3D-mpi-io/include_CLANG.mk deleted file mode 100644 index a1e37e0..0000000 --- a/BasicSolver/3D-mpi-io/include_CLANG.mk +++ /dev/null @@ -1,21 +0,0 @@ -ifeq ($(ENABLE_MPI),true) -CC = mpicc -DEFINES = -D_MPI -else -CC = cc -endif - -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 = -Ofast -std=c17 -LFLAGS = $(OPENMP) -lm -DEFINES += -D_GNU_SOURCE# -DDEBUG -INCLUDES = diff --git a/BasicSolver/3D-mpi-io/include_GCC.mk b/BasicSolver/3D-mpi-io/include_GCC.mk deleted file mode 100644 index 427e798..0000000 --- a/BasicSolver/3D-mpi-io/include_GCC.mk +++ /dev/null @@ -1,14 +0,0 @@ -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/BasicSolver/3D-mpi-io/include_ICC.mk b/BasicSolver/3D-mpi-io/include_ICC.mk deleted file mode 100644 index 94b8e20..0000000 --- a/BasicSolver/3D-mpi-io/include_ICC.mk +++ /dev/null @@ -1,14 +0,0 @@ -CC = icc -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/BasicSolver/3D-mpi-io/src/allocate.c b/BasicSolver/3D-mpi-io/src/allocate.c deleted file mode 100644 index 6af6c7f..0000000 --- a/BasicSolver/3D-mpi-io/src/allocate.c +++ /dev/null @@ -1,38 +0,0 @@ -/* - * Copyright (C) 2022 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 - -#include "allocate.h" - -void* allocate(size_t 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/BasicSolver/3D-mpi-io/src/allocate.h b/BasicSolver/3D-mpi-io/src/allocate.h deleted file mode 100644 index 1537eb2..0000000 --- a/BasicSolver/3D-mpi-io/src/allocate.h +++ /dev/null @@ -1,13 +0,0 @@ -/* - * Copyright (C) 2022 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(size_t alignment, size_t bytesize); - -#endif diff --git a/BasicSolver/3D-mpi-io/src/comm.c b/BasicSolver/3D-mpi-io/src/comm.c deleted file mode 100644 index c03688d..0000000 --- a/BasicSolver/3D-mpi-io/src/comm.c +++ /dev/null @@ -1,359 +0,0 @@ -/* - * Copyright (C) 2022 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. - */ -#if defined(_MPI) -#include -#endif -#include -#include -#include - -#include "allocate.h" -#include "comm.h" - -#if defined(_MPI) -// subroutines local to this module -static int sizeOfRank(int rank, int size, int N) -{ - return N / size + ((N % size > rank) ? 1 : 0); -} - -static void setupCommunication(Comm* c, Direction direction, int layer) -{ - int imaxLocal = c->imaxLocal; - int jmaxLocal = c->jmaxLocal; - int kmaxLocal = c->kmaxLocal; - - size_t dblsize = sizeof(double); - int sizes[NDIMS]; - int subSizes[NDIMS]; - int starts[NDIMS]; - int offset = 0; - - sizes[IDIM] = imaxLocal + 2; - sizes[JDIM] = jmaxLocal + 2; - sizes[KDIM] = kmaxLocal + 2; - - if (layer == HALO) { - offset = 1; - } - - switch (direction) { - case LEFT: - subSizes[IDIM] = 1; - subSizes[JDIM] = jmaxLocal; - subSizes[KDIM] = kmaxLocal; - starts[IDIM] = 1 - offset; - starts[JDIM] = 1; - starts[KDIM] = 1; - break; - case RIGHT: - subSizes[IDIM] = 1; - subSizes[JDIM] = jmaxLocal; - subSizes[KDIM] = kmaxLocal; - starts[IDIM] = imaxLocal + offset; - starts[JDIM] = 1; - starts[KDIM] = 1; - break; - case BOTTOM: - subSizes[IDIM] = imaxLocal; - subSizes[JDIM] = 1; - subSizes[KDIM] = kmaxLocal; - starts[IDIM] = 1; - starts[JDIM] = 1 - offset; - starts[KDIM] = 1; - break; - case TOP: - subSizes[IDIM] = imaxLocal; - subSizes[JDIM] = 1; - subSizes[KDIM] = kmaxLocal; - starts[IDIM] = 1; - starts[JDIM] = jmaxLocal + offset; - starts[KDIM] = 1; - break; - case FRONT: - subSizes[IDIM] = imaxLocal; - subSizes[JDIM] = jmaxLocal; - subSizes[KDIM] = 1; - starts[IDIM] = 1; - starts[JDIM] = 1; - starts[KDIM] = 1 - offset; - break; - case BACK: - subSizes[IDIM] = imaxLocal; - subSizes[JDIM] = jmaxLocal; - subSizes[KDIM] = 1; - starts[IDIM] = 1; - starts[JDIM] = 1; - starts[KDIM] = kmaxLocal + offset; - break; - case NDIRS: - printf("ERROR!\n"); - break; - } - - if (layer == HALO) { - MPI_Type_create_subarray(NDIMS, - sizes, - subSizes, - starts, - MPI_ORDER_C, - MPI_DOUBLE, - &c->rbufferTypes[direction]); - MPI_Type_commit(&c->rbufferTypes[direction]); - } else if (layer == BULK) { - MPI_Type_create_subarray(NDIMS, - sizes, - subSizes, - starts, - MPI_ORDER_C, - MPI_DOUBLE, - &c->sbufferTypes[direction]); - MPI_Type_commit(&c->sbufferTypes[direction]); - } -} - -static int sum(int* sizes, int position) -{ - int sum = 0; - - for (int i = 0; i < position; i++) { - sum += sizes[i]; - } - - return sum; -} -#endif - -// exported subroutines -void commReduction(double* v, int op) -{ -#if defined(_MPI) - if (op == MAX) { - MPI_Allreduce(MPI_IN_PLACE, v, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); - } else if (op == SUM) { - MPI_Allreduce(MPI_IN_PLACE, v, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); - } -#endif -} - -int commIsBoundary(Comm* c, Direction direction) -{ -#if defined(_MPI) - switch (direction) { - case LEFT: - return c->coords[ICORD] == 0; - break; - case RIGHT: - return c->coords[ICORD] == (c->dims[ICORD] - 1); - break; - case BOTTOM: - return c->coords[JCORD] == 0; - break; - case TOP: - return c->coords[JCORD] == (c->dims[JCORD] - 1); - break; - case FRONT: - return c->coords[KCORD] == 0; - break; - case BACK: - return c->coords[KCORD] == (c->dims[KCORD] - 1); - break; - case NDIRS: - printf("ERROR!\n"); - break; - } -#endif - - return 1; -} - -void commExchange(Comm* c, double* grid) -{ -#if defined(_MPI) - int counts[6] = { 1, 1, 1, 1, 1, 1 }; - MPI_Aint displs[6] = { 0, 0, 0, 0, 0, 0 }; - - MPI_Neighbor_alltoallw(grid, - counts, - displs, - c->sbufferTypes, - grid, - counts, - displs, - c->rbufferTypes, - c->comm); -#endif -} - -void commShift(Comm* c, double* f, double* g, double* h) -{ -#if defined(_MPI) - MPI_Request requests[6] = { MPI_REQUEST_NULL, - MPI_REQUEST_NULL, - MPI_REQUEST_NULL, - MPI_REQUEST_NULL, - MPI_REQUEST_NULL, - MPI_REQUEST_NULL }; - - /* shift G */ - /* receive ghost cells from bottom neighbor */ - MPI_Irecv(g, - 1, - c->rbufferTypes[BOTTOM], - c->neighbours[BOTTOM], - 0, - c->comm, - &requests[0]); - - /* send ghost cells to top neighbor */ - MPI_Isend(g, 1, c->sbufferTypes[TOP], c->neighbours[TOP], 0, c->comm, &requests[1]); - - /* shift F */ - /* receive ghost cells from left neighbor */ - MPI_Irecv(f, 1, c->rbufferTypes[LEFT], c->neighbours[LEFT], 1, c->comm, &requests[2]); - - /* send ghost cells to right neighbor */ - MPI_Isend(f, - 1, - c->sbufferTypes[RIGHT], - c->neighbours[RIGHT], - 1, - c->comm, - &requests[3]); - - /* shift H */ - /* receive ghost cells from front neighbor */ - MPI_Irecv(h, - 1, - c->rbufferTypes[FRONT], - c->neighbours[FRONT], - 2, - c->comm, - &requests[4]); - - /* send ghost cells to back neighbor */ - MPI_Isend(h, 1, c->sbufferTypes[BACK], c->neighbours[BACK], 2, c->comm, &requests[5]); - - MPI_Waitall(6, requests, MPI_STATUSES_IGNORE); -#endif -} - -void commGetOffsets(Comm* c, int offsets[], int kmax, int jmax, int imax) -{ -#if defined(_MPI) - int sum = 0; - - for (int i = 0; i < c->coords[ICORD]; i++) { - sum += sizeOfRank(i, c->dims[ICORD], imax); - } - offsets[IDIM] = sum; - sum = 0; - - for (int i = 0; i < c->coords[JCORD]; i++) { - sum += sizeOfRank(i, c->dims[JCORD], jmax); - } - offsets[JDIM] = sum; - sum = 0; - - for (int i = 0; i < c->coords[KCORD]; i++) { - sum += sizeOfRank(i, c->dims[KCORD], kmax); - } - offsets[KDIM] = sum; -#endif -} - -void commPrintConfig(Comm* c) -{ -#if defined(_MPI) - fflush(stdout); - MPI_Barrier(MPI_COMM_WORLD); - if (commIsMaster(c)) { - printf("Communication setup:\n"); - } - - for (int i = 0; i < c->size; i++) { - if (i == c->rank) { - printf("\tRank %d of %d\n", c->rank, c->size); - printf("\tNeighbours (front, back, bottom, top, left, right): %d, %d, %d, " - "%d, %d, %d\n", - c->neighbours[FRONT], - c->neighbours[BACK], - c->neighbours[BOTTOM], - c->neighbours[TOP], - c->neighbours[LEFT], - c->neighbours[RIGHT]); - printf("\tCoordinates (k,j,i) %d %d %d\n", - c->coords[KCORD], - c->coords[JCORD], - c->coords[ICORD]); - printf("\tLocal domain size (k,j,i) %dx%dx%d\n", - c->kmaxLocal, - c->jmaxLocal, - c->imaxLocal); - fflush(stdout); - } - } - MPI_Barrier(MPI_COMM_WORLD); -#endif -} - -void commInit(Comm* c, int argc, char** argv) -{ -#if defined(_MPI) - MPI_Init(&argc, &argv); - MPI_Comm_rank(MPI_COMM_WORLD, &(c->rank)); - MPI_Comm_size(MPI_COMM_WORLD, &(c->size)); -#endif -} - -void commPartition(Comm* c, int kmax, int jmax, int imax) -{ -#if defined(_MPI) - int dims[NDIMS] = { 0, 0, 0 }; - int periods[NDIMS] = { 0, 0, 0 }; - MPI_Dims_create(c->size, NDIMS, dims); - MPI_Cart_create(MPI_COMM_WORLD, NCORDS, dims, periods, 0, &c->comm); - MPI_Cart_shift(c->comm, ICORD, 1, &c->neighbours[LEFT], &c->neighbours[RIGHT]); - MPI_Cart_shift(c->comm, JCORD, 1, &c->neighbours[BOTTOM], &c->neighbours[TOP]); - MPI_Cart_shift(c->comm, KCORD, 1, &c->neighbours[FRONT], &c->neighbours[BACK]); - MPI_Cart_get(c->comm, NCORDS, c->dims, periods, c->coords); - - c->imaxLocal = sizeOfRank(c->rank, dims[ICORD], imax); - c->jmaxLocal = sizeOfRank(c->rank, dims[JCORD], jmax); - c->kmaxLocal = sizeOfRank(c->rank, dims[KCORD], kmax); - - // setup buffer types for communication - setupCommunication(c, LEFT, BULK); - setupCommunication(c, LEFT, HALO); - setupCommunication(c, RIGHT, BULK); - setupCommunication(c, RIGHT, HALO); - setupCommunication(c, BOTTOM, BULK); - setupCommunication(c, BOTTOM, HALO); - setupCommunication(c, TOP, BULK); - setupCommunication(c, TOP, HALO); - setupCommunication(c, FRONT, BULK); - setupCommunication(c, FRONT, HALO); - setupCommunication(c, BACK, BULK); - setupCommunication(c, BACK, HALO); -#else - c->imaxLocal = imax; - c->jmaxLocal = jmax; - c->kmaxLocal = kmax; -#endif -} - -void commFinalize(Comm* c) -{ -#if defined(_MPI) - for (int i = 0; i < NDIRS; i++) { - MPI_Type_free(&c->sbufferTypes[i]); - MPI_Type_free(&c->rbufferTypes[i]); - } - - MPI_Finalize(); -#endif -} diff --git a/BasicSolver/3D-mpi-io/src/comm.h b/BasicSolver/3D-mpi-io/src/comm.h deleted file mode 100644 index 017555c..0000000 --- a/BasicSolver/3D-mpi-io/src/comm.h +++ /dev/null @@ -1,51 +0,0 @@ -/* - * Copyright (C) 2024 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 __COMM_H_ -#define __COMM_H_ -#if defined(_MPI) -#include -#endif -/* - * Spatial directions: - * ICORD (0) from 0 (LEFT) to imax (RIGHT) - * JCORD (1) from 0 (BOTTOM) to jmax (TOP) - * KCORD (2) from 0 (FRONT) to kmax (BACK) - * All derived Subarray types are in C ordering - * with indices KDIM (0), JDIM(1) and IDIM(2) - * */ -typedef enum direction { LEFT = 0, RIGHT, BOTTOM, TOP, FRONT, BACK, NDIRS } Direction; -typedef enum coordinates { ICORD = 0, JCORD, KCORD, NCORDS } Coordinates; -typedef enum dimension { KDIM = 0, JDIM, IDIM, NDIMS } Dimension; -enum layer { HALO = 0, BULK }; -enum op { MAX = 0, SUM }; - -typedef struct { - int rank; - int size; -#if defined(_MPI) - MPI_Comm comm; - MPI_Datatype sbufferTypes[NDIRS]; - MPI_Datatype rbufferTypes[NDIRS]; -#endif - int neighbours[NDIRS]; - int coords[NDIMS], dims[NDIMS]; - int imaxLocal, jmaxLocal, kmaxLocal; - MPI_File fh; -} Comm; - -extern void commInit(Comm* c, int argc, char** argv); -extern void commPartition(Comm* c, int kmax, int jmax, int imax); -extern void commFinalize(Comm* comm); -extern void commPrintConfig(Comm*); -extern void commExchange(Comm*, double*); -extern void commShift(Comm* c, double* f, double* g, double* h); -extern void commReduction(double* v, int op); -extern int commIsBoundary(Comm* c, Direction direction); -extern void commGetOffsets(Comm* c, int offsets[], int kmax, int jmax, int imax); - -static inline int commIsMaster(Comm* c) { return c->rank == 0; } -#endif // __COMM_H_ diff --git a/BasicSolver/3D-mpi-io/src/grid.h b/BasicSolver/3D-mpi-io/src/grid.h deleted file mode 100644 index c963429..0000000 --- a/BasicSolver/3D-mpi-io/src/grid.h +++ /dev/null @@ -1,16 +0,0 @@ -/* - * Copyright (C) 2022 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, dz; - int imax, jmax, kmax; - double xlength, ylength, zlength; -} Grid; - -#endif // __GRID_H_ diff --git a/BasicSolver/3D-mpi-io/src/likwid-marker.h b/BasicSolver/3D-mpi-io/src/likwid-marker.h deleted file mode 100644 index c3770c0..0000000 --- a/BasicSolver/3D-mpi-io/src/likwid-marker.h +++ /dev/null @@ -1,54 +0,0 @@ -/* - * ======================================================================================= - * - * 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/BasicSolver/3D-mpi-io/src/main.c b/BasicSolver/3D-mpi-io/src/main.c deleted file mode 100644 index 8a1c88a..0000000 --- a/BasicSolver/3D-mpi-io/src/main.c +++ /dev/null @@ -1,84 +0,0 @@ -/* - * Copyright (C) 2024 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 -#include - -#include "allocate.h" -#include "comm.h" -#include "parameter.h" -#include "progress.h" -#include "solver.h" -#include "timing.h" -#include "vtkWriter.h" - -int main(int argc, char** argv) -{ - double timeStart, timeStop; - Parameter p; - Solver s; - - commInit(&s.comm, argc, argv); - initParameter(&p); - - if (argc != 2) { - printf("Usage: %s \n", argv[0]); - exit(EXIT_SUCCESS); - } - - readParameter(&p, argv[1]); - commPartition(&s.comm, p.kmax, p.jmax, p.imax); - if (commIsMaster(&s.comm)) { - printParameter(&p); - } - initSolver(&s, &p); -#ifndef VERBOSE - initProgress(s.te); -#endif - - double tau = s.tau; - double te = s.te; - double t = 0.0; - - timeStart = getTimeStamp(); - while (t <= te) { - if (tau > 0.0) computeTimestep(&s); - setBoundaryConditions(&s); - setSpecialBoundaryCondition(&s); - computeFG(&s); - computeRHS(&s); - solve(&s); - adaptUV(&s); - t += s.dt; - -#ifdef VERBOSE - if (commIsMaster(&s.comm)) { - printf("TIME %f , TIMESTEP %f\n", t, s.dt); - } -#else - printProgress(t); -#endif - } - timeStop = getTimeStamp(); -#ifndef VERBOSE - stopProgress(); -#endif - if (commIsMaster(&s.comm)) { - printf("Solution took %.2fs\n", timeStop - timeStart); - } - - VtkOptions opts = { .grid = s.grid, .comm = s.comm }; - vtkOpen(&opts, s.problem); - vtkScalar(&opts, "pressure", s.p); - vtkVector(&opts, "velocity", (VtkVector) { s.u, s.v, s.w }); - vtkClose(&opts); - - commFinalize(&s.comm); - return EXIT_SUCCESS; -} diff --git a/BasicSolver/3D-mpi-io/src/parameter.c b/BasicSolver/3D-mpi-io/src/parameter.c deleted file mode 100644 index d4ea5b6..0000000 --- a/BasicSolver/3D-mpi-io/src/parameter.c +++ /dev/null @@ -1,126 +0,0 @@ -/* - * Copyright (C) 2022 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->zlength = 1.0; - param->imax = 100; - param->jmax = 100; - param->kmax = 100; - param->itermax = 1000; - param->eps = 0.0001; - param->omg = 1.7; - param->re = 100.0; - param->gamma = 0.9; - param->tau = 0.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_REAL(zlength); - PARSE_INT(imax); - PARSE_INT(jmax); - PARSE_INT(kmax); - PARSE_INT(itermax); - 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_REAL(gz); - PARSE_STRING(name); - PARSE_INT(bcLeft); - PARSE_INT(bcRight); - PARSE_INT(bcBottom); - PARSE_INT(bcTop); - PARSE_INT(bcFront); - PARSE_INT(bcBack); - PARSE_REAL(u_init); - PARSE_REAL(v_init); - PARSE_REAL(w_init); - PARSE_REAL(p_init); - } - } - - fclose(fp); -} - -void printParameter(Parameter* param) -{ - printf("Parameters for %s\n", param->name); - printf("Boundary conditions Left:%d Right:%d Bottom:%d Top:%d Front:%d " - "Back:%d\n", - param->bcLeft, - param->bcRight, - param->bcBottom, - param->bcTop, - param->bcFront, - param->bcBack); - printf("\tReynolds number: %.2f\n", param->re); - printf("\tInit arrays: U:%.2f V:%.2f W:%.2f P:%.2f\n", - param->u_init, - param->v_init, - param->w_init, - param->p_init); - printf("Geometry data:\n"); - printf("\tDomain box size (x, y, z): %.2f, %.2f, %.2f\n", - param->xlength, - param->ylength, - param->zlength); - printf("\tCells (x, y, z): %d, %d, %d\n", param->imax, param->jmax, param->kmax); - 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); -} diff --git a/BasicSolver/3D-mpi-io/src/parameter.h b/BasicSolver/3D-mpi-io/src/parameter.h deleted file mode 100644 index 6dddf5f..0000000 --- a/BasicSolver/3D-mpi-io/src/parameter.h +++ /dev/null @@ -1,26 +0,0 @@ -/* - * Copyright (C) 2022 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 { - int imax, jmax, kmax; - double xlength, ylength, zlength; - int itermax; - double eps, omg; - double re, tau, gamma; - double te, dt; - double gx, gy, gz; - char* name; - int bcLeft, bcRight, bcBottom, bcTop, bcFront, bcBack; - double u_init, v_init, w_init, p_init; -} Parameter; - -void initParameter(Parameter*); -void readParameter(Parameter*, const char*); -void printParameter(Parameter*); -#endif diff --git a/BasicSolver/3D-mpi-io/src/progress.c b/BasicSolver/3D-mpi-io/src/progress.c deleted file mode 100644 index a9b82bd..0000000 --- a/BasicSolver/3D-mpi-io/src/progress.c +++ /dev/null @@ -1,51 +0,0 @@ -/* - * Copyright (C) 2022 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/BasicSolver/3D-mpi-io/src/progress.h b/BasicSolver/3D-mpi-io/src/progress.h deleted file mode 100644 index 4d02cdb..0000000 --- a/BasicSolver/3D-mpi-io/src/progress.h +++ /dev/null @@ -1,14 +0,0 @@ -/* - * Copyright (C) 2022 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(void); - -#endif diff --git a/BasicSolver/3D-mpi-io/src/solver.c b/BasicSolver/3D-mpi-io/src/solver.c deleted file mode 100644 index 1ba1042..0000000 --- a/BasicSolver/3D-mpi-io/src/solver.c +++ /dev/null @@ -1,853 +0,0 @@ -/* - * Copyright (C) 2022 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 "comm.h" -#include "parameter.h" -#include "solver.h" -#include "util.h" - -#define P(i, j, k) \ - p[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)] -#define F(i, j, k) \ - f[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)] -#define G(i, j, k) \ - g[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)] -#define H(i, j, k) \ - h[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)] -#define U(i, j, k) \ - u[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)] -#define V(i, j, k) \ - v[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)] -#define W(i, j, k) \ - w[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)] -#define RHS(i, j, k) \ - rhs[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)] - -static void printConfig(Solver* s) -{ - if (commIsMaster(&s->comm)) { - printf("Parameters for #%s#\n", s->problem); - printf("BC Left:%d Right:%d Bottom:%d Top:%d Front:%d Back:%d\n", - s->bcLeft, - s->bcRight, - s->bcBottom, - s->bcTop, - s->bcFront, - s->bcBack); - printf("\tReynolds number: %.2f\n", s->re); - printf("\tGx Gy: %.2f %.2f %.2f\n", s->gx, s->gy, s->gz); - printf("Geometry data:\n"); - printf("\tDomain box size (x, y, z): %.2f, %.2f, %.2f\n", - s->grid.xlength, - s->grid.ylength, - s->grid.zlength); - printf("\tCells (x, y, z): %d, %d, %d\n", - s->grid.imax, - s->grid.jmax, - s->grid.kmax); - printf("\tCell size (dx, dy, dz): %f, %f, %f\n", - s->grid.dx, - s->grid.dy, - s->grid.dz); - printf("Timestep parameters:\n"); - printf("\tDefault stepsize: %.2f, Final time %.2f\n", s->dt, s->te); - printf("\tdt bound: %.6f\n", s->dtBound); - printf("\tTau factor: %.2f\n", s->tau); - printf("Iterative parameters:\n"); - printf("\tMax iterations: %d\n", s->itermax); - printf("\tepsilon (stopping tolerance) : %f\n", s->eps); - printf("\tgamma factor: %f\n", s->gamma); - printf("\tomega (SOR relaxation): %f\n", s->omega); - } - commPrintConfig(&s->comm); -} - -void initSolver(Solver* s, Parameter* params) -{ - s->problem = params->name; - s->bcLeft = params->bcLeft; - s->bcRight = params->bcRight; - s->bcBottom = params->bcBottom; - s->bcTop = params->bcTop; - s->bcFront = params->bcFront; - s->bcBack = params->bcBack; - - s->grid.imax = params->imax; - s->grid.jmax = params->jmax; - s->grid.kmax = params->kmax; - s->grid.xlength = params->xlength; - s->grid.ylength = params->ylength; - s->grid.zlength = params->zlength; - s->grid.dx = params->xlength / params->imax; - s->grid.dy = params->ylength / params->jmax; - s->grid.dz = params->zlength / params->kmax; - - s->eps = params->eps; - s->omega = params->omg; - s->itermax = params->itermax; - s->re = params->re; - s->gx = params->gx; - s->gy = params->gy; - s->gz = params->gz; - s->dt = params->dt; - s->te = params->te; - s->tau = params->tau; - s->gamma = params->gamma; - - /* allocate arrays */ - int imaxLocal = s->comm.imaxLocal; - int jmaxLocal = s->comm.jmaxLocal; - int kmaxLocal = s->comm.kmaxLocal; - size_t size = (imaxLocal + 2) * (jmaxLocal + 2) * (kmaxLocal + 2); - - s->u = allocate(64, size * sizeof(double)); - s->v = allocate(64, size * sizeof(double)); - s->w = allocate(64, size * sizeof(double)); - s->p = allocate(64, size * sizeof(double)); - s->rhs = allocate(64, size * sizeof(double)); - s->f = allocate(64, size * sizeof(double)); - s->g = allocate(64, size * sizeof(double)); - s->h = allocate(64, size * sizeof(double)); - - for (int i = 0; i < size; i++) { - s->u[i] = params->u_init; - s->v[i] = params->v_init; - s->w[i] = params->w_init; - s->p[i] = params->p_init; - s->rhs[i] = 0.0; - s->f[i] = 0.0; - s->g[i] = 0.0; - s->h[i] = 0.0; - } - - double dx = s->grid.dx; - double dy = s->grid.dy; - double dz = s->grid.dz; - - double invSqrSum = 1.0 / (dx * dx) + 1.0 / (dy * dy) + 1.0 / (dz * dz); - s->dtBound = 0.5 * s->re * 1.0 / invSqrSum; - -#ifdef VERBOSE - printConfig(s); -#endif /* VERBOSE */ -} - -void computeRHS(Solver* s) -{ - int imaxLocal = s->comm.imaxLocal; - int jmaxLocal = s->comm.jmaxLocal; - int kmaxLocal = s->comm.kmaxLocal; - - double idx = 1.0 / s->grid.dx; - double idy = 1.0 / s->grid.dy; - double idz = 1.0 / s->grid.dz; - double idt = 1.0 / s->dt; - - double* rhs = s->rhs; - double* f = s->f; - double* g = s->g; - double* h = s->h; - - commShift(&s->comm, f, g, h); - - for (int k = 1; k < kmaxLocal + 1; k++) { - for (int j = 1; j < jmaxLocal + 1; j++) { - for (int i = 1; i < imaxLocal + 1; i++) { - RHS(i, j, k) = ((F(i, j, k) - F(i - 1, j, k)) * idx + - (G(i, j, k) - G(i, j - 1, k)) * idy + - (H(i, j, k) - H(i, j, k - 1)) * idz) * - idt; - } - } - } -} - -void solve(Solver* s) -{ - int imaxLocal = s->comm.imaxLocal; - int jmaxLocal = s->comm.jmaxLocal; - int kmaxLocal = s->comm.kmaxLocal; - - int imax = s->grid.imax; - int jmax = s->grid.jmax; - int kmax = s->grid.kmax; - - double eps = s->eps; - int itermax = s->itermax; - double dx2 = s->grid.dx * s->grid.dx; - double dy2 = s->grid.dy * s->grid.dy; - double dz2 = s->grid.dz * s->grid.dz; - double idx2 = 1.0 / dx2; - double idy2 = 1.0 / dy2; - double idz2 = 1.0 / dz2; - - double factor = s->omega * 0.5 * (dx2 * dy2 * dz2) / - (dy2 * dz2 + dx2 * dz2 + dx2 * dy2); - double* p = s->p; - double* rhs = s->rhs; - double epssq = eps * eps; - int it = 0; - double res = 1.0; - int pass, ksw, jsw, isw; - - while ((res >= epssq) && (it < itermax)) { - ksw = 1; - - for (pass = 0; pass < 2; pass++) { - jsw = ksw; - commExchange(&s->comm, p); - - for (int k = 1; k < kmaxLocal + 1; k++) { - isw = jsw; - for (int j = 1; j < jmaxLocal + 1; j++) { - for (int i = isw; i < imaxLocal + 1; i += 2) { - - double r = - RHS(i, j, k) - - ((P(i + 1, j, k) - 2.0 * P(i, j, k) + P(i - 1, j, k)) * idx2 + - (P(i, j + 1, k) - 2.0 * P(i, j, k) + P(i, j - 1, k)) * - idy2 + - (P(i, j, k + 1) - 2.0 * P(i, j, k) + P(i, j, k - 1)) * - idz2); - - P(i, j, k) -= (factor * r); - res += (r * r); - } - isw = 3 - isw; - } - jsw = 3 - jsw; - } - ksw = 3 - ksw; - } - - if (commIsBoundary(&s->comm, FRONT)) { - for (int j = 1; j < jmaxLocal + 1; j++) { - for (int i = 1; i < imaxLocal + 1; i++) { - P(i, j, 0) = P(i, j, 1); - } - } - } - - if (commIsBoundary(&s->comm, BACK)) { - for (int j = 1; j < jmaxLocal + 1; j++) { - for (int i = 1; i < imaxLocal + 1; i++) { - P(i, j, kmaxLocal + 1) = P(i, j, kmaxLocal); - } - } - } - - if (commIsBoundary(&s->comm, BOTTOM)) { - for (int k = 1; k < kmaxLocal + 1; k++) { - for (int i = 1; i < imaxLocal + 1; i++) { - P(i, 0, k) = P(i, 1, k); - } - } - } - - if (commIsBoundary(&s->comm, TOP)) { - for (int k = 1; k < kmaxLocal + 1; k++) { - for (int i = 1; i < imaxLocal + 1; i++) { - P(i, jmaxLocal + 1, k) = P(i, jmaxLocal, k); - } - } - } - - if (commIsBoundary(&s->comm, LEFT)) { - for (int k = 1; k < kmaxLocal + 1; k++) { - for (int j = 1; j < jmaxLocal + 1; j++) { - P(0, j, k) = P(1, j, k); - } - } - } - - if (commIsBoundary(&s->comm, RIGHT)) { - for (int k = 1; k < kmaxLocal + 1; k++) { - for (int j = 1; j < jmaxLocal + 1; j++) { - P(imaxLocal + 1, j, k) = P(imaxLocal, j, k); - } - } - } - - commReduction(&res, SUM); - res = res / (double)(imax * jmax * kmax); -#ifdef DEBUG - if (commIsMaster(&s->comm)) { - printf("%d Residuum: %e\n", it, res); - } -#endif - commExchange(&s->comm, p); - it++; - } - -#ifdef VERBOSE - if (commIsMaster(&s->comm)) { - printf("Solver took %d iterations to reach %f\n", it, sqrt(res)); - } -#endif -} - -static double maxElement(Solver* s, double* m) -{ - int size = (s->comm.imaxLocal + 2) * (s->comm.jmaxLocal + 2) * - (s->comm.kmaxLocal + 2); - double maxval = DBL_MIN; - - for (int i = 0; i < size; i++) { - maxval = MAX(maxval, fabs(m[i])); - } - commReduction(&maxval, MAX); - return maxval; -} - -void normalizePressure(Solver* s) -{ - int imaxLocal = s->comm.imaxLocal; - int jmaxLocal = s->comm.jmaxLocal; - int kmaxLocal = s->comm.kmaxLocal; - - double* p = s->p; - double avgP = 0.0; - - for (int k = 1; k < kmaxLocal + 1; k++) { - for (int j = 1; j < jmaxLocal + 1; j++) { - for (int i = 1; i < imaxLocal + 1; i++) { - avgP += P(i, j, k); - } - } - } - commReduction(&avgP, SUM); - avgP /= (s->grid.imax * s->grid.jmax * s->grid.kmax); - - for (int k = 1; k < kmaxLocal + 1; k++) { - for (int j = 1; j < jmaxLocal + 1; j++) { - for (int i = 1; i < imaxLocal + 1; i++) { - P(i, j, k) = P(i, j, k) - avgP; - } - } - } -} - -void computeTimestep(Solver* s) -{ - double dt = s->dtBound; - double dx = s->grid.dx; - double dy = s->grid.dy; - double dz = s->grid.dz; - - double umax = maxElement(s, s->u); - double vmax = maxElement(s, s->v); - double wmax = maxElement(s, s->w); - - if (umax > 0) { - dt = (dt > dx / umax) ? dx / umax : dt; - } - if (vmax > 0) { - dt = (dt > dy / vmax) ? dy / vmax : dt; - } - if (wmax > 0) { - dt = (dt > dz / wmax) ? dz / wmax : dt; - } - - s->dt = dt * s->tau; -} - -void setBoundaryConditions(Solver* s) -{ - int imaxLocal = s->comm.imaxLocal; - int jmaxLocal = s->comm.jmaxLocal; - int kmaxLocal = s->comm.kmaxLocal; - - double* u = s->u; - double* v = s->v; - double* w = s->w; - - if (commIsBoundary(&s->comm, TOP)) { - switch (s->bcTop) { - case NOSLIP: - for (int k = 1; k < kmaxLocal + 1; k++) { - for (int i = 1; i < imaxLocal + 1; i++) { - U(i, jmaxLocal + 1, k) = -U(i, jmaxLocal, k); - V(i, jmaxLocal, k) = 0.0; - W(i, jmaxLocal + 1, k) = -W(i, jmaxLocal, k); - } - } - break; - case SLIP: - for (int k = 1; k < kmaxLocal + 1; k++) { - for (int i = 1; i < imaxLocal + 1; i++) { - U(i, jmaxLocal + 1, k) = U(i, jmaxLocal, k); - V(i, jmaxLocal, k) = 0.0; - W(i, jmaxLocal + 1, k) = W(i, jmaxLocal, k); - } - } - break; - case OUTFLOW: - for (int k = 1; k < kmaxLocal + 1; k++) { - for (int i = 1; i < imaxLocal + 1; i++) { - U(i, jmaxLocal + 1, k) = U(i, jmaxLocal, k); - V(i, jmaxLocal, k) = V(i, jmaxLocal - 1, k); - W(i, jmaxLocal + 1, k) = W(i, jmaxLocal, k); - } - } - break; - case PERIODIC: - break; - } - } - - if (commIsBoundary(&s->comm, BOTTOM)) { - switch (s->bcBottom) { - case NOSLIP: - for (int k = 1; k < kmaxLocal + 1; k++) { - for (int i = 1; i < imaxLocal + 1; i++) { - U(i, 0, k) = -U(i, 1, k); - V(i, 0, k) = 0.0; - W(i, 0, k) = -W(i, 1, k); - } - } - break; - case SLIP: - for (int k = 1; k < kmaxLocal + 1; k++) { - for (int i = 1; i < imaxLocal + 1; i++) { - U(i, 0, k) = U(i, 1, k); - V(i, 0, k) = 0.0; - W(i, 0, k) = W(i, 1, k); - } - } - break; - case OUTFLOW: - for (int k = 1; k < kmaxLocal + 1; k++) { - for (int i = 1; i < imaxLocal + 1; i++) { - U(i, 0, k) = U(i, 1, k); - V(i, 0, k) = V(i, 1, k); - W(i, 0, k) = W(i, 1, k); - } - } - break; - case PERIODIC: - break; - } - } - - if (commIsBoundary(&s->comm, LEFT)) { - switch (s->bcLeft) { - case NOSLIP: - for (int k = 1; k < kmaxLocal + 1; k++) { - for (int j = 1; j < jmaxLocal + 1; j++) { - U(0, j, k) = 0.0; - V(0, j, k) = -V(1, j, k); - W(0, j, k) = -W(1, j, k); - } - } - break; - case SLIP: - for (int k = 1; k < kmaxLocal + 1; k++) { - for (int j = 1; j < jmaxLocal + 1; j++) { - U(0, j, k) = 0.0; - V(0, j, k) = V(1, j, k); - W(0, j, k) = W(1, j, k); - } - } - break; - case OUTFLOW: - for (int k = 1; k < kmaxLocal + 1; k++) { - for (int j = 1; j < jmaxLocal + 1; j++) { - U(0, j, k) = U(1, j, k); - V(0, j, k) = V(1, j, k); - W(0, j, k) = W(1, j, k); - } - } - break; - case PERIODIC: - break; - } - } - - if (commIsBoundary(&s->comm, RIGHT)) { - switch (s->bcRight) { - case NOSLIP: - for (int k = 1; k < kmaxLocal + 1; k++) { - for (int j = 1; j < jmaxLocal + 1; j++) { - U(imaxLocal, j, k) = 0.0; - V(imaxLocal + 1, j, k) = -V(imaxLocal, j, k); - W(imaxLocal + 1, j, k) = -W(imaxLocal, j, k); - } - } - break; - case SLIP: - for (int k = 1; k < kmaxLocal + 1; k++) { - for (int j = 1; j < jmaxLocal + 1; j++) { - U(imaxLocal, j, k) = 0.0; - V(imaxLocal + 1, j, k) = V(imaxLocal, j, k); - W(imaxLocal + 1, j, k) = W(imaxLocal, j, k); - } - } - break; - case OUTFLOW: - for (int k = 1; k < kmaxLocal + 1; k++) { - for (int j = 1; j < jmaxLocal + 1; j++) { - U(imaxLocal, j, k) = U(imaxLocal - 1, j, k); - V(imaxLocal + 1, j, k) = V(imaxLocal, j, k); - W(imaxLocal + 1, j, k) = W(imaxLocal, j, k); - } - } - break; - case PERIODIC: - break; - } - } - - if (commIsBoundary(&s->comm, FRONT)) { - switch (s->bcFront) { - case NOSLIP: - for (int j = 1; j < jmaxLocal + 1; j++) { - for (int i = 1; i < imaxLocal + 1; i++) { - U(i, j, 0) = -U(i, j, 1); - V(i, j, 0) = -V(i, j, 1); - W(i, j, 0) = 0.0; - } - } - break; - case SLIP: - for (int j = 1; j < jmaxLocal + 1; j++) { - for (int i = 1; i < imaxLocal + 1; i++) { - U(i, j, 0) = U(i, j, 1); - V(i, j, 0) = V(i, j, 1); - W(i, j, 0) = 0.0; - } - } - break; - case OUTFLOW: - for (int j = 1; j < jmaxLocal + 1; j++) { - for (int i = 1; i < imaxLocal + 1; i++) { - U(i, j, 0) = U(i, j, 1); - V(i, j, 0) = V(i, j, 1); - W(i, j, 0) = W(i, j, 1); - } - } - break; - case PERIODIC: - break; - } - } - - if (commIsBoundary(&s->comm, BACK)) { - switch (s->bcBack) { - case NOSLIP: - for (int j = 1; j < jmaxLocal + 1; j++) { - for (int i = 1; i < imaxLocal + 1; i++) { - U(i, j, kmaxLocal + 1) = -U(i, j, kmaxLocal); - V(i, j, kmaxLocal + 1) = -V(i, j, kmaxLocal); - W(i, j, kmaxLocal) = 0.0; - } - } - break; - case SLIP: - for (int j = 1; j < jmaxLocal + 1; j++) { - for (int i = 1; i < imaxLocal + 1; i++) { - U(i, j, kmaxLocal + 1) = U(i, j, kmaxLocal); - V(i, j, kmaxLocal + 1) = V(i, j, kmaxLocal); - W(i, j, kmaxLocal) = 0.0; - } - } - break; - case OUTFLOW: - for (int j = 1; j < jmaxLocal + 1; j++) { - for (int i = 1; i < imaxLocal + 1; i++) { - U(i, j, kmaxLocal + 1) = U(i, j, kmaxLocal); - V(i, j, kmaxLocal + 1) = V(i, j, kmaxLocal); - W(i, j, kmaxLocal) = W(i, j, kmaxLocal - 1); - } - } - break; - case PERIODIC: - break; - } - } -} - -void setSpecialBoundaryCondition(Solver* s) -{ - int imaxLocal = s->comm.imaxLocal; - int jmaxLocal = s->comm.jmaxLocal; - int kmaxLocal = s->comm.kmaxLocal; - - double* u = s->u; - - if (strcmp(s->problem, "dcavity") == 0) { - if (commIsBoundary(&s->comm, TOP)) { - for (int k = 1; k < kmaxLocal; k++) { - for (int i = 1; i < imaxLocal; i++) { - U(i, jmaxLocal + 1, k) = 2.0 - U(i, jmaxLocal, k); - } - } - } - } else if (strcmp(s->problem, "canal") == 0) { - if (commIsBoundary(&s->comm, LEFT)) { - for (int k = 1; k < kmaxLocal + 1; k++) { - for (int j = 1; j < jmaxLocal + 1; j++) { - U(0, j, k) = 2.0; - } - } - } - } -} - -void computeFG(Solver* s) -{ - int imaxLocal = s->comm.imaxLocal; - int jmaxLocal = s->comm.jmaxLocal; - int kmaxLocal = s->comm.kmaxLocal; - - double* u = s->u; - double* v = s->v; - double* w = s->w; - double* f = s->f; - double* g = s->g; - double* h = s->h; - - double gx = s->gx; - double gy = s->gy; - double gz = s->gz; - double dt = s->dt; - - double gamma = s->gamma; - double inverseRe = 1.0 / s->re; - double inverseDx = 1.0 / s->grid.dx; - double inverseDy = 1.0 / s->grid.dy; - double inverseDz = 1.0 / s->grid.dz; - double du2dx, dv2dy, dw2dz; - double duvdx, duwdx, duvdy, dvwdy, duwdz, dvwdz; - double du2dx2, du2dy2, du2dz2; - double dv2dx2, dv2dy2, dv2dz2; - double dw2dx2, dw2dy2, dw2dz2; - - commExchange(&s->comm, u); - commExchange(&s->comm, v); - commExchange(&s->comm, w); - - for (int k = 1; k < kmaxLocal + 1; k++) { - for (int j = 1; j < jmaxLocal + 1; j++) { - for (int i = 1; i < imaxLocal + 1; i++) { - du2dx = inverseDx * 0.25 * - ((U(i, j, k) + U(i + 1, j, k)) * - (U(i, j, k) + U(i + 1, j, k)) - - (U(i, j, k) + U(i - 1, j, k)) * - (U(i, j, k) + U(i - 1, j, k))) + - gamma * inverseDx * 0.25 * - (fabs(U(i, j, k) + U(i + 1, j, k)) * - (U(i, j, k) - U(i + 1, j, k)) + - fabs(U(i, j, k) + U(i - 1, j, k)) * - (U(i, j, k) - U(i - 1, j, k))); - - duvdy = inverseDy * 0.25 * - ((V(i, j, k) + V(i + 1, j, k)) * - (U(i, j, k) + U(i, j + 1, k)) - - (V(i, j - 1, k) + V(i + 1, j - 1, k)) * - (U(i, j, k) + U(i, j - 1, k))) + - gamma * inverseDy * 0.25 * - (fabs(V(i, j, k) + V(i + 1, j, k)) * - (U(i, j, k) - U(i, j + 1, k)) + - fabs(V(i, j - 1, k) + V(i + 1, j - 1, k)) * - (U(i, j, k) - U(i, j - 1, k))); - - duwdz = inverseDz * 0.25 * - ((W(i, j, k) + W(i + 1, j, k)) * - (U(i, j, k) + U(i, j, k + 1)) - - (W(i, j, k - 1) + W(i + 1, j, k - 1)) * - (U(i, j, k) + U(i, j, k - 1))) + - gamma * inverseDz * 0.25 * - (fabs(W(i, j, k) + W(i + 1, j, k)) * - (U(i, j, k) - U(i, j, k + 1)) + - fabs(W(i, j, k - 1) + W(i + 1, j, k - 1)) * - (U(i, j, k) - U(i, j, k - 1))); - - du2dx2 = inverseDx * inverseDx * - (U(i + 1, j, k) - 2.0 * U(i, j, k) + U(i - 1, j, k)); - du2dy2 = inverseDy * inverseDy * - (U(i, j + 1, k) - 2.0 * U(i, j, k) + U(i, j - 1, k)); - du2dz2 = inverseDz * inverseDz * - (U(i, j, k + 1) - 2.0 * U(i, j, k) + U(i, j, k - 1)); - F(i, j, k) = U(i, j, k) + dt * (inverseRe * (du2dx2 + du2dy2 + du2dz2) - - du2dx - duvdy - duwdz + gx); - - duvdx = inverseDx * 0.25 * - ((U(i, j, k) + U(i, j + 1, k)) * - (V(i, j, k) + V(i + 1, j, k)) - - (U(i - 1, j, k) + U(i - 1, j + 1, k)) * - (V(i, j, k) + V(i - 1, j, k))) + - gamma * inverseDx * 0.25 * - (fabs(U(i, j, k) + U(i, j + 1, k)) * - (V(i, j, k) - V(i + 1, j, k)) + - fabs(U(i - 1, j, k) + U(i - 1, j + 1, k)) * - (V(i, j, k) - V(i - 1, j, k))); - - dv2dy = inverseDy * 0.25 * - ((V(i, j, k) + V(i, j + 1, k)) * - (V(i, j, k) + V(i, j + 1, k)) - - (V(i, j, k) + V(i, j - 1, k)) * - (V(i, j, k) + V(i, j - 1, k))) + - gamma * inverseDy * 0.25 * - (fabs(V(i, j, k) + V(i, j + 1, k)) * - (V(i, j, k) - V(i, j + 1, k)) + - fabs(V(i, j, k) + V(i, j - 1, k)) * - (V(i, j, k) - V(i, j - 1, k))); - - dvwdz = inverseDz * 0.25 * - ((W(i, j, k) + W(i, j + 1, k)) * - (V(i, j, k) + V(i, j, k + 1)) - - (W(i, j, k - 1) + W(i, j + 1, k - 1)) * - (V(i, j, k) + V(i, j, k + 1))) + - gamma * inverseDz * 0.25 * - (fabs(W(i, j, k) + W(i, j + 1, k)) * - (V(i, j, k) - V(i, j, k + 1)) + - fabs(W(i, j, k - 1) + W(i, j + 1, k - 1)) * - (V(i, j, k) - V(i, j, k + 1))); - - dv2dx2 = inverseDx * inverseDx * - (V(i + 1, j, k) - 2.0 * V(i, j, k) + V(i - 1, j, k)); - dv2dy2 = inverseDy * inverseDy * - (V(i, j + 1, k) - 2.0 * V(i, j, k) + V(i, j - 1, k)); - dv2dz2 = inverseDz * inverseDz * - (V(i, j, k + 1) - 2.0 * V(i, j, k) + V(i, j, k - 1)); - G(i, j, k) = V(i, j, k) + dt * (inverseRe * (dv2dx2 + dv2dy2 + dv2dz2) - - duvdx - dv2dy - dvwdz + gy); - - duwdx = inverseDx * 0.25 * - ((U(i, j, k) + U(i, j, k + 1)) * - (W(i, j, k) + W(i + 1, j, k)) - - (U(i - 1, j, k) + U(i - 1, j, k + 1)) * - (W(i, j, k) + W(i - 1, j, k))) + - gamma * inverseDx * 0.25 * - (fabs(U(i, j, k) + U(i, j, k + 1)) * - (W(i, j, k) - W(i + 1, j, k)) + - fabs(U(i - 1, j, k) + U(i - 1, j, k + 1)) * - (W(i, j, k) - W(i - 1, j, k))); - - dvwdy = inverseDy * 0.25 * - ((V(i, j, k) + V(i, j, k + 1)) * - (W(i, j, k) + W(i, j + 1, k)) - - (V(i, j - 1, k + 1) + V(i, j - 1, k)) * - (W(i, j, k) + W(i, j - 1, k))) + - gamma * inverseDy * 0.25 * - (fabs(V(i, j, k) + V(i, j, k + 1)) * - (W(i, j, k) - W(i, j + 1, k)) + - fabs(V(i, j - 1, k + 1) + V(i, j - 1, k)) * - (W(i, j, k) - W(i, j - 1, k))); - - dw2dz = inverseDz * 0.25 * - ((W(i, j, k) + W(i, j, k + 1)) * - (W(i, j, k) + W(i, j, k + 1)) - - (W(i, j, k) + W(i, j, k - 1)) * - (W(i, j, k) + W(i, j, k - 1))) + - gamma * inverseDz * 0.25 * - (fabs(W(i, j, k) + W(i, j, k + 1)) * - (W(i, j, k) - W(i, j, k + 1)) + - fabs(W(i, j, k) + W(i, j, k - 1)) * - (W(i, j, k) - W(i, j, k - 1))); - - dw2dx2 = inverseDx * inverseDx * - (W(i + 1, j, k) - 2.0 * W(i, j, k) + W(i - 1, j, k)); - dw2dy2 = inverseDy * inverseDy * - (W(i, j + 1, k) - 2.0 * W(i, j, k) + W(i, j - 1, k)); - dw2dz2 = inverseDz * inverseDz * - (W(i, j, k + 1) - 2.0 * W(i, j, k) + W(i, j, k - 1)); - H(i, j, k) = W(i, j, k) + dt * (inverseRe * (dw2dx2 + dw2dy2 + dw2dz2) - - duwdx - dvwdy - dw2dz + gz); - } - } - } - - /* ----------------------------- boundary of F --------------------------- - */ - if (commIsBoundary(&s->comm, LEFT)) { - for (int k = 1; k < kmaxLocal + 1; k++) { - for (int j = 1; j < jmaxLocal + 1; j++) { - F(0, j, k) = U(0, j, k); - } - } - } - - if (commIsBoundary(&s->comm, RIGHT)) { - for (int k = 1; k < kmaxLocal + 1; k++) { - for (int j = 1; j < jmaxLocal + 1; j++) { - F(imaxLocal, j, k) = U(imaxLocal, j, k); - } - } - } - - /* ----------------------------- boundary of G --------------------------- - */ - if (commIsBoundary(&s->comm, BOTTOM)) { - for (int k = 1; k < kmaxLocal + 1; k++) { - for (int i = 1; i < imaxLocal + 1; i++) { - G(i, 0, k) = V(i, 0, k); - } - } - } - - if (commIsBoundary(&s->comm, TOP)) { - for (int k = 1; k < kmaxLocal + 1; k++) { - for (int i = 1; i < imaxLocal + 1; i++) { - G(i, jmaxLocal, k) = V(i, jmaxLocal, k); - } - } - } - - /* ----------------------------- boundary of H --------------------------- - */ - if (commIsBoundary(&s->comm, FRONT)) { - for (int j = 1; j < jmaxLocal + 1; j++) { - for (int i = 1; i < imaxLocal + 1; i++) { - H(i, j, 0) = W(i, j, 0); - } - } - } - - if (commIsBoundary(&s->comm, BACK)) { - for (int j = 1; j < jmaxLocal + 1; j++) { - for (int i = 1; i < imaxLocal + 1; i++) { - H(i, j, kmaxLocal) = W(i, j, kmaxLocal); - } - } - } -} - -void adaptUV(Solver* s) -{ - int imaxLocal = s->comm.imaxLocal; - int jmaxLocal = s->comm.jmaxLocal; - int kmaxLocal = s->comm.kmaxLocal; - - double* p = s->p; - double* u = s->u; - double* v = s->v; - double* w = s->w; - double* f = s->f; - double* g = s->g; - double* h = s->h; - - double factorX = s->dt / s->grid.dx; - double factorY = s->dt / s->grid.dy; - double factorZ = s->dt / s->grid.dz; - - for (int k = 1; k < kmaxLocal + 1; k++) { - for (int j = 1; j < jmaxLocal + 1; j++) { - for (int i = 1; i < imaxLocal + 1; i++) { - U(i, j, k) = F(i, j, k) - (P(i + 1, j, k) - P(i, j, k)) * factorX; - V(i, j, k) = G(i, j, k) - (P(i, j + 1, k) - P(i, j, k)) * factorY; - W(i, j, k) = H(i, j, k) - (P(i, j, k + 1) - P(i, j, k)) * factorZ; - } - } - } -} diff --git a/BasicSolver/3D-mpi-io/src/solver.h b/BasicSolver/3D-mpi-io/src/solver.h deleted file mode 100644 index f951ec8..0000000 --- a/BasicSolver/3D-mpi-io/src/solver.h +++ /dev/null @@ -1,45 +0,0 @@ -/* - * Copyright (C) 2022 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 "comm.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, *h; - double *u, *v, *w; - /* parameters */ - double eps, omega; - double re, tau, gamma; - double gx, gy, gz; - /* time stepping */ - int itermax; - double dt, te; - double dtBound; - char* problem; - int bcLeft, bcRight, bcBottom, bcTop, bcFront, bcBack; - /* communication */ - Comm comm; -} Solver; - -extern void initSolver(Solver*, Parameter*); -extern void computeRHS(Solver*); -extern void solve(Solver*); -extern void normalizePressure(Solver*); -extern void computeTimestep(Solver*); -extern void setBoundaryConditions(Solver*); -extern void setSpecialBoundaryCondition(Solver*); -extern void computeFG(Solver*); -extern void adaptUV(Solver*); -#endif diff --git a/BasicSolver/3D-mpi-io/src/test.c b/BasicSolver/3D-mpi-io/src/test.c deleted file mode 100644 index 5f61b9f..0000000 --- a/BasicSolver/3D-mpi-io/src/test.c +++ /dev/null @@ -1,129 +0,0 @@ -/* - * Copyright (C) 2022 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 "test.h" - -#define G(v, i, j, k) \ - v[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)] - -void testInit(Solver* s) -{ - int imaxLocal = s->comm.imaxLocal; - int jmaxLocal = s->comm.jmaxLocal; - int kmaxLocal = s->comm.kmaxLocal; - int myrank = s->comm.rank; - double* p = s->p; - double* f = s->f; - double* g = s->g; - double* h = s->h; - - for (int k = 0; k < kmaxLocal + 2; k++) { - for (int j = 0; j < jmaxLocal + 2; j++) { - for (int i = 0; i < imaxLocal + 2; i++) { - G(p, i, j, k) = 10.0; - G(f, i, j, k) = myrank + 1.0; - G(g, i, j, k) = myrank + 1.0; - G(h, i, j, k) = myrank + 1.0; - } - } - } - - for (int k = 1; k < kmaxLocal + 1; k++) { - for (int j = 1; j < jmaxLocal + 1; j++) { - for (int i = 1; i < imaxLocal + 1; i++) { - G(p, i, j, k) = myrank + 1.0; - G(f, i, j, k) = myrank + 1.0; - G(g, i, j, k) = myrank + 1.0; - G(h, i, j, k) = myrank + 1.0; - } - } - } -} - -static char* direction2String(Direction dir) -{ - switch (dir) { - case LEFT: - return "left"; - break; - case RIGHT: - return "right"; - break; - case BOTTOM: - return "bottom"; - break; - case TOP: - return "top"; - break; - case FRONT: - return "front"; - break; - case BACK: - return "back"; - break; - case NDIRS: - return "ERROR"; - default: - return "ERROR"; - } -} - -static void printPlane(Solver* s, double* a, int ymax, int xmax, Direction dir) -{ - int imaxLocal = s->comm.imaxLocal; - int jmaxLocal = s->comm.jmaxLocal; - int kmaxLocal = s->comm.kmaxLocal; - char filename[50]; - snprintf(filename, 50, "halo-%s-r%d.txt", direction2String(dir), s->comm.rank); - FILE* fh = fopen(filename, "w"); - - for (int y = 0; y < ymax; y++) { - for (int x = 0; x < xmax; x++) { - switch (dir) { - case LEFT: - fprintf(fh, "%12.8f ", G(a, 0, x, y)); - break; - case RIGHT: - fprintf(fh, "%12.8f ", G(a, imaxLocal + 1, x, y)); - break; - case BOTTOM: - fprintf(fh, "%12.8f ", G(a, x, 0, y)); - break; - case TOP: - fprintf(fh, "%12.8f ", G(a, x, jmaxLocal + 1, y)); - break; - case FRONT: - fprintf(fh, "%12.8f ", G(a, x, y, 0)); - break; - case BACK: - fprintf(fh, "%12.8f ", G(a, x, y, kmaxLocal + 1)); - break; - case NDIRS: - printf("ERROR\n"); - break; - } - } - fprintf(fh, "\n"); - } - fclose(fh); -} - -void testPrintHalo(Solver* s, double* a) -{ - int imaxLocal = s->comm.imaxLocal; - int jmaxLocal = s->comm.jmaxLocal; - int kmaxLocal = s->comm.kmaxLocal; - - printPlane(s, a, kmaxLocal + 2, imaxLocal + 2, BOTTOM); - printPlane(s, a, kmaxLocal + 2, imaxLocal + 2, TOP); - printPlane(s, a, kmaxLocal + 2, jmaxLocal + 2, LEFT); - printPlane(s, a, kmaxLocal + 2, jmaxLocal + 2, RIGHT); - printPlane(s, a, jmaxLocal + 2, imaxLocal + 2, FRONT); - printPlane(s, a, jmaxLocal + 2, imaxLocal + 2, BACK); -} diff --git a/BasicSolver/3D-mpi-io/src/test.h b/BasicSolver/3D-mpi-io/src/test.h deleted file mode 100644 index ecefe4e..0000000 --- a/BasicSolver/3D-mpi-io/src/test.h +++ /dev/null @@ -1,13 +0,0 @@ -/* - * Copyright (C) 2022 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 __TEST_H_ -#define __TEST_H_ -#include "solver.h" - -extern void testInit(Solver* s); -extern void testPrintHalo(Solver* s, double* a); -#endif // __TEST_H_ diff --git a/BasicSolver/3D-mpi-io/src/timing.c b/BasicSolver/3D-mpi-io/src/timing.c deleted file mode 100644 index e578f18..0000000 --- a/BasicSolver/3D-mpi-io/src/timing.c +++ /dev/null @@ -1,22 +0,0 @@ -/* - * Copyright (C) 2022 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(void) -{ - struct timespec ts; - clock_gettime(CLOCK_MONOTONIC, &ts); - return (double)ts.tv_sec + (double)ts.tv_nsec * 1.e-9; -} - -double getTimeResolution(void) -{ - struct timespec ts; - clock_getres(CLOCK_MONOTONIC, &ts); - return (double)ts.tv_sec + (double)ts.tv_nsec * 1.e-9; -} diff --git a/BasicSolver/3D-mpi-io/src/timing.h b/BasicSolver/3D-mpi-io/src/timing.h deleted file mode 100644 index 58fb5ac..0000000 --- a/BasicSolver/3D-mpi-io/src/timing.h +++ /dev/null @@ -1,13 +0,0 @@ -/* - * Copyright (C) 2022 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(void); -extern double getTimeResolution(void); - -#endif // __TIMING_H_ diff --git a/BasicSolver/3D-mpi-io/src/util.h b/BasicSolver/3D-mpi-io/src/util.h deleted file mode 100644 index 657b009..0000000 --- a/BasicSolver/3D-mpi-io/src/util.h +++ /dev/null @@ -1,22 +0,0 @@ -/* - * Copyright (C) 2022 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/BasicSolver/3D-mpi-io/src/vtkWriter.h b/BasicSolver/3D-mpi-io/src/vtkWriter.h deleted file mode 100644 index 7756826..0000000 --- a/BasicSolver/3D-mpi-io/src/vtkWriter.h +++ /dev/null @@ -1,29 +0,0 @@ -/* - * Copyright (C) 2022 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 - -#include "comm.h" -#include "grid.h" - -typedef struct VtkOptions { - Grid grid; - MPI_File fh; - Comm comm; -} VtkOptions; - -typedef struct VtkVector { - double *u, *v, *w; -} VtkVector; - -extern void vtkOpen(VtkOptions* opts, char* filename); -extern void vtkVector(VtkOptions* opts, char* name, VtkVector vec); -extern void vtkScalar(VtkOptions* opts, char* name, double* p); -extern void vtkClose(VtkOptions* opts); -#endif // __VTKWRITER_H_ diff --git a/BasicSolver/3D-mpi/Makefile b/BasicSolver/3D-mpi/Makefile index 57f99f4..f61416c 100644 --- a/BasicSolver/3D-mpi/Makefile +++ b/BasicSolver/3D-mpi/Makefile @@ -1,5 +1,5 @@ #======================================================================================= -# Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg. +# Copyright (C) NHR@FAU, University Erlangen-Nuremberg. # All rights reserved. # Use of this source code is governed by a MIT-style # license that can be found in the LICENSE file. @@ -18,13 +18,18 @@ 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)/vtkWriter-$(VTK_OUTPUT_FMT).o SOURCES = $(SRC) $(wildcard $(SRC_DIR)/*.h) +ifeq ($(VTK_OUTPUT_FMT),mpi) +DEFINES += -D_VTK_WRITER_MPI +endif + CPPFLAGS := $(CPPFLAGS) $(DEFINES) $(OPTIONS) $(INCLUDES) -${TARGET}: $(BUILD_DIR) $(OBJ) +${TARGET}: sanity-checks $(BUILD_DIR) $(OBJ) $(info ===> LINKING $(TARGET)) $(Q)${LINKER} ${LFLAGS} -o $(TARGET) $(OBJ) $(LIBS) @@ -65,6 +70,14 @@ format: done @echo "Done" +sanity-checks: +ifeq ($(VTK_OUTPUT_FMT),mpi) +ifeq ($(ENABLE_MPI),false) + $(error VTK_OUTPUT_FMT mpi only supported for ENABLE_MPI true!) +endif +endif + + $(BUILD_DIR): @mkdir $(BUILD_DIR) diff --git a/BasicSolver/3D-mpi/config.mk b/BasicSolver/3D-mpi/config.mk index c163666..52154f4 100644 --- a/BasicSolver/3D-mpi/config.mk +++ b/BasicSolver/3D-mpi/config.mk @@ -2,6 +2,8 @@ TAG ?= CLANG ENABLE_MPI ?= true ENABLE_OPENMP ?= false +# Supported: seq, mpi +VTK_OUTPUT_FMT ?= seq #Feature options OPTIONS += -DARRAY_ALIGNMENT=64 diff --git a/BasicSolver/3D-mpi/include_CLANG.mk b/BasicSolver/3D-mpi/include_CLANG.mk index 5b036ab..3641cad 100644 --- a/BasicSolver/3D-mpi/include_CLANG.mk +++ b/BasicSolver/3D-mpi/include_CLANG.mk @@ -2,7 +2,7 @@ ifeq ($(ENABLE_MPI),true) CC = mpicc DEFINES = -D_MPI else -CC = cc +CC = cc endif GCC = cc @@ -18,4 +18,4 @@ VERSION = --version CFLAGS = -Ofast -std=c17 LFLAGS = $(OPENMP) -lm DEFINES += -D_GNU_SOURCE# -DDEBUG -INCLUDES = +INCLUDES = -I/opt/homebrew/include diff --git a/BasicSolver/3D-mpi/include_GCC.mk b/BasicSolver/3D-mpi/include_GCC.mk index 427e798..cf1e2aa 100644 --- a/BasicSolver/3D-mpi/include_GCC.mk +++ b/BasicSolver/3D-mpi/include_GCC.mk @@ -1,4 +1,10 @@ +ifeq ($(ENABLE_MPI),true) +CC = mpicc +DEFINES = -D_MPI +else CC = gcc +endif + GCC = gcc LINKER = $(CC) @@ -9,6 +15,6 @@ endif VERSION = --version CFLAGS = -Ofast -ffreestanding -std=c99 $(OPENMP) LFLAGS = $(OPENMP) -DEFINES = -D_GNU_SOURCE +DEFINES += -D_GNU_SOURCE INCLUDES = LIBS = diff --git a/BasicSolver/3D-mpi/include_ICC.mk b/BasicSolver/3D-mpi/include_ICC.mk index 94b8e20..02d7076 100644 --- a/BasicSolver/3D-mpi/include_ICC.mk +++ b/BasicSolver/3D-mpi/include_ICC.mk @@ -1,4 +1,10 @@ +ifeq ($(ENABLE_MPI),true) +CC = mpiicc +DEFINES = -D_MPI +else CC = icc +endif + GCC = gcc LINKER = $(CC) @@ -9,6 +15,6 @@ endif VERSION = --version CFLAGS = -O3 -xHost -qopt-zmm-usage=high -std=c99 $(OPENMP) LFLAGS = $(OPENMP) -DEFINES = -D_GNU_SOURCE +DEFINES += -D_GNU_SOURCE INCLUDES = LIBS = diff --git a/BasicSolver/3D-mpi/src/allocate.c b/BasicSolver/3D-mpi/src/allocate.c index 6af6c7f..cf2efd6 100644 --- a/BasicSolver/3D-mpi/src/allocate.c +++ b/BasicSolver/3D-mpi/src/allocate.c @@ -1,5 +1,5 @@ /* - * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg. + * Copyright (C) NHR@FAU, University Erlangen-Nuremberg. * All rights reserved. * Use of this source code is governed by a MIT-style * license that can be found in the LICENSE file. diff --git a/BasicSolver/3D-mpi/src/allocate.h b/BasicSolver/3D-mpi/src/allocate.h index 1537eb2..77f4ba0 100644 --- a/BasicSolver/3D-mpi/src/allocate.h +++ b/BasicSolver/3D-mpi/src/allocate.h @@ -1,5 +1,5 @@ /* - * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg. + * Copyright (C) NHR@FAU, University Erlangen-Nuremberg. * All rights reserved. * Use of this source code is governed by a MIT-style * license that can be found in the LICENSE file. diff --git a/BasicSolver/3D-mpi/src/comm.c b/BasicSolver/3D-mpi/src/comm.c index fb8f519..8a0c870 100644 --- a/BasicSolver/3D-mpi/src/comm.c +++ b/BasicSolver/3D-mpi/src/comm.c @@ -1,5 +1,5 @@ /* - * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg. + * 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. @@ -296,6 +296,30 @@ void commShift(Comm* c, double* f, double* g, double* h) #endif } +void commGetOffsets(Comm* c, int offsets[], int kmax, int jmax, int imax) +{ +#if defined(_MPI) + int sum = 0; + + for (int i = 0; i < c->coords[ICORD]; i++) { + sum += sizeOfRank(i, c->dims[ICORD], imax); + } + offsets[IDIM] = sum; + sum = 0; + + for (int i = 0; i < c->coords[JCORD]; i++) { + sum += sizeOfRank(i, c->dims[JCORD], jmax); + } + offsets[JDIM] = sum; + sum = 0; + + for (int i = 0; i < c->coords[KCORD]; i++) { + sum += sizeOfRank(i, c->dims[KCORD], kmax); + } + offsets[KDIM] = sum; +#endif +} + #define G(v, i, j, k) \ v[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)] @@ -445,7 +469,7 @@ void commCollectResult(Comm* c, for (int k = 1; k < kmaxLocal + 1; k++) { for (int j = 1; j < jmaxLocal + 1; j++) { for (int i = 1; i < imaxLocal + 1; i++) { - pg[idx++] = G(s->p, i, j, k); + pg[idx++] = G(p, i, j, k); } } } @@ -455,7 +479,7 @@ void commCollectResult(Comm* c, for (int k = 1; k < kmaxLocal + 1; k++) { for (int j = 1; j < jmaxLocal + 1; j++) { for (int i = 1; i < imaxLocal + 1; i++) { - ug[idx++] = (G(s->u, i, j, k) + G(s->u, i - 1, j, k)) / 2.0; + ug[idx++] = (G(u, i, j, k) + G(u, i - 1, j, k)) / 2.0; } } } @@ -465,7 +489,7 @@ void commCollectResult(Comm* c, for (int k = 1; k < kmaxLocal + 1; k++) { for (int j = 1; j < jmaxLocal + 1; j++) { for (int i = 1; i < imaxLocal + 1; i++) { - vg[idx++] = (G(s->v, i, j, k) + G(s->v, i, j - 1, k)) / 2.0; + vg[idx++] = (G(v, i, j, k) + G(v, i, j - 1, k)) / 2.0; } } } @@ -475,7 +499,7 @@ void commCollectResult(Comm* c, for (int k = 1; k < kmaxLocal + 1; k++) { for (int j = 1; j < jmaxLocal + 1; j++) { for (int i = 1; i < imaxLocal + 1; i++) { - wg[idx++] = (G(s->w, i, j, k) + G(s->w, i, j, k - 1)) / 2.0; + wg[idx++] = (G(w, i, j, k) + G(w, i, j, k - 1)) / 2.0; } } } @@ -524,8 +548,8 @@ void commInit(Comm* c, int argc, char** argv) MPI_Comm_rank(MPI_COMM_WORLD, &(c->rank)); MPI_Comm_size(MPI_COMM_WORLD, &(c->size)); #else - c->rank = 0; - c->size = 1; + c->rank = 0; + c->size = 1; #endif } diff --git a/BasicSolver/3D-mpi/src/comm.h b/BasicSolver/3D-mpi/src/comm.h index 64f7082..48102fb 100644 --- a/BasicSolver/3D-mpi/src/comm.h +++ b/BasicSolver/3D-mpi/src/comm.h @@ -1,5 +1,5 @@ /* - * Copyright (C) 2024 NHR@FAU, University Erlangen-Nuremberg. + * 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. @@ -44,6 +44,7 @@ extern void commExchange(Comm*, double*); extern void commShift(Comm* c, double* f, double* g, double* h); extern void commReduction(double* v, int op); extern int commIsBoundary(Comm* c, Direction direction); +extern void commGetOffsets(Comm* c, int offsets[], int kmax, int jmax, int imax); extern void commCollectResult(Comm* c, double* ug, double* vg, diff --git a/BasicSolver/3D-mpi/src/grid.h b/BasicSolver/3D-mpi/src/grid.h index c963429..0689ebc 100644 --- a/BasicSolver/3D-mpi/src/grid.h +++ b/BasicSolver/3D-mpi/src/grid.h @@ -1,5 +1,5 @@ /* - * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg. + * 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. diff --git a/BasicSolver/3D-mpi/src/main.c b/BasicSolver/3D-mpi/src/main.c index cad1f8e..9562307 100644 --- a/BasicSolver/3D-mpi/src/main.c +++ b/BasicSolver/3D-mpi/src/main.c @@ -1,11 +1,9 @@ /* - * Copyright (C) 2024 NHR@FAU, University Erlangen-Nuremberg. + * Copyright (C) NHR@FAU, University Erlangen-Nuremberg. * All rights reserved. * Use of this source code is governed by a MIT-style * license that can be found in the LICENSE file. */ -#include -#include #include #include #include @@ -73,6 +71,13 @@ int main(int argc, char** argv) printf("Solution took %.2fs\n", timeStop - timeStart); } +#ifdef _VTK_WRITER_MPI + VtkOptions opts = { .grid = s.grid, .comm = s.comm }; + vtkOpen(&opts, s.problem); + vtkScalar(&opts, "pressure", s.p); + vtkVector(&opts, "velocity", (VtkVector) { s.u, s.v, s.w }); + vtkClose(&opts); +#else double *pg, *ug, *vg, *wg; if (commIsMaster(&s.comm)) { @@ -104,6 +109,7 @@ int main(int argc, char** argv) vtkVector(&opts, "velocity", (VtkVector) { ug, vg, wg }); vtkClose(&opts); } +#endif commFinalize(&s.comm); return EXIT_SUCCESS; diff --git a/BasicSolver/3D-mpi/src/parameter.c b/BasicSolver/3D-mpi/src/parameter.c index d4ea5b6..2128c2a 100644 --- a/BasicSolver/3D-mpi/src/parameter.c +++ b/BasicSolver/3D-mpi/src/parameter.c @@ -1,5 +1,5 @@ /* - * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg. + * 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. diff --git a/BasicSolver/3D-mpi/src/parameter.h b/BasicSolver/3D-mpi/src/parameter.h index 6dddf5f..d946d5c 100644 --- a/BasicSolver/3D-mpi/src/parameter.h +++ b/BasicSolver/3D-mpi/src/parameter.h @@ -1,5 +1,5 @@ /* - * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg. + * 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. diff --git a/BasicSolver/3D-mpi/src/progress.c b/BasicSolver/3D-mpi/src/progress.c index a9b82bd..d5393c4 100644 --- a/BasicSolver/3D-mpi/src/progress.c +++ b/BasicSolver/3D-mpi/src/progress.c @@ -1,5 +1,5 @@ /* - * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg. + * 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. diff --git a/BasicSolver/3D-mpi/src/progress.h b/BasicSolver/3D-mpi/src/progress.h index 4d02cdb..314e921 100644 --- a/BasicSolver/3D-mpi/src/progress.h +++ b/BasicSolver/3D-mpi/src/progress.h @@ -1,5 +1,5 @@ /* - * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg. + * Copyright (C) NHR@FAU, University Erlangen-Nuremberg. * All rights reserved. * Use of this source code is governed by a MIT-style * license that can be found in the LICENSE file. diff --git a/BasicSolver/3D-mpi/src/solver.c b/BasicSolver/3D-mpi/src/solver.c index 1a6282f..121b4dc 100644 --- a/BasicSolver/3D-mpi/src/solver.c +++ b/BasicSolver/3D-mpi/src/solver.c @@ -1,5 +1,5 @@ /* - * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg. + * 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. @@ -207,24 +207,24 @@ void solve(Solver* s) jsw = ksw; commExchange(&s->comm, p); - for (int k = 1; k < kmaxLocal + 1; k++) { + for (int k = 1; k < kmaxLocal + 1; k++) { isw = jsw; - for (int j = 1; j < jmaxLocal + 1; j++) { + for (int j = 1; j < jmaxLocal + 1; j++) { for (int i = isw; i < imaxLocal + 1; i += 2) { double r = RHS(i, j, k) - ((P(i + 1, j, k) - 2.0 * P(i, j, k) + P(i - 1, j, k)) * idx2 + - (P(i, j + 1, k) - 2.0 * P(i, j, k) + P(i, j - 1, k)) * - idy2 + - (P(i, j, k + 1) - 2.0 * P(i, j, k) + P(i, j, k - 1)) * - idz2); + (P(i, j + 1, k) - 2.0 * P(i, j, k) + P(i, j - 1, k)) * + idy2 + + (P(i, j, k + 1) - 2.0 * P(i, j, k) + P(i, j, k - 1)) * + idz2); - P(i, j, k) -= (factor * r); - res += (r * r); - } - isw = 3 - isw; + P(i, j, k) -= (factor * r); + res += (r * r); } + isw = 3 - isw; + } jsw = 3 - jsw; } ksw = 3 - ksw; diff --git a/BasicSolver/3D-mpi/src/solver.h b/BasicSolver/3D-mpi/src/solver.h index f951ec8..ae734f7 100644 --- a/BasicSolver/3D-mpi/src/solver.h +++ b/BasicSolver/3D-mpi/src/solver.h @@ -1,5 +1,5 @@ /* - * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg. + * 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. diff --git a/BasicSolver/3D-mpi/src/test.c b/BasicSolver/3D-mpi/src/test.c index 1665ce6..4f27776 100644 --- a/BasicSolver/3D-mpi/src/test.c +++ b/BasicSolver/3D-mpi/src/test.c @@ -1,5 +1,5 @@ /* - * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg. + * Copyright (C) NHR@FAU, University Erlangen-Nuremberg. * All rights reserved. * Use of this source code is governed by a MIT-style * license that can be found in the LICENSE file. @@ -26,10 +26,21 @@ void testInit(Solver* s) for (int k = 0; k < kmaxLocal + 2; k++) { for (int j = 0; j < jmaxLocal + 2; j++) { for (int i = 0; i < imaxLocal + 2; i++) { - G(p, i, j, k) = myrank; - G(f, i, j, k) = myrank; - G(g, i, j, k) = myrank; - G(h, i, j, k) = myrank; + G(p, i, j, k) = 10.0; + G(f, i, j, k) = myrank + 1.0; + G(g, i, j, k) = myrank + 1.0; + G(h, i, j, k) = myrank + 1.0; + } + } + } + + for (int k = 1; k < kmaxLocal + 1; k++) { + for (int j = 1; j < jmaxLocal + 1; j++) { + for (int i = 1; i < imaxLocal + 1; i++) { + G(p, i, j, k) = myrank + 1.0; + G(f, i, j, k) = myrank + 1.0; + G(g, i, j, k) = myrank + 1.0; + G(h, i, j, k) = myrank + 1.0; } } } diff --git a/BasicSolver/3D-mpi/src/test.h b/BasicSolver/3D-mpi/src/test.h index ecefe4e..779b509 100644 --- a/BasicSolver/3D-mpi/src/test.h +++ b/BasicSolver/3D-mpi/src/test.h @@ -1,5 +1,5 @@ /* - * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg. + * 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. diff --git a/BasicSolver/3D-mpi/src/timing.c b/BasicSolver/3D-mpi/src/timing.c index e578f18..78b01c4 100644 --- a/BasicSolver/3D-mpi/src/timing.c +++ b/BasicSolver/3D-mpi/src/timing.c @@ -1,5 +1,5 @@ /* - * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg. + * Copyright (C) NHR@FAU, University Erlangen-Nuremberg. * All rights reserved. * Use of this source code is governed by a MIT-style * license that can be found in the LICENSE file. diff --git a/BasicSolver/3D-mpi/src/timing.h b/BasicSolver/3D-mpi/src/timing.h index 58fb5ac..ed05a8c 100644 --- a/BasicSolver/3D-mpi/src/timing.h +++ b/BasicSolver/3D-mpi/src/timing.h @@ -1,5 +1,5 @@ /* - * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg. + * Copyright (C) NHR@FAU, University Erlangen-Nuremberg. * All rights reserved. * Use of this source code is governed by a MIT-style * license that can be found in the LICENSE file. diff --git a/BasicSolver/3D-mpi/src/util.h b/BasicSolver/3D-mpi/src/util.h index 657b009..c27b2ba 100644 --- a/BasicSolver/3D-mpi/src/util.h +++ b/BasicSolver/3D-mpi/src/util.h @@ -1,5 +1,5 @@ /* - * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg. + * Copyright (C) NHR@FAU, University Erlangen-Nuremberg. * All rights reserved. * Use of this source code is governed by a MIT-style * license that can be found in the LICENSE file. diff --git a/BasicSolver/3D-mpi-io/src/vtkWriter.c b/BasicSolver/3D-mpi/src/vtkWriter-mpi.c similarity index 84% rename from BasicSolver/3D-mpi-io/src/vtkWriter.c rename to BasicSolver/3D-mpi/src/vtkWriter-mpi.c index fe9dc3f..7ff124d 100644 --- a/BasicSolver/3D-mpi-io/src/vtkWriter.c +++ b/BasicSolver/3D-mpi/src/vtkWriter-mpi.c @@ -1,5 +1,5 @@ /* - * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg. + * 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. @@ -14,24 +14,32 @@ #include "comm.h" #include "vtkWriter.h" +// reset fileview for output of string headers static void resetFileview(VtkOptions* o) { - // reset fileview for output of string header MPI_Offset disp; + MPI_File_sync(o->fh); MPI_Barrier(o->comm.comm); MPI_File_get_size(o->fh, &disp); MPI_File_set_view(o->fh, disp, MPI_CHAR, MPI_CHAR, "native", MPI_INFO_NULL); - // printf("Rank %d disp %lld\n", o->comm.rank, disp); +} + +static void writeVersion(VtkOptions* o) +{ + char header[50] = "# vtk DataFile Version 3.0\n"; + // always overwrite exiting files + MPI_File_set_view(o->fh, 0, MPI_CHAR, MPI_CHAR, "native", MPI_INFO_NULL); + + if (commIsMaster(&o->comm)) { + MPI_File_write(o->fh, header, (int)strlen(header), MPI_CHAR, MPI_STATUS_IGNORE); + } } static void writeHeader(VtkOptions* o) { - const size_t MAX_HEADER = 200; - - char* header = (char*)malloc(MAX_HEADER); + char header[400]; char* cursor = header; - cursor += sprintf(cursor, "# vtk DataFile Version 3.0\n"); cursor += sprintf(cursor, "PAMPI cfd solver output\n"); cursor += sprintf(cursor, "BINARY\n"); @@ -52,13 +60,10 @@ static void writeHeader(VtkOptions* o) o->grid.imax * o->grid.jmax * o->grid.kmax); if (commIsMaster(&o->comm)) { - MPI_File_write(o->fh, header, strlen(header), MPI_CHAR, MPI_STATUS_IGNORE); + MPI_File_write(o->fh, header, (int)strlen(header), MPI_CHAR, MPI_STATUS_IGNORE); } - - free(header); } -// TODO Check inputs for length void vtkOpen(VtkOptions* o, char* problem) { char filename[50]; @@ -66,7 +71,7 @@ void vtkOpen(VtkOptions* o, char* problem) snprintf(filename, 50, "%s-p%d.vtk", problem, o->comm.size); MPI_File_open(o->comm.comm, filename, - MPI_MODE_WRONLY | MPI_MODE_CREATE | MPI_MODE_EXCL, + MPI_MODE_WRONLY | MPI_MODE_CREATE, MPI_INFO_NULL, &o->fh); @@ -74,35 +79,32 @@ void vtkOpen(VtkOptions* o, char* problem) printf("Writing VTK output for %s\n", problem); } - resetFileview(o); + writeVersion(o); writeHeader(o); } void vtkScalar(VtkOptions* o, char* name, double* s) { resetFileview(o); - if (commIsMaster(&o->comm)) printf("Register scalar %s\n", name); - const size_t MAX_HEADER = 100; - char* header = (char*)malloc(MAX_HEADER); + char header[100]; char* cursor = header; cursor += sprintf(cursor, "SCALARS %s double 1\n", name); cursor += sprintf(cursor, "LOOKUP_TABLE default\n"); if (commIsMaster(&o->comm)) { - MPI_File_write(o->fh, header, strlen(header), MPI_CHAR, MPI_STATUS_IGNORE); + MPI_File_write(o->fh, header, (int)strlen(header), MPI_CHAR, MPI_STATUS_IGNORE); } - free(header); - int offsets[NDIMS]; - commGetOffsets(&o->comm, offsets, o->grid.imax, o->grid.jmax, o->grid.kmax); + commGetOffsets(&o->comm, offsets, o->grid.kmax, o->grid.jmax, o->grid.imax); // set global view in file MPI_Offset disp; MPI_Datatype fileViewType; + MPI_File_sync(o->fh); MPI_Barrier(o->comm.comm); MPI_File_get_size(o->fh, &disp); @@ -172,15 +174,16 @@ void vtkVector(VtkOptions* o, char* name, VtkVector vec) resetFileview(o); if (commIsMaster(&o->comm)) { - MPI_File_write(o->fh, header, strlen(header), MPI_CHAR, MPI_STATUS_IGNORE); + MPI_File_write(o->fh, header, (int)strlen(header), MPI_CHAR, MPI_STATUS_IGNORE); } int offsets[NDIMS]; - commGetOffsets(&o->comm, offsets, o->grid.imax, o->grid.jmax, o->grid.kmax); + commGetOffsets(&o->comm, offsets, o->grid.kmax, o->grid.jmax, o->grid.imax); // set global view in file MPI_Offset disp; MPI_Datatype fileViewType, vectorType; + MPI_File_sync(o->fh); MPI_Barrier(o->comm.comm); MPI_File_get_size(o->fh, &disp); diff --git a/BasicSolver/3D-mpi/src/vtkWriter.c b/BasicSolver/3D-mpi/src/vtkWriter-seq.c similarity index 98% rename from BasicSolver/3D-mpi/src/vtkWriter.c rename to BasicSolver/3D-mpi/src/vtkWriter-seq.c index 3433b3a..7bd8a14 100644 --- a/BasicSolver/3D-mpi/src/vtkWriter.c +++ b/BasicSolver/3D-mpi/src/vtkWriter-seq.c @@ -1,5 +1,5 @@ /* - * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg. + * 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. diff --git a/BasicSolver/3D-mpi/src/vtkWriter.h b/BasicSolver/3D-mpi/src/vtkWriter.h index d166f62..4b4c2fe 100644 --- a/BasicSolver/3D-mpi/src/vtkWriter.h +++ b/BasicSolver/3D-mpi/src/vtkWriter.h @@ -1,5 +1,5 @@ /* - * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg. + * 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. @@ -8,14 +8,20 @@ #define __VTKWRITER_H_ #include +#include "comm.h" #include "grid.h" typedef enum VtkFormat { ASCII = 0, BINARY } VtkFormat; typedef struct VtkOptions { - VtkFormat fmt; Grid grid; +#ifdef _VTK_WRITER_MPI + MPI_File fh; +#else FILE* fh; + VtkFormat fmt; +#endif // _VTK_WRITER_MPI + Comm comm; } VtkOptions; typedef struct VtkVector {