Compare commits
No commits in common. "cc2aef66b456ed87251e2213ef0731892f79afa8" and "4e8535f4a10f03090651c887cf9cf47850b84b7e" have entirely different histories.
cc2aef66b4
...
4e8535f4a1
71
BasicSolver/3D-mpi-io/Makefile
Normal file
71
BasicSolver/3D-mpi-io/Makefile
Normal file
@ -0,0 +1,71 @@
|
||||
#=======================================================================================
|
||||
# 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)
|
||||
78
BasicSolver/3D-mpi-io/README.md
Normal file
78
BasicSolver/3D-mpi-io/README.md
Normal file
@ -0,0 +1,78 @@
|
||||
# C source skeleton
|
||||
|
||||
## Build
|
||||
|
||||
1. Configure the toolchain and additional options in `config.mk`:
|
||||
```
|
||||
# Supported: GCC, CLANG, ICC
|
||||
TAG ?= GCC
|
||||
ENABLE_OPENMP ?= false
|
||||
|
||||
OPTIONS += -DARRAY_ALIGNMENT=64
|
||||
#OPTIONS += -DVERBOSE
|
||||
#OPTIONS += -DVERBOSE_AFFINITY
|
||||
#OPTIONS += -DVERBOSE_DATASIZE
|
||||
#OPTIONS += -DVERBOSE_TIMER
|
||||
```
|
||||
|
||||
The verbosity options enable detailed output about solver, affinity settings, allocation sizes and timer resolution.
|
||||
For debugging you may want to set the VERBOSE option:
|
||||
```
|
||||
# Supported: GCC, CLANG, ICC
|
||||
TAG ?= GCC
|
||||
ENABLE_OPENMP ?= false
|
||||
|
||||
OPTIONS += -DARRAY_ALIGNMENT=64
|
||||
OPTIONS += -DVERBOSE
|
||||
#OPTIONS += -DVERBOSE_AFFINITY
|
||||
#OPTIONS += -DVERBOSE_DATASIZE
|
||||
#OPTIONS += -DVERBOSE_TIMER
|
||||
`
|
||||
|
||||
2. Build with:
|
||||
```
|
||||
make
|
||||
```
|
||||
|
||||
You can build multiple toolchains in the same directory, but notice that the Makefile is only acting on the one currently set.
|
||||
Intermediate build results are located in the `<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
|
||||
```
|
||||
52
BasicSolver/3D-mpi-io/canal.par
Normal file
52
BasicSolver/3D-mpi-io/canal.par
Normal file
@ -0,0 +1,52 @@
|
||||
#==============================================================================
|
||||
# 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
|
||||
#===============================================================================
|
||||
13
BasicSolver/3D-mpi-io/config.mk
Normal file
13
BasicSolver/3D-mpi-io/config.mk
Normal file
@ -0,0 +1,13 @@
|
||||
# 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
|
||||
52
BasicSolver/3D-mpi-io/dcavity.par
Normal file
52
BasicSolver/3D-mpi-io/dcavity.par
Normal file
@ -0,0 +1,52 @@
|
||||
#==============================================================================
|
||||
# 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
|
||||
#===============================================================================
|
||||
21
BasicSolver/3D-mpi-io/include_CLANG.mk
Normal file
21
BasicSolver/3D-mpi-io/include_CLANG.mk
Normal file
@ -0,0 +1,21 @@
|
||||
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 =
|
||||
14
BasicSolver/3D-mpi-io/include_GCC.mk
Normal file
14
BasicSolver/3D-mpi-io/include_GCC.mk
Normal file
@ -0,0 +1,14 @@
|
||||
CC = gcc
|
||||
GCC = gcc
|
||||
LINKER = $(CC)
|
||||
|
||||
ifeq ($(ENABLE_OPENMP),true)
|
||||
OPENMP = -fopenmp
|
||||
endif
|
||||
|
||||
VERSION = --version
|
||||
CFLAGS = -Ofast -ffreestanding -std=c99 $(OPENMP)
|
||||
LFLAGS = $(OPENMP)
|
||||
DEFINES = -D_GNU_SOURCE
|
||||
INCLUDES =
|
||||
LIBS =
|
||||
14
BasicSolver/3D-mpi-io/include_ICC.mk
Normal file
14
BasicSolver/3D-mpi-io/include_ICC.mk
Normal file
@ -0,0 +1,14 @@
|
||||
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 =
|
||||
38
BasicSolver/3D-mpi-io/src/allocate.c
Normal file
38
BasicSolver/3D-mpi-io/src/allocate.c
Normal file
@ -0,0 +1,38 @@
|
||||
/*
|
||||
* 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;
|
||||
}
|
||||
13
BasicSolver/3D-mpi-io/src/allocate.h
Normal file
13
BasicSolver/3D-mpi-io/src/allocate.h
Normal file
@ -0,0 +1,13 @@
|
||||
/*
|
||||
* 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
|
||||
359
BasicSolver/3D-mpi-io/src/comm.c
Normal file
359
BasicSolver/3D-mpi-io/src/comm.c
Normal file
@ -0,0 +1,359 @@
|
||||
/*
|
||||
* 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
|
||||
}
|
||||
51
BasicSolver/3D-mpi-io/src/comm.h
Normal file
51
BasicSolver/3D-mpi-io/src/comm.h
Normal file
@ -0,0 +1,51 @@
|
||||
/*
|
||||
* 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_
|
||||
16
BasicSolver/3D-mpi-io/src/grid.h
Normal file
16
BasicSolver/3D-mpi-io/src/grid.h
Normal file
@ -0,0 +1,16 @@
|
||||
/*
|
||||
* 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_
|
||||
54
BasicSolver/3D-mpi-io/src/likwid-marker.h
Normal file
54
BasicSolver/3D-mpi-io/src/likwid-marker.h
Normal file
@ -0,0 +1,54 @@
|
||||
/*
|
||||
* =======================================================================================
|
||||
*
|
||||
* Author: Jan Eitzinger (je), jan.eitzinger@fau.de
|
||||
* Copyright (c) 2020 RRZE, University Erlangen-Nuremberg
|
||||
*
|
||||
* Permission is hereby granted, free of charge, to any person obtaining a copy
|
||||
* of this software and associated documentation files (the "Software"), to deal
|
||||
* in the Software without restriction, including without limitation the rights
|
||||
* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
|
||||
* copies of the Software, and to permit persons to whom the Software is
|
||||
* furnished to do so, subject to the following conditions:
|
||||
*
|
||||
* The above copyright notice and this permission notice shall be included in all
|
||||
* copies or substantial portions of the Software.
|
||||
*
|
||||
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
|
||||
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
|
||||
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
|
||||
* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
|
||||
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
|
||||
* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
|
||||
* SOFTWARE.
|
||||
*
|
||||
* =======================================================================================
|
||||
*/
|
||||
#ifndef LIKWID_MARKERS_H
|
||||
#define LIKWID_MARKERS_H
|
||||
|
||||
#ifdef LIKWID_PERFMON
|
||||
#include <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*/
|
||||
126
BasicSolver/3D-mpi-io/src/parameter.c
Normal file
126
BasicSolver/3D-mpi-io/src/parameter.c
Normal file
@ -0,0 +1,126 @@
|
||||
/*
|
||||
* 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);
|
||||
}
|
||||
26
BasicSolver/3D-mpi-io/src/parameter.h
Normal file
26
BasicSolver/3D-mpi-io/src/parameter.h
Normal file
@ -0,0 +1,26 @@
|
||||
/*
|
||||
* 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
|
||||
51
BasicSolver/3D-mpi-io/src/progress.c
Normal file
51
BasicSolver/3D-mpi-io/src/progress.c
Normal file
@ -0,0 +1,51 @@
|
||||
/*
|
||||
* 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);
|
||||
}
|
||||
14
BasicSolver/3D-mpi-io/src/progress.h
Normal file
14
BasicSolver/3D-mpi-io/src/progress.h
Normal file
@ -0,0 +1,14 @@
|
||||
/*
|
||||
* 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
|
||||
853
BasicSolver/3D-mpi-io/src/solver.c
Normal file
853
BasicSolver/3D-mpi-io/src/solver.c
Normal file
@ -0,0 +1,853 @@
|
||||
/*
|
||||
* 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;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
45
BasicSolver/3D-mpi-io/src/solver.h
Normal file
45
BasicSolver/3D-mpi-io/src/solver.h
Normal file
@ -0,0 +1,45 @@
|
||||
/*
|
||||
* 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
|
||||
129
BasicSolver/3D-mpi-io/src/test.c
Normal file
129
BasicSolver/3D-mpi-io/src/test.c
Normal file
@ -0,0 +1,129 @@
|
||||
/*
|
||||
* 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);
|
||||
}
|
||||
13
BasicSolver/3D-mpi-io/src/test.h
Normal file
13
BasicSolver/3D-mpi-io/src/test.h
Normal file
@ -0,0 +1,13 @@
|
||||
/*
|
||||
* 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_
|
||||
22
BasicSolver/3D-mpi-io/src/timing.c
Normal file
22
BasicSolver/3D-mpi-io/src/timing.c
Normal file
@ -0,0 +1,22 @@
|
||||
/*
|
||||
* 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;
|
||||
}
|
||||
13
BasicSolver/3D-mpi-io/src/timing.h
Normal file
13
BasicSolver/3D-mpi-io/src/timing.h
Normal file
@ -0,0 +1,13 @@
|
||||
/*
|
||||
* 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_
|
||||
22
BasicSolver/3D-mpi-io/src/util.h
Normal file
22
BasicSolver/3D-mpi-io/src/util.h
Normal file
@ -0,0 +1,22 @@
|
||||
/*
|
||||
* 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,5 +1,5 @@
|
||||
/*
|
||||
* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
|
||||
* 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.
|
||||
@ -14,32 +14,24 @@
|
||||
#include "comm.h"
|
||||
#include "vtkWriter.h"
|
||||
|
||||
// reset fileview for output of string headers
|
||||
static void resetFileview(VtkOptions* o)
|
||||
{
|
||||
// reset fileview for output of string header
|
||||
MPI_Offset disp;
|
||||
MPI_File_sync(o->fh);
|
||||
MPI_Barrier(o->comm.comm);
|
||||
MPI_File_get_size(o->fh, &disp);
|
||||
MPI_File_set_view(o->fh, disp, MPI_CHAR, MPI_CHAR, "native", MPI_INFO_NULL);
|
||||
}
|
||||
|
||||
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);
|
||||
}
|
||||
// printf("Rank %d disp %lld\n", o->comm.rank, disp);
|
||||
}
|
||||
|
||||
static void writeHeader(VtkOptions* o)
|
||||
{
|
||||
char header[400];
|
||||
const size_t MAX_HEADER = 200;
|
||||
|
||||
char* header = (char*)malloc(MAX_HEADER);
|
||||
char* cursor = header;
|
||||
|
||||
cursor += sprintf(cursor, "# vtk DataFile Version 3.0\n");
|
||||
cursor += sprintf(cursor, "PAMPI cfd solver output\n");
|
||||
cursor += sprintf(cursor, "BINARY\n");
|
||||
|
||||
@ -60,10 +52,13 @@ static void writeHeader(VtkOptions* o)
|
||||
o->grid.imax * o->grid.jmax * o->grid.kmax);
|
||||
|
||||
if (commIsMaster(&o->comm)) {
|
||||
MPI_File_write(o->fh, header, (int)strlen(header), MPI_CHAR, MPI_STATUS_IGNORE);
|
||||
MPI_File_write(o->fh, header, strlen(header), MPI_CHAR, MPI_STATUS_IGNORE);
|
||||
}
|
||||
|
||||
free(header);
|
||||
}
|
||||
|
||||
// TODO Check inputs for length
|
||||
void vtkOpen(VtkOptions* o, char* problem)
|
||||
{
|
||||
char filename[50];
|
||||
@ -71,7 +66,7 @@ void vtkOpen(VtkOptions* o, char* problem)
|
||||
snprintf(filename, 50, "%s-p%d.vtk", problem, o->comm.size);
|
||||
MPI_File_open(o->comm.comm,
|
||||
filename,
|
||||
MPI_MODE_WRONLY | MPI_MODE_CREATE,
|
||||
MPI_MODE_WRONLY | MPI_MODE_CREATE | MPI_MODE_EXCL,
|
||||
MPI_INFO_NULL,
|
||||
&o->fh);
|
||||
|
||||
@ -79,31 +74,34 @@ void vtkOpen(VtkOptions* o, char* problem)
|
||||
printf("Writing VTK output for %s\n", problem);
|
||||
}
|
||||
|
||||
writeVersion(o);
|
||||
resetFileview(o);
|
||||
writeHeader(o);
|
||||
}
|
||||
|
||||
void vtkScalar(VtkOptions* o, char* name, double* s)
|
||||
{
|
||||
resetFileview(o);
|
||||
if (commIsMaster(&o->comm)) printf("Register scalar %s\n", name);
|
||||
|
||||
char header[100];
|
||||
if (commIsMaster(&o->comm)) printf("Register scalar %s\n", name);
|
||||
const size_t MAX_HEADER = 100;
|
||||
|
||||
char* header = (char*)malloc(MAX_HEADER);
|
||||
char* cursor = header;
|
||||
|
||||
cursor += sprintf(cursor, "SCALARS %s double\n", name);
|
||||
|
||||
if (commIsMaster(&o->comm)) {
|
||||
MPI_File_write(o->fh, header, (int)strlen(header), MPI_CHAR, MPI_STATUS_IGNORE);
|
||||
MPI_File_write(o->fh, header, strlen(header), MPI_CHAR, MPI_STATUS_IGNORE);
|
||||
}
|
||||
|
||||
free(header);
|
||||
|
||||
int offsets[NDIMS];
|
||||
commGetOffsets(&o->comm, offsets, o->grid.kmax, o->grid.jmax, o->grid.imax);
|
||||
commGetOffsets(&o->comm, offsets, o->grid.imax, o->grid.jmax, o->grid.kmax);
|
||||
|
||||
// set global view in file
|
||||
MPI_Offset disp;
|
||||
MPI_Datatype fileViewType;
|
||||
MPI_File_sync(o->fh);
|
||||
MPI_Barrier(o->comm.comm);
|
||||
MPI_File_get_size(o->fh, &disp);
|
||||
|
||||
@ -173,16 +171,15 @@ void vtkVector(VtkOptions* o, char* name, VtkVector vec)
|
||||
|
||||
resetFileview(o);
|
||||
if (commIsMaster(&o->comm)) {
|
||||
MPI_File_write(o->fh, header, (int)strlen(header), MPI_CHAR, MPI_STATUS_IGNORE);
|
||||
MPI_File_write(o->fh, header, strlen(header), MPI_CHAR, MPI_STATUS_IGNORE);
|
||||
}
|
||||
|
||||
int offsets[NDIMS];
|
||||
commGetOffsets(&o->comm, offsets, o->grid.kmax, o->grid.jmax, o->grid.imax);
|
||||
commGetOffsets(&o->comm, offsets, o->grid.imax, o->grid.jmax, o->grid.kmax);
|
||||
|
||||
// set global view in file
|
||||
MPI_Offset disp;
|
||||
MPI_Datatype fileViewType, vectorType;
|
||||
MPI_File_sync(o->fh);
|
||||
MPI_Barrier(o->comm.comm);
|
||||
MPI_File_get_size(o->fh, &disp);
|
||||
|
||||
29
BasicSolver/3D-mpi-io/src/vtkWriter.h
Normal file
29
BasicSolver/3D-mpi-io/src/vtkWriter.h
Normal file
@ -0,0 +1,29 @@
|
||||
/*
|
||||
* 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) NHR@FAU, University Erlangen-Nuremberg.
|
||||
# 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.
|
||||
@ -18,18 +18,13 @@ include $(MAKE_DIR)/include_$(TAG).mk
|
||||
INCLUDES += -I$(SRC_DIR) -I$(BUILD_DIR)
|
||||
|
||||
VPATH = $(SRC_DIR)
|
||||
SRC = $(filter-out $(wildcard $(SRC_DIR)/*-*.c),$(wildcard $(SRC_DIR)/*.c))
|
||||
SRC = $(wildcard $(SRC_DIR)/*.c)
|
||||
ASM = $(patsubst $(SRC_DIR)/%.c, $(BUILD_DIR)/%.s, $(SRC))
|
||||
OBJ = $(patsubst $(SRC_DIR)/%.c, $(BUILD_DIR)/%.o, $(SRC))
|
||||
OBJ += $(BUILD_DIR)/vtkWriter-$(VTK_OUTPUT_FMT).o
|
||||
SOURCES = $(SRC) $(wildcard $(SRC_DIR)/*.h)
|
||||
ifeq ($(VTK_OUTPUT_FMT),mpi)
|
||||
DEFINES += -D_VTK_WRITER_MPI
|
||||
endif
|
||||
|
||||
CPPFLAGS := $(CPPFLAGS) $(DEFINES) $(OPTIONS) $(INCLUDES)
|
||||
|
||||
${TARGET}: sanity-checks $(BUILD_DIR) $(OBJ)
|
||||
${TARGET}: $(BUILD_DIR) $(OBJ)
|
||||
$(info ===> LINKING $(TARGET))
|
||||
$(Q)${LINKER} ${LFLAGS} -o $(TARGET) $(OBJ) $(LIBS)
|
||||
|
||||
@ -70,14 +65,6 @@ format:
|
||||
done
|
||||
@echo "Done"
|
||||
|
||||
sanity-checks:
|
||||
ifeq ($(VTK_OUTPUT_FMT),mpi)
|
||||
ifeq ($(ENABLE_MPI),false)
|
||||
$(error VTK_OUTPUT_FMT mpi only supported for ENABLE_MPI true!)
|
||||
endif
|
||||
endif
|
||||
|
||||
|
||||
$(BUILD_DIR):
|
||||
@mkdir $(BUILD_DIR)
|
||||
|
||||
|
||||
@ -2,8 +2,6 @@
|
||||
TAG ?= CLANG
|
||||
ENABLE_MPI ?= true
|
||||
ENABLE_OPENMP ?= false
|
||||
# Supported: seq, mpi
|
||||
VTK_OUTPUT_FMT ?= seq
|
||||
|
||||
#Feature options
|
||||
OPTIONS += -DARRAY_ALIGNMENT=64
|
||||
|
||||
@ -18,4 +18,4 @@ VERSION = --version
|
||||
CFLAGS = -Ofast -std=c17
|
||||
LFLAGS = $(OPENMP) -lm
|
||||
DEFINES += -D_GNU_SOURCE# -DDEBUG
|
||||
INCLUDES = -I/opt/homebrew/include
|
||||
INCLUDES =
|
||||
|
||||
@ -1,10 +1,4 @@
|
||||
ifeq ($(ENABLE_MPI),true)
|
||||
CC = mpicc
|
||||
DEFINES = -D_MPI
|
||||
else
|
||||
CC = gcc
|
||||
endif
|
||||
|
||||
GCC = gcc
|
||||
LINKER = $(CC)
|
||||
|
||||
@ -15,6 +9,6 @@ endif
|
||||
VERSION = --version
|
||||
CFLAGS = -Ofast -ffreestanding -std=c99 $(OPENMP)
|
||||
LFLAGS = $(OPENMP)
|
||||
DEFINES += -D_GNU_SOURCE
|
||||
DEFINES = -D_GNU_SOURCE
|
||||
INCLUDES =
|
||||
LIBS =
|
||||
|
||||
@ -1,5 +1,5 @@
|
||||
/*
|
||||
* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
|
||||
* 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.
|
||||
|
||||
@ -1,5 +1,5 @@
|
||||
/*
|
||||
* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
|
||||
* 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.
|
||||
|
||||
@ -1,5 +1,5 @@
|
||||
/*
|
||||
* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
|
||||
* 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.
|
||||
@ -296,30 +296,6 @@ void commShift(Comm* c, double* f, double* g, double* h)
|
||||
#endif
|
||||
}
|
||||
|
||||
void commGetOffsets(Comm* c, int offsets[], int kmax, int jmax, int imax)
|
||||
{
|
||||
#if defined(_MPI)
|
||||
int sum = 0;
|
||||
|
||||
for (int i = 0; i < c->coords[ICORD]; i++) {
|
||||
sum += sizeOfRank(i, c->dims[ICORD], imax);
|
||||
}
|
||||
offsets[IDIM] = sum;
|
||||
sum = 0;
|
||||
|
||||
for (int i = 0; i < c->coords[JCORD]; i++) {
|
||||
sum += sizeOfRank(i, c->dims[JCORD], jmax);
|
||||
}
|
||||
offsets[JDIM] = sum;
|
||||
sum = 0;
|
||||
|
||||
for (int i = 0; i < c->coords[KCORD]; i++) {
|
||||
sum += sizeOfRank(i, c->dims[KCORD], kmax);
|
||||
}
|
||||
offsets[KDIM] = sum;
|
||||
#endif
|
||||
}
|
||||
|
||||
#define G(v, i, j, k) \
|
||||
v[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)]
|
||||
|
||||
|
||||
@ -1,5 +1,5 @@
|
||||
/*
|
||||
* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
|
||||
* 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.
|
||||
@ -44,7 +44,6 @@ extern void commExchange(Comm*, double*);
|
||||
extern void commShift(Comm* c, double* f, double* g, double* h);
|
||||
extern void commReduction(double* v, int op);
|
||||
extern int commIsBoundary(Comm* c, Direction direction);
|
||||
extern void commGetOffsets(Comm* c, int offsets[], int kmax, int jmax, int imax);
|
||||
extern void commCollectResult(Comm* c,
|
||||
double* ug,
|
||||
double* vg,
|
||||
|
||||
@ -1,5 +1,5 @@
|
||||
/*
|
||||
* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
|
||||
* 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.
|
||||
|
||||
@ -1,9 +1,11 @@
|
||||
/*
|
||||
* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
|
||||
* 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>
|
||||
@ -72,13 +74,6 @@ int main(int argc, char** argv)
|
||||
}
|
||||
|
||||
timeStart = getTimeStamp();
|
||||
#ifdef _VTK_WRITER_MPI
|
||||
VtkOptions opts = { .grid = s.grid, .comm = s.comm };
|
||||
vtkOpen(&opts, s.problem);
|
||||
vtkScalar(&opts, "pressure", s.p);
|
||||
vtkVector(&opts, "velocity", (VtkVector) { s.u, s.v, s.w });
|
||||
vtkClose(&opts);
|
||||
#else
|
||||
double *pg, *ug, *vg, *wg;
|
||||
|
||||
if (commIsMaster(&s.comm)) {
|
||||
@ -110,8 +105,6 @@ int main(int argc, char** argv)
|
||||
vtkVector(&opts, "velocity", (VtkVector) { ug, vg, wg });
|
||||
vtkClose(&opts);
|
||||
}
|
||||
#endif
|
||||
|
||||
timeStop = getTimeStamp();
|
||||
|
||||
if (commIsMaster(&s.comm)) {
|
||||
|
||||
@ -1,5 +1,5 @@
|
||||
/*
|
||||
* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
|
||||
* 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.
|
||||
|
||||
@ -1,5 +1,5 @@
|
||||
/*
|
||||
* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
|
||||
* 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.
|
||||
|
||||
@ -1,5 +1,5 @@
|
||||
/*
|
||||
* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
|
||||
* 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.
|
||||
|
||||
@ -1,5 +1,5 @@
|
||||
/*
|
||||
* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
|
||||
* 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.
|
||||
|
||||
@ -1,5 +1,5 @@
|
||||
/*
|
||||
* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
|
||||
* 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.
|
||||
|
||||
@ -1,5 +1,5 @@
|
||||
/*
|
||||
* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
|
||||
* 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.
|
||||
|
||||
@ -1,5 +1,5 @@
|
||||
/*
|
||||
* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
|
||||
* 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.
|
||||
@ -26,21 +26,10 @@ void testInit(Solver* s)
|
||||
for (int k = 0; k < kmaxLocal + 2; k++) {
|
||||
for (int j = 0; j < jmaxLocal + 2; j++) {
|
||||
for (int i = 0; i < imaxLocal + 2; i++) {
|
||||
G(p, i, j, k) = 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;
|
||||
G(p, i, j, k) = myrank;
|
||||
G(f, i, j, k) = myrank;
|
||||
G(g, i, j, k) = myrank;
|
||||
G(h, i, j, k) = myrank;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@ -1,5 +1,5 @@
|
||||
/*
|
||||
* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
|
||||
* 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.
|
||||
|
||||
@ -1,5 +1,5 @@
|
||||
/*
|
||||
* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
|
||||
* 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.
|
||||
|
||||
@ -1,5 +1,5 @@
|
||||
/*
|
||||
* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
|
||||
* 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.
|
||||
|
||||
@ -1,5 +1,5 @@
|
||||
/*
|
||||
* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
|
||||
* 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.
|
||||
|
||||
@ -1,5 +1,5 @@
|
||||
/*
|
||||
* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
|
||||
* 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.
|
||||
@ -1,5 +1,5 @@
|
||||
/*
|
||||
* Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
|
||||
* 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.
|
||||
@ -8,20 +8,14 @@
|
||||
#define __VTKWRITER_H_
|
||||
#include <stdio.h>
|
||||
|
||||
#include "comm.h"
|
||||
#include "grid.h"
|
||||
|
||||
typedef enum VtkFormat { ASCII = 0, BINARY } VtkFormat;
|
||||
|
||||
typedef struct VtkOptions {
|
||||
Grid grid;
|
||||
#ifdef _VTK_WRITER_MPI
|
||||
MPI_File fh;
|
||||
#else
|
||||
FILE* fh;
|
||||
VtkFormat fmt;
|
||||
#endif // _VTK_WRITER_MPI
|
||||
Comm comm;
|
||||
Grid grid;
|
||||
FILE* fh;
|
||||
} VtkOptions;
|
||||
|
||||
typedef struct VtkVector {
|
||||
|
||||
Loading…
Reference in New Issue
Block a user