2D-MPI Multigrid complete and tested

This commit is contained in:
Aditya Ujeniya 2024-03-05 21:45:17 +01:00
parent 9c8ee975c5
commit b4bef94066
51 changed files with 1300314 additions and 40197 deletions

View File

@ -0,0 +1,76 @@
#=======================================================================================
# 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
distclean: clean
$(info ===> DIST CLEAN)
@rm -f $(TARGET)
info:
$(info $(CFLAGS))
$(Q)$(CC) $(VERSION)
asm: $(BUILD_DIR) $(ASM)
vis:
$(info ===> REMOVING VIZUALISATION FILES)
@rm -f vtk_files/particle*.vtk
@rm -f vis_files/particle*.dat
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)

View File

@ -0,0 +1,48 @@
# 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_AFFINITY
#OPTIONS += -DVERBOSE_DATASIZE
#OPTIONS += -DVERBOSE_TIMER
```
The verbosity options enable detailed output about affinity settings, allocation sizes and timer resolution.
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.

View File

@ -0,0 +1,73 @@
#==============================================================================
# Laminar Canal Flow
#==============================================================================
# Problem specific Data:
# ---------------------
name backstep # name of flow setup
bcTop 1 # flags for boundary conditions
bcBottom 1 # 1 = no-slip 3 = outflow
bcLeft 3 # 2 = free-slip 4 = periodic
bcRight 3 #
gx 0.0 # Body forces (e.g. gravity)
gy 0.0 #
re 36500.0 # Reynolds number
u_init 1.0 # initial value for velocity in x-direction
v_init 0.0 # initial value for velocity in y-direction
p_init 1.0 # initial value for pressure
# Geometry Data:
# -------------
xlength 7.0 # domain size in x-direction
ylength 1.5 # domain size in y-direction
imax 210 # number of interior cells in x-direction
jmax 45 # number of interior cells in y-direction
# Time Data:
# ---------
te 60.0 # final time
dt 0.02 # time stepsize
tau 0.5 # safety factor for time stepsize control (<0 constant delt)
# Pressure Iteration Data:
# -----------------------
itermax 500 # maximal number of pressure iteration in one time step
eps 0.0001 # stopping tolerance for pressure iteration
rho 0.52
omg 1.8 # relaxation parameter for SOR iteration
gamma 0.9 # upwind differencing factor gamma
levels 5 # Multigrid levels
# Particle Tracing Data:
# -----------------------
numberOfParticles 200
startTime 0
injectTimePeriod 1.0
writeTimePeriod 0.5
x1 0.0
y1 0.5
x2 0.0
y2 1.5
# Obstacle Geometry Data:
# -----------------------
# Shape 0 disable, 1 Rectangle/Square, 2 Circle
shape 1
xCenter 0.0
yCenter 0.0
xRectLength 2.0
yRectLength 1.0
circleRadius 1.0
#===============================================================================

View File

@ -0,0 +1,73 @@
#==============================================================================
# Laminar Canal Flow
#==============================================================================
# Problem specific Data:
# ---------------------
name canal # name of flow setup
bcTop 1 # flags for boundary conditions
bcBottom 1 # 1 = no-slip 3 = outflow
bcLeft 3 # 2 = free-slip 4 = periodic
bcRight 3 #
gx 0.0 # Body forces (e.g. gravity)
gy 0.0 #
re 100.0 # Reynolds number
u_init 1.0 # initial value for velocity in x-direction
v_init 0.0 # initial value for velocity in y-direction
p_init 1.0 # initial value for pressure
# Geometry Data:
# -------------
xlength 30.0 # domain size in x-direction
ylength 4.0 # domain size in y-direction
imax 1200 # number of interior cells in x-direction
jmax 160 # number of interior cells in y-direction
# Time Data:
# ---------
te 60.0 # final time
dt 0.02 # time stepsize
tau 0.5 # safety factor for time stepsize control (<0 constant delt)
# Pressure Iteration Data:
# -----------------------
itermax 500 # maximal number of pressure iteration in one time step
eps 0.0001 # stopping tolerance for pressure iteration
rho 0.52
omg 1.8 # relaxation parameter for SOR iteration
gamma 0.9 # upwind differencing factor gamma
levels 5 # Multigrid levels
# Particle Tracing Data:
# -----------------------
numberOfParticles 60
startTime 500
injectTimePeriod 2.0
writeTimePeriod 0.5
x1 1.0
y1 0.0
x2 1.0
y2 4.0
# Obstacle Geometry Data:
# -----------------------
# Shape 0 disable, 1 Rectangle/Square, 2 Circle
shape 0
xCenter 10.0
yCenter 2
xRectLength 6.0
yRectLength 1.0
circleRadius 1.0
#===============================================================================

View File

@ -0,0 +1,10 @@
# Supported: GCC, CLANG, ICC
TAG ?= ICC
ENABLE_OPENMP ?= false
#Feature options
OPTIONS += -DARRAY_ALIGNMENT=64
OPTIONS += -DVERBOSE
#OPTIONS += -DVERBOSE_AFFINITY
#OPTIONS += -DVERBOSE_DATASIZE
#OPTIONS += -DVERBOSE_TIMER

View File

@ -0,0 +1,61 @@
#==============================================================================
# Driven Cavity
#==============================================================================
# Problem specific Data:
# ---------------------
name dcavity # name of flow setup
bcTop 1 # flags for boundary conditions
bcBottom 1 # 1 = no-slip 3 = outflow
bcLeft 1 # 2 = free-slip 4 = periodic
bcRight 1 #
gx 0.0 # Body forces (e.g. gravity)
gy 0.0 #
re 10.0 # Reynolds number
u_init 0.0 # initial value for velocity in x-direction
v_init 0.0 # initial value for velocity in y-direction
p_init 0.0 # initial value for pressure
# Geometry Data:
# -------------
xlength 1.0 # domain size in x-direction
ylength 1.0 # domain size in y-direction
imax 100 # number of interior cells in x-direction
jmax 100 # number of interior cells in y-direction
# Time Data:
# ---------
te 10.0 # final time
dt 0.02 # time stepsize
tau 0.5 # safety factor for time stepsize control (<0 constant delt)
# Pressure Iteration Data:
# -----------------------
itermax 1000 # maximal number of pressure iteration in one time step
eps 0.001 # stopping tolerance for pressure iteration
rho 0.5 #
omg 1.7 # relaxation parameter for SOR iteration
gamma 0.9 # upwind differencing factor gamma
levels 5 # Multigrid levels
# Particle Tracing Data:
# -----------------------
numberOfParticles 30
startTime 2
injectTimePeriod 2.0
writeTimePeriod 0.1
x1 0.1
y1 0.5
x2 0.9
y2 0.5
#===============================================================================

View File

@ -0,0 +1,16 @@
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 = -Ofast -std=c99 $(OPENMP)
#CFLAGS = -Ofast -fnt-store=aggressive -std=c99 $(OPENMP) #AMD CLANG
LFLAGS = $(OPENMP)
DEFINES = -D_GNU_SOURCE# -DDEBUG
INCLUDES = -I/usr/local/include

View 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 =

View File

@ -0,0 +1,14 @@
CC = mpiicc
GCC = gcc
LINKER = $(CC)
ifeq ($(ENABLE_OPENMP),true)
OPENMP = -qopenmp
endif
VERSION = --version
CFLAGS = -O3 -xHost -qopt-zmm-usage=high -std=c99 $(OPENMP)
LFLAGS = $(OPENMP)
DEFINES = -D_GNU_SOURCE
INCLUDES =
LIBS =

View File

@ -0,0 +1,73 @@
#==============================================================================
# Laminar Canal Flow
#==============================================================================
# Problem specific Data:
# ---------------------
name karman # name of flow setup
bcTop 1 # flags for boundary conditions
bcBottom 1 # 1 = no-slip 3 = outflow
bcLeft 3 # 2 = free-slip 4 = periodic
bcRight 3 #
gx 0.0 # Body forces (e.g. gravity)
gy 0.0 #
re 5050.0 # Reynolds number
u_init 1.0 # initial value for velocity in x-direction
v_init 0.0 # initial value for velocity in y-direction
p_init 0.0 # initial value for pressure
# Geometry Data:
# -------------
xlength 30.0 # domain size in x-direction
ylength 8.0 # domain size in y-direction
imax 1400 # number of interior cells in x-direction
jmax 320 # number of interior cells in y-direction
# Time Data:
# ---------
te 150.0 # final time
dt 0.02 # time stepsize
tau 0.5 # safety factor for time stepsize control (<0 constant delt)
# Pressure Iteration Data:
# -----------------------
itermax 200 # maximal number of pressure iteration in one time step
eps 0.001 # stopping tolerance for pressure iteration
rho 0.52
omg 1.75 # relaxation parameter for SOR iteration
gamma 0.9 # upwind differencing factor gamma
levels 5 # Multigrid levels
# Particle Tracing Data:
# -----------------------
numberOfParticles 200
startTime 150
injectTimePeriod 1.0
writeTimePeriod 0.5
x1 0.0
y1 3.8
x2 0.0
y2 4.1
# Obstacle Geometry Data:
# -----------------------
# Shape 0 disable, 1 Rectangle/Square, 2 Circle
shape 2
xCenter 5.0
yCenter 4.0
xRectLength 2.0
yRectLength 1.0
circleRadius 1.0
#===============================================================================

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,17 @@
#!/bin/bash
rm data.csv
echo "Rank, SOR, RB" >> data.csv
for i in {1..72}
do
res=$(mpirun -n $i ./exe-ICC dcavity.par 1 | grep "Solution took" | cut -c 14- )
res2=$(mpirun -n $i ./exe-ICC dcavity.par 2 | grep "Solution took" | cut -c 14- )
s="$i, $res, $res2"
s=${s//"s"/}
echo $s >> data.csv
done

View File

@ -0,0 +1,61 @@
/*
* 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.
*/
#ifdef __linux__
#ifdef _OPENMP
#include <pthread.h>
#include <sched.h>
#include <stdio.h>
#include <stdlib.h>
#include <sys/syscall.h>
#include <sys/types.h>
#include <unistd.h>
#define MAX_NUM_THREADS 128
#define gettid() syscall(SYS_gettid)
static int getProcessorID(cpu_set_t* cpu_set)
{
int processorId;
for (processorId = 0; processorId < MAX_NUM_THREADS; processorId++) {
if (CPU_ISSET(processorId, cpu_set)) {
break;
}
}
return processorId;
}
int affinity_getProcessorId()
{
cpu_set_t cpu_set;
CPU_ZERO(&cpu_set);
sched_getaffinity(gettid(), sizeof(cpu_set_t), &cpu_set);
return getProcessorID(&cpu_set);
}
void affinity_pinThread(int processorId)
{
cpu_set_t cpuset;
pthread_t thread;
thread = pthread_self();
CPU_ZERO(&cpuset);
CPU_SET(processorId, &cpuset);
pthread_setaffinity_np(thread, sizeof(cpu_set_t), &cpuset);
}
void affinity_pinProcess(int processorId)
{
cpu_set_t cpuset;
CPU_ZERO(&cpuset);
CPU_SET(processorId, &cpuset);
sched_setaffinity(0, sizeof(cpu_set_t), &cpuset);
}
#endif /*_OPENMP*/
#endif /*__linux__*/

View 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 AFFINITY_H
#define AFFINITY_H
extern int affinity_getProcessorId();
extern void affinity_pinProcess(int);
extern void affinity_pinThread(int);
#endif /*AFFINITY_H*/

View 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;
}

View 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

View 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*/

View File

@ -0,0 +1,92 @@
/*
* 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 <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include "parameter.h"
#include "particletracing.h"
#include "progress.h"
#include "timing.h"
#include <mpi.h>
int main(int argc, char** argv)
{
int rank;
double S, E;
Parameter params;
Solver solver;
ParticleTracer particletracer;
MPI_Init(&argc, &argv);
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
initParameter(&params);
if (argc < 2) {
printf("Usage: %s <configFile>\n", argv[0]);
exit(EXIT_SUCCESS);
}
readParameter(&params, argv[1]);
initSolver(&solver, &params);
initParticleTracer(&particletracer, &params, &solver);
if (rank == 0) {
printParameter(&params);
}
printParticleTracerParameters(&particletracer);
initProgress(solver.te);
double tau = solver.tau;
double te = solver.te;
double t = 0.0;
S = getTimeStamp();
while (t <= te) {
if (tau > 0.0) {
computeTimestep(&solver);
}
setBoundaryConditions(&solver);
setSpecialBoundaryCondition(&solver);
setObjectBoundaryCondition(&solver);
computeFG(&solver);
computeRHS(&solver);
multiGrid(&solver);
adaptUV(&solver);
trace(&particletracer, solver.u, solver.v, solver.s, t);
t += solver.dt;
#ifdef VERBOSE
if (rank == 0) {
printf("TIME %f , TIMESTEP %f\n", t, solver.dt);
}
#else
printProgress(t);
#endif
}
E = getTimeStamp();
stopProgress();
if (rank == 0) {
printf("Solution took %.2fs\n", E - S);
}
collectResult(&solver);
freeParticles(&particletracer);
MPI_Finalize();
return EXIT_SUCCESS;
}

View File

@ -0,0 +1,135 @@
/*
* 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->imax = 100;
param->jmax = 100;
param->itermax = 1000;
param->eps = 0.0001;
param->omg = 1.7;
param->re = 100.0;
param->gamma = 0.9;
param->tau = 0.5;
param->rho = 0.99;
param->levels = 5;
}
void readParameter(Parameter* param, const char* filename)
{
FILE* fp = fopen(filename, "r");
char line[MAXLINE];
int i;
if (!fp) {
fprintf(stderr, "Could not open parameter file: %s\n", filename);
exit(EXIT_FAILURE);
}
while (!feof(fp)) {
line[0] = '\0';
fgets(line, MAXLINE, fp);
for (i = 0; line[i] != '\0' && line[i] != '#'; i++)
;
line[i] = '\0';
char* tok = strtok(line, " ");
char* val = strtok(NULL, " ");
#define PARSE_PARAM(p, f) \
if (strncmp(tok, #p, sizeof(#p) / sizeof(#p[0]) - 1) == 0) { \
param->p = f(val); \
}
#define PARSE_STRING(p) PARSE_PARAM(p, strdup)
#define PARSE_INT(p) PARSE_PARAM(p, atoi)
#define PARSE_REAL(p) PARSE_PARAM(p, atof)
if (tok != NULL && val != NULL) {
PARSE_REAL(xlength);
PARSE_REAL(ylength);
PARSE_INT(imax);
PARSE_INT(jmax);
PARSE_INT(itermax);
PARSE_INT(levels);
PARSE_REAL(eps);
PARSE_REAL(omg);
PARSE_REAL(re);
PARSE_REAL(tau);
PARSE_REAL(gamma);
PARSE_REAL(dt);
PARSE_REAL(te);
PARSE_REAL(gx);
PARSE_REAL(gy);
PARSE_STRING(name);
PARSE_INT(bcLeft);
PARSE_INT(bcRight);
PARSE_INT(bcBottom);
PARSE_INT(bcTop);
PARSE_REAL(u_init);
PARSE_REAL(v_init);
PARSE_REAL(p_init);
PARSE_REAL(rho);
/* Added new particle tracing parameters */
PARSE_INT(numberOfParticles);
PARSE_REAL(startTime);
PARSE_REAL(injectTimePeriod);
PARSE_REAL(writeTimePeriod);
PARSE_REAL(x1);
PARSE_REAL(y1);
PARSE_REAL(x2);
PARSE_REAL(y2);
/* Added obstacle geometry parameters */
PARSE_INT(shape);
PARSE_REAL(xCenter);
PARSE_REAL(yCenter);
PARSE_REAL(xRectLength);
PARSE_REAL(yRectLength);
PARSE_REAL(circleRadius);
}
}
fclose(fp);
}
void printParameter(Parameter* param)
{
printf("Parameters for %s\n", param->name);
printf("Boundary conditions Left:%d Right:%d Bottom:%d Top:%d\n",
param->bcLeft,
param->bcRight,
param->bcBottom,
param->bcTop);
printf("\tReynolds number: %.2f\n", param->re);
printf("\tInit arrays: U:%.2f V:%.2f P:%.2f\n",
param->u_init,
param->v_init,
param->p_init);
printf("Geometry data:\n");
printf("\tDomain box size (x, y): %.2f, %.2f\n", param->xlength, param->ylength);
printf("\tCells (x, y): %d, %d\n", param->imax, param->jmax);
printf("Timestep parameters:\n");
printf("\tDefault stepsize: %.2f, Final time %.2f\n", param->dt, param->te);
printf("\tTau factor: %.2f\n", param->tau);
printf("Iterative solver parameters:\n");
printf("\tMax iterations: %d\n", param->itermax);
printf("\tepsilon (stopping tolerance) : %f\n", param->eps);
printf("\tgamma (stopping tolerance) : %f\n", param->gamma);
printf("\tomega (SOR relaxation): %f\n", param->omg);
}

View File

@ -0,0 +1,34 @@
/*
* 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 {
double xlength, ylength;
int imax, jmax;
int itermax, levels;
double eps, omg, rho;
double re, tau, gamma;
double te, dt;
double gx, gy;
char* name;
int bcLeft, bcRight, bcBottom, bcTop;
double u_init, v_init, p_init;
int numberOfParticles;
double startTime, injectTimePeriod, writeTimePeriod;
double x1, y1, x2, y2;
int shape;
double xCenter, yCenter, xRectLength, yRectLength, circleRadius;
} Parameter;
void initParameter(Parameter*);
void readParameter(Parameter*, const char*);
void printParameter(Parameter*);
#endif

View File

@ -0,0 +1,629 @@
/*
* 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"
#define U(i, j) u[(j) * (imaxLocal + 2) + (i)]
#define V(i, j) v[(j) * (imaxLocal + 2) + (i)]
#define S(i, j) s[(j) * (imaxLocal + 2) + (i)]
static int ts = 0;
#define IDIM 0
#define JDIM 1
#define XOFFSET 0
#define YOFFSET 1
#define XOFFSETEND 2
#define YOFFSETEND 3
static double sum(int* sizes, int size)
{
double sum = 0;
for (int i = 0; i < size; ++i) {
sum += sizes[i];
}
return sum;
}
void printUV(ParticleTracer* particletracer, double* u, double* v)
{
int imaxLocal = particletracer->imaxLocal;
for (int i = 0; i < particletracer->size; i++) {
if (i == particletracer->rank) {
printf(
"\n### RANK %d #######################################################\n",
particletracer->rank);
printf("\nGrid U : \n");
for (int j = 0; j < particletracer->jmaxLocal + 2; j++) {
printf("%02d: ", j);
for (int i = 0; i < particletracer->imaxLocal + 2; i++) {
printf("%4.2f ", u[j * (imaxLocal + 2) + i]);
}
printf("\n");
}
fflush(stdout);
printf("\nGrid V : \n");
for (int j = 0; j < particletracer->jmaxLocal + 2; j++) {
printf("%02d: ", j);
for (int i = 0; i < particletracer->imaxLocal + 2; i++) {
printf("%4.2f ", v[j * (imaxLocal + 2) + i]);
}
printf("\n");
}
fflush(stdout);
}
MPI_Barrier(MPI_COMM_WORLD);
}
}
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)
{
double x, y;
for (int i = 0; i < particletracer->numberOfParticles; ++i) {
x = particletracer->linSpaceLine[i].x;
y = particletracer->linSpaceLine[i].y;
if (x >= particletracer->xOffset && y >= particletracer->yOffset &&
x <= particletracer->xOffsetEnd && y <= particletracer->yOffsetEnd) {
// 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].flag = true;
++(particletracer->pointer);
++(particletracer->totalParticles);
}
}
}
void advanceParticles(ParticleTracer* particletracer,
double* restrict u,
double* restrict v,
int* restrict s,
double time)
{
int imax = particletracer->imax;
int jmax = particletracer->jmax;
int imaxLocal = particletracer->imaxLocal;
int jmaxLocal = particletracer->jmaxLocal;
double dx = particletracer->dx;
double dy = particletracer->dy;
double xlength = particletracer->xlength;
double ylength = particletracer->ylength;
Particle buff[particletracer->size][100];
for (int i = 0; i < particletracer->size; ++i) {
for (int j = 0; j < 100; ++j) {
buff[i][j].x = 0.0;
buff[i][j].y = 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 x = xTemp - particletracer->xOffset;
double y = yTemp - particletracer->yOffset;
int iCoord = (int)(x / dx) + 1;
int jCoord = (int)((y + 0.5 * dy) / dy) + 1;
double x1 = (double)(iCoord - 1) * dx;
double y1 = ((double)(jCoord - 1) - 0.5) * dy;
double x2 = (double)iCoord * dx;
double y2 = ((double)jCoord - 0.5) * dy;
double u_n = (1.0 / (dx * dy)) *
((x2 - x) * (y2 - y) * U(iCoord - 1, jCoord - 1) +
(x - x1) * (y2 - y) * U(iCoord, jCoord - 1) +
(x2 - x) * (y - y1) * U(iCoord - 1, jCoord) +
(x - x1) * (y - y1) * U(iCoord, jCoord));
double new_x = (x + 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;
x1 = ((double)(iCoord - 1) - 0.5) * dx;
y1 = (double)(jCoord - 1) * dy;
x2 = ((double)iCoord - 0.5) * dx;
y2 = (double)jCoord * dy;
double v_n = (1.0 / (dx * dy)) *
((x2 - x) * (y2 - y) * V(iCoord - 1, jCoord - 1) +
(x - x1) * (y2 - y) * V(iCoord, jCoord - 1) +
(x2 - x) * (y - y1) * V(iCoord - 1, jCoord) +
(x - x1) * (y - y1) * V(iCoord, jCoord));
double new_y = (y + particletracer->yOffset) + particletracer->dt * v_n;
particletracer->particlePool[i].y = new_y;
// 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 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]) &&
i != particletracer->rank) {
buff[i][particleBufIndex[i]].x = new_x;
buff[i][particleBufIndex[i]].y = new_y;
buff[i][particleBufIndex[i]].flag = true;
++particleBufIndex[i];
}
}
particletracer->particlePool[i].flag = false;
}
int i_new = new_x / dx, j_new = new_y / dy;
int iOffset = particletracer->xOffset / dx,
jOffset = particletracer->yOffset / dy;
if (S(i_new - iOffset, j_new - jOffset) != NONE) {
particletracer->particlePool[i].flag = false;
}
}
}
for (int i = 0; i < particletracer->size; ++i) {
if (i != particletracer->rank) {
MPI_Send(buff[i],
100,
particletracer->mpi_particle,
i,
0,
particletracer->comm);
}
}
for (int i = 0; i < particletracer->size; ++i) {
if (i != particletracer->rank) {
MPI_Recv(buff[i],
100,
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 < 100; ++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].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].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, "vis_files/particles_%d.dat", 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;
fprintf(fp, "%f %f\n", x, y);
// 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;
fprintf(fp, "%f %f\n", x, y);
}
// 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)
{
int dims[NDIMS] = { 0, 0 };
int periods[NDIMS] = { 0, 0 };
/* 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->xlength = params->xlength;
particletracer->ylength = params->ylength;
particletracer->x1 = params->x1;
particletracer->y1 = params->y1;
particletracer->x2 = params->x2;
particletracer->y2 = params->y2;
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->imaxLocal = solver->imaxLocal;
particletracer->jmaxLocal = solver->jmaxLocal;
particletracer->estimatedNumParticles = (particletracer->imaxLocal *
particletracer->jmaxLocal);
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].flag = false;
}
particletracer->linSpaceLine = malloc(
sizeof(Particle) * particletracer->numberOfParticles);
/* duplicating communication from solver */
MPI_Comm_dup(solver->comm, &particletracer->comm);
particletracer->rank = solver->rank;
particletracer->size = solver->size;
particletracer->offset = (double*)malloc(sizeof(double) * 4 * particletracer->size);
memcpy(particletracer->dims, solver->dims, sizeof(solver->dims));
memcpy(particletracer->coords, solver->coords, sizeof(solver->coords));
memcpy(particletracer->iNeighbours, solver->iNeighbours, sizeof(solver->iNeighbours));
memcpy(particletracer->jNeighbours, solver->jNeighbours, sizeof(solver->jNeighbours));
particletracer->xLocal = particletracer->imaxLocal * particletracer->dx;
particletracer->yLocal = particletracer->jmaxLocal * particletracer->dy;
double offset[4][particletracer->size];
particletracer->xOffset = solver->xOffset;
particletracer->yOffset = solver->yOffset;
particletracer->xOffsetEnd = particletracer->xOffset + particletracer->xLocal;
particletracer->yOffsetEnd = particletracer->yOffset + particletracer->yLocal;
printf("Rank : %d, xOffset : %.2f, yOffset : %.2f, xOffsetEnd : %.2f, yOffsetEnd : "
"%.2f\n",
particletracer->rank,
particletracer->xOffset,
particletracer->yOffset,
particletracer->xOffsetEnd,
particletracer->yOffsetEnd);
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->xOffsetEnd,
1,
MPI_DOUBLE,
offset[2],
1,
MPI_DOUBLE,
particletracer->comm);
MPI_Allgather(&particletracer->yOffsetEnd,
1,
MPI_DOUBLE,
offset[3],
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 and its xOffset : %.2f, yOffset : %.2f, xOffsetEnd :
// %.2f, yOffsetEnd : %.2f\n", i, particletracer->offset[i +
// particletracer->size * XOFFSET], particletracer->offset[i +
// particletracer->size * YOFFSET], particletracer->offset[i +
// particletracer->size * XOFFSETEND], particletracer->offset[i +
// particletracer->size * YOFFSETEND]);
// }
// }
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].flag = true;
// if(particletracer->rank == 1) printf("\nRank : %d, with linspace X : %.2f and
// linspace Y : %.2f\n", particletracer->rank, particletracer->linSpaceLine[i].x ,
// particletracer->linSpaceLine[i].y);
}
// Create the mpi_particle datatype
MPI_Datatype mpi_particle;
int lengths[3] = { 1, 1, 1 };
MPI_Aint displacements[3];
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.flag, &displacements[2]);
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);
MPI_Datatype types[3] = { MPI_DOUBLE, MPI_DOUBLE, MPI_C_BOOL };
MPI_Type_create_struct(3,
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, x2 : %.2f, y2 : %.2f\n",
particletracer->x1,
particletracer->y1,
particletracer->x2,
particletracer->y2);
printf("\tPointer : %d, TotalParticles : %d\n",
particletracer->pointer,
particletracer->totalParticles);
printf("\tdt : %.2f, dx : %.2f, dy : %.2f\n",
particletracer->dt,
particletracer->dx,
particletracer->dy);
printf("\tcoord[0] : %d, coord[1] : %d\n",
particletracer->coords[IDIM],
particletracer->coords[JDIM]);
printf("\txOffset : %.2f, yOffset : %.2f\n",
particletracer->xOffset,
particletracer->yOffset);
printf("\txOffsetEnd : %.2f, yOffsetEnd : %.2f\n",
particletracer->xOffsetEnd,
particletracer->yOffsetEnd);
printf("\txLocal : %.2f, yLocal : %.2f\n",
particletracer->xLocal,
particletracer->yLocal);
}
void trace(ParticleTracer* particletracer,
double* restrict u,
double* restrict v,
int* restrict s,
double time)
{
if (time >= particletracer->startTime) {
if ((time - particletracer->lastInjectTime) >= particletracer->injectTimePeriod) {
injectParticles(particletracer);
particletracer->lastInjectTime = time;
}
if ((time - particletracer->lastWriteTime) >= particletracer->writeTimePeriod) {
writeParticles(particletracer);
particletracer->lastWriteTime = time;
}
advanceParticles(particletracer, u, v, s, 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].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].flag = memPool[i].flag;
++totalParticles;
}
}
particletracer->totalParticles = totalParticles;
particletracer->pointer = totalParticles;
memcpy(particletracer->particlePool, tempPool, totalParticles * sizeof(Particle));
}

View File

@ -0,0 +1,65 @@
/*
* 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 2
typedef enum COORD { X = 0, Y, NCOORD } COORD;
typedef struct {
double x, y;
bool flag;
} Particle;
typedef struct {
int numberOfParticles, totalParticles;
double startTime, injectTimePeriod, writeTimePeriod, lastInjectTime, lastUpdateTime,
lastWriteTime;
int estimatedNumParticles;
double dx, dy, dt;
Particle* linSpaceLine;
Particle* particlePool;
int pointer;
double imax, jmax, xlength, ylength, imaxLocal, jmaxLocal;
double x1, y1, x2, y2;
MPI_Comm comm;
MPI_Datatype mpi_particle;
int rank, size;
int iNeighbours[NDIMS], jNeighbours[NDIMS];
int coords[NDIMS], dims[NDIMS];
double xLocal, yLocal, xOffset, yOffset, xOffsetEnd, yOffsetEnd;
double* offset;
} ParticleTracer;
extern void initParticleTracer(ParticleTracer*, Parameter*, Solver*);
extern void injectParticles(ParticleTracer*);
extern void advanceParticles(ParticleTracer*, double*, double*, int*, double);
extern void freeParticles(ParticleTracer*);
extern void writeParticles(ParticleTracer*);
extern void printParticleTracerParameters(ParticleTracer*);
extern void printParticles(ParticleTracer*);
extern void trace(ParticleTracer*, double*, double*, int*, double);
extern void compress(ParticleTracer*);
#endif

View File

@ -0,0 +1,60 @@
/*
* 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 <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "progress.h"
static double _end;
static int _current;
static int _rank = -1;
void initProgress(double end)
{
MPI_Comm_rank(MPI_COMM_WORLD, &_rank);
_end = end;
_current = 0;
if (_rank == 0) {
printf("[ ]");
fflush(stdout);
}
}
void printProgress(double current)
{
if (_rank == 0) {
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()
{
if (_rank == 0) {
printf("\n");
fflush(stdout);
}
}

View 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

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,91 @@
/*
* 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 "parameter.h"
#include <mpi.h>
#define NDIMS 2
enum OBJECTBOUNDARY {
NONE = 0,
TOP,
BOTTOM,
LEFT,
RIGHT,
TOPLEFT,
BOTTOMLEFT,
TOPRIGHT,
BOTTOMRIGHT,
LOCAL
};
enum BC { NOSLIP = 1, SLIP, OUTFLOW, PERIODIC };
/// @brief
enum SHAPE { NOSHAPE = 0, RECT, CIRCLE };
typedef struct {
/* geometry and grid information */
double dx, dy;
int imax, jmax;
double xlength, ylength;
/* arrays */
double *p, *rhs, **r, **e;
;
double *f, *g;
double *u, *v;
int* s;
/* parameters */
double eps, omega, rho;
double re, tau, gamma;
double gx, gy;
/* time stepping */
int itermax, levels, currentlevel;
double dt, te;
double dtBound;
char* problem;
int bcLeft, bcRight, bcBottom, bcTop;
/* mpi */
int rank;
int size;
MPI_Comm comm;
MPI_Datatype bufferTypes[NDIMS * 2];
MPI_Aint sdispls[NDIMS * 2], rdispls[NDIMS * 2];
MPI_Datatype bufferIntTypes[NDIMS * 2];
MPI_Aint sdisplsInt[NDIMS * 2], rdisplsInt[NDIMS * 2];
int iNeighbours[NDIMS], jNeighbours[NDIMS];
int coords[NDIMS], dims[NDIMS];
int imaxLocal, jmaxLocal;
double xLocal, yLocal, xOffset, yOffset, xOffsetEnd, yOffsetEnd;
} 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 computeTimestep(Solver*);
extern void setBoundaryConditions(Solver*);
extern void setSpecialBoundaryCondition(Solver*);
extern void setObjectBoundaryCondition(Solver*);
extern void setPressureBoundaryCondition(Solver*);
extern void computeFG(Solver*);
extern void adaptUV(Solver*);
extern void collectResult(Solver*);
extern void writeResult(Solver*, double*, double*, double*);
extern void debugExchange(Solver*, double*);
extern void print(Solver*, double*);
extern void printInt(Solver*, int*);
#endif

View 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(); }

View 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_

View 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_

View File

@ -0,0 +1,7 @@
set terminal png size 1024,768 enhanced font ,12
set output 'p.png'
set datafile separator whitespace
set grid
set hidden3d
splot 'pressure.dat' using 1:2:3 with lines

View File

@ -0,0 +1,6 @@
set terminal png size 3600,768 enhanced font ,28
set output 'velocity.png'
set size ratio -1
set datafile separator whitespace
plot 'velocity.dat' using 1:2:3:4:5 with vectors filled head size 0.01,20,60 lc palette

File diff suppressed because it is too large Load Diff

Binary file not shown.

After

Width:  |  Height:  |  Size: 338 KiB

View File

@ -0,0 +1,10 @@
unset border; unset tics; unset key;
set term gif animate delay 50
set output "trace.gif"
set xrange [0:30]
set yrange [0:4]
do for [ts=0:120] {
plot "particles_".ts.".dat" with points pointtype 7
}
unset output

View File

@ -0,0 +1,12 @@
unset border; unset tics; unset key;
set term gif animate delay 10
set output "trace.gif"
set xrange [0:7]
set yrange [0:1.5]
set object 1 rect from 0.0,0.0 to 1.0,0.5 lw 2
do for [ts=0:300] {
plot "particles_".ts.".dat" with points pointtype 7
}
unset output

View File

@ -0,0 +1,13 @@
unset border; unset tics; unset key;
set term gif animate delay 10
set output "trace.gif"
set xrange [0:30]
set yrange [0:8]
set size ratio -1
set object 1 circle front at 5.0,4.0 size 1.0 fillcolor rgb "black" lw 2
do for [ts=0:500] {
plot "particles_".ts.".dat" with points pointtype 7 pointsize 0.3
}
unset output

Binary file not shown.

After

Width:  |  Height:  |  Size: 751 KiB

View File

@ -44,6 +44,7 @@ eps 0.0001 # stopping tolerance for pressure iteration
rho 0.52 rho 0.52
omg 1.8 # relaxation parameter for SOR iteration omg 1.8 # relaxation parameter for SOR iteration
gamma 0.9 # upwind differencing factor gamma gamma 0.9 # upwind differencing factor gamma
levels 5 # Multigrid levels
# Particle Tracing Data: # Particle Tracing Data:
# ----------------------- # -----------------------

View File

@ -26,8 +26,8 @@ p_init 1.0 # initial value for pressure
xlength 30.0 # domain size in x-direction xlength 30.0 # domain size in x-direction
ylength 4.0 # domain size in y-direction ylength 4.0 # domain size in y-direction
imax 200 # number of interior cells in x-direction imax 1200 # number of interior cells in x-direction
jmax 50 # number of interior cells in y-direction jmax 160 # number of interior cells in y-direction
# Time Data: # Time Data:
# --------- # ---------

File diff suppressed because it is too large Load Diff

View File

@ -680,9 +680,10 @@ double smoothRB(Solver* solver)
void multiGrid(Solver* solver) void multiGrid(Solver* solver)
{ {
double res = 0.0; double res = 0.0;
int imax = solver->imax; int imax = solver->imax;
int jmax = solver->jmax; int jmax = solver->jmax;
if (solver->currentlevel == (solver->levels - 1)) { if (solver->currentlevel == (solver->levels - 1) ||
(solver->imax <= 40 || solver->jmax <= 40)) {
for (int i = 0; i < 5; i++) { for (int i = 0; i < 5; i++) {
smoothRB(solver); smoothRB(solver);
} }
@ -692,18 +693,7 @@ void multiGrid(Solver* solver)
for (int i = 0; i < 5; i++) { for (int i = 0; i < 5; i++) {
smoothRB(solver); smoothRB(solver);
if (solver->currentlevel == 0) { if (solver->currentlevel == 0) {
setPressureBoundaryCondition(solver);
double* p = solver->p;
for (int i = 1; i < imax + 1; i++) {
P(i, 0) = P(i, 1);
P(i, jmax + 1) = P(i, jmax);
}
for (int j = 1; j < jmax + 1; j++) {
P(0, j) = P(1, j);
P(imax + 1, j) = P(imax, j);
}
} }
} }
Solver coarseSolver = copySolver(solver); Solver coarseSolver = copySolver(solver);
@ -726,46 +716,26 @@ void multiGrid(Solver* solver)
correct(solver); correct(solver);
if (solver->currentlevel == 0) { if (solver->currentlevel == 0) {
setPressureBoundaryCondition(solver);
double* p = solver->p;
for (int i = 1; i < imax + 1; i++) {
P(i, 0) = P(i, 1);
P(i, jmax + 1) = P(i, jmax);
}
for (int j = 1; j < jmax + 1; j++) {
P(0, j) = P(1, j);
P(imax + 1, j) = P(imax, j);
}
} }
for (int i = 0; i < 5; i++) { for (int i = 0; i < 5; i++) {
res = smoothRB(solver); res = smoothRB(solver);
if (solver->currentlevel == 0) { if (solver->currentlevel == 0) {
setPressureBoundaryCondition(solver);
double* p = solver->p;
for (int i = 1; i < imax + 1; i++) {
P(i, 0) = P(i, 1);
P(i, jmax + 1) = P(i, jmax);
}
for (int j = 1; j < jmax + 1; j++) {
P(0, j) = P(1, j);
P(imax + 1, j) = P(imax, j);
}
} }
} }
#ifdef VERBOSE #ifdef VERBOSE
if (solver->currentlevel == 0) {printf("Residuum: %.6f\n", res); } if (solver->currentlevel == 0) {
printf("Residuum: %.6f\n", res);
}
#endif #endif
} }
void restrictMG(Solver* solver) void restrictMG(Solver* solver)
{ {
int imax = solver->imax; int imax = solver->imax / 2;
int jmax = solver->jmax; int jmax = solver->jmax / 2;
double* r = solver->r[solver->currentlevel + 1]; double* r = solver->r[solver->currentlevel + 1];
double* oldr = solver->r[solver->currentlevel]; double* oldr = solver->r[solver->currentlevel];
@ -892,3 +862,19 @@ void writeResult(Solver* solver)
fclose(fp); fclose(fp);
} }
void setPressureBoundaryCondition(Solver* solver)
{
double* p = solver->p;
int imax = solver->imax;
int jmax = solver->jmax;
for (int i = 1; i < imax + 1; i++) {
P(i, 0) = P(i, 1);
P(i, jmax + 1) = P(i, jmax);
}
for (int j = 1; j < jmax + 1; j++) {
P(0, j) = P(1, j);
P(imax + 1, j) = P(imax, j);
}
}

View File

@ -59,6 +59,7 @@ extern void computeTimestep(Solver*);
extern void setBoundaryConditions(Solver*); extern void setBoundaryConditions(Solver*);
extern void setSpecialBoundaryCondition(Solver*); extern void setSpecialBoundaryCondition(Solver*);
extern void setObjectBoundaryCondition(Solver*); extern void setObjectBoundaryCondition(Solver*);
extern void setPressureBoundaryCondition(Solver*);
extern void computeFG(Solver*); extern void computeFG(Solver*);
extern void adaptUV(Solver*); extern void adaptUV(Solver*);
extern void writeResult(Solver*); extern void writeResult(Solver*);

File diff suppressed because it is too large Load Diff

Binary file not shown.

Before

Width:  |  Height:  |  Size: 117 KiB

After

Width:  |  Height:  |  Size: 17 KiB

View File

@ -49,7 +49,7 @@ gamma 0.9 # upwind differencing factor gamma
# ----------------------- # -----------------------
numberOfParticles 60 numberOfParticles 60
startTime 0 startTime 100
injectTimePeriod 2.0 injectTimePeriod 2.0
writeTimePeriod 0.5 writeTimePeriod 0.5

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

Binary file not shown.

Before

Width:  |  Height:  |  Size: 117 KiB

After

Width:  |  Height:  |  Size: 34 KiB

View File

@ -1,5 +1,5 @@
# Supported: GCC, CLANG, ICC # Supported: GCC, CLANG, ICC
TAG ?= CLANG TAG ?= ICC
ENABLE_OPENMP ?= false ENABLE_OPENMP ?= false
#Feature options #Feature options

Binary file not shown.

View File

@ -1,4 +1,4 @@
CC = icc CC = mpiicc
GCC = gcc GCC = gcc
LINKER = $(CC) LINKER = $(CC)