3D Seq MultiGrid complete

BasicSolver/2D-seq-multigrid/backstep.par

@@ -15,7 +15,7 @@ bcRight    3			#
 gx    0.0      # Body forces (e.g. gravity)
 gy    0.0      #
 
-re    7500.0	   # Reynolds number
+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
@@ -32,7 +32,7 @@ jmax       45	   # number of interior cells in y-direction
 # Time Data:
 # ---------
 
-te      100.0   # final time
+te      60.0   # final time
 dt      0.02    # time stepsize
 tau     0.5     # safety factor for time stepsize control (<0 constant delt)
 
@@ -44,7 +44,6 @@ 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:
 # -----------------------
@@ -66,8 +65,8 @@ y2                  1.5
 shape               1 
 xCenter             0.0
 yCenter             0.0
-xRectLength         1.0
-yRectLength         2.0
+xRectLength         2.0
+yRectLength         1.0
 circleRadius        1.0
 
 #===============================================================================

BasicSolver/2D-seq-multigrid/dcavity.par

@@ -50,7 +50,7 @@ levels   5         # Multigrid levels
 # -----------------------
 
 numberOfParticles   200
-startTime           0
+startTime           100
 injectTimePeriod    0.5
 writeTimePeriod     0.2
 

BasicSolver/2D-seq-multigrid/src/solver.c

