WIP: Pull Request for a complete Solver package #1
75
BasicSolver/3D-mpi-multigrid/Makefile
Normal file
75
BasicSolver/3D-mpi-multigrid/Makefile
Normal file
@ -0,0 +1,75 @@
|
||||
#=======================================================================================
|
||||
# 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: vis
|
||||
$(info ===> CLEAN)
|
||||
@rm -rf $(BUILD_DIR)
|
||||
@rm -f tags
|
||||
|
||||
vis:
|
||||
$(info ===> REMOVING VIZUALISATION FILES)
|
||||
@rm -f vtk_files/particle*.vtk
|
||||
|
||||
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-multigrid/README.md
Normal file
78
BasicSolver/3D-mpi-multigrid/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
|
||||
```
|
||||
82
BasicSolver/3D-mpi-multigrid/backstep.par
Normal file
82
BasicSolver/3D-mpi-multigrid/backstep.par
Normal file
@ -0,0 +1,82 @@
|
||||
#==============================================================================
|
||||
# Laminar Canal Flow
|
||||
#==============================================================================
|
||||
|
||||
# Problem specific Data:
|
||||
# ---------------------
|
||||
|
||||
name backstep # 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 5000.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 1.0 # initial value for pressure
|
||||
|
||||
# Geometry Data:
|
||||
# -------------
|
||||
|
||||
xlength 7.0 # domain size in x-direction
|
||||
ylength 1.5 # domain size in y-direction
|
||||
zlength 1.0 # domain size in z-direction
|
||||
imax 70 # number of interior cells in x-direction
|
||||
jmax 15 # number of interior cells in y-direction
|
||||
kmax 10 # 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
|
||||
rho 0.52
|
||||
omg 1.7 # relaxation parameter for SOR iteration
|
||||
gamma 0.9 # upwind differencing factor gamma
|
||||
|
||||
# Particle Tracing Data:
|
||||
# -----------------------
|
||||
|
||||
numberOfParticles 500
|
||||
startTime 30
|
||||
injectTimePeriod 1.0
|
||||
writeTimePeriod 0.2
|
||||
|
||||
x1 0.0
|
||||
y1 0.5
|
||||
z1 0.0
|
||||
x2 0.0
|
||||
y2 1.45
|
||||
z2 1.0
|
||||
|
||||
# Obstacle Geometry Data:
|
||||
# -----------------------
|
||||
# Shape 0 disable, 1 Rectangle/Square, 2 Circle
|
||||
|
||||
shape 1
|
||||
xCenter 0.0
|
||||
yCenter 0.0
|
||||
zCenter 0.0
|
||||
xRectLength 2.0
|
||||
yRectLength 1.0
|
||||
zRectLength 2.0
|
||||
circleRadius 1.0
|
||||
|
||||
#===============================================================================
|
||||
File diff suppressed because it is too large
Load Diff
81
BasicSolver/3D-mpi-multigrid/canal.par
Normal file
81
BasicSolver/3D-mpi-multigrid/canal.par
Normal file
@ -0,0 +1,81 @@
|
||||
#==============================================================================
|
||||
# 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 5000.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 30.0 # domain size in x-direction
|
||||
ylength 4.0 # domain size in y-direction
|
||||
zlength 4.0 # domain size in z-direction
|
||||
imax 120 # number of interior cells in x-direction
|
||||
jmax 40 # number of interior cells in y-direction
|
||||
kmax 40 # number of interior cells in z-direction
|
||||
|
||||
# Time Data:
|
||||
# ---------
|
||||
|
||||
te 60.0 # final time
|
||||
dt 0.02 # time stepsize
|
||||
tau 0.5 # safety factor for time stepsize control (<0 constant delt)
|
||||
|
||||
# Pressure Iteration Data:
|
||||
# -----------------------
|
||||
|
||||
itermax 500 # maximal number of pressure iteration in one time step
|
||||
eps 0.0001 # stopping tolerance for pressure iteration
|
||||
omg 0.52
|
||||
omg 1.7 # relaxation parameter for SOR iteration
|
||||
gamma 0.9 # upwind differencing factor gamma
|
||||
|
||||
# Particle Tracing Data:
|
||||
# -----------------------
|
||||
|
||||
numberOfParticles 200
|
||||
startTime 100
|
||||
injectTimePeriod 2.0
|
||||
writeTimePeriod 1.0
|
||||
|
||||
x1 1.0
|
||||
y1 0.0
|
||||
z1 1.0
|
||||
x2 1.0
|
||||
y2 4.0
|
||||
z2 1.0
|
||||
|
||||
# Obstacle Geometry Data:
|
||||
# -----------------------
|
||||
# Shape 0 disable, 1 Rectangle/Square, 2 Circle
|
||||
|
||||
shape 0
|
||||
xCenter 10.0
|
||||
yCenter 2.0
|
||||
zCenter 2.0
|
||||
xRectLength 8.0
|
||||
yRectLength 2.0
|
||||
zRectLength 2.0
|
||||
circleRadius 1.0
|
||||
#===============================================================================
|
||||
12
BasicSolver/3D-mpi-multigrid/config.mk
Normal file
12
BasicSolver/3D-mpi-multigrid/config.mk
Normal file
@ -0,0 +1,12 @@
|
||||
# Supported: GCC, CLANG, ICC
|
||||
TAG ?= ICC
|
||||
ENABLE_OPENMP ?= false
|
||||
|
||||
#Feature options
|
||||
OPTIONS += -DARRAY_ALIGNMENT=64
|
||||
OPTIONS += -DVERBOSE
|
||||
#OPTIONS += -DDEBUG
|
||||
#OPTIONS += -DBOUNDCHECK
|
||||
#OPTIONS += -DVERBOSE_AFFINITY
|
||||
#OPTIONS += -DVERBOSE_DATASIZE
|
||||
#OPTIONS += -DVERBOSE_TIMER
|
||||
81
BasicSolver/3D-mpi-multigrid/dcavity.par
Normal file
81
BasicSolver/3D-mpi-multigrid/dcavity.par
Normal file
@ -0,0 +1,81 @@
|
||||
#==============================================================================
|
||||
# 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 40 # number of interior cells in x-direction
|
||||
jmax 40 # number of interior cells in y-direction
|
||||
kmax 40 # number of interior cells in z-direction
|
||||
|
||||
# Time Data:
|
||||
# ---------
|
||||
|
||||
te 50.0 # final time
|
||||
dt 0.02 # time stepsize
|
||||
tau 0.5 # safety factor for time stepsize control (<0 constant delt)
|
||||
|
||||
# Pressure Iteration Data:
|
||||
# -----------------------
|
||||
|
||||
itermax 1000 # maximal number of pressure iteration in one time step
|
||||
eps 0.001 # stopping tolerance for pressure iteration
|
||||
rho 0.5
|
||||
omg 1.7 # relaxation parameter for SOR iteration
|
||||
gamma 0.9 # upwind differencing factor gamma
|
||||
|
||||
# Particle Tracing Data:
|
||||
# -----------------------
|
||||
|
||||
numberOfParticles 30
|
||||
startTime 10
|
||||
injectTimePeriod 3.0
|
||||
writeTimePeriod 1.0
|
||||
|
||||
x1 0.1
|
||||
y1 0.0
|
||||
z1 1.0
|
||||
x2 1.0
|
||||
y2 4.0
|
||||
z2 1.0
|
||||
|
||||
# Obstacle Geometry Data:
|
||||
# -----------------------
|
||||
# Shape 0 disable, 1 Rectangle/Square, 2 Circle
|
||||
|
||||
shape 1
|
||||
xCenter 10.0
|
||||
yCenter 2.0
|
||||
zCenter 2.0
|
||||
xRectLength 8.0
|
||||
yRectLength 2.0
|
||||
zRectLength 2.0
|
||||
circleRadius 1.0
|
||||
#===============================================================================
|
||||
17
BasicSolver/3D-mpi-multigrid/include_CLANG.mk
Normal file
17
BasicSolver/3D-mpi-multigrid/include_CLANG.mk
Normal file
@ -0,0 +1,17 @@
|
||||
CC = mpicc
|
||||
GCC = cc
|
||||
LINKER = $(CC)
|
||||
|
||||
ifeq ($(ENABLE_OPENMP),true)
|
||||
OPENMP = -fopenmp
|
||||
#OPENMP = -Xpreprocessor -fopenmp #required on Macos with homebrew libomp
|
||||
LIBS = # -lomp
|
||||
endif
|
||||
|
||||
VERSION = --version
|
||||
# CFLAGS = -O3 -std=c17 $(OPENMP)
|
||||
CFLAGS = -Ofast -std=c17
|
||||
#CFLAGS = -Ofast -fnt-store=aggressive -std=c99 $(OPENMP) #AMD CLANG
|
||||
LFLAGS = $(OPENMP) -lm
|
||||
DEFINES = -D_GNU_SOURCE# -DDEBUG
|
||||
INCLUDES =
|
||||
14
BasicSolver/3D-mpi-multigrid/include_GCC.mk
Normal file
14
BasicSolver/3D-mpi-multigrid/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-multigrid/include_ICC.mk
Normal file
14
BasicSolver/3D-mpi-multigrid/include_ICC.mk
Normal file
@ -0,0 +1,14 @@
|
||||
CC = mpiicc
|
||||
GCC = icc
|
||||
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 =
|
||||
82
BasicSolver/3D-mpi-multigrid/karman.par
Normal file
82
BasicSolver/3D-mpi-multigrid/karman.par
Normal file
@ -0,0 +1,82 @@
|
||||
#==============================================================================
|
||||
# Turbulent Canal Flow
|
||||
#==============================================================================
|
||||
|
||||
# Problem specific Data:
|
||||
# ---------------------
|
||||
|
||||
name karman # 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 5050.0 # Reynolds number
|
||||
|
||||
u_init 1.0 # initial value for velocity in x-direction
|
||||
v_init 0.0 # initial value for velocity in y-direction
|
||||
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 8.0 # domain size in y-direction
|
||||
zlength 8.0 # domain size in z-direction
|
||||
imax 200 # number of interior cells in x-direction
|
||||
jmax 80 # number of interior cells in y-direction
|
||||
kmax 80 # number of interior cells in z-direction
|
||||
|
||||
# Time Data:
|
||||
# ---------
|
||||
|
||||
te 250.0 # final time
|
||||
dt 0.02 # time stepsize
|
||||
tau 0.5 # safety factor for time stepsize control (<0 constant delt)
|
||||
|
||||
# Pressure Iteration Data:
|
||||
# -----------------------
|
||||
|
||||
itermax 200 # maximal number of pressure iteration in one time step
|
||||
eps 0.001 # stopping tolerance for pressure iteration
|
||||
rho 0.52
|
||||
omg 1.75 # relaxation parameter for SOR iteration
|
||||
gamma 0.9 # upwind differencing factor gamma
|
||||
|
||||
# Particle Tracing Data:
|
||||
# -----------------------
|
||||
|
||||
numberOfParticles 4000
|
||||
startTime 50
|
||||
injectTimePeriod 1.0
|
||||
writeTimePeriod 5.0
|
||||
|
||||
x1 0.0
|
||||
y1 3.6
|
||||
z1 3.6
|
||||
x2 0.0
|
||||
y2 4.7
|
||||
z2 4.7
|
||||
|
||||
# Obstacle Geometry Data:
|
||||
# -----------------------
|
||||
# Shape 0 disable, 1 Rectangle/Square, 2 Circle
|
||||
|
||||
shape 2
|
||||
xCenter 5.0
|
||||
yCenter 4.0
|
||||
zCenter 4.0
|
||||
xRectLength 8.0
|
||||
yRectLength 2.0
|
||||
zRectLength 2.0
|
||||
circleRadius 2.0
|
||||
|
||||
#===============================================================================
|
||||
35
BasicSolver/3D-mpi-multigrid/src/allocate.c
Normal file
35
BasicSolver/3D-mpi-multigrid/src/allocate.c
Normal file
@ -0,0 +1,35 @@
|
||||
/*
|
||||
* 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 <stdio.h>
|
||||
#include <stdlib.h>
|
||||
|
||||
void* allocate(int alignment, size_t bytesize)
|
||||
{
|
||||
int errorCode;
|
||||
void* ptr;
|
||||
|
||||
errorCode = posix_memalign(&ptr, alignment, bytesize);
|
||||
|
||||
if (errorCode) {
|
||||
if (errorCode == EINVAL) {
|
||||
fprintf(stderr, "Error: Alignment parameter is not a power of two\n");
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
if (errorCode == ENOMEM) {
|
||||
fprintf(stderr, "Error: Insufficient memory to fulfill the request\n");
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
}
|
||||
|
||||
if (ptr == NULL) {
|
||||
fprintf(stderr, "Error: posix_memalign failed!\n");
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
|
||||
return ptr;
|
||||
}
|
||||
13
BasicSolver/3D-mpi-multigrid/src/allocate.h
Normal file
13
BasicSolver/3D-mpi-multigrid/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(int alignment, size_t bytesize);
|
||||
|
||||
#endif
|
||||
544
BasicSolver/3D-mpi-multigrid/src/comm.c
Normal file
544
BasicSolver/3D-mpi-multigrid/src/comm.c
Normal file
@ -0,0 +1,544 @@
|
||||
/*
|
||||
* 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 <mpi.h>
|
||||
#include <stddef.h>
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
|
||||
#include "allocate.h"
|
||||
#include "comm.h"
|
||||
|
||||
// subroutines local to this module
|
||||
static int sizeOfRank(int rank, int size, int N)
|
||||
{
|
||||
return N / size + ((N % size > rank) ? 1 : 0);
|
||||
}
|
||||
|
||||
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]);
|
||||
MPI_Type_create_subarray(NDIMS,
|
||||
sizes,
|
||||
subSizes,
|
||||
starts,
|
||||
MPI_ORDER_C,
|
||||
MPI_INT,
|
||||
&c->rbufferTypesInt[direction]);
|
||||
MPI_Type_commit(&c->rbufferTypesInt[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]);
|
||||
MPI_Type_create_subarray(NDIMS,
|
||||
sizes,
|
||||
subSizes,
|
||||
starts,
|
||||
MPI_ORDER_C,
|
||||
MPI_INT,
|
||||
&c->sbufferTypesInt[direction]);
|
||||
MPI_Type_commit(&c->sbufferTypesInt[direction]);
|
||||
}
|
||||
}
|
||||
|
||||
static void assembleResult(Comm* c,
|
||||
double* src,
|
||||
double* dst,
|
||||
int imaxLocal[],
|
||||
int jmaxLocal[],
|
||||
int kmaxLocal[],
|
||||
int offset[],
|
||||
int kmax,
|
||||
int jmax,
|
||||
int imax)
|
||||
{
|
||||
int numRequests = 1;
|
||||
|
||||
if (c->rank == 0) {
|
||||
numRequests = c->size + 1;
|
||||
}
|
||||
|
||||
MPI_Request requests[numRequests];
|
||||
|
||||
/* all ranks send their interpolated bulk array */
|
||||
MPI_Isend(src,
|
||||
c->imaxLocal * c->jmaxLocal * c->kmaxLocal,
|
||||
MPI_DOUBLE,
|
||||
0,
|
||||
0,
|
||||
c->comm,
|
||||
&requests[0]);
|
||||
|
||||
/* rank 0 assembles the subdomains */
|
||||
if (c->rank == 0) {
|
||||
for (int i = 0; i < c->size; i++) {
|
||||
MPI_Datatype domainType;
|
||||
int oldSizes[NDIMS] = { kmax, jmax, imax };
|
||||
int newSizes[NDIMS] = { kmaxLocal[i], jmaxLocal[i], imaxLocal[i] };
|
||||
int starts[NDIMS] = { offset[i * NDIMS + KDIM],
|
||||
offset[i * NDIMS + JDIM],
|
||||
offset[i * NDIMS + IDIM] };
|
||||
MPI_Type_create_subarray(NDIMS,
|
||||
oldSizes,
|
||||
newSizes,
|
||||
starts,
|
||||
MPI_ORDER_C,
|
||||
MPI_DOUBLE,
|
||||
&domainType);
|
||||
MPI_Type_commit(&domainType);
|
||||
|
||||
MPI_Irecv(dst, 1, domainType, i, 0, c->comm, &requests[i + 1]);
|
||||
MPI_Type_free(&domainType);
|
||||
}
|
||||
}
|
||||
|
||||
MPI_Waitall(numRequests, requests, MPI_STATUSES_IGNORE);
|
||||
}
|
||||
|
||||
static int sum(int* sizes, int position)
|
||||
{
|
||||
int sum = 0;
|
||||
|
||||
for (int i = 0; i < position; i++) {
|
||||
sum += sizes[i];
|
||||
}
|
||||
|
||||
return sum;
|
||||
}
|
||||
|
||||
// exported subroutines
|
||||
void commReduction(double* v, int op)
|
||||
{
|
||||
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);
|
||||
}
|
||||
}
|
||||
|
||||
int commIsBoundary(Comm* c, Direction direction)
|
||||
{
|
||||
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;
|
||||
}
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
void commExchange(Comm* c, double* grid)
|
||||
{
|
||||
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);
|
||||
}
|
||||
|
||||
void commExchangeInt(Comm* c, int* grid)
|
||||
{
|
||||
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->sbufferTypesInt,
|
||||
grid,
|
||||
counts,
|
||||
displs,
|
||||
c->rbufferTypesInt,
|
||||
c->comm);
|
||||
}
|
||||
|
||||
void commShift(Comm* c, double* f, double* g, double* h)
|
||||
{
|
||||
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);
|
||||
}
|
||||
|
||||
#define G(v, i, j, k) \
|
||||
v[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)]
|
||||
|
||||
void commCollectResult(Comm* c,
|
||||
double* ug,
|
||||
double* vg,
|
||||
double* wg,
|
||||
double* pg,
|
||||
double* u,
|
||||
double* v,
|
||||
double* w,
|
||||
double* p,
|
||||
int kmax,
|
||||
int jmax,
|
||||
int imax)
|
||||
{
|
||||
int imaxLocal = c->imaxLocal;
|
||||
int jmaxLocal = c->jmaxLocal;
|
||||
int kmaxLocal = c->kmaxLocal;
|
||||
|
||||
int offset[c->size * NDIMS];
|
||||
int imaxLocalAll[c->size];
|
||||
int jmaxLocalAll[c->size];
|
||||
int kmaxLocalAll[c->size];
|
||||
|
||||
MPI_Gather(&imaxLocal, 1, MPI_INT, imaxLocalAll, 1, MPI_INT, 0, MPI_COMM_WORLD);
|
||||
MPI_Gather(&jmaxLocal, 1, MPI_INT, jmaxLocalAll, 1, MPI_INT, 0, MPI_COMM_WORLD);
|
||||
MPI_Gather(&kmaxLocal, 1, MPI_INT, kmaxLocalAll, 1, MPI_INT, 0, MPI_COMM_WORLD);
|
||||
|
||||
if (c->rank == 0) {
|
||||
for (int i = 0; i < c->size; i++) {
|
||||
int coords[NCORDS];
|
||||
MPI_Cart_coords(c->comm, i, NDIMS, coords);
|
||||
offset[i * NDIMS + IDIM] = sum(imaxLocalAll, coords[ICORD]);
|
||||
offset[i * NDIMS + JDIM] = sum(jmaxLocalAll, coords[JCORD]);
|
||||
offset[i * NDIMS + KDIM] = sum(kmaxLocalAll, coords[KCORD]);
|
||||
printf("Rank: %d, Coords(k,j,i): %d %d %d, Size(k,j,i): %d %d %d, "
|
||||
"Offset(k,j,i): %d %d %d\n",
|
||||
i,
|
||||
coords[KCORD],
|
||||
coords[JCORD],
|
||||
coords[ICORD],
|
||||
kmaxLocalAll[i],
|
||||
jmaxLocalAll[i],
|
||||
imaxLocalAll[i],
|
||||
offset[i * NDIMS + KDIM],
|
||||
offset[i * NDIMS + JDIM],
|
||||
offset[i * NDIMS + IDIM]);
|
||||
}
|
||||
}
|
||||
|
||||
size_t bytesize = imaxLocal * jmaxLocal * kmaxLocal * sizeof(double);
|
||||
double* tmp = allocate(64, bytesize);
|
||||
int idx = 0;
|
||||
|
||||
/* collect P */
|
||||
for (int k = 1; k < kmaxLocal + 1; k++) {
|
||||
for (int j = 1; j < jmaxLocal + 1; j++) {
|
||||
for (int i = 1; i < imaxLocal + 1; i++) {
|
||||
tmp[idx++] = G(p, i, j, k);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
assembleResult(c,
|
||||
tmp,
|
||||
pg,
|
||||
imaxLocalAll,
|
||||
jmaxLocalAll,
|
||||
kmaxLocalAll,
|
||||
offset,
|
||||
kmax,
|
||||
jmax,
|
||||
imax);
|
||||
|
||||
/* collect U */
|
||||
idx = 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++) {
|
||||
tmp[idx++] = (G(u, i, j, k) + G(u, i - 1, j, k)) / 2.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
assembleResult(c,
|
||||
tmp,
|
||||
ug,
|
||||
imaxLocalAll,
|
||||
jmaxLocalAll,
|
||||
kmaxLocalAll,
|
||||
offset,
|
||||
kmax,
|
||||
jmax,
|
||||
imax);
|
||||
|
||||
/* collect V */
|
||||
idx = 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++) {
|
||||
tmp[idx++] = (G(v, i, j, k) + G(v, i, j - 1, k)) / 2.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
assembleResult(c,
|
||||
tmp,
|
||||
vg,
|
||||
imaxLocalAll,
|
||||
jmaxLocalAll,
|
||||
kmaxLocalAll,
|
||||
offset,
|
||||
kmax,
|
||||
jmax,
|
||||
imax);
|
||||
|
||||
/* collect W */
|
||||
idx = 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++) {
|
||||
tmp[idx++] = (G(w, i, j, k) + G(w, i, j, k - 1)) / 2.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
assembleResult(c,
|
||||
tmp,
|
||||
wg,
|
||||
imaxLocalAll,
|
||||
jmaxLocalAll,
|
||||
kmaxLocalAll,
|
||||
offset,
|
||||
kmax,
|
||||
jmax,
|
||||
imax);
|
||||
|
||||
free(tmp);
|
||||
}
|
||||
|
||||
void commPrintConfig(Comm* c)
|
||||
{
|
||||
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);
|
||||
}
|
||||
|
||||
void commInit(Comm* c, int kmax, int jmax, int imax)
|
||||
{
|
||||
/* setup communication */
|
||||
MPI_Comm_rank(MPI_COMM_WORLD, &(c->rank));
|
||||
MPI_Comm_size(MPI_COMM_WORLD, &(c->size));
|
||||
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->coords[KDIM], dims[ICORD], imax);
|
||||
// printf("\nRank : %d\nimaxLocal : %d -> c->coords[IDIM] : %d , dims[ICORD] : %d,
|
||||
// imax : %d\n", c->rank, c->imaxLocal, c->coords[IDIM], dims[ICORD], imax);
|
||||
c->jmaxLocal = sizeOfRank(c->coords[JDIM], dims[JCORD], jmax);
|
||||
// printf("jmaxLocal : %d -> c->coords[JDIM] : %d , dims[JCORD] : %d, imax : %d\n",
|
||||
// c->jmaxLocal, c->coords[JDIM], dims[JCORD], jmax);
|
||||
c->kmaxLocal = sizeOfRank(c->coords[IDIM], dims[KCORD], kmax);
|
||||
// printf("kmaxLocal : %d -> c->coords[KDIM] : %d , dims[KCORD] : %d, imax : %d\n",
|
||||
// c->kmaxLocal, c->coords[KDIM], dims[KCORD], kmax);
|
||||
|
||||
// sizeOfRank(int rank, int size, int N)
|
||||
// { return N / size + ((N % size > rank) ? 1 : 0); }
|
||||
|
||||
// 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);
|
||||
}
|
||||
|
||||
void commFree(Comm* c)
|
||||
{
|
||||
for (int i = 0; i < NDIRS; i++) {
|
||||
MPI_Type_free(&c->sbufferTypes[i]);
|
||||
MPI_Type_free(&c->rbufferTypes[i]);
|
||||
}
|
||||
}
|
||||
63
BasicSolver/3D-mpi-multigrid/src/comm.h
Normal file
63
BasicSolver/3D-mpi-multigrid/src/comm.h
Normal file
@ -0,0 +1,63 @@
|
||||
/*
|
||||
* 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 __COMM_H_
|
||||
#define __COMM_H_
|
||||
#include <mpi.h>
|
||||
/*
|
||||
* 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;
|
||||
MPI_Comm comm;
|
||||
MPI_Datatype sbufferTypes[NDIRS];
|
||||
MPI_Datatype rbufferTypes[NDIRS];
|
||||
MPI_Datatype sbufferTypesInt[NDIRS];
|
||||
MPI_Datatype rbufferTypesInt[NDIRS];
|
||||
int neighbours[NDIRS];
|
||||
int coords[NDIMS], dims[NDIMS];
|
||||
int imaxLocal, jmaxLocal, kmaxLocal;
|
||||
MPI_File fh;
|
||||
} Comm;
|
||||
|
||||
extern void setupCommunication(Comm*, Direction, int);
|
||||
|
||||
extern void commInit(Comm* comm, int kmax, int jmax, int imax);
|
||||
extern void commFree(Comm* comm);
|
||||
extern void commPrintConfig(Comm*);
|
||||
extern void commExchange(Comm*, double*);
|
||||
extern void commExchangeInt(Comm*, int*);
|
||||
|
||||
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 commCollectResult(Comm* c,
|
||||
double* ug,
|
||||
double* vg,
|
||||
double* wg,
|
||||
double* pg,
|
||||
double* u,
|
||||
double* v,
|
||||
double* w,
|
||||
double* p,
|
||||
int kmax,
|
||||
int jmax,
|
||||
int imax);
|
||||
|
||||
static inline int commIsMaster(Comm* c) { return c->rank == 0; }
|
||||
#endif // __COMM_H_
|
||||
16
BasicSolver/3D-mpi-multigrid/src/grid.h
Normal file
16
BasicSolver/3D-mpi-multigrid/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-multigrid/src/likwid-marker.h
Normal file
54
BasicSolver/3D-mpi-multigrid/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*/
|
||||
140
BasicSolver/3D-mpi-multigrid/src/main.c
Normal file
140
BasicSolver/3D-mpi-multigrid/src/main.c
Normal file
@ -0,0 +1,140 @@
|
||||
/*
|
||||
* 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 <float.h>
|
||||
#include <limits.h>
|
||||
#include <mpi.h>
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <unistd.h>
|
||||
|
||||
#include "allocate.h"
|
||||
#include "parameter.h"
|
||||
#include "particletracing.h"
|
||||
#include "progress.h"
|
||||
#include "solver.h"
|
||||
#include "test.h"
|
||||
#include "timing.h"
|
||||
#include "vtkWriter.h"
|
||||
|
||||
|
||||
int main(int argc, char** argv)
|
||||
{
|
||||
int rank;
|
||||
|
||||
double timeStart, timeStop;
|
||||
Parameter params;
|
||||
ParticleTracer particletracer;
|
||||
Solver solver;
|
||||
|
||||
MPI_Init(&argc, &argv);
|
||||
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
|
||||
initParameter(¶ms);
|
||||
|
||||
if (argc < 2) {
|
||||
printf("Usage: %s <configFile>\n", argv[0]);
|
||||
exit(EXIT_SUCCESS);
|
||||
}
|
||||
|
||||
readParameter(¶ms, argv[1]);
|
||||
|
||||
if (commIsMaster(&solver.comm)) {
|
||||
printParameter(¶ms);
|
||||
}
|
||||
|
||||
initSolver(&solver, ¶ms);
|
||||
|
||||
initParticleTracer(&particletracer, ¶ms, &solver);
|
||||
|
||||
|
||||
#ifndef VERBOSE
|
||||
initProgress(solver.te);
|
||||
#endif
|
||||
|
||||
// printf("Rank : %d, imaxLocal : %d, jmaxLocal : %d, kmaxLocal : %d\n", solver.comm.rank, solver.comm.imaxLocal, solver.comm.jmaxLocal, solver.comm.kmaxLocal);
|
||||
// if(solver.comm.rank == 1) printGrid(&solver, solver.seg);
|
||||
// exit(0);
|
||||
|
||||
// printParticleTracerParameters(&particletracer);
|
||||
|
||||
double tau = solver.tau;
|
||||
double te = solver.te;
|
||||
double t = 0.0;
|
||||
int nt = 0;
|
||||
|
||||
timeStart = getTimeStamp();
|
||||
|
||||
while (t <= te) {
|
||||
if (tau > 0.0) computeTimestep(&solver);
|
||||
setBoundaryConditions(&solver);
|
||||
setSpecialBoundaryCondition(&solver);
|
||||
setObjectBoundaryCondition(&solver);
|
||||
|
||||
computeFG(&solver);
|
||||
computeRHS(&solver);
|
||||
// if (nt % 100 == 0) normalizePressure(&solver);
|
||||
multiGrid(&solver);
|
||||
adaptUV(&solver);
|
||||
|
||||
trace(&particletracer, solver.u, solver.v, solver.w, solver.seg, t);
|
||||
|
||||
t += solver.dt;
|
||||
nt++;
|
||||
|
||||
#ifdef VERBOSE
|
||||
if (rank == 0) {
|
||||
printf("TIME %f , TIMESTEP %f\n", t, solver.dt);
|
||||
}
|
||||
#else
|
||||
printProgress(t);
|
||||
#endif
|
||||
}
|
||||
|
||||
timeStop = getTimeStamp();
|
||||
stopProgress();
|
||||
if (commIsMaster(&solver.comm)) {
|
||||
printf("Solution took %.2fs\n", timeStop - timeStart);
|
||||
}
|
||||
|
||||
// testInit(&solver);
|
||||
// commExchange(&solver.comm, solver.p);
|
||||
// testPrintHalo(&solver, solver.p);
|
||||
|
||||
double *pg, *ug, *vg, *wg;
|
||||
|
||||
if (commIsMaster(&solver.comm)) {
|
||||
size_t bytesize = solver.grid.imax * solver.grid.jmax * solver.grid.kmax *
|
||||
sizeof(double);
|
||||
|
||||
pg = allocate(64, bytesize);
|
||||
ug = allocate(64, bytesize);
|
||||
vg = allocate(64, bytesize);
|
||||
wg = allocate(64, bytesize);
|
||||
}
|
||||
|
||||
commCollectResult(&solver.comm,
|
||||
ug,
|
||||
vg,
|
||||
wg,
|
||||
pg,
|
||||
solver.u,
|
||||
solver.v,
|
||||
solver.w,
|
||||
solver.p,
|
||||
solver.grid.kmax,
|
||||
solver.grid.jmax,
|
||||
solver.grid.imax);
|
||||
|
||||
VtkOptions opts = { .grid = solver.grid, .comm = solver.comm };
|
||||
vtkOpen(&opts, solver.problem);
|
||||
vtkScalar(&opts, "pressure", pg);
|
||||
vtkVector(&opts, "velocity", (VtkVector) { ug, vg, wg });
|
||||
vtkClose(&opts);
|
||||
|
||||
commFree(&solver.comm);
|
||||
MPI_Finalize();
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
152
BasicSolver/3D-mpi-multigrid/src/parameter.c
Normal file
152
BasicSolver/3D-mpi-multigrid/src/parameter.c
Normal file
@ -0,0 +1,152 @@
|
||||
/*
|
||||
* 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;
|
||||
param->rho = 0.99;
|
||||
param->levels = 5;
|
||||
}
|
||||
|
||||
void readParameter(Parameter* param, const char* filename)
|
||||
{
|
||||
FILE* fp = fopen(filename, "r");
|
||||
char line[MAXLINE];
|
||||
int i;
|
||||
|
||||
if (!fp) {
|
||||
fprintf(stderr, "Could not open parameter file: %s\n", filename);
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
|
||||
while (!feof(fp)) {
|
||||
line[0] = '\0';
|
||||
fgets(line, MAXLINE, fp);
|
||||
for (i = 0; line[i] != '\0' && line[i] != '#'; i++)
|
||||
;
|
||||
line[i] = '\0';
|
||||
|
||||
char* tok = strtok(line, " ");
|
||||
char* val = strtok(NULL, " ");
|
||||
|
||||
#define PARSE_PARAM(p, f) \
|
||||
if (strncmp(tok, #p, sizeof(#p) / sizeof(#p[0]) - 1) == 0) { \
|
||||
param->p = f(val); \
|
||||
}
|
||||
#define PARSE_STRING(p) PARSE_PARAM(p, strdup)
|
||||
#define PARSE_INT(p) PARSE_PARAM(p, atoi)
|
||||
#define PARSE_REAL(p) PARSE_PARAM(p, atof)
|
||||
|
||||
if (tok != NULL && val != NULL) {
|
||||
PARSE_REAL(xlength);
|
||||
PARSE_REAL(ylength);
|
||||
PARSE_REAL(zlength);
|
||||
PARSE_INT(imax);
|
||||
PARSE_INT(jmax);
|
||||
PARSE_INT(kmax);
|
||||
PARSE_INT(levels);
|
||||
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);
|
||||
PARSE_REAL(rho);
|
||||
|
||||
/* Added new particle tracing parameters */
|
||||
PARSE_INT(numberOfParticles);
|
||||
PARSE_REAL(startTime);
|
||||
PARSE_REAL(injectTimePeriod);
|
||||
PARSE_REAL(writeTimePeriod);
|
||||
PARSE_REAL(x1);
|
||||
PARSE_REAL(y1);
|
||||
PARSE_REAL(z1);
|
||||
PARSE_REAL(x2);
|
||||
PARSE_REAL(y2);
|
||||
PARSE_REAL(z2);
|
||||
|
||||
/* Added obstacle geometry parameters */
|
||||
PARSE_INT(shape);
|
||||
PARSE_REAL(xCenter);
|
||||
PARSE_REAL(yCenter);
|
||||
PARSE_REAL(zCenter);
|
||||
PARSE_REAL(xRectLength);
|
||||
PARSE_REAL(yRectLength);
|
||||
PARSE_REAL(zRectLength);
|
||||
PARSE_REAL(circleRadius);
|
||||
}
|
||||
}
|
||||
|
||||
fclose(fp);
|
||||
}
|
||||
|
||||
void printParameter(Parameter* param)
|
||||
{
|
||||
printf("Parameters for %s\n", param->name);
|
||||
printf("Boundary conditions Left:%d Right:%d Bottom:%d Top:%d 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);
|
||||
}
|
||||
35
BasicSolver/3D-mpi-multigrid/src/parameter.h
Normal file
35
BasicSolver/3D-mpi-multigrid/src/parameter.h
Normal file
@ -0,0 +1,35 @@
|
||||
/*
|
||||
* 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, levels;
|
||||
|
||||
double eps, omg, rho;
|
||||
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;
|
||||
|
||||
int numberOfParticles;
|
||||
double startTime, injectTimePeriod, writeTimePeriod;
|
||||
|
||||
double x1, y1, z1, x2, y2, z2;
|
||||
|
||||
int shape;
|
||||
double xCenter, yCenter, zCenter, xRectLength, yRectLength, zRectLength, circleRadius;
|
||||
} Parameter;
|
||||
|
||||
void initParameter(Parameter*);
|
||||
void readParameter(Parameter*, const char*);
|
||||
void printParameter(Parameter*);
|
||||
#endif
|
||||
729
BasicSolver/3D-mpi-multigrid/src/particletracing.c
Normal file
729
BasicSolver/3D-mpi-multigrid/src/particletracing.c
Normal file
@ -0,0 +1,729 @@
|
||||
/*
|
||||
* 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 "particletracing.h"
|
||||
#include "solver.h"
|
||||
|
||||
#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 S(i, j, k) \
|
||||
seg[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)]
|
||||
|
||||
static int ts = 0;
|
||||
|
||||
#define XOFFSET 0
|
||||
#define YOFFSET 1
|
||||
#define ZOFFSET 2
|
||||
#define XOFFSETEND 3
|
||||
#define YOFFSETEND 4
|
||||
#define ZOFFSETEND 5
|
||||
|
||||
static double sum(int* sizes, int size)
|
||||
{
|
||||
double sum = 0;
|
||||
|
||||
for (int i = 0; i < size; ++i) {
|
||||
sum += sizes[i];
|
||||
}
|
||||
|
||||
return sum;
|
||||
}
|
||||
|
||||
static double sumOffset(double* sizes, int init, int offset, int coord)
|
||||
{
|
||||
double sum = 0;
|
||||
|
||||
for (int i = init - offset; coord > 0; i -= offset, --coord) {
|
||||
sum += sizes[i];
|
||||
}
|
||||
|
||||
return sum;
|
||||
}
|
||||
|
||||
void printParticles(ParticleTracer* particletracer)
|
||||
{
|
||||
for (int i = 0; i < particletracer->totalParticles; ++i) {
|
||||
printf("Rank : %d Particle position X : %.2f, Y : %.2f, flag : %d, total pt : "
|
||||
"%d, pointer : %d, xOffset : %.2f, yOffset : %.2f, xOffsetEnd : %.2f, "
|
||||
"yOffsetEnd : %.2f\n",
|
||||
particletracer->rank,
|
||||
particletracer->particlePool[i].x,
|
||||
particletracer->particlePool[i].y,
|
||||
particletracer->particlePool[i].flag,
|
||||
particletracer->totalParticles,
|
||||
particletracer->pointer,
|
||||
particletracer->xOffset,
|
||||
particletracer->yOffset,
|
||||
particletracer->xOffsetEnd,
|
||||
particletracer->yOffsetEnd);
|
||||
}
|
||||
}
|
||||
void injectParticles(ParticleTracer* particletracer, int* restrict seg)
|
||||
{
|
||||
double x, y, z;
|
||||
int imaxLocal = particletracer->imaxLocal;
|
||||
int jmaxLocal = particletracer->jmaxLocal;
|
||||
int kmaxLocal = particletracer->kmaxLocal;
|
||||
for (int i = 0; i < particletracer->numberOfParticles; ++i) {
|
||||
x = particletracer->linSpaceLine[i].x;
|
||||
y = particletracer->linSpaceLine[i].y;
|
||||
z = particletracer->linSpaceLine[i].z;
|
||||
if (x >= particletracer->xOffset && y >= particletracer->yOffset &&
|
||||
z >= particletracer->zOffset && x <= particletracer->xOffsetEnd &&
|
||||
y <= particletracer->yOffsetEnd && y <= particletracer->zOffsetEnd) {
|
||||
// printf("\nRank : %d\n", particletracer->rank);
|
||||
// printf("\t%.2f >= %.2f && %.2f >= %.2f && %.2f <= %.2f && %.2f <= %.2f\n",x
|
||||
// , particletracer->xOffset ,y , particletracer->yOffset, x ,
|
||||
// particletracer->xOffsetEnd ,y , particletracer->yOffsetEnd);
|
||||
|
||||
particletracer->particlePool[particletracer->pointer].x = x;
|
||||
particletracer->particlePool[particletracer->pointer].y = y;
|
||||
particletracer->particlePool[particletracer->pointer].z = z;
|
||||
|
||||
int i = particletracer->particlePool[particletracer->pointer].x /
|
||||
particletracer->dx;
|
||||
int j = particletracer->particlePool[particletracer->pointer].y /
|
||||
particletracer->dy;
|
||||
int k = particletracer->particlePool[particletracer->pointer].z /
|
||||
particletracer->dz;
|
||||
if (S(i, j, k) == NONE) {
|
||||
particletracer->particlePool[particletracer->pointer].flag = true;
|
||||
++(particletracer->pointer);
|
||||
++(particletracer->totalParticles);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void advanceParticles(ParticleTracer* particletracer,
|
||||
double* restrict u,
|
||||
double* restrict v,
|
||||
double* restrict w,
|
||||
int* restrict seg,
|
||||
double time)
|
||||
{
|
||||
int imax = particletracer->imax;
|
||||
int jmax = particletracer->jmax;
|
||||
int kmax = particletracer->kmax;
|
||||
int imaxLocal = particletracer->imaxLocal;
|
||||
int jmaxLocal = particletracer->jmaxLocal;
|
||||
int kmaxLocal = particletracer->kmaxLocal;
|
||||
|
||||
double dx = particletracer->dx;
|
||||
double dy = particletracer->dy;
|
||||
double dz = particletracer->dz;
|
||||
|
||||
double xlength = particletracer->xlength;
|
||||
double ylength = particletracer->ylength;
|
||||
double zlength = particletracer->zlength;
|
||||
|
||||
Particle buff[particletracer->size][4000];
|
||||
|
||||
for (int i = 0; i < particletracer->size; ++i) {
|
||||
for (int j = 0; j < 4000; ++j) {
|
||||
buff[i][j].x = 0.0;
|
||||
buff[i][j].y = 0.0;
|
||||
buff[i][j].z = 0.0;
|
||||
buff[i][j].flag = false;
|
||||
}
|
||||
}
|
||||
int particleBufIndex[particletracer->size];
|
||||
|
||||
memset(particleBufIndex, 0, sizeof(particleBufIndex));
|
||||
|
||||
for (int i = 0; i < particletracer->totalParticles; ++i) {
|
||||
if (particletracer->particlePool[i].flag == true) {
|
||||
double xTemp = particletracer->particlePool[i].x;
|
||||
double yTemp = particletracer->particlePool[i].y;
|
||||
double zTemp = particletracer->particlePool[i].z;
|
||||
|
||||
double x = xTemp - particletracer->xOffset;
|
||||
double y = yTemp - particletracer->yOffset;
|
||||
double z = zTemp - particletracer->zOffset;
|
||||
|
||||
int iCoord = (int)(x / dx) + 1;
|
||||
int jCoord = (int)((y + 0.5 * dy) / dy) + 1;
|
||||
int kCoord = (int)((z + 0.5 * dz) / dz) + 1;
|
||||
|
||||
double x1 = (double)(iCoord - 1) * dx;
|
||||
double y1 = ((double)(jCoord - 1) - 0.5) * dy;
|
||||
double z1 = ((double)(kCoord - 1) - 0.5) * dz;
|
||||
|
||||
double x2 = (double)iCoord * dx;
|
||||
double y2 = ((double)jCoord - 0.5) * dy;
|
||||
double z2 = ((double)kCoord - 0.5) * dz;
|
||||
|
||||
double u_n =
|
||||
(1.0 / (dx * dy * dz)) *
|
||||
((x2 - x) * (y2 - y) * (z2 - z) * U(iCoord - 1, jCoord - 1, kCoord - 1) +
|
||||
(x - x1) * (y2 - y) * (z2 - z) * U(iCoord, jCoord - 1, kCoord - 1) +
|
||||
(x2 - x) * (y - y1) * (z2 - z) * U(iCoord - 1, jCoord, kCoord - 1) +
|
||||
(x - x1) * (y - y1) * (z2 - z) * U(iCoord, jCoord, kCoord - 1) +
|
||||
(x2 - x) * (y2 - y) * (z - z1) * U(iCoord - 1, jCoord - 1, kCoord) +
|
||||
(x - x1) * (y2 - y) * (z - z1) * U(iCoord, jCoord - 1, kCoord) +
|
||||
(x2 - x) * (y - y1) * (z - z1) * U(iCoord - 1, jCoord, kCoord) +
|
||||
(x - x1) * (y - y1) * (z - z1) * U(iCoord, jCoord, kCoord));
|
||||
|
||||
double new_x = (x + particletracer->xOffset) + particletracer->dt * u_n;
|
||||
particletracer->particlePool[i].x = new_x;
|
||||
|
||||
iCoord = (int)((x + 0.5 * dx) / dx) + 1;
|
||||
jCoord = (int)(y / dy) + 1;
|
||||
kCoord = (int)((z + 0.5 * dz) / dz) + 1;
|
||||
|
||||
x1 = ((double)(iCoord - 1) - 0.5) * dx;
|
||||
y1 = (double)(jCoord - 1) * dy;
|
||||
z1 = ((double)(kCoord - 1) - 0.5) * dz;
|
||||
|
||||
x2 = ((double)iCoord - 0.5) * dx;
|
||||
y2 = (double)jCoord * dy;
|
||||
z2 = ((double)kCoord - 0.5) * dz;
|
||||
|
||||
double v_n =
|
||||
(1.0 / (dx * dy * dz)) *
|
||||
((x2 - x) * (y2 - y) * (z2 - z) * V(iCoord - 1, jCoord - 1, kCoord - 1) +
|
||||
(x - x1) * (y2 - y) * (z2 - z) * V(iCoord, jCoord - 1, kCoord - 1) +
|
||||
(x2 - x) * (y - y1) * (z2 - z) * V(iCoord - 1, jCoord, kCoord - 1) +
|
||||
(x - x1) * (y - y1) * (z2 - z) * V(iCoord, jCoord, kCoord - 1) +
|
||||
(x2 - x) * (y2 - y) * (z - z1) * V(iCoord - 1, jCoord - 1, kCoord) +
|
||||
(x - x1) * (y2 - y) * (z - z1) * V(iCoord, jCoord - 1, kCoord) +
|
||||
(x2 - x) * (y - y1) * (z - z1) * V(iCoord - 1, jCoord, kCoord) +
|
||||
(x - x1) * (y - y1) * (z - z1) * V(iCoord, jCoord, kCoord));
|
||||
|
||||
double new_y = (y + particletracer->yOffset) + particletracer->dt * v_n;
|
||||
particletracer->particlePool[i].y = new_y;
|
||||
|
||||
iCoord = (int)((x + 0.5 * dx) / dx) + 1;
|
||||
jCoord = (int)((y + 0.5 * dy) / dy) + 1;
|
||||
kCoord = (int)(z / dz) + 1;
|
||||
|
||||
x1 = ((double)(iCoord - 1) - 0.5) * dx;
|
||||
y1 = ((double)(jCoord - 1) - 0.5) * dy;
|
||||
z1 = (double)(kCoord - 1) * dz;
|
||||
|
||||
x2 = ((double)iCoord - 0.5) * dx;
|
||||
y2 = ((double)jCoord - 0.5) * dy;
|
||||
z2 = (double)kCoord * dz;
|
||||
|
||||
double w_n =
|
||||
(1.0 / (dx * dy * dz)) *
|
||||
((x2 - x) * (y2 - y) * (z2 - z) * W(iCoord - 1, jCoord - 1, kCoord - 1) +
|
||||
(x - x1) * (y2 - y) * (z2 - z) * W(iCoord, jCoord - 1, kCoord - 1) +
|
||||
(x2 - x) * (y - y1) * (z2 - z) * W(iCoord - 1, jCoord, kCoord - 1) +
|
||||
(x - x1) * (y - y1) * (z2 - z) * W(iCoord, jCoord, kCoord - 1) +
|
||||
(x2 - x) * (y2 - y) * (z - z1) * W(iCoord - 1, jCoord - 1, kCoord) +
|
||||
(x - x1) * (y2 - y) * (z - z1) * W(iCoord, jCoord - 1, kCoord) +
|
||||
(x2 - x) * (y - y1) * (z - z1) * W(iCoord - 1, jCoord, kCoord) +
|
||||
(x - x1) * (y - y1) * (z - z1) * W(iCoord, jCoord, kCoord));
|
||||
|
||||
double new_z = (z + particletracer->zOffset) + particletracer->dt * w_n;
|
||||
particletracer->particlePool[i].z = new_z;
|
||||
|
||||
// printf("Rank : %d\n", particletracer->rank);
|
||||
// printf("\tOld X : %.2f, translated X : %.2f, xOffset : %.2f, New X : %.2f,
|
||||
// iCoord : %d\n\tOld Y : %.2f, translated X : %.2f, yOffset : %.2f, New Y :
|
||||
// %.2f, jCoord : %d\n\n",xTemp, x, particletracer->xOffset, new_x, iCoord,
|
||||
// yTemp, y, particletracer->yOffset , new_y, jCoord); printf("\tU(iCoord - 1,
|
||||
// jCoord - 1) : %.2f, U(iCoord, jCoord - 1) : %.2f, U(iCoord - 1, jCoord) :
|
||||
// %.2f, U(iCoord, jCoord) : %.2f\n", U(iCoord - 1, jCoord - 1), U(iCoord,
|
||||
// jCoord - 1), U(iCoord - 1, jCoord), U(iCoord, jCoord)); printf("\tV(iCoord
|
||||
// - 1, jCoord - 1) : %.2f, V(iCoord, jCoord - 1) : %.2f, V(iCoord - 1,
|
||||
// jCoord) : %.2f, V(iCoord, jCoord) : %.2f\n\n", V(iCoord - 1, jCoord - 1),
|
||||
// V(iCoord, jCoord - 1), V(iCoord - 1, jCoord), V(iCoord, jCoord));
|
||||
// printf("\t U N : %.2f, V N : %.2f\n\n", u_n, v_n);
|
||||
|
||||
if (((new_x < particletracer->xOffset) ||
|
||||
(new_x >= particletracer->xOffsetEnd) ||
|
||||
(new_y < particletracer->yOffset) ||
|
||||
(new_y >= particletracer->yOffsetEnd) ||
|
||||
(new_z < particletracer->zOffset) ||
|
||||
(new_z >= particletracer->zOffsetEnd))) {
|
||||
// New logic to transfer particles to neighbouring ranks or discard the
|
||||
// particle.
|
||||
|
||||
for (int i = 0; i < particletracer->size; ++i) {
|
||||
if ((new_x >=
|
||||
particletracer->offset[i + particletracer->size * XOFFSET]) &&
|
||||
(new_x <= particletracer
|
||||
->offset[i + particletracer->size * XOFFSETEND]) &&
|
||||
(new_y >=
|
||||
particletracer->offset[i + particletracer->size * YOFFSET]) &&
|
||||
(new_y <= particletracer
|
||||
->offset[i + particletracer->size * YOFFSETEND]) &&
|
||||
(new_z >=
|
||||
particletracer->offset[i + particletracer->size * ZOFFSET]) &&
|
||||
(new_z <= particletracer
|
||||
->offset[i + particletracer->size * ZOFFSETEND]) &&
|
||||
i != particletracer->rank) {
|
||||
buff[i][particleBufIndex[i]].x = new_x;
|
||||
buff[i][particleBufIndex[i]].y = new_y;
|
||||
buff[i][particleBufIndex[i]].z = new_z;
|
||||
buff[i][particleBufIndex[i]].flag = true;
|
||||
++particleBufIndex[i];
|
||||
}
|
||||
}
|
||||
particletracer->particlePool[i].flag = false;
|
||||
|
||||
int i_new = new_x / dx, j_new = new_y / dy, k_new = new_z / dz;
|
||||
|
||||
int iOffset = particletracer->xOffset / dx,
|
||||
jOffset = particletracer->yOffset / dy,
|
||||
kOffset = particletracer->zOffset / dz;
|
||||
if (S(i_new - iOffset, j_new - jOffset, k_new - kOffset) != NONE) {
|
||||
particletracer->particlePool[i].flag = false;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
for (int i = 0; i < particletracer->size; ++i) {
|
||||
if (i != particletracer->rank) {
|
||||
MPI_Send(buff[i],
|
||||
4000,
|
||||
particletracer->mpi_particle,
|
||||
i,
|
||||
0,
|
||||
particletracer->comm);
|
||||
}
|
||||
}
|
||||
for (int i = 0; i < particletracer->size; ++i) {
|
||||
if (i != particletracer->rank) {
|
||||
MPI_Recv(buff[i],
|
||||
4000,
|
||||
particletracer->mpi_particle,
|
||||
i,
|
||||
0,
|
||||
particletracer->comm,
|
||||
MPI_STATUS_IGNORE);
|
||||
}
|
||||
}
|
||||
for (int i = 0; i < particletracer->size; ++i) {
|
||||
if (i != particletracer->rank) {
|
||||
for (int j = 0; j < 4000; ++j) {
|
||||
if (buff[i][j].flag == true) {
|
||||
particletracer->particlePool[particletracer->pointer].x = buff[i][j]
|
||||
.x;
|
||||
particletracer->particlePool[particletracer->pointer].y = buff[i][j]
|
||||
.y;
|
||||
particletracer->particlePool[particletracer->pointer].z = buff[i][j]
|
||||
.z;
|
||||
particletracer->particlePool[particletracer->pointer].flag = true;
|
||||
++(particletracer->pointer);
|
||||
++(particletracer->totalParticles);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void freeParticles(ParticleTracer* particletracer)
|
||||
{
|
||||
free(particletracer->particlePool);
|
||||
free(particletracer->linSpaceLine);
|
||||
free(particletracer->offset);
|
||||
}
|
||||
|
||||
void writeParticles(ParticleTracer* particletracer)
|
||||
{
|
||||
int collectedBuffIndex[particletracer->size];
|
||||
|
||||
MPI_Gather(&particletracer->totalParticles,
|
||||
1,
|
||||
MPI_INT,
|
||||
collectedBuffIndex,
|
||||
1,
|
||||
MPI_INT,
|
||||
0,
|
||||
particletracer->comm);
|
||||
|
||||
if (particletracer->rank != 0) {
|
||||
Particle buff[particletracer->totalParticles];
|
||||
for (int i = 0; i < particletracer->totalParticles; ++i) {
|
||||
buff[i].x = particletracer->particlePool[i].x;
|
||||
buff[i].y = particletracer->particlePool[i].y;
|
||||
buff[i].z = particletracer->particlePool[i].z;
|
||||
buff[i].flag = particletracer->particlePool[i].flag;
|
||||
// printf("Rank : %d sending to rank 0 X : %.2f, Y : %.2f with totalpt :
|
||||
// %d\n", particletracer->rank, buff[i].x, buff[i].y,
|
||||
// particletracer->totalParticles);
|
||||
}
|
||||
MPI_Send(buff,
|
||||
particletracer->totalParticles,
|
||||
particletracer->mpi_particle,
|
||||
0,
|
||||
1,
|
||||
particletracer->comm);
|
||||
}
|
||||
if (particletracer->rank == 0) {
|
||||
char filename[50];
|
||||
FILE* fp;
|
||||
|
||||
snprintf(filename, 50, "vtk_files/particles_%d.vtk", ts);
|
||||
fp = fopen(filename, "w");
|
||||
|
||||
if (fp == NULL) {
|
||||
printf("Error!\n");
|
||||
exit(EXIT_FAILURE);
|
||||
}
|
||||
fprintf(fp, "# vtk DataFile Version 3.0\n");
|
||||
fprintf(fp, "PAMPI cfd solver particle tracing file\n");
|
||||
fprintf(fp, "ASCII\n");
|
||||
|
||||
fprintf(fp, "DATASET UNSTRUCTURED_GRID\n");
|
||||
fprintf(fp, "FIELD FieldData 2\n");
|
||||
fprintf(fp, "TIME 1 1 double\n");
|
||||
fprintf(fp, "%d\n", ts);
|
||||
fprintf(fp, "CYCLE 1 1 int\n");
|
||||
fprintf(fp, "1\n");
|
||||
|
||||
int overallTotalParticles = sum(collectedBuffIndex, particletracer->size);
|
||||
|
||||
fprintf(fp, "POINTS %d float\n", overallTotalParticles);
|
||||
|
||||
printf("Total particles : %d\n", overallTotalParticles);
|
||||
|
||||
for (int i = 1; i < particletracer->size; ++i) {
|
||||
Particle recvBuff[collectedBuffIndex[i]];
|
||||
MPI_Recv(&recvBuff,
|
||||
collectedBuffIndex[i],
|
||||
particletracer->mpi_particle,
|
||||
i,
|
||||
1,
|
||||
particletracer->comm,
|
||||
MPI_STATUS_IGNORE);
|
||||
|
||||
for (int j = 0; j < collectedBuffIndex[i]; ++j) {
|
||||
double x = recvBuff[j].x;
|
||||
double y = recvBuff[j].y;
|
||||
double z = recvBuff[j].z;
|
||||
fprintf(fp, "%f %f %f\n", x, y, z);
|
||||
// printf("Rank : 0 receiving from rank %d X : %.2f, Y : %.2f with totalpt
|
||||
// : %d\n", i, x, y, particletracer->totalParticles);
|
||||
}
|
||||
}
|
||||
for (int i = 0; i < particletracer->totalParticles; ++i) {
|
||||
double x = particletracer->particlePool[i].x;
|
||||
double y = particletracer->particlePool[i].y;
|
||||
double z = particletracer->particlePool[i].z;
|
||||
fprintf(fp, "%f %f %f\n", x, y, z);
|
||||
}
|
||||
|
||||
fprintf(fp, "CELLS %d %d\n", overallTotalParticles, 2 * overallTotalParticles);
|
||||
|
||||
for (int i = 0; i < overallTotalParticles; ++i) {
|
||||
fprintf(fp, "1 %d\n", i);
|
||||
}
|
||||
|
||||
fprintf(fp, "CELL_TYPES %d\n", overallTotalParticles);
|
||||
|
||||
for (int i = 0; i < overallTotalParticles; ++i) {
|
||||
fprintf(fp, "1\n");
|
||||
}
|
||||
|
||||
fclose(fp);
|
||||
}
|
||||
|
||||
++ts;
|
||||
}
|
||||
|
||||
void initParticleTracer(ParticleTracer* particletracer, Parameter* params, Solver* solver)
|
||||
{
|
||||
/* initializing local properties from params */
|
||||
|
||||
particletracer->numberOfParticles = params->numberOfParticles;
|
||||
particletracer->startTime = params->startTime;
|
||||
particletracer->injectTimePeriod = params->injectTimePeriod;
|
||||
particletracer->writeTimePeriod = params->writeTimePeriod;
|
||||
|
||||
particletracer->dt = params->dt;
|
||||
particletracer->dx = params->xlength / params->imax;
|
||||
particletracer->dy = params->ylength / params->jmax;
|
||||
particletracer->dz = params->zlength / params->kmax;
|
||||
|
||||
particletracer->xlength = params->xlength;
|
||||
particletracer->ylength = params->ylength;
|
||||
particletracer->zlength = params->zlength;
|
||||
|
||||
particletracer->x1 = params->x1;
|
||||
particletracer->y1 = params->y1;
|
||||
particletracer->z1 = params->z1;
|
||||
particletracer->x2 = params->x2;
|
||||
particletracer->y2 = params->y2;
|
||||
particletracer->z2 = params->z2;
|
||||
|
||||
particletracer->lastInjectTime = params->startTime;
|
||||
particletracer->lastUpdateTime = params->startTime;
|
||||
particletracer->lastWriteTime = params->startTime;
|
||||
|
||||
particletracer->pointer = 0;
|
||||
particletracer->totalParticles = 0;
|
||||
|
||||
particletracer->imax = params->imax;
|
||||
particletracer->jmax = params->jmax;
|
||||
particletracer->kmax = params->kmax;
|
||||
|
||||
particletracer->imaxLocal = solver->comm.imaxLocal;
|
||||
particletracer->jmaxLocal = solver->comm.jmaxLocal;
|
||||
particletracer->kmaxLocal = solver->comm.kmaxLocal;
|
||||
|
||||
particletracer->estimatedNumParticles = (particletracer->imaxLocal *
|
||||
particletracer->jmaxLocal *
|
||||
particletracer->kmaxLocal);
|
||||
|
||||
particletracer->particlePool = malloc(
|
||||
sizeof(Particle) * particletracer->estimatedNumParticles);
|
||||
|
||||
for (int i = 0; i < particletracer->estimatedNumParticles; ++i) {
|
||||
particletracer->particlePool[i].x = 0.0;
|
||||
particletracer->particlePool[i].y = 0.0;
|
||||
particletracer->particlePool[i].z = 0.0;
|
||||
particletracer->particlePool[i].flag = false;
|
||||
}
|
||||
|
||||
particletracer->linSpaceLine = malloc(
|
||||
sizeof(Particle) * particletracer->numberOfParticles);
|
||||
|
||||
/* duplicating communication from solver */
|
||||
MPI_Comm_dup(solver->comm.comm, &particletracer->comm);
|
||||
|
||||
particletracer->rank = solver->comm.rank;
|
||||
particletracer->size = solver->comm.size;
|
||||
|
||||
particletracer->offset = (double*)malloc(sizeof(double) * 4 * particletracer->size);
|
||||
|
||||
memcpy(particletracer->dims, solver->comm.dims, sizeof(solver->comm.dims));
|
||||
|
||||
memcpy(particletracer->coords, solver->comm.coords, sizeof(solver->comm.coords));
|
||||
|
||||
double offset[6][particletracer->size];
|
||||
|
||||
particletracer->xOffset = solver->xOffset;
|
||||
particletracer->yOffset = solver->yOffset;
|
||||
particletracer->zOffset = solver->zOffset;
|
||||
particletracer->xOffsetEnd = solver->xOffsetEnd;
|
||||
particletracer->yOffsetEnd = solver->yOffsetEnd;
|
||||
particletracer->zOffsetEnd = solver->zOffsetEnd;
|
||||
|
||||
// printf("Rank : %d, xOffset : %.2f, yOffset : %.2f, zOffset : %.2f, xOffsetEnd :
|
||||
// %.2f, yOffsetEnd : %.2f, zOffsetEnd : %.2f\n", particletracer->rank,
|
||||
// particletracer->xOffset, particletracer->yOffset, particletracer->zOffset,
|
||||
// particletracer->xOffsetEnd, particletracer->yOffsetEnd,
|
||||
// particletracer->zOffsetEnd);
|
||||
|
||||
MPI_Allgather(&particletracer->xOffset,
|
||||
1,
|
||||
MPI_DOUBLE,
|
||||
offset[0],
|
||||
1,
|
||||
MPI_DOUBLE,
|
||||
particletracer->comm);
|
||||
MPI_Allgather(&particletracer->yOffset,
|
||||
1,
|
||||
MPI_DOUBLE,
|
||||
offset[1],
|
||||
1,
|
||||
MPI_DOUBLE,
|
||||
particletracer->comm);
|
||||
MPI_Allgather(&particletracer->zOffset,
|
||||
1,
|
||||
MPI_DOUBLE,
|
||||
offset[2],
|
||||
1,
|
||||
MPI_DOUBLE,
|
||||
particletracer->comm);
|
||||
MPI_Allgather(&particletracer->xOffsetEnd,
|
||||
1,
|
||||
MPI_DOUBLE,
|
||||
offset[3],
|
||||
1,
|
||||
MPI_DOUBLE,
|
||||
particletracer->comm);
|
||||
MPI_Allgather(&particletracer->yOffsetEnd,
|
||||
1,
|
||||
MPI_DOUBLE,
|
||||
offset[4],
|
||||
1,
|
||||
MPI_DOUBLE,
|
||||
particletracer->comm);
|
||||
MPI_Allgather(&particletracer->zOffsetEnd,
|
||||
1,
|
||||
MPI_DOUBLE,
|
||||
offset[5],
|
||||
1,
|
||||
MPI_DOUBLE,
|
||||
particletracer->comm);
|
||||
|
||||
memcpy(particletracer->offset, offset, sizeof(offset));
|
||||
|
||||
// if(particletracer->rank == 0)
|
||||
// {
|
||||
// for(int i = 0;i < particletracer->size; ++i)
|
||||
// {
|
||||
// printf("Rank : %d, xLocal : %.2f, yLocal : %.2f, zLocal : %.2f\n", i,
|
||||
// xLocal[i], yLocal[i], zLocal[i]);
|
||||
// }
|
||||
// for(int i = 0;i < particletracer->size; ++i)
|
||||
// {
|
||||
// printf("Rank : %d and its xOffset : %.2f, yOffset : %.2f, zOffset :
|
||||
// %.2f,xOffsetEnd : %.2f, yOffsetEnd : %.2f, zOffsetEnd : %.2f\n", i,
|
||||
// particletracer->offset[i + particletracer->size * XOFFSET],
|
||||
// particletracer->offset[i + particletracer->size * YOFFSET],
|
||||
// particletracer->offset[i + particletracer->size * ZOFFSET] ,
|
||||
// particletracer->offset[i + particletracer->size * XOFFSETEND],
|
||||
// particletracer->offset[i + particletracer->size * YOFFSETEND],
|
||||
// particletracer->offset[i + particletracer->size * ZOFFSETEND]);
|
||||
// }
|
||||
// }
|
||||
|
||||
for (int i = 0; i < particletracer->numberOfParticles; ++i) {
|
||||
// double spacing = (double)i /
|
||||
// (double)(particletracer->numberOfParticles - 1);
|
||||
// particletracer->linSpaceLine[i].x = spacing * particletracer->x1 + (1.0 -
|
||||
// spacing) * particletracer->x2; particletracer->linSpaceLine[i].y = spacing *
|
||||
// particletracer->y1 + (1.0 - spacing) * particletracer->y2;
|
||||
// particletracer->linSpaceLine[i].z = spacing * particletracer->z1 + (1.0 -
|
||||
// spacing) * particletracer->z2;
|
||||
|
||||
particletracer->linSpaceLine[i].x = particletracer->x1;
|
||||
|
||||
particletracer->linSpaceLine[i].y = (((double)rand() /
|
||||
RAND_MAX) *
|
||||
(particletracer->y2 - particletracer->y1)) +
|
||||
particletracer->y1;
|
||||
particletracer->linSpaceLine[i].z = (((double)rand() /
|
||||
RAND_MAX) *
|
||||
(particletracer->z2 - particletracer->z1)) +
|
||||
particletracer->z1;
|
||||
particletracer->linSpaceLine[i].flag = true;
|
||||
}
|
||||
|
||||
// Create the mpi_particle datatype
|
||||
MPI_Datatype mpi_particle;
|
||||
int lengths[4] = { 1, 1, 1, 1 };
|
||||
|
||||
MPI_Aint displacements[4];
|
||||
Particle dummy_particle;
|
||||
MPI_Aint base_address;
|
||||
MPI_Get_address(&dummy_particle, &base_address);
|
||||
MPI_Get_address(&dummy_particle.x, &displacements[0]);
|
||||
MPI_Get_address(&dummy_particle.y, &displacements[1]);
|
||||
MPI_Get_address(&dummy_particle.z, &displacements[2]);
|
||||
MPI_Get_address(&dummy_particle.flag, &displacements[3]);
|
||||
displacements[0] = MPI_Aint_diff(displacements[0], base_address);
|
||||
displacements[1] = MPI_Aint_diff(displacements[1], base_address);
|
||||
displacements[2] = MPI_Aint_diff(displacements[2], base_address);
|
||||
displacements[3] = MPI_Aint_diff(displacements[3], base_address);
|
||||
MPI_Datatype types[4] = { MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_C_BOOL };
|
||||
MPI_Type_create_struct(4,
|
||||
lengths,
|
||||
displacements,
|
||||
types,
|
||||
&particletracer->mpi_particle);
|
||||
MPI_Type_commit(&particletracer->mpi_particle);
|
||||
}
|
||||
|
||||
void printParticleTracerParameters(ParticleTracer* particletracer)
|
||||
{
|
||||
printf("Particle Tracing data:\n");
|
||||
printf("Rank : %d\n", particletracer->rank);
|
||||
printf("\tNumber of particles : %d being injected for every period of %.2f\n",
|
||||
particletracer->numberOfParticles,
|
||||
particletracer->injectTimePeriod);
|
||||
printf("\tstartTime : %.2f\n", particletracer->startTime);
|
||||
printf("\t(Line along which the particles are to be injected) \n\tx1 : %.2f, y1 : "
|
||||
"%.2f, z1 : %.2f, x2 : %.2f, y2 : %.2f, z2 : %.2f\n",
|
||||
particletracer->x1,
|
||||
particletracer->y1,
|
||||
particletracer->z1,
|
||||
particletracer->x2,
|
||||
particletracer->y2,
|
||||
particletracer->z2);
|
||||
printf("\tPointer : %d, TotalParticles : %d\n",
|
||||
particletracer->pointer,
|
||||
particletracer->totalParticles);
|
||||
printf("\tdt : %.2f, dx : %.2f, dy : %.2f, dz : %.2f\n",
|
||||
particletracer->dt,
|
||||
particletracer->dx,
|
||||
particletracer->dy,
|
||||
particletracer->dz);
|
||||
printf("\tcoord[0] : %d, coord[1] : %d, coord[2] : %d\n",
|
||||
particletracer->coords[IDIM],
|
||||
particletracer->coords[JDIM],
|
||||
particletracer->coords[KDIM]);
|
||||
printf("\txOffset : %.2f, yOffset : %.2f, zOffset : %.2f\n",
|
||||
particletracer->xOffset,
|
||||
particletracer->yOffset,
|
||||
particletracer->zOffset);
|
||||
printf("\txOffsetEnd : %.2f, yOffsetEnd : %.2f, zOffsetEnd : %.2f\n",
|
||||
particletracer->xOffsetEnd,
|
||||
particletracer->yOffsetEnd,
|
||||
particletracer->zOffsetEnd);
|
||||
printf("\txLocal : %.2f, yLocal : %.2f, zLocal : %.2f\n",
|
||||
particletracer->xLocal,
|
||||
particletracer->yLocal,
|
||||
particletracer->zLocal);
|
||||
}
|
||||
|
||||
void trace(ParticleTracer* particletracer,
|
||||
double* restrict u,
|
||||
double* restrict v,
|
||||
double* restrict w,
|
||||
int* restrict seg,
|
||||
double time)
|
||||
{
|
||||
if (time >= particletracer->startTime) {
|
||||
if ((time - particletracer->lastInjectTime) >= particletracer->injectTimePeriod) {
|
||||
injectParticles(particletracer, seg);
|
||||
particletracer->lastInjectTime = time;
|
||||
}
|
||||
if ((time - particletracer->lastWriteTime) >= particletracer->writeTimePeriod) {
|
||||
writeParticles(particletracer);
|
||||
particletracer->lastWriteTime = time;
|
||||
}
|
||||
advanceParticles(particletracer, u, v, w, seg, time);
|
||||
compress(particletracer);
|
||||
particletracer->lastUpdateTime = time;
|
||||
}
|
||||
}
|
||||
|
||||
void compress(ParticleTracer* particletracer)
|
||||
{
|
||||
Particle* memPool = particletracer->particlePool;
|
||||
Particle tempPool[particletracer->totalParticles];
|
||||
|
||||
for (int i = 0; i < particletracer->totalParticles; ++i) {
|
||||
tempPool[i].x = 0.0;
|
||||
tempPool[i].y = 0.0;
|
||||
tempPool[i].z = 0.0;
|
||||
tempPool[i].flag = false;
|
||||
}
|
||||
int totalParticles = 0;
|
||||
|
||||
for (int i = 0; i < particletracer->totalParticles; ++i) {
|
||||
if (memPool[i].flag == true) {
|
||||
tempPool[totalParticles].x = memPool[i].x;
|
||||
tempPool[totalParticles].y = memPool[i].y;
|
||||
tempPool[totalParticles].z = memPool[i].z;
|
||||
tempPool[totalParticles].flag = memPool[i].flag;
|
||||
++totalParticles;
|
||||
}
|
||||
}
|
||||
|
||||
particletracer->totalParticles = totalParticles;
|
||||
particletracer->pointer = totalParticles;
|
||||
|
||||
memcpy(particletracer->particlePool, tempPool, totalParticles * sizeof(Particle));
|
||||
}
|
||||
62
BasicSolver/3D-mpi-multigrid/src/particletracing.h
Normal file
62
BasicSolver/3D-mpi-multigrid/src/particletracing.h
Normal file
@ -0,0 +1,62 @@
|
||||
/*
|
||||
* 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 __PARTICLETRACING_H_
|
||||
#define __PARTICLETRACING_H_
|
||||
#include "allocate.h"
|
||||
#include "parameter.h"
|
||||
#include "solver.h"
|
||||
#include <mpi.h>
|
||||
#include <stdbool.h>
|
||||
|
||||
#define NDIMS 3
|
||||
|
||||
typedef enum COORD { X = 0, Y, NCOORD } COORD;
|
||||
typedef struct {
|
||||
double x, y, z;
|
||||
bool flag;
|
||||
} Particle;
|
||||
|
||||
typedef struct {
|
||||
int numberOfParticles, totalParticles;
|
||||
double startTime, injectTimePeriod, writeTimePeriod, lastInjectTime, lastUpdateTime,
|
||||
lastWriteTime;
|
||||
|
||||
int estimatedNumParticles;
|
||||
|
||||
double dx, dy, dz, dt;
|
||||
Particle* linSpaceLine;
|
||||
Particle* particlePool;
|
||||
|
||||
int pointer;
|
||||
|
||||
double imax, jmax, kmax, xlength, ylength, zlength, imaxLocal, jmaxLocal, kmaxLocal;
|
||||
|
||||
double x1, y1, z1, x2, y2, z2;
|
||||
|
||||
MPI_Comm comm;
|
||||
MPI_Datatype mpi_particle;
|
||||
|
||||
int rank, size;
|
||||
int coords[NDIMS], dims[NDIMS];
|
||||
|
||||
double xLocal, yLocal, zLocal, xOffset, yOffset, zOffset, xOffsetEnd, yOffsetEnd,
|
||||
zOffsetEnd;
|
||||
|
||||
double* offset;
|
||||
|
||||
} ParticleTracer;
|
||||
|
||||
void initParticleTracer(ParticleTracer*, Parameter*, Solver*);
|
||||
void injectParticles(ParticleTracer*, int*);
|
||||
void advanceParticles(ParticleTracer*, double*, double*, double*, int*, double);
|
||||
void freeParticles(ParticleTracer*);
|
||||
void writeParticles(ParticleTracer*);
|
||||
void printParticleTracerParameters(ParticleTracer*);
|
||||
void printParticles(ParticleTracer*);
|
||||
void trace(ParticleTracer*, double*, double*, double*, int*, double);
|
||||
void compress(ParticleTracer*);
|
||||
#endif
|
||||
51
BasicSolver/3D-mpi-multigrid/src/progress.c
Normal file
51
BasicSolver/3D-mpi-multigrid/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-multigrid/src/progress.h
Normal file
14
BasicSolver/3D-mpi-multigrid/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();
|
||||
|
||||
#endif
|
||||
1818
BasicSolver/3D-mpi-multigrid/src/solver.c
Normal file
1818
BasicSolver/3D-mpi-multigrid/src/solver.c
Normal file
File diff suppressed because it is too large
Load Diff
96
BasicSolver/3D-mpi-multigrid/src/solver.h
Normal file
96
BasicSolver/3D-mpi-multigrid/src/solver.h
Normal file
@ -0,0 +1,96 @@
|
||||
/*
|
||||
* 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 OBJECTBOUNDARY {
|
||||
NONE = 0,
|
||||
/* Front Corners */
|
||||
FRONTTOPLEFTCORNER,
|
||||
FRONTTOPRIGHTCORNER,
|
||||
FRONTBOTTOMLEFTCORNER,
|
||||
FRONTBOTTOMRIGHTCORNER,
|
||||
/* Back Corners */
|
||||
BACKTOPLEFTCORNER,
|
||||
BACKTOPRIGHTCORNER,
|
||||
BACKBOTTOMLEFTCORNER,
|
||||
BACKBOTTOMRIGHTCORNER,
|
||||
/* Faces */
|
||||
FRONTFACE,
|
||||
BACKFACE,
|
||||
LEFTFACE,
|
||||
RIGHTFACE,
|
||||
TOPFACE,
|
||||
BOTTOMFACE,
|
||||
/* Front Lines remaining after Corners and Faces */
|
||||
FRONTLEFTLINE,
|
||||
FRONTRIGHTLINE,
|
||||
FRONTTOPLINE,
|
||||
FRONTBOTTOMLINE,
|
||||
/* Bottom Lines remaining after Corners and Faces */
|
||||
BACKLEFTLINE,
|
||||
BACKRIGHTLINE,
|
||||
BACKTOPLINE,
|
||||
BACKBOTTOMLINE,
|
||||
/* Mid Lines remaining after Corners and Faces */
|
||||
MIDTOPLEFTLINE,
|
||||
MIDTOPRIGHTLINE,
|
||||
MIDBOTTOMLEFTLINE,
|
||||
MIDBOTTOMRIGHTLINE,
|
||||
/* Local where its an object but not a boundary */
|
||||
LOCAL
|
||||
};
|
||||
|
||||
enum BC { NOSLIP = 1, SLIP, OUTFLOW, PERIODIC };
|
||||
enum SHAPE { NOSHAPE = 0, RECT, CIRCLE };
|
||||
|
||||
typedef struct {
|
||||
/* geometry and grid information */
|
||||
Grid grid;
|
||||
/* arrays */
|
||||
double *p, *rhs, **r, **e;
|
||||
double *f, *g, *h;
|
||||
double *u, *v, *w;
|
||||
int* seg;
|
||||
double xLocal, yLocal, zLocal, xOffset, yOffset, zOffset, xOffsetEnd, yOffsetEnd,
|
||||
zOffsetEnd;
|
||||
/* parameters */
|
||||
double eps, omega, rho;
|
||||
double re, tau, gamma;
|
||||
double gx, gy, gz;
|
||||
/* time stepping */
|
||||
int itermax, levels, currentlevel;
|
||||
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 double smoothRB(Solver*);
|
||||
extern void restrictMG(Solver*);
|
||||
extern void prolongate(Solver*);
|
||||
extern void correct(Solver*);
|
||||
extern Solver copySolver(Solver*);
|
||||
extern void multiGrid(Solver*);
|
||||
extern void normalizePressure(Solver*);
|
||||
extern void computeTimestep(Solver*);
|
||||
extern void setBoundaryConditions(Solver*);
|
||||
extern void setSpecialBoundaryCondition(Solver*);
|
||||
extern void setPressureBoundaryCondition(Solver*);
|
||||
extern void setObjectBoundaryCondition(Solver*);
|
||||
extern void setObjectPBoundaryCondition(Solver*);
|
||||
extern void printGrid(Solver*, int*);
|
||||
extern void computeFG(Solver*);
|
||||
extern void adaptUV(Solver*);
|
||||
#endif
|
||||
118
BasicSolver/3D-mpi-multigrid/src/test.c
Normal file
118
BasicSolver/3D-mpi-multigrid/src/test.c
Normal file
@ -0,0 +1,118 @@
|
||||
/*
|
||||
* 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) = myrank;
|
||||
G(f, i, j, k) = myrank;
|
||||
G(g, i, j, k) = myrank;
|
||||
G(h, i, j, k) = myrank;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
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-multigrid/src/test.h
Normal file
13
BasicSolver/3D-mpi-multigrid/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_
|
||||
24
BasicSolver/3D-mpi-multigrid/src/timing.c
Normal file
24
BasicSolver/3D-mpi-multigrid/src/timing.c
Normal file
@ -0,0 +1,24 @@
|
||||
/*
|
||||
* 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()
|
||||
{
|
||||
struct timespec ts;
|
||||
clock_gettime(CLOCK_MONOTONIC, &ts);
|
||||
return (double)ts.tv_sec + (double)ts.tv_nsec * 1.e-9;
|
||||
}
|
||||
|
||||
double getTimeResolution()
|
||||
{
|
||||
struct timespec ts;
|
||||
clock_getres(CLOCK_MONOTONIC, &ts);
|
||||
return (double)ts.tv_sec + (double)ts.tv_nsec * 1.e-9;
|
||||
}
|
||||
|
||||
double getTimeStamp_() { return getTimeStamp(); }
|
||||
14
BasicSolver/3D-mpi-multigrid/src/timing.h
Normal file
14
BasicSolver/3D-mpi-multigrid/src/timing.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 __TIMING_H_
|
||||
#define __TIMING_H_
|
||||
|
||||
extern double getTimeStamp();
|
||||
extern double getTimeResolution();
|
||||
extern double getTimeStamp_();
|
||||
|
||||
#endif // __TIMING_H_
|
||||
22
BasicSolver/3D-mpi-multigrid/src/util.h
Normal file
22
BasicSolver/3D-mpi-multigrid/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_
|
||||
169
BasicSolver/3D-mpi-multigrid/src/vtkWriter.c
Normal file
169
BasicSolver/3D-mpi-multigrid/src/vtkWriter.c
Normal file
@ -0,0 +1,169 @@
|
||||
/*
|
||||
* 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 <stdbool.h>
|
||||
#include <stdio.h>
|
||||
#include <string.h>
|
||||
|
||||
#include "vtkWriter.h"
|
||||
#define G(v, i, j, k) v[(k) * imax * jmax + (j) * imax + (i)]
|
||||
|
||||
static double floatSwap(double f)
|
||||
{
|
||||
union {
|
||||
double f;
|
||||
char b[8];
|
||||
} dat1, dat2;
|
||||
|
||||
dat1.f = f;
|
||||
dat2.b[0] = dat1.b[7];
|
||||
dat2.b[1] = dat1.b[6];
|
||||
dat2.b[2] = dat1.b[5];
|
||||
dat2.b[3] = dat1.b[4];
|
||||
dat2.b[4] = dat1.b[3];
|
||||
dat2.b[5] = dat1.b[2];
|
||||
dat2.b[6] = dat1.b[1];
|
||||
dat2.b[7] = dat1.b[0];
|
||||
return dat2.f;
|
||||
}
|
||||
|
||||
static void writeHeader(VtkOptions* o)
|
||||
{
|
||||
fprintf(o->fh, "# vtk DataFile Version 3.0\n");
|
||||
fprintf(o->fh, "PAMPI cfd solver output\n");
|
||||
if (o->fmt == ASCII) {
|
||||
fprintf(o->fh, "ASCII\n");
|
||||
} else if (o->fmt == BINARY) {
|
||||
fprintf(o->fh, "BINARY\n");
|
||||
}
|
||||
|
||||
fprintf(o->fh, "DATASET STRUCTURED_POINTS\n");
|
||||
fprintf(o->fh, "DIMENSIONS %d %d %d\n", o->grid.imax, o->grid.jmax, o->grid.kmax);
|
||||
fprintf(o->fh,
|
||||
"ORIGIN %f %f %f\n",
|
||||
o->grid.dx * 0.5,
|
||||
o->grid.dy * 0.5,
|
||||
o->grid.dz * 0.5);
|
||||
fprintf(o->fh, "SPACING %f %f %f\n", o->grid.dx, o->grid.dy, o->grid.dz);
|
||||
fprintf(o->fh, "POINT_DATA %d\n", o->grid.imax * o->grid.jmax * o->grid.kmax);
|
||||
}
|
||||
|
||||
void vtkOpen(VtkOptions* o, char* problem)
|
||||
{
|
||||
char filename[50];
|
||||
|
||||
if (o->mode == UNIX) {
|
||||
if (commIsMaster(&o->comm)) {
|
||||
snprintf(filename, 50, "%s-p%d.vtk", problem, o->comm.size);
|
||||
o->fh = fopen(filename, "w");
|
||||
writeHeader(o);
|
||||
}
|
||||
} else if (o->mode == MPI) {
|
||||
}
|
||||
|
||||
if (commIsMaster(&o->comm)) printf("Writing VTK output for %s\n", problem);
|
||||
}
|
||||
|
||||
static void writeScalar(VtkOptions* o, double* s)
|
||||
{
|
||||
int imax = o->grid.imax;
|
||||
int jmax = o->grid.jmax;
|
||||
int kmax = o->grid.kmax;
|
||||
|
||||
for (int k = 0; k < kmax; k++) {
|
||||
for (int j = 0; j < jmax; j++) {
|
||||
for (int i = 0; i < imax; i++) {
|
||||
if (o->fmt == ASCII) {
|
||||
fprintf(o->fh, "%f\n", G(s, i, j, k));
|
||||
} else if (o->fmt == BINARY) {
|
||||
fwrite((double[1]) { floatSwap(G(s, i, j, k)) },
|
||||
sizeof(double),
|
||||
1,
|
||||
o->fh);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
if (o->fmt == BINARY) fprintf(o->fh, "\n");
|
||||
}
|
||||
|
||||
static bool isInitialized(FILE* ptr)
|
||||
{
|
||||
if (ptr == NULL) {
|
||||
printf("vtkWriter not initialize! Call vtkOpen first!\n");
|
||||
return false;
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
void vtkScalar(VtkOptions* o, char* name, double* s)
|
||||
{
|
||||
if (commIsMaster(&o->comm)) printf("Register scalar %s\n", name);
|
||||
|
||||
if (o->mode == UNIX) {
|
||||
if (commIsMaster(&o->comm)) {
|
||||
if (!isInitialized(o->fh)) return;
|
||||
fprintf(o->fh, "SCALARS %s double 1\n", name);
|
||||
fprintf(o->fh, "LOOKUP_TABLE default\n");
|
||||
writeScalar(o, s);
|
||||
}
|
||||
} else if (o->mode == MPI) {
|
||||
}
|
||||
}
|
||||
|
||||
static void writeVector(VtkOptions* o, VtkVector vec)
|
||||
{
|
||||
int imax = o->grid.imax;
|
||||
int jmax = o->grid.jmax;
|
||||
int kmax = o->grid.kmax;
|
||||
|
||||
for (int k = 0; k < kmax; k++) {
|
||||
for (int j = 0; j < jmax; j++) {
|
||||
for (int i = 0; i < imax; i++) {
|
||||
if (o->fmt == ASCII) {
|
||||
fprintf(o->fh,
|
||||
"%f %f %f\n",
|
||||
G(vec.u, i, j, k),
|
||||
G(vec.v, i, j, k),
|
||||
G(vec.w, i, j, k));
|
||||
} else if (o->fmt == BINARY) {
|
||||
fwrite((double[3]) { floatSwap(G(vec.u, i, j, k)),
|
||||
floatSwap(G(vec.v, i, j, k)),
|
||||
floatSwap(G(vec.w, i, j, k)) },
|
||||
sizeof(double),
|
||||
3,
|
||||
o->fh);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
if (o->fmt == BINARY) fprintf(o->fh, "\n");
|
||||
}
|
||||
|
||||
void vtkVector(VtkOptions* o, char* name, VtkVector vec)
|
||||
{
|
||||
if (commIsMaster(&o->comm)) printf("Register vector %s\n", name);
|
||||
|
||||
if (o->mode == UNIX) {
|
||||
if (commIsMaster(&o->comm)) {
|
||||
if (!isInitialized(o->fh)) return;
|
||||
fprintf(o->fh, "VECTORS %s double\n", name);
|
||||
writeVector(o, vec);
|
||||
}
|
||||
} else if (o->mode == MPI) {
|
||||
}
|
||||
}
|
||||
|
||||
void vtkClose(VtkOptions* o)
|
||||
{
|
||||
if (o->mode == UNIX) {
|
||||
if (commIsMaster(&o->comm)) {
|
||||
fclose(o->fh);
|
||||
o->fh = NULL;
|
||||
}
|
||||
} else if (o->mode == MPI) {
|
||||
}
|
||||
}
|
||||
33
BasicSolver/3D-mpi-multigrid/src/vtkWriter.h
Normal file
33
BasicSolver/3D-mpi-multigrid/src/vtkWriter.h
Normal file
@ -0,0 +1,33 @@
|
||||
/*
|
||||
* 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 <stdio.h>
|
||||
|
||||
#include "comm.h"
|
||||
#include "grid.h"
|
||||
|
||||
typedef enum VtkFormat { ASCII = 0, BINARY } VtkFormat;
|
||||
typedef enum VtkMode { UNIX = 0, MPI } VtkMode;
|
||||
|
||||
typedef struct VtkOptions {
|
||||
VtkFormat fmt;
|
||||
VtkMode mode;
|
||||
Grid grid;
|
||||
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_
|
||||
2560011
BasicSolver/3D-mpi/karman-p4.vtk
2560011
BasicSolver/3D-mpi/karman-p4.vtk
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
@ -1416,9 +1416,9 @@ void multiGrid(Solver* solver)
|
||||
|
||||
void restrictMG(Solver* solver)
|
||||
{
|
||||
int imax = solver->grid.imax;
|
||||
int jmax = solver->grid.jmax;
|
||||
int kmax = solver->grid.kmax;
|
||||
int imax = solver->grid.imax / 2;
|
||||
int jmax = solver->grid.jmax / 2;
|
||||
int kmax = solver->grid.kmax / 2;
|
||||
|
||||
double* r = solver->r[solver->currentlevel + 1];
|
||||
double* oldr = solver->r[solver->currentlevel];
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user