Merge 3D mpi versions
This commit is contained in:
parent
d53da89b95
commit
91cbb29f6e
@ -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)
|
|
||||||
@ -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 `<TOOLCHAIN>` 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 `<TOOLCHAIN>` 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
|
|
||||||
```
|
|
||||||
@ -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
|
|
||||||
#===============================================================================
|
|
||||||
@ -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
|
|
||||||
@ -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
|
|
||||||
#===============================================================================
|
|
||||||
@ -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 =
|
|
||||||
@ -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 =
|
|
||||||
@ -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 =
|
|
||||||
@ -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 <errno.h>
|
|
||||||
#include <stddef.h>
|
|
||||||
#include <stdio.h>
|
|
||||||
#include <stdlib.h>
|
|
||||||
|
|
||||||
#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;
|
|
||||||
}
|
|
||||||
@ -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 <stdlib.h>
|
|
||||||
|
|
||||||
extern void* allocate(size_t alignment, size_t bytesize);
|
|
||||||
|
|
||||||
#endif
|
|
||||||
@ -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 <mpi.h>
|
|
||||||
#endif
|
|
||||||
#include <stddef.h>
|
|
||||||
#include <stdio.h>
|
|
||||||
#include <stdlib.h>
|
|
||||||
|
|
||||||
#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
|
|
||||||
}
|
|
||||||
@ -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 <mpi.h>
|
|
||||||
#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_
|
|
||||||
@ -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_
|
|
||||||
@ -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 <likwid.h>
|
|
||||||
#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*/
|
|
||||||
@ -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 <float.h>
|
|
||||||
#include <limits.h>
|
|
||||||
#include <stdio.h>
|
|
||||||
#include <stdlib.h>
|
|
||||||
#include <unistd.h>
|
|
||||||
|
|
||||||
#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 <configFile>\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;
|
|
||||||
}
|
|
||||||
@ -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 <stdio.h>
|
|
||||||
#include <stdlib.h>
|
|
||||||
#include <string.h>
|
|
||||||
|
|
||||||
#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);
|
|
||||||
}
|
|
||||||
@ -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
|
|
||||||
@ -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 <math.h>
|
|
||||||
#include <stdio.h>
|
|
||||||
#include <stdlib.h>
|
|
||||||
#include <string.h>
|
|
||||||
|
|
||||||
#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);
|
|
||||||
}
|
|
||||||
@ -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
|
|
||||||
@ -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 <float.h>
|
|
||||||
#include <math.h>
|
|
||||||
#include <stdio.h>
|
|
||||||
#include <stdlib.h>
|
|
||||||
#include <string.h>
|
|
||||||
|
|
||||||
#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;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
@ -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
|
|
||||||
@ -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 <stdio.h>
|
|
||||||
#include <string.h>
|
|
||||||
|
|
||||||
#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);
|
|
||||||
}
|
|
||||||
@ -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_
|
|
||||||
@ -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 <stdlib.h>
|
|
||||||
#include <time.h>
|
|
||||||
|
|
||||||
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;
|
|
||||||
}
|
|
||||||
@ -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_
|
|
||||||
@ -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_
|
|
||||||
@ -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 <mpi.h>
|
|
||||||
#include <stdio.h>
|
|
||||||
|
|
||||||
#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_
|
|
||||||
@ -1,5 +1,5 @@
|
|||||||
#=======================================================================================
|
#=======================================================================================
|
||||||
# Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
|
# Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
|
||||||
# All rights reserved.
|
# All rights reserved.
|
||||||
# Use of this source code is governed by a MIT-style
|
# Use of this source code is governed by a MIT-style
|
||||||
# license that can be found in the LICENSE file.
|
# 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)
|
INCLUDES += -I$(SRC_DIR) -I$(BUILD_DIR)
|
||||||
|
|
||||||
VPATH = $(SRC_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))
|
ASM = $(patsubst $(SRC_DIR)/%.c, $(BUILD_DIR)/%.s, $(SRC))
|
||||||
OBJ = $(patsubst $(SRC_DIR)/%.c, $(BUILD_DIR)/%.o, $(SRC))
|
OBJ = $(patsubst $(SRC_DIR)/%.c, $(BUILD_DIR)/%.o, $(SRC))
|
||||||
|
OBJ += $(BUILD_DIR)/vtkWriter-$(VTK_OUTPUT_FMT).o
|
||||||
SOURCES = $(SRC) $(wildcard $(SRC_DIR)/*.h)
|
SOURCES = $(SRC) $(wildcard $(SRC_DIR)/*.h)
|
||||||
|
ifeq ($(VTK_OUTPUT_FMT),mpi)
|
||||||
|
DEFINES += -D_VTK_WRITER_MPI
|
||||||
|
endif
|
||||||
|
|
||||||
CPPFLAGS := $(CPPFLAGS) $(DEFINES) $(OPTIONS) $(INCLUDES)
|
CPPFLAGS := $(CPPFLAGS) $(DEFINES) $(OPTIONS) $(INCLUDES)
|
||||||
|
|
||||||
${TARGET}: $(BUILD_DIR) $(OBJ)
|
${TARGET}: sanity-checks $(BUILD_DIR) $(OBJ)
|
||||||
$(info ===> LINKING $(TARGET))
|
$(info ===> LINKING $(TARGET))
|
||||||
$(Q)${LINKER} ${LFLAGS} -o $(TARGET) $(OBJ) $(LIBS)
|
$(Q)${LINKER} ${LFLAGS} -o $(TARGET) $(OBJ) $(LIBS)
|
||||||
|
|
||||||
@ -65,6 +70,14 @@ format:
|
|||||||
done
|
done
|
||||||
@echo "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):
|
$(BUILD_DIR):
|
||||||
@mkdir $(BUILD_DIR)
|
@mkdir $(BUILD_DIR)
|
||||||
|
|
||||||
|
|||||||
@ -2,6 +2,8 @@
|
|||||||
TAG ?= CLANG
|
TAG ?= CLANG
|
||||||
ENABLE_MPI ?= true
|
ENABLE_MPI ?= true
|
||||||
ENABLE_OPENMP ?= false
|
ENABLE_OPENMP ?= false
|
||||||
|
# Supported: seq, mpi
|
||||||
|
VTK_OUTPUT_FMT ?= seq
|
||||||
|
|
||||||
#Feature options
|
#Feature options
|
||||||
OPTIONS += -DARRAY_ALIGNMENT=64
|
OPTIONS += -DARRAY_ALIGNMENT=64
|
||||||
|
|||||||
@ -18,4 +18,4 @@ VERSION = --version
|
|||||||
CFLAGS = -Ofast -std=c17
|
CFLAGS = -Ofast -std=c17
|
||||||
LFLAGS = $(OPENMP) -lm
|
LFLAGS = $(OPENMP) -lm
|
||||||
DEFINES += -D_GNU_SOURCE# -DDEBUG
|
DEFINES += -D_GNU_SOURCE# -DDEBUG
|
||||||
INCLUDES =
|
INCLUDES = -I/opt/homebrew/include
|
||||||
|
|||||||
@ -1,4 +1,10 @@
|
|||||||
|
ifeq ($(ENABLE_MPI),true)
|
||||||
|
CC = mpicc
|
||||||
|
DEFINES = -D_MPI
|
||||||
|
else
|
||||||
CC = gcc
|
CC = gcc
|
||||||
|
endif
|
||||||
|
|
||||||
GCC = gcc
|
GCC = gcc
|
||||||
LINKER = $(CC)
|
LINKER = $(CC)
|
||||||
|
|
||||||
@ -9,6 +15,6 @@ endif
|
|||||||
VERSION = --version
|
VERSION = --version
|
||||||
CFLAGS = -Ofast -ffreestanding -std=c99 $(OPENMP)
|
CFLAGS = -Ofast -ffreestanding -std=c99 $(OPENMP)
|
||||||
LFLAGS = $(OPENMP)
|
LFLAGS = $(OPENMP)
|
||||||
DEFINES = -D_GNU_SOURCE
|
DEFINES += -D_GNU_SOURCE
|
||||||
INCLUDES =
|
INCLUDES =
|
||||||
LIBS =
|
LIBS =
|
||||||
|
|||||||
@ -1,4 +1,10 @@
|
|||||||
|
ifeq ($(ENABLE_MPI),true)
|
||||||
|
CC = mpiicc
|
||||||
|
DEFINES = -D_MPI
|
||||||
|
else
|
||||||
CC = icc
|
CC = icc
|
||||||
|
endif
|
||||||
|
|
||||||
GCC = gcc
|
GCC = gcc
|
||||||
LINKER = $(CC)
|
LINKER = $(CC)
|
||||||
|
|
||||||
@ -9,6 +15,6 @@ endif
|
|||||||
VERSION = --version
|
VERSION = --version
|
||||||
CFLAGS = -O3 -xHost -qopt-zmm-usage=high -std=c99 $(OPENMP)
|
CFLAGS = -O3 -xHost -qopt-zmm-usage=high -std=c99 $(OPENMP)
|
||||||
LFLAGS = $(OPENMP)
|
LFLAGS = $(OPENMP)
|
||||||
DEFINES = -D_GNU_SOURCE
|
DEFINES += -D_GNU_SOURCE
|
||||||
INCLUDES =
|
INCLUDES =
|
||||||
LIBS =
|
LIBS =
|
||||||
|
|||||||
@ -1,5 +1,5 @@
|
|||||||
/*
|
/*
|
||||||
* Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
|
* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
|
||||||
* All rights reserved.
|
* All rights reserved.
|
||||||
* Use of this source code is governed by a MIT-style
|
* Use of this source code is governed by a MIT-style
|
||||||
* license that can be found in the LICENSE file.
|
* license that can be found in the LICENSE file.
|
||||||
|
|||||||
@ -1,5 +1,5 @@
|
|||||||
/*
|
/*
|
||||||
* Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
|
* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
|
||||||
* All rights reserved.
|
* All rights reserved.
|
||||||
* Use of this source code is governed by a MIT-style
|
* Use of this source code is governed by a MIT-style
|
||||||
* license that can be found in the LICENSE file.
|
* license that can be found in the LICENSE file.
|
||||||
|
|||||||
@ -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.
|
* All rights reserved. This file is part of nusif-solver.
|
||||||
* Use of this source code is governed by a MIT style
|
* Use of this source code is governed by a MIT style
|
||||||
* license that can be found in the LICENSE file.
|
* license that can be found in the LICENSE file.
|
||||||
@ -296,6 +296,30 @@ void commShift(Comm* c, double* f, double* g, double* h)
|
|||||||
#endif
|
#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) \
|
#define G(v, i, j, k) \
|
||||||
v[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)]
|
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 k = 1; k < kmaxLocal + 1; k++) {
|
||||||
for (int j = 1; j < jmaxLocal + 1; j++) {
|
for (int j = 1; j < jmaxLocal + 1; j++) {
|
||||||
for (int i = 1; i < imaxLocal + 1; i++) {
|
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 k = 1; k < kmaxLocal + 1; k++) {
|
||||||
for (int j = 1; j < jmaxLocal + 1; j++) {
|
for (int j = 1; j < jmaxLocal + 1; j++) {
|
||||||
for (int i = 1; i < imaxLocal + 1; i++) {
|
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 k = 1; k < kmaxLocal + 1; k++) {
|
||||||
for (int j = 1; j < jmaxLocal + 1; j++) {
|
for (int j = 1; j < jmaxLocal + 1; j++) {
|
||||||
for (int i = 1; i < imaxLocal + 1; i++) {
|
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 k = 1; k < kmaxLocal + 1; k++) {
|
||||||
for (int j = 1; j < jmaxLocal + 1; j++) {
|
for (int j = 1; j < jmaxLocal + 1; j++) {
|
||||||
for (int i = 1; i < imaxLocal + 1; i++) {
|
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;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|||||||
@ -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.
|
* All rights reserved. This file is part of nusif-solver.
|
||||||
* Use of this source code is governed by a MIT style
|
* Use of this source code is governed by a MIT style
|
||||||
* license that can be found in the LICENSE file.
|
* 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 commShift(Comm* c, double* f, double* g, double* h);
|
||||||
extern void commReduction(double* v, int op);
|
extern void commReduction(double* v, int op);
|
||||||
extern int commIsBoundary(Comm* c, Direction direction);
|
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,
|
extern void commCollectResult(Comm* c,
|
||||||
double* ug,
|
double* ug,
|
||||||
double* vg,
|
double* vg,
|
||||||
|
|||||||
@ -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.
|
* All rights reserved. This file is part of nusif-solver.
|
||||||
* Use of this source code is governed by a MIT style
|
* Use of this source code is governed by a MIT style
|
||||||
* license that can be found in the LICENSE file.
|
* license that can be found in the LICENSE file.
|
||||||
|
|||||||
@ -1,11 +1,9 @@
|
|||||||
/*
|
/*
|
||||||
* Copyright (C) 2024 NHR@FAU, University Erlangen-Nuremberg.
|
* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
|
||||||
* All rights reserved.
|
* All rights reserved.
|
||||||
* Use of this source code is governed by a MIT-style
|
* Use of this source code is governed by a MIT-style
|
||||||
* license that can be found in the LICENSE file.
|
* license that can be found in the LICENSE file.
|
||||||
*/
|
*/
|
||||||
#include <float.h>
|
|
||||||
#include <limits.h>
|
|
||||||
#include <stdio.h>
|
#include <stdio.h>
|
||||||
#include <stdlib.h>
|
#include <stdlib.h>
|
||||||
#include <unistd.h>
|
#include <unistd.h>
|
||||||
@ -73,6 +71,13 @@ int main(int argc, char** argv)
|
|||||||
printf("Solution took %.2fs\n", timeStop - timeStart);
|
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;
|
double *pg, *ug, *vg, *wg;
|
||||||
|
|
||||||
if (commIsMaster(&s.comm)) {
|
if (commIsMaster(&s.comm)) {
|
||||||
@ -104,6 +109,7 @@ int main(int argc, char** argv)
|
|||||||
vtkVector(&opts, "velocity", (VtkVector) { ug, vg, wg });
|
vtkVector(&opts, "velocity", (VtkVector) { ug, vg, wg });
|
||||||
vtkClose(&opts);
|
vtkClose(&opts);
|
||||||
}
|
}
|
||||||
|
#endif
|
||||||
|
|
||||||
commFinalize(&s.comm);
|
commFinalize(&s.comm);
|
||||||
return EXIT_SUCCESS;
|
return EXIT_SUCCESS;
|
||||||
|
|||||||
@ -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.
|
* All rights reserved. This file is part of nusif-solver.
|
||||||
* Use of this source code is governed by a MIT style
|
* Use of this source code is governed by a MIT style
|
||||||
* license that can be found in the LICENSE file.
|
* license that can be found in the LICENSE file.
|
||||||
|
|||||||
@ -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.
|
* All rights reserved. This file is part of nusif-solver.
|
||||||
* Use of this source code is governed by a MIT style
|
* Use of this source code is governed by a MIT style
|
||||||
* license that can be found in the LICENSE file.
|
* license that can be found in the LICENSE file.
|
||||||
|
|||||||
@ -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.
|
* All rights reserved. This file is part of nusif-solver.
|
||||||
* Use of this source code is governed by a MIT style
|
* Use of this source code is governed by a MIT style
|
||||||
* license that can be found in the LICENSE file.
|
* license that can be found in the LICENSE file.
|
||||||
|
|||||||
@ -1,5 +1,5 @@
|
|||||||
/*
|
/*
|
||||||
* Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
|
* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
|
||||||
* All rights reserved.
|
* All rights reserved.
|
||||||
* Use of this source code is governed by a MIT-style
|
* Use of this source code is governed by a MIT-style
|
||||||
* license that can be found in the LICENSE file.
|
* license that can be found in the LICENSE file.
|
||||||
|
|||||||
@ -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.
|
* All rights reserved. This file is part of nusif-solver.
|
||||||
* Use of this source code is governed by a MIT style
|
* Use of this source code is governed by a MIT style
|
||||||
* license that can be found in the LICENSE file.
|
* license that can be found in the LICENSE file.
|
||||||
|
|||||||
@ -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.
|
* All rights reserved. This file is part of nusif-solver.
|
||||||
* Use of this source code is governed by a MIT style
|
* Use of this source code is governed by a MIT style
|
||||||
* license that can be found in the LICENSE file.
|
* license that can be found in the LICENSE file.
|
||||||
|
|||||||
@ -1,5 +1,5 @@
|
|||||||
/*
|
/*
|
||||||
* Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
|
* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
|
||||||
* All rights reserved.
|
* All rights reserved.
|
||||||
* Use of this source code is governed by a MIT-style
|
* Use of this source code is governed by a MIT-style
|
||||||
* license that can be found in the LICENSE file.
|
* 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 k = 0; k < kmaxLocal + 2; k++) {
|
||||||
for (int j = 0; j < jmaxLocal + 2; j++) {
|
for (int j = 0; j < jmaxLocal + 2; j++) {
|
||||||
for (int i = 0; i < imaxLocal + 2; i++) {
|
for (int i = 0; i < imaxLocal + 2; i++) {
|
||||||
G(p, i, j, k) = myrank;
|
G(p, i, j, k) = 10.0;
|
||||||
G(f, i, j, k) = myrank;
|
G(f, i, j, k) = myrank + 1.0;
|
||||||
G(g, i, j, k) = myrank;
|
G(g, i, j, k) = myrank + 1.0;
|
||||||
G(h, i, j, k) = myrank;
|
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;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|||||||
@ -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.
|
* All rights reserved. This file is part of nusif-solver.
|
||||||
* Use of this source code is governed by a MIT style
|
* Use of this source code is governed by a MIT style
|
||||||
* license that can be found in the LICENSE file.
|
* license that can be found in the LICENSE file.
|
||||||
|
|||||||
@ -1,5 +1,5 @@
|
|||||||
/*
|
/*
|
||||||
* Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
|
* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
|
||||||
* All rights reserved.
|
* All rights reserved.
|
||||||
* Use of this source code is governed by a MIT-style
|
* Use of this source code is governed by a MIT-style
|
||||||
* license that can be found in the LICENSE file.
|
* license that can be found in the LICENSE file.
|
||||||
|
|||||||
@ -1,5 +1,5 @@
|
|||||||
/*
|
/*
|
||||||
* Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
|
* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
|
||||||
* All rights reserved.
|
* All rights reserved.
|
||||||
* Use of this source code is governed by a MIT-style
|
* Use of this source code is governed by a MIT-style
|
||||||
* license that can be found in the LICENSE file.
|
* license that can be found in the LICENSE file.
|
||||||
|
|||||||
@ -1,5 +1,5 @@
|
|||||||
/*
|
/*
|
||||||
* Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
|
* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
|
||||||
* All rights reserved.
|
* All rights reserved.
|
||||||
* Use of this source code is governed by a MIT-style
|
* Use of this source code is governed by a MIT-style
|
||||||
* license that can be found in the LICENSE file.
|
* license that can be found in the LICENSE file.
|
||||||
|
|||||||
@ -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.
|
* All rights reserved. This file is part of nusif-solver.
|
||||||
* Use of this source code is governed by a MIT style
|
* Use of this source code is governed by a MIT style
|
||||||
* license that can be found in the LICENSE file.
|
* license that can be found in the LICENSE file.
|
||||||
@ -14,24 +14,32 @@
|
|||||||
#include "comm.h"
|
#include "comm.h"
|
||||||
#include "vtkWriter.h"
|
#include "vtkWriter.h"
|
||||||
|
|
||||||
|
// reset fileview for output of string headers
|
||||||
static void resetFileview(VtkOptions* o)
|
static void resetFileview(VtkOptions* o)
|
||||||
{
|
{
|
||||||
// reset fileview for output of string header
|
|
||||||
MPI_Offset disp;
|
MPI_Offset disp;
|
||||||
|
MPI_File_sync(o->fh);
|
||||||
MPI_Barrier(o->comm.comm);
|
MPI_Barrier(o->comm.comm);
|
||||||
MPI_File_get_size(o->fh, &disp);
|
MPI_File_get_size(o->fh, &disp);
|
||||||
MPI_File_set_view(o->fh, disp, MPI_CHAR, MPI_CHAR, "native", MPI_INFO_NULL);
|
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)
|
static void writeHeader(VtkOptions* o)
|
||||||
{
|
{
|
||||||
const size_t MAX_HEADER = 200;
|
char header[400];
|
||||||
|
|
||||||
char* header = (char*)malloc(MAX_HEADER);
|
|
||||||
char* cursor = header;
|
char* cursor = header;
|
||||||
|
|
||||||
cursor += sprintf(cursor, "# vtk DataFile Version 3.0\n");
|
|
||||||
cursor += sprintf(cursor, "PAMPI cfd solver output\n");
|
cursor += sprintf(cursor, "PAMPI cfd solver output\n");
|
||||||
cursor += sprintf(cursor, "BINARY\n");
|
cursor += sprintf(cursor, "BINARY\n");
|
||||||
|
|
||||||
@ -52,13 +60,10 @@ static void writeHeader(VtkOptions* o)
|
|||||||
o->grid.imax * o->grid.jmax * o->grid.kmax);
|
o->grid.imax * o->grid.jmax * o->grid.kmax);
|
||||||
|
|
||||||
if (commIsMaster(&o->comm)) {
|
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)
|
void vtkOpen(VtkOptions* o, char* problem)
|
||||||
{
|
{
|
||||||
char filename[50];
|
char filename[50];
|
||||||
@ -66,7 +71,7 @@ void vtkOpen(VtkOptions* o, char* problem)
|
|||||||
snprintf(filename, 50, "%s-p%d.vtk", problem, o->comm.size);
|
snprintf(filename, 50, "%s-p%d.vtk", problem, o->comm.size);
|
||||||
MPI_File_open(o->comm.comm,
|
MPI_File_open(o->comm.comm,
|
||||||
filename,
|
filename,
|
||||||
MPI_MODE_WRONLY | MPI_MODE_CREATE | MPI_MODE_EXCL,
|
MPI_MODE_WRONLY | MPI_MODE_CREATE,
|
||||||
MPI_INFO_NULL,
|
MPI_INFO_NULL,
|
||||||
&o->fh);
|
&o->fh);
|
||||||
|
|
||||||
@ -74,35 +79,32 @@ void vtkOpen(VtkOptions* o, char* problem)
|
|||||||
printf("Writing VTK output for %s\n", problem);
|
printf("Writing VTK output for %s\n", problem);
|
||||||
}
|
}
|
||||||
|
|
||||||
resetFileview(o);
|
writeVersion(o);
|
||||||
writeHeader(o);
|
writeHeader(o);
|
||||||
}
|
}
|
||||||
|
|
||||||
void vtkScalar(VtkOptions* o, char* name, double* s)
|
void vtkScalar(VtkOptions* o, char* name, double* s)
|
||||||
{
|
{
|
||||||
resetFileview(o);
|
resetFileview(o);
|
||||||
|
|
||||||
if (commIsMaster(&o->comm)) printf("Register scalar %s\n", name);
|
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;
|
char* cursor = header;
|
||||||
|
|
||||||
cursor += sprintf(cursor, "SCALARS %s double 1\n", name);
|
cursor += sprintf(cursor, "SCALARS %s double 1\n", name);
|
||||||
cursor += sprintf(cursor, "LOOKUP_TABLE default\n");
|
cursor += sprintf(cursor, "LOOKUP_TABLE default\n");
|
||||||
|
|
||||||
if (commIsMaster(&o->comm)) {
|
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];
|
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
|
// set global view in file
|
||||||
MPI_Offset disp;
|
MPI_Offset disp;
|
||||||
MPI_Datatype fileViewType;
|
MPI_Datatype fileViewType;
|
||||||
|
MPI_File_sync(o->fh);
|
||||||
MPI_Barrier(o->comm.comm);
|
MPI_Barrier(o->comm.comm);
|
||||||
MPI_File_get_size(o->fh, &disp);
|
MPI_File_get_size(o->fh, &disp);
|
||||||
|
|
||||||
@ -172,15 +174,16 @@ void vtkVector(VtkOptions* o, char* name, VtkVector vec)
|
|||||||
|
|
||||||
resetFileview(o);
|
resetFileview(o);
|
||||||
if (commIsMaster(&o->comm)) {
|
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];
|
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
|
// set global view in file
|
||||||
MPI_Offset disp;
|
MPI_Offset disp;
|
||||||
MPI_Datatype fileViewType, vectorType;
|
MPI_Datatype fileViewType, vectorType;
|
||||||
|
MPI_File_sync(o->fh);
|
||||||
MPI_Barrier(o->comm.comm);
|
MPI_Barrier(o->comm.comm);
|
||||||
MPI_File_get_size(o->fh, &disp);
|
MPI_File_get_size(o->fh, &disp);
|
||||||
|
|
||||||
@ -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.
|
* All rights reserved. This file is part of nusif-solver.
|
||||||
* Use of this source code is governed by a MIT style
|
* Use of this source code is governed by a MIT style
|
||||||
* license that can be found in the LICENSE file.
|
* license that can be found in the LICENSE file.
|
||||||
@ -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.
|
* All rights reserved. This file is part of nusif-solver.
|
||||||
* Use of this source code is governed by a MIT style
|
* Use of this source code is governed by a MIT style
|
||||||
* license that can be found in the LICENSE file.
|
* license that can be found in the LICENSE file.
|
||||||
@ -8,14 +8,20 @@
|
|||||||
#define __VTKWRITER_H_
|
#define __VTKWRITER_H_
|
||||||
#include <stdio.h>
|
#include <stdio.h>
|
||||||
|
|
||||||
|
#include "comm.h"
|
||||||
#include "grid.h"
|
#include "grid.h"
|
||||||
|
|
||||||
typedef enum VtkFormat { ASCII = 0, BINARY } VtkFormat;
|
typedef enum VtkFormat { ASCII = 0, BINARY } VtkFormat;
|
||||||
|
|
||||||
typedef struct VtkOptions {
|
typedef struct VtkOptions {
|
||||||
VtkFormat fmt;
|
|
||||||
Grid grid;
|
Grid grid;
|
||||||
|
#ifdef _VTK_WRITER_MPI
|
||||||
|
MPI_File fh;
|
||||||
|
#else
|
||||||
FILE* fh;
|
FILE* fh;
|
||||||
|
VtkFormat fmt;
|
||||||
|
#endif // _VTK_WRITER_MPI
|
||||||
|
Comm comm;
|
||||||
} VtkOptions;
|
} VtkOptions;
|
||||||
|
|
||||||
typedef struct VtkVector {
|
typedef struct VtkVector {
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user