@@ -201,23 +201,6 @@ void initSolver(Solver* solver, Parameter* params)
     if (params->shape != NOSHAPE) {
         for (int j = 1; j < jmax + 1; j++) {
             for (int i = 1; i < imax + 1; i++) {
-                // if( S(i+1,j+1) == NONE && S(i+1,j) == NONE && S(i,j+1) == NONE &&
-                // S(i-1,j-1) == LOCAL && S(i,j) == LOCAL ) S(i,j) = BOTTOMRIGHT; else if(
-                // S(i-1,j+1) == NONE && S(i-1,j) == NONE && S(i,j+1) == NONE &&
-                // S(i+1,j-1) == LOCAL && S(i,j) == LOCAL ) S(i,j) = BOTTOMLEFT; else if(
-                // S(i+1,j-1) == NONE
-                // && S(i,j-1) == NONE && S(i+1,j) == NONE && S(i-1,j+1) == LOCAL &&
-                // S(i,j) == LOCAL ) S(i,j) = TOPRIGHT; else if( S(i-1,j-1) == NONE &&
-                // S(i,j-1) == NONE
-                // && S(i-1,j) == NONE && S(i+1,j+1) == LOCAL && S(i,j) == LOCAL ) S(i,j)
-                // = TOPLEFT; else if( S(i+1,j) == NONE && S(i-1,j) == LOCAL && S(i,j) ==
-                // LOCAL ) S(i,j) = RIGHT; else if( S(i,j+1) == NONE && S(i,j-1) == LOCAL
-                // && S(i,j)
-                // == LOCAL ) S(i,j) = BOTTOM; else if( S(i-1,j) == NONE && S(i+1,j) ==
-                // LOCAL
-                // && S(i,j) == LOCAL ) S(i,j) = LEFT; else if( S(i,j-1) == NONE &&
-                // S(i,j+1)
-                // == LOCAL && S(i,j) == LOCAL ) S(i,j) = TOP;
 
                 if (S(i, j - 1) == NONE && S(i, j + 1) == LOCAL && S(i, j) == LOCAL)
                     S(i, j) = BOTTOM; // TOP

BasicSolver/2D-seq-multigrid/src/solver.h

@@ -49,7 +49,6 @@ typedef struct {
 extern void initSolver(Solver*, Parameter*);
 extern void computeRHS(Solver*);
 extern double smoothRB(Solver*);
-extern void residualRB(Solver*);
 extern void restrictMG(Solver*);
 extern void prolongate(Solver*);
 extern void correct(Solver*);

BasicSolver/2D-seq/backstep.par

@@ -15,7 +15,7 @@ bcRight    3			#
 gx    0.0      # Body forces (e.g. gravity)
 gy    0.0      #
 
-re    7500.0	   # Reynolds number
+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
@@ -49,7 +49,7 @@ gamma    0.9       # upwind differencing factor gamma
 # -----------------------
 
 numberOfParticles   200
-startTime           0
+startTime           100
 injectTimePeriod    1.0
 writeTimePeriod     0.5
 
@@ -65,8 +65,8 @@ y2                  1.5
 shape               1 
 xCenter             0.0
 yCenter             0.0
-xRectLength         1.0
-yRectLength         2.0
+xRectLength         2.0
+yRectLength         1.0
 circleRadius        1.0
 
 #===============================================================================

BasicSolver/2D-seq/dcavity.par

@@ -44,12 +44,13 @@ eps      0.001		# stopping tolerance for pressure iteration
 rho      0.5
 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
+startTime           100
 injectTimePeriod    0.5
 writeTimePeriod     0.2
 

BasicSolver/2D-seq/karman.par

@@ -44,12 +44,13 @@ 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           50
+startTime           201
 injectTimePeriod    1.0
 writeTimePeriod     0.5
 

BasicSolver/2D-seq/vector.plot

@@ -1,7 +1,6 @@
 set terminal png size 3600,768 enhanced font ,28
 set output 'velocity.png'
-set xrange[0:7]
-set yrange[0:1.5]
+
 set size ratio -1
 set datafile separator whitespace
 

BasicSolver/3D-seq-multigrid/Makefile

@@ -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)

BasicSolver/3D-seq-multigrid/README.md

@@ -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
+```

BasicSolver/3D-seq-multigrid/backstep.par

@@ -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           105
+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
+
+#===============================================================================

BasicSolver/3D-seq-multigrid/canal.par

@@ -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            100.0	   # Reynolds number
+
+u_init        1.0      # initial value for velocity in x-direction
+v_init        0.0      # initial value for velocity in y-direction
+w_init        0.0      # initial value for velocity in z-direction
+p_init        0.0      # initial value for pressure
+
+# Geometry Data:
+# -------------
+
+xlength       30.0     # domain size in x-direction
+ylength       4.0	   # domain size in y-direction
+zlength       4.0	   # domain size in z-direction
+imax          100      # 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
+#===============================================================================

BasicSolver/3D-seq-multigrid/config.mk

@@ -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

BasicSolver/3D-seq-multigrid/dcavity.par

@@ -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       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
+
+# Particle Tracing Data:
+# -----------------------
+
+numberOfParticles   30
+startTime           100
+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               0 
+xCenter             10.0
+yCenter             2.0
+zCenter             2.0
+xRectLength         8.0
+yRectLength         2.0
+zRectLength         2.0
+circleRadius        1.0
+#===============================================================================

BasicSolver/3D-seq-multigrid/include_CLANG.mk

@@ -0,0 +1,17 @@
+CC   = clang
+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 -Weverything
+#CFLAGS   = -Ofast -fnt-store=aggressive  -std=c99 $(OPENMP) #AMD CLANG
+LFLAGS   = $(OPENMP) -lm
+DEFINES  = -D_GNU_SOURCE# -DDEBUG
+INCLUDES =

BasicSolver/3D-seq-multigrid/include_GCC.mk

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

BasicSolver/3D-seq-multigrid/include_ICC.mk

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

BasicSolver/3D-seq-multigrid/karman.par

@@ -0,0 +1,82 @@
+#==============================================================================
+#                            Laminar 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           500
+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
+
+#===============================================================================

BasicSolver/3D-seq-multigrid/src/allocate.c

@@ -0,0 +1,38 @@
+/*
+ * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
+ * All rights reserved.
+ * Use of this source code is governed by a MIT-style
+ * license that can be found in the LICENSE file.
+ */
+#include <errno.h>
+#include <stddef.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "allocate.h"
+
+void* allocate(size_t alignment, size_t bytesize)
+{
+    int errorCode;
+    void* ptr;
+
+    errorCode = posix_memalign(&ptr, alignment, bytesize);
+
+    if (errorCode) {
+        if (errorCode == EINVAL) {
+            fprintf(stderr, "Error: Alignment parameter is not a power of two\n");
+            exit(EXIT_FAILURE);
+        }
+        if (errorCode == ENOMEM) {
+            fprintf(stderr, "Error: Insufficient memory to fulfill the request\n");
+            exit(EXIT_FAILURE);
+        }
+    }
+
+    if (ptr == NULL) {
+        fprintf(stderr, "Error: posix_memalign failed!\n");
+        exit(EXIT_FAILURE);
+    }
+
+    return ptr;
+}

BasicSolver/3D-seq-multigrid/src/allocate.h

@@ -0,0 +1,13 @@
+/*
+ * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
+ * All rights reserved.
+ * Use of this source code is governed by a MIT-style
+ * license that can be found in the LICENSE file.
+ */
+#ifndef __ALLOCATE_H_
+#define __ALLOCATE_H_
+#include <stdlib.h>
+
+extern void* allocate(size_t alignment, size_t bytesize);
+
+#endif

BasicSolver/3D-seq-multigrid/src/grid.h

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

BasicSolver/3D-seq-multigrid/src/likwid-marker.h

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

BasicSolver/3D-seq-multigrid/src/main.c

@@ -0,0 +1,153 @@
+/*
+ * 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 "allocate.h"
+#include "parameter.h"
+#include "particletracing.h"
+#include "progress.h"
+#include "solver.h"
+#include "timing.h"
+#include "vtkWriter.h"
+
+#define G(v, i, j, k) v[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)]
+
+enum VARIANT { SOR = 1, RB, RBA };
+
+static void createBulkArrays(Solver* s, double* pg, double* ug, double* vg, double* wg)
+{
+    int imax = s->grid.imax;
+    int jmax = s->grid.jmax;
+    int kmax = s->grid.kmax;
+    int idx  = 0;
+
+    for (int k = 1; k < kmax + 1; k++) {
+        for (int j = 1; j < jmax + 1; j++) {
+            for (int i = 1; i < imax + 1; i++) {
+                pg[idx++] = G(s->p, i, j, k);
+            }
+        }
+    }
+
+    idx = 0;
+
+    for (int k = 1; k < kmax + 1; k++) {
+        for (int j = 1; j < jmax + 1; j++) {
+            for (int i = 1; i < imax + 1; i++) {
+                ug[idx++] = (G(s->u, i, j, k) + G(s->u, i - 1, j, k)) / 2.0;
+            }
+        }
+    }
+
+    idx = 0;
+
+    for (int k = 1; k < kmax + 1; k++) {
+        for (int j = 1; j < jmax + 1; j++) {
+            for (int i = 1; i < imax + 1; i++) {
+                vg[idx++] = (G(s->v, i, j, k) + G(s->v, i, j - 1, k)) / 2.0;
+            }
+        }
+    }
+
+    idx = 0;
+
+    for (int k = 1; k < kmax + 1; k++) {
+        for (int j = 1; j < jmax + 1; j++) {
+            for (int i = 1; i < imax + 1; i++) {
+                wg[idx++] = (G(s->w, i, j, k) + G(s->w, i, j, k - 1)) / 2.0;
+            }
+        }
+    }
+}
+
+int main(int argc, char** argv)
+{
+    int rank;
+
+    double timeStart, timeStop;
+    Parameter params;
+    Solver s;
+    ParticleTracer particletracer;
+    initParameter(&params);
+
+    if (argc < 2) {
+        printf("Usage: %s <configFile>\n", argv[0]);
+        exit(EXIT_SUCCESS);
+    }
+
+    readParameter(&params, argv[1]);
+
+
+    printParameter(&params);
+    initSolver(&s, &params);
+    initParticleTracer(&particletracer, &params);
+    printParticleTracerParameters(&particletracer);
+#ifndef VERBOSE
+    initProgress(s.te);
+#endif
+
+    double tau = s.tau;
+    double te  = s.te;
+    double t   = 0.0;
+
+    // printGrid(&s, s.seg);
+    // exit(0)
+
+    timeStart = getTimeStamp();
+
+    while (t <= te) {
+        if (tau > 0.0) computeTimestep(&s);
+        setBoundaryConditions(&s);
+        setSpecialBoundaryCondition(&s);
+        setObjectBoundaryCondition(&s);
+
+        computeFG(&s);
+        computeRHS(&s);
+
+        multiGrid(&s);
+
+        adaptUV(&s);
+
+        trace(&particletracer, s.u, s.v, s.w, s.seg, t);
+
+        t += s.dt;
+
+#ifdef VERBOSE
+        if (rank == 0) {
+            printf("TIME %f , TIMESTEP %f\n", t, s.dt);
+        }
+#else
+        printProgress(t);
+#endif
+    }
+    timeStop = getTimeStamp();
+#ifndef VERBOSE
+    stopProgress();
+#endif
+    printf("Solution took %.2fs\n", timeStop - timeStart);
+
+    double *pg, *ug, *vg, *wg;
+
+    size_t bytesize = (size_t)(s.grid.imax * s.grid.jmax * s.grid.kmax) * sizeof(double);
+
+    pg = allocate(64, bytesize);
+    ug = allocate(64, bytesize);
+    vg = allocate(64, bytesize);
+    wg = allocate(64, bytesize);
+
+    createBulkArrays(&s, pg, ug, vg, wg);
+    VtkOptions opts = { .grid = s.grid };
+    vtkOpen(&opts, s.problem);
+    vtkScalar(&opts, "pressure", pg);
+    vtkVector(&opts, "velocity", (VtkVector) { ug, vg, wg });
+    vtkClose(&opts);
+    return EXIT_SUCCESS;
+}

BasicSolver/3D-seq-multigrid/src/parameter.c

@@ -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(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_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);
+}

BasicSolver/3D-seq-multigrid/src/parameter.h

@@ -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 {
+    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;
+
+extern void initParameter(Parameter*);
+extern void readParameter(Parameter*, const char*);
+extern void printParameter(Parameter*);
+#endif

BasicSolver/3D-seq-multigrid/src/particletracing.c

@@ -0,0 +1,325 @@
+/*
+ * 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 "vtkWriter.h"
+
+#define U(i, j, k) u[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)]
+#define V(i, j, k) v[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)]
+#define W(i, j, k) w[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)]
+#define S(i, j, k) seg[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)]
+
+static int ts     = 0;
+unsigned int seed = 32767;
+void printParticles(ParticleTracer* particletracer)
+{
+    for (int i = 0; i < particletracer->totalParticles; ++i) {
+        printf("Particle position X : %.2f, Y : %.2f, flag : %d\n",
+            particletracer->particlePool[i].x,
+            particletracer->particlePool[i].y,
+            particletracer->particlePool[i].flag);
+    }
+}
+void injectParticles(ParticleTracer* particletracer, int* seg)
+{
+
+    int imax = particletracer->imax;
+    int jmax = particletracer->jmax;
+    int kmax = particletracer->kmax;
+
+    for (int i = 0; i < particletracer->numberOfParticles; ++i) {
+
+        particletracer->particlePool[particletracer->pointer].x = particletracer->x1;
+        particletracer->particlePool[particletracer->pointer].y = (((double)rand() /
+                                                                  RAND_MAX) *
+                                                                  (particletracer->y2 - particletracer->y1)) +
+                                                                  particletracer->y1;
+        particletracer->particlePool[particletracer->pointer].z = (((double)rand() /
+                                                                  RAND_MAX) *
+                                                                  (particletracer->z2 - particletracer->z1)) +
+                                                                  particletracer->z1;
+
+        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+1, j, k) == NONE)
+        {
+            particletracer->particlePool[particletracer->pointer].flag = true;
+            ++(particletracer->pointer);
+            ++(particletracer->totalParticles);
+        }
+        else   particletracer->particlePool[particletracer->pointer].flag = false;
+    }
+}
+
+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;
+
+    double dx = particletracer->dx;
+    double dy = particletracer->dy;
+    double dz = particletracer->dz;
+
+    double xlength = particletracer->xlength;
+    double ylength = particletracer->ylength;
+    double zlength = particletracer->zlength;
+
+    for (int i = 0; i < particletracer->totalParticles; ++i) {
+        if (particletracer->particlePool[i].flag == true) {
+            double x = particletracer->particlePool[i].x;
+            double y = particletracer->particlePool[i].y;
+            double z = particletracer->particlePool[i].z;
+
+            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->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->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->dt * w_n;
+            particletracer->particlePool[i].z = new_z;
+
+            // printf("\tOld X : %.2f, New X : %.2f, iCoord : %d\n\tOld Y : %.2f, New Y :
+            // %.2f, jCoord : %d\n\n", x, new_x, iCoord, y, 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);
+            // printf("\t j-1 * (imax + 2) + i-1 = %d with element from U : %.2f", (jCoord
+            // - 1) * (200 + 2) + (iCoord - 1), u[(jCoord - 1) * (imax + 2) + (iCoord -
+            // 1)]); printf("\nimax : %d, jmax : %d\n", imax, jmax);
+
+            if (((new_x < 0.0) || (new_x > xlength) || (new_y < 0.0) ||
+                    (new_y > ylength) || (new_z < 0.0) || (new_z > zlength))) {
+                particletracer->particlePool[i].flag = false;
+            }
+            int i_new = new_x / dx, j_new = new_y / dy, k_new = new_z / dz;
+
+            if (S(i_new, j_new, k_new) != NONE) {
+                particletracer->particlePool[i].flag = false;
+            }
+        }
+    }
+}
+
+void freeParticles(ParticleTracer* particletracer)
+{
+    free(particletracer->particlePool);
+    free(particletracer->linSpaceLine);
+}
+
+void writeParticles(ParticleTracer* particletracer)
+{
+    VtkOptions opts = { .particletracer = particletracer };
+
+    char filename[50];
+    snprintf(filename, 50, "vtk_files/particles%d.vtk", ts);
+    vtkOpenPT(&opts, filename, ts);
+    vtkParticle(&opts, "particle");
+    vtkClose(&opts);
+    ++ts;
+}
+
+void initParticleTracer(ParticleTracer* particletracer, Parameter* 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->estimatedNumParticles = ((params->te - params->startTime) + 2) *
+                                            params->numberOfParticles;
+
+    particletracer->particlePool = malloc(
+        sizeof(Particle) * particletracer->estimatedNumParticles);
+    particletracer->linSpaceLine = malloc(
+        sizeof(Particle) * particletracer->numberOfParticles);
+}
+
+void printParticleTracerParameters(ParticleTracer* particletracer)
+{
+    printf("Particle Tracing data:\n");
+    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);
+}
+
+void trace(ParticleTracer* particletracer,
+    double* u,
+    double* v,
+    double* w,
+    int* seg,
+    double time)
+{
+    if (time >= particletracer->startTime) {
+        // printParticles(particletracer);
+        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];
+    int totalParticles = 0;
+
+    for (int i = 0; i < particletracer->totalParticles; ++i) {
+        if (memPool[i].flag == 1) {
+            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 + 1;
+
+    memcpy(particletracer->particlePool, tempPool, totalParticles * sizeof(Particle));
+}

BasicSolver/3D-seq-multigrid/src/particletracing.h

@@ -0,0 +1,48 @@
+/*
+ * 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 "particletracing.h"
+#include "solver.h"
+#include <stdbool.h>
+
+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, activeParticles;
+
+    double dx, dy, dz, dt;
+    Particle* linSpaceLine;
+    Particle* particlePool;
+
+    int pointer;
+
+    double imax, jmax, kmax, xlength, ylength, zlength;
+
+    double x1, y1, x2, y2, z1, z2;
+} ParticleTracer;
+
+extern void initParticleTracer(ParticleTracer*, Parameter*);
+extern void injectParticles(ParticleTracer*, int* seg);
+extern void advanceParticles(ParticleTracer*, double*, double*, double*, int*, double);
+extern void freeParticles(ParticleTracer*);
+extern void writeParticles(ParticleTracer*);
+extern void printParticleTracerParameters(ParticleTracer*);
+extern void printParticles(ParticleTracer*);
+extern void compress(ParticleTracer*);
+extern void trace(ParticleTracer*, double*, double*, double*, int*, double);
+#endif

BasicSolver/3D-seq-multigrid/src/progress.c

@@ -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);
+}

BasicSolver/3D-seq-multigrid/src/progress.h

@@ -0,0 +1,14 @@
+/*
+ * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
+ * All rights reserved.
+ * Use of this source code is governed by a MIT-style
+ * license that can be found in the LICENSE file.
+ */
+#ifndef __PROGRESS_H_
+#define __PROGRESS_H_
+
+extern void initProgress(double);
+extern void printProgress(double);
+extern void stopProgress(void);
+
+#endif

BasicSolver/3D-seq-multigrid/src/solver.h

@@ -0,0 +1,95 @@
+/*
+ * 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 "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,
+    /*Ghost cells boundary */
+    OUTSIDEBOUNDARY
+};
+enum BC { NOSLIP = 1, SLIP, OUTFLOW, PERIODIC };
+/// @brief
+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;
+    /* 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;
+} 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 setObjectBoundaryCondition(Solver*);
+extern void setPressureBoundaryCondition(Solver*);
+extern void setSpecialBoundaryCondition(Solver*);
+extern void computeFG(Solver*);
+extern void adaptUV(Solver*);
+extern void writeResult(Solver*);
+extern void printGrid(Solver*, int*);
+
+#endif

BasicSolver/3D-seq-multigrid/src/timing.c

@@ -0,0 +1,22 @@
+/*
+ * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
+ * All rights reserved.
+ * Use of this source code is governed by a MIT-style
+ * license that can be found in the LICENSE file.
+ */
+#include <stdlib.h>
+#include <time.h>
+
+double getTimeStamp(void)
+{
+    struct timespec ts;
+    clock_gettime(CLOCK_MONOTONIC, &ts);
+    return (double)ts.tv_sec + (double)ts.tv_nsec * 1.e-9;
+}
+
+double getTimeResolution(void)
+{
+    struct timespec ts;
+    clock_getres(CLOCK_MONOTONIC, &ts);
+    return (double)ts.tv_sec + (double)ts.tv_nsec * 1.e-9;
+}

BasicSolver/3D-seq-multigrid/src/timing.h

@@ -0,0 +1,13 @@
+/*
+ * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
+ * All rights reserved.
+ * Use of this source code is governed by a MIT-style
+ * license that can be found in the LICENSE file.
+ */
+#ifndef __TIMING_H_
+#define __TIMING_H_
+
+extern double getTimeStamp(void);
+extern double getTimeResolution(void);
+
+#endif // __TIMING_H_

BasicSolver/3D-seq-multigrid/src/util.h

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

BasicSolver/3D-seq-multigrid/src/vtkWriter.c

@@ -0,0 +1,198 @@
+/*
+ * 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 "vtkWriter.h"
+#define G(v, i, j, k) v[(k) * imax * jmax + (j) * imax + (i)]
+
+static float floatSwap(float f)
+{
+    union {
+        float f;
+        char b[4];
+    } dat1, dat2;
+
+    dat1.f    = f;
+    dat2.b[0] = dat1.b[3];
+    dat2.b[1] = dat1.b[2];
+    dat2.b[2] = dat1.b[1];
+    dat2.b[3] = 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];
+    snprintf(filename, 50, "%s.vtk", problem);
+    o->fh = fopen(filename, "w");
+    writeHeader(o);
+
+    printf("Writing VTK output for %s\n", problem);
+}
+
+void vtkScalar(VtkOptions* o, char* name, double* s)
+{
+    int imax = o->grid.imax;
+    int jmax = o->grid.jmax;
+    int kmax = o->grid.kmax;
+
+    printf("Register scalar %s\n", name);
+
+    if (o->fh == NULL) {
+        printf("vtkWriter not initialize! Call vtkOpen first!\n");
+        exit(EXIT_FAILURE);
+    }
+    fprintf(o->fh, "SCALARS %s float 1\n", name);
+    fprintf(o->fh, "LOOKUP_TABLE default\n");
+
+    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((float[1]) { floatSwap(G(s, i, j, k)) },
+                        sizeof(float),
+                        1,
+                        o->fh);
+                }
+            }
+        }
+    }
+    if (o->fmt == BINARY) fprintf(o->fh, "\n");
+}
+
+void vtkVector(VtkOptions* o, char* name, VtkVector vec)
+{
+    int imax = o->grid.imax;
+    int jmax = o->grid.jmax;
+    int kmax = o->grid.kmax;
+
+    if (o->fh == NULL) {
+        printf("vtkWriter not initialize! Call vtkOpen first!\n");
+        exit(EXIT_FAILURE);
+    }
+
+    fprintf(o->fh, "VECTORS %s float\n", name);
+
+    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 /*&& k >= 20*/) {
+                    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((float[3]) { floatSwap(G(vec.u, i, j, k)),
+                               floatSwap(G(vec.v, i, j, k)),
+                               floatSwap(G(vec.w, i, j, k)) },
+                        sizeof(float),
+                        3,
+                        o->fh);
+                }
+            }
+        }
+    }
+    if (o->fmt == BINARY) fprintf(o->fh, "\n");
+}
+
+void vtkClose(VtkOptions* o)
+{
+    fclose(o->fh);
+    o->fh = NULL;
+}
+
+static void writeHeaderPT(VtkOptions* o, int ts)
+{
+    fprintf(o->fh, "# vtk DataFile Version 3.0\n");
+    fprintf(o->fh, "PAMPI cfd solver particle tracing file\n");
+    if (o->fmt == ASCII) {
+        fprintf(o->fh, "ASCII\n");
+    } else if (o->fmt == BINARY) {
+        fprintf(o->fh, "BINARY\n");
+    }
+
+    fprintf(o->fh, "DATASET UNSTRUCTURED_GRID\n");
+    fprintf(o->fh, "FIELD FieldData 2\n");
+    fprintf(o->fh, "TIME 1 1 double\n");
+    fprintf(o->fh, "%d\n", ts);
+    fprintf(o->fh, "CYCLE 1 1 int\n");
+    fprintf(o->fh, "1\n");
+}
+
+void vtkOpenPT(VtkOptions* o, char* problem, int ts)
+{
+    o->fh = fopen(problem, "w");
+
+    if (o->fh == NULL) {
+        printf("vtkWriter not initialize! Call vtkOpen first!\n");
+        exit(EXIT_FAILURE);
+    }
+
+    writeHeaderPT(o, ts);
+
+    printf("Writing VTK output for %s\n", problem);
+}
+
+void vtkParticle(VtkOptions* o, char* name)
+{
+    Particle* particlePool = o->particletracer->particlePool;
+
+    if (o->fh == NULL) {
+        printf("vtkWriter not initialize! Call vtkOpen first!\n");
+        exit(EXIT_FAILURE);
+    }
+
+    fprintf(o->fh, "POINTS %d float\n", o->particletracer->totalParticles);
+
+    for (int i = 0; i < o->particletracer->totalParticles; ++i) {
+        double x = particlePool[i].x;
+        double y = particlePool[i].y;
+        double z = particlePool[i].z;
+        fprintf(o->fh, "%.2f %.2f %.2f\n", x, y, z);
+    }
+
+    fprintf(o->fh,
+        "CELLS %d %d\n",
+        o->particletracer->totalParticles,
+        2 * o->particletracer->totalParticles);
+
+    for (int i = 0; i < o->particletracer->totalParticles; ++i) {
+        fprintf(o->fh, "1 %d\n", i);
+    }
+
+    fprintf(o->fh, "CELL_TYPES %d\n", o->particletracer->totalParticles);
+
+    for (int i = 0; i < o->particletracer->totalParticles; ++i) {
+        fprintf(o->fh, "1\n");
+    }
+}

BasicSolver/3D-seq-multigrid/src/vtkWriter.h

@@ -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 "particletracing.h"
+#include "grid.h"
+
+typedef enum VtkFormat { ASCII = 0, BINARY } VtkFormat;
+
+typedef struct VtkOptions {
+    VtkFormat fmt;
+    Grid grid;
+    FILE* fh;
+    ParticleTracer* particletracer;
+} 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);
+
+extern void vtkOpenPT(VtkOptions* opts, char* filename, int ts);
+extern void vtkParticle(VtkOptions* opts, char* name);
+#endif // __VTKWRITER_H_

BasicSolver/3D-seq/canal.par

@@ -18,7 +18,7 @@ gx     0.0      # Body forces (e.g. gravity)
 gy     0.0      #
 gz     0.0      #
 
-re            5000.0	   # Reynolds number
+re            100.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
@@ -38,7 +38,7 @@ kmax          40	   # number of interior cells in z-direction
 # Time Data:
 # ---------
 
-te       200.0   # final time
+te       60.0   # final time
 dt       0.02    # time stepsize
 tau      0.5     # safety factor for time stepsize control (<0 constant delt)
 
@@ -55,7 +55,7 @@ gamma         0.9       # upwind differencing factor gamma
 # -----------------------
 
 numberOfParticles   200
-startTime           0
+startTime           100
 injectTimePeriod    2.0
 writeTimePeriod     1.0
 
@@ -70,7 +70,7 @@ z2                  1.0
 # -----------------------
 # Shape 0 disable, 1 Rectangle/Square, 2 Circle
 
-shape               1 
+shape               0 
 xCenter             10.0
 yCenter             2.0
 zCenter             2.0