Fix bug. Add output for compression.

BasicSolver/2D-seq/src/affinity.c

@@ -1,61 +0,0 @@
-/*
- * Copyright (C)  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__*/

BasicSolver/2D-seq/src/affinity.h

@@ -1,14 +0,0 @@
-/*
- * Copyright (C)  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*/

BasicSolver/2D-seq/src/main.c

@@ -20,6 +20,7 @@ int main(int argc, char** argv)
     Parameter p;
     Discretization d;
     Solver s;
+    
     initParameter(&p);
 
     if (argc != 2) {
@@ -31,6 +32,7 @@ int main(int argc, char** argv)
     printParameter(&p);
     initDiscretization(&d, &p);
     initSolver(&s, &d, &p);
+    
 #ifndef VERBOSE
     initProgress(d.te);
 #endif

BasicSolver/2D-seq/src/parameter.c

@@ -9,22 +9,20 @@
 #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->imax    = 128;
-    param->jmax    = 100;
+    param->jmax    = 128;
     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;
 }
 
@@ -81,7 +79,6 @@ void readParameter(Parameter* param, const char* filename)
             PARSE_REAL(u_init);
             PARSE_REAL(v_init);
             PARSE_REAL(p_init);
-            PARSE_REAL(rho);
         }
     }
 
@@ -112,6 +109,5 @@ void printParameter(Parameter* param)
     printf("\tepsilon (stopping tolerance) : %f\n", param->eps);
     printf("\tgamma (stopping tolerance) : %f\n", param->gamma);
     printf("\tomega (SOR relaxation): %f\n", param->omg);
-    printf("\trho (SOR relaxation): %f\n", param->rho);
     printf("\tMultiGrid levels : %d\n", param->levels);
 }

BasicSolver/2D-seq/src/progress.h

@@ -9,6 +9,6 @@
 
 extern void initProgress(double);
 extern void printProgress(double);
-extern void stopProgress();
+extern void stopProgress(void);
 
 #endif

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

@@ -8,12 +8,11 @@
 #include <stdlib.h>
 
 #include "allocate.h"
-#include "discretization.h"
-#include "parameter.h"
 #include "solver.h"
 #include "util.h"
 
 #define FINEST_LEVEL   0
+#define COARSEST_LEVEL (s->levels - 1)
 #define S(i, j)        s[(j) * (imax + 2) + (i)]
 #define E(i, j)        e[(j) * (imax + 2) + (i)]
 #define R(i, j)        r[(j) * (imax + 2) + (i)]
@@ -51,6 +50,7 @@ static void prolongate(Solver* s, int level, int imax, int jmax)
 static void correct(Solver* s, double* p, int level, int imax, int jmax)
 {
     double* e = s->e[level];
+
     for (int j = 1; j < jmax + 1; ++j) {
         for (int i = 1; i < imax + 1; ++i) {
             P(i, j) += E(i, j);
@@ -106,12 +106,52 @@ static double smooth(Solver* s, double* p, double* rhs, int level, int imax, int
     return res;
 }
 
+static double multiGrid(Solver* s, double* p, double* rhs, int level, int imax, int jmax)
+{
+    double res = 0.0;
+
+    // coarsest level TODO: Use direct solver?
+    if (level == COARSEST_LEVEL) {
+        for (int i = 0; i < 5; i++) {
+            smooth(s, p, rhs, level, imax, jmax);
+        }
+        return res;
+    }
+
+    // pre-smoothing TODO: Make smoothing steps configurable?
+    for (int i = 0; i < 5; i++) {
+        smooth(s, p, rhs, level, imax, jmax);
+        if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax);
+    }
+
+    // restrict
+    restrictMG(s, level, imax, jmax);
+
+    // MGSolver on residual and error.
+    // TODO: What if there is a rest?
+    multiGrid(s, s->e[level + 1], s->r[level + 1], level + 1, imax / 2, jmax / 2);
+
+    // prolongate
+    prolongate(s, level, imax, jmax);
+
+    // correct p on finer level using residual
+    correct(s, p, level, imax, jmax);
+    if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax);
+
+    // post-smoothing
+    for (int i = 0; i < 5; i++) {
+        res = smooth(s, p, rhs, level, imax, jmax);
+        if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax);
+    }
+
+    return res;
+}
+
 void initSolver(Solver* s, Discretization* d, Parameter* p)
 {
     s->eps     = p->eps;
     s->omega   = p->omg;
     s->itermax = p->itermax;
-    s->rho     = p->rho;
     s->levels  = p->levels;
     s->grid    = &d->grid;
 
@@ -136,52 +176,6 @@ void initSolver(Solver* s, Discretization* d, Parameter* p)
     }
 }
 
-double multiGrid(Solver* solver, double* p, double* rhs, int level, int imax, int jmax)
-{
-    double res = 0.0;
-
-    // coarsest level TODO: Use direct solver?
-    if (level == (solver->levels - 1)) {
-        for (int i = 0; i < 5; i++) {
-            smooth(solver, p, rhs, level, imax, jmax);
-        }
-        return res;
-    }
-
-    // pre-smoothing TODO: Make smoothing steps configurable?
-    for (int i = 0; i < 5; i++) {
-        smooth(solver, p, rhs, level, imax, jmax);
-        if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax);
-    }
-
-    // restrict
-    restrictMG(solver, level, imax, jmax);
-
-    // MGSolver on residual and error.
-    // TODO: What if there is a rest?
-    multiGrid(solver,
-        solver->e[level + 1],
-        solver->r[level],
-        level + 1,
-        imax / 2,
-        jmax / 2);
-
-    // prolongate
-    prolongate(solver, level, imax, jmax);
-
-    // correct p on finer level using residual
-    correct(solver, p, level, imax, jmax);
-    if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax);
-
-    // post-smoothing
-    for (int i = 0; i < 5; i++) {
-        res = smooth(solver, p, rhs, level, imax, jmax);
-        if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax);
-    }
-
-    return res;
-}
-
 void solve(Solver* s, double* p, double* rhs)
 {
     double res = multiGrid(s, p, rhs, 0, s->grid->imax, s->grid->jmax);

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

@@ -4,7 +4,6 @@
  * Use of this source code is governed by a MIT style
  * license that can be found in the LICENSE file.
  */
-#include "discretization.h"
 #include "solver.h"
 #include "util.h"
 

BasicSolver/3D-seq/Makefile

@@ -18,9 +18,10 @@ include $(MAKE_DIR)/include_$(TAG).mk
 INCLUDES  += -I$(SRC_DIR) -I$(BUILD_DIR)
 
 VPATH     = $(SRC_DIR)
-SRC       = $(wildcard $(SRC_DIR)/*.c)
+SRC       = $(filter-out $(wildcard $(SRC_DIR)/*-*.c),$(wildcard $(SRC_DIR)/*.c))
 ASM       = $(patsubst $(SRC_DIR)/%.c, $(BUILD_DIR)/%.s, $(SRC))
 OBJ       = $(patsubst $(SRC_DIR)/%.c, $(BUILD_DIR)/%.o, $(SRC))
+OBJ      += $(BUILD_DIR)/solver-$(SOLVER).o
 SOURCES   = $(SRC) $(wildcard $(SRC_DIR)/*.h)
 CPPFLAGS := $(CPPFLAGS) $(DEFINES) $(OPTIONS) $(INCLUDES)
 
@@ -47,6 +48,8 @@ clean:
 distclean: clean
 	$(info ===>  DIST CLEAN)
 	@rm -f $(TARGET)
+	@rm -f *.dat
+	@rm -f *.png
 
 info:
 	$(info $(CFLAGS))

BasicSolver/3D-seq/config.mk

@@ -1,12 +1,12 @@
 # Supported: GCC, CLANG, ICC
 TAG ?= CLANG
 ENABLE_OPENMP ?= false
+# Supported: sor, mg
+SOLVER ?= mg
+# Run in debug settings
+DEBUG ?= false
 
 #Feature options
 OPTIONS +=  -DARRAY_ALIGNMENT=64
-OPTIONS +=  -DVERBOSE
+#OPTIONS +=  -DVERBOSE
 #OPTIONS +=  -DDEBUG
-#OPTIONS +=  -DBOUNDCHECK
-#OPTIONS +=  -DVERBOSE_AFFINITY
-#OPTIONS +=  -DVERBOSE_DATASIZE
-#OPTIONS +=  -DVERBOSE_TIMER

BasicSolver/3D-seq/dcavity.par

@@ -47,6 +47,8 @@ tau      0.5		# safety factor for time stepsize control (<0 constant delt)
 
 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
 #===============================================================================

BasicSolver/3D-seq/include_CLANG.mk

@@ -2,16 +2,18 @@ CC   = clang
 GCC  = cc
 LINKER = $(CC)
 
-ifeq ($(ENABLE_OPENMP),true)
+ifeq ($(strip $(ENABLE_OPENMP)),true)
 OPENMP   = -fopenmp
 #OPENMP   = -Xpreprocessor -fopenmp #required on Macos with homebrew libomp
 LIBS     = # -lomp
 endif
+ifeq ($(strip $(DEBUG)),true)
+CFLAGS   = -O0 -g -std=c17
+else
+CFLAGS   = -O3 -std=c17 $(OPENMP)
+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
+DEFINES  = -D_GNU_SOURCE
 INCLUDES =

BasicSolver/3D-seq/src/discretization.c

@@ -10,131 +10,120 @@
 #include <string.h>
 
 #include "allocate.h"
+#include "discretization.h"
 #include "parameter.h"
-#include "solver.h"
 #include "util.h"
 
-#define P(i, j, k)   p[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)]
+static void printConfig(Discretization* d)
-#define F(i, j, k)   f[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)]
-#define G(i, j, k)   g[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)]
-#define H(i, j, k)   h[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)]
-#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 RHS(i, j, k) rhs[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)]
-
-static void printConfig(Solver* s)
 {
-    printf("Parameters for #%s#\n", s->problem);
+    printf("Parameters for #%s#\n", d->problem);
     printf("BC Left:%d Right:%d Bottom:%d Top:%d Front:%d Back:%d\n",
-        s->bcLeft,
+        d->bcLeft,
-        s->bcRight,
+        d->bcRight,
-        s->bcBottom,
+        d->bcBottom,
-        s->bcTop,
+        d->bcTop,
-        s->bcFront,
+        d->bcFront,
-        s->bcBack);
+        d->bcBack);
-    printf("\tReynolds number: %.2f\n", s->re);
+    printf("\tReynolds number: %.2f\n", d->re);
-    printf("\tGx Gy: %.2f %.2f %.2f\n", s->gx, s->gy, s->gz);
+    printf("\tGx Gy: %.2f %.2f %.2f\n", d->gx, d->gy, d->gz);
     printf("Geometry data:\n");
     printf("\tDomain box size (x, y, z): %.2f, %.2f, %.2f\n",
-        s->grid.xlength,
+        d->grid.xlength,
-        s->grid.ylength,
+        d->grid.ylength,
-        s->grid.zlength);
+        d->grid.zlength);
-    printf("\tCells (x, y, z): %d, %d, %d\n", s->grid.imax, s->grid.jmax, s->grid.kmax);
+    printf("\tCells (x, y, z): %d, %d, %d\n", d->grid.imax, d->grid.jmax, d->grid.kmax);
-    printf("\tCell size (dx, dy, dz): %f, %f, %f\n", s->grid.dx, s->grid.dy, s->grid.dz);
+    printf("\tCell size (dx, dy, dz): %f, %f, %f\n", d->grid.dx, d->grid.dy, d->grid.dz);
     printf("Timestep parameters:\n");
-    printf("\tDefault stepsize: %.2f, Final time %.2f\n", s->dt, s->te);
+    printf("\tDefault stepsize: %.2f, Final time %.2f\n", d->dt, d->te);
-    printf("\tdt bound: %.6f\n", s->dtBound);
+    printf("\tdt bound: %.6f\n", d->dtBound);
-    printf("\tTau factor: %.2f\n", s->tau);
+    printf("\tTau factor: %.2f\n", d->tau);
     printf("Iterative parameters:\n");
-    printf("\tMax iterations: %d\n", s->itermax);
+    printf("\tepsilon (stopping tolerance) : %f\n", d->eps);
-    printf("\tepsilon (stopping tolerance) : %f\n", s->eps);
+    printf("\tgamma factor: %f\n", d->gamma);
-    printf("\tgamma factor: %f\n", s->gamma);
+    printf("\tomega (SOR relaxation): %f\n", d->omega);
-    printf("\tomega (SOR relaxation): %f\n", s->omega);
 }
 
-void initSolver(Solver* s, Parameter* params)
+void initDiscretization(Discretization* d, Parameter* p)
 {
-    s->problem  = params->name;
+    d->problem  = p->name;
-    s->bcLeft   = params->bcLeft;
+    d->bcLeft   = p->bcLeft;
-    s->bcRight  = params->bcRight;
+    d->bcRight  = p->bcRight;
-    s->bcBottom = params->bcBottom;
+    d->bcBottom = p->bcBottom;
-    s->bcTop    = params->bcTop;
+    d->bcTop    = p->bcTop;
-    s->bcFront  = params->bcFront;
+    d->bcFront  = p->bcFront;
-    s->bcBack   = params->bcBack;
+    d->bcBack   = p->bcBack;
 
-    s->grid.imax    = params->imax;
+    d->grid.imax    = p->imax;
-    s->grid.jmax    = params->jmax;
+    d->grid.jmax    = p->jmax;
-    s->grid.kmax    = params->kmax;
+    d->grid.kmax    = p->kmax;
-    s->grid.xlength = params->xlength;
+    d->grid.xlength = p->xlength;
-    s->grid.ylength = params->ylength;
+    d->grid.ylength = p->ylength;
-    s->grid.zlength = params->zlength;
+    d->grid.zlength = p->zlength;
-    s->grid.dx      = params->xlength / params->imax;
+    d->grid.dx      = p->xlength / p->imax;
-    s->grid.dy      = params->ylength / params->jmax;
+    d->grid.dy      = p->ylength / p->jmax;
-    s->grid.dz      = params->zlength / params->kmax;
+    d->grid.dz      = p->zlength / p->kmax;
 
-    s->eps     = params->eps;
+    d->eps   = p->eps;
-    s->omega   = params->omg;
+    d->omega = p->omg;
-    s->itermax = params->itermax;
+    d->re    = p->re;
-    s->re      = params->re;
+    d->gx    = p->gx;
-    s->gx      = params->gx;
+    d->gy    = p->gy;
-    s->gy      = params->gy;
+    d->gz    = p->gz;
-    s->gz      = params->gz;
+    d->dt    = p->dt;
-    s->dt      = params->dt;
+    d->te    = p->te;
-    s->te      = params->te;
+    d->tau   = p->tau;
-    s->tau     = params->tau;
+    d->gamma = p->gamma;
-    s->gamma   = params->gamma;
 
-    int imax        = s->grid.imax;
+    int imax        = d->grid.imax;
-    int jmax        = s->grid.jmax;
+    int jmax        = d->grid.jmax;
-    int kmax        = s->grid.kmax;
+    int kmax        = d->grid.kmax;
     size_t bytesize = (imax + 2) * (jmax + 2) * (kmax + 2) * sizeof(double);
-    s->u            = allocate(64, bytesize);
+    d->u            = allocate(64, bytesize);
-    s->v            = allocate(64, bytesize);
+    d->v            = allocate(64, bytesize);
-    s->w            = allocate(64, bytesize);
+    d->w            = allocate(64, bytesize);
-    s->p            = allocate(64, bytesize);
+    d->p            = allocate(64, bytesize);
-    s->rhs          = allocate(64, bytesize);
+    d->rhs          = allocate(64, bytesize);
-    s->f            = allocate(64, bytesize);
+    d->f            = allocate(64, bytesize);
-    s->g            = allocate(64, bytesize);
+    d->g            = allocate(64, bytesize);
-    s->h            = allocate(64, bytesize);
+    d->h            = allocate(64, bytesize);
 
     for (int i = 0; i < (imax + 2) * (jmax + 2) * (kmax + 2); i++) {
-        s->u[i]   = params->u_init;
+        d->u[i]   = p->u_init;
-        s->v[i]   = params->v_init;
+        d->v[i]   = p->v_init;
-        s->w[i]   = params->w_init;
+        d->w[i]   = p->w_init;
-        s->p[i]   = params->p_init;
+        d->p[i]   = p->p_init;
-        s->rhs[i] = 0.0;
+        d->rhs[i] = 0.0;
-        s->f[i]   = 0.0;
+        d->f[i]   = 0.0;
-        s->g[i]   = 0.0;
+        d->g[i]   = 0.0;
-        s->h[i]   = 0.0;
+        d->h[i]   = 0.0;
     }
 
-    double dx = s->grid.dx;
+    double dx = d->grid.dx;
-    double dy = s->grid.dy;
+    double dy = d->grid.dy;
-    double dz = s->grid.dz;
+    double dz = d->grid.dz;
 
     double invSqrSum = 1.0 / (dx * dx) + 1.0 / (dy * dy) + 1.0 / (dz * dz);
-    s->dtBound       = 0.5 * s->re * 1.0 / invSqrSum;
+    d->dtBound       = 0.5 * d->re * 1.0 / invSqrSum;
 
 #ifdef VERBOSE
     printConfig(s);
 #endif /* VERBOSE */
 }
 
-void computeRHS(Solver* s)
+void computeRHS(Discretization* d)
 {
-    int imax   = s->grid.imax;
+    int imax   = d->grid.imax;
-    int jmax   = s->grid.jmax;
+    int jmax   = d->grid.jmax;
-    int kmax   = s->grid.kmax;
+    int kmax   = d->grid.kmax;
-    double idx = 1.0 / s->grid.dx;
+    double idx = 1.0 / d->grid.dx;
-    double idy = 1.0 / s->grid.dy;
+    double idy = 1.0 / d->grid.dy;
-    double idz = 1.0 / s->grid.dz;
+    double idz = 1.0 / d->grid.dz;
-    double idt = 1.0 / s->dt;
+    double idt = 1.0 / d->dt;
 
-    double* rhs = s->rhs;
+    double* rhs = d->rhs;
-    double* f   = s->f;
+    double* f   = d->f;
-    double* g   = s->g;
+    double* g   = d->g;
-    double* h   = s->h;
+    double* h   = d->h;
 
     for (int k = 1; k < kmax + 1; k++) {
         for (int j = 1; j < jmax + 1; j++) {
@@ -148,94 +137,9 @@ void computeRHS(Solver* s)
     }
 }
 
-void solve(Solver* s)
+static double maxElement(Discretization* d, double* m)
 {
-    int imax      = s->grid.imax;
+    int size      = (d->grid.imax + 2) * (d->grid.jmax + 2) * (d->grid.kmax + 2);
-    int jmax      = s->grid.jmax;
-    int kmax      = s->grid.kmax;
-    double eps    = s->eps;
-    int itermax   = s->itermax;
-    double dx2    = s->grid.dx * s->grid.dx;
-    double dy2    = s->grid.dy * s->grid.dy;
-    double dz2    = s->grid.dz * s->grid.dz;
-    double idx2   = 1.0 / dx2;
-    double idy2   = 1.0 / dy2;
-    double idz2   = 1.0 / dz2;
-    double factor = s->omega * 0.5 * (dx2 * dy2 * dz2) /
-                    (dy2 * dz2 + dx2 * dz2 + dx2 * dy2);
-    double* p    = s->p;
-    double* rhs  = s->rhs;
-    double epssq = eps * eps;
-    int it       = 0;
-    double res   = 1.0;
-    int pass, ksw, jsw, isw;
-
-    while ((res >= epssq) && (it < itermax)) {
-        res = 0.0;
-        ksw = 1;
-
-        for (pass = 0; pass < 2; pass++) {
-            jsw = ksw;
-
-            for (int k = 1; k < kmax + 1; k++) {
-                isw = jsw;
-                for (int j = 1; j < jmax + 1; j++) {
-                    for (int i = isw; i < imax + 1; i += 2) {
-
-                        double r =
-                            RHS(i, j, k) -
-                            ((P(i + 1, j, k) - 2.0 * P(i, j, k) + P(i - 1, j, k)) * idx2 +
-                                (P(i, j + 1, k) - 2.0 * P(i, j, k) + P(i, j - 1, k)) *
-                                    idy2 +
-                                (P(i, j, k + 1) - 2.0 * P(i, j, k) + P(i, j, k - 1)) *
-                                    idz2);
-
-                        P(i, j, k) -= (factor * r);
-                        res += (r * r);
-                    }
-                    isw = 3 - isw;
-                }
-                jsw = 3 - jsw;
-            }
-            ksw = 3 - ksw;
-        }
-
-        for (int j = 1; j < jmax + 1; j++) {
-            for (int i = 1; i < imax + 1; i++) {
-                P(i, j, 0)        = P(i, j, 1);
-                P(i, j, kmax + 1) = P(i, j, kmax);
-            }
-        }
-
-        for (int k = 1; k < kmax + 1; k++) {
-            for (int i = 1; i < imax + 1; i++) {
-                P(i, 0, k)        = P(i, 1, k);
-                P(i, jmax + 1, k) = P(i, jmax, k);
-            }
-        }
-
-        for (int k = 1; k < kmax + 1; k++) {
-            for (int j = 1; j < jmax + 1; j++) {
-                P(0, j, k)        = P(1, j, k);
-                P(imax + 1, j, k) = P(imax, j, k);
-            }
-        }
-
-        res = res / (double)(imax * jmax * kmax);
-#ifdef DEBUG
-        printf("%d Residuum: %e\n", it, res);
-#endif
-        it++;
-    }
-
-#ifdef VERBOSE
-    printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
-#endif
-}
-
-static double maxElement(Solver* s, double* m)
-{
-    int size      = (s->grid.imax + 2) * (s->grid.jmax + 2) * (s->grid.kmax + 2);
     double maxval = DBL_MIN;
 
     for (int i = 0; i < size; i++) {
@@ -245,10 +149,10 @@ static double maxElement(Solver* s, double* m)
     return maxval;
 }
 
-void normalizePressure(Solver* s)
+void normalizePressure(Discretization* d)
 {
-    int size    = (s->grid.imax + 2) * (s->grid.jmax + 2) * (s->grid.kmax + 2);
+    int size    = (d->grid.imax + 2) * (d->grid.jmax + 2) * (d->grid.kmax + 2);
-    double* p   = s->p;
+    double* p   = d->p;
     double avgP = 0.0;
 
     for (int i = 0; i < size; i++) {
@@ -261,16 +165,16 @@ void normalizePressure(Solver* s)
     }
 }
 
-void computeTimestep(Solver* s)
+void computeTimestep(Discretization* d)
 {
-    double dt = s->dtBound;
+    double dt = d->dtBound;
-    double dx = s->grid.dx;
+    double dx = d->grid.dx;
-    double dy = s->grid.dy;
+    double dy = d->grid.dy;
-    double dz = s->grid.dz;
+    double dz = d->grid.dz;
 
-    double umax = maxElement(s, s->u);
+    double umax = maxElement(d, d->u);
-    double vmax = maxElement(s, s->v);
+    double vmax = maxElement(d, d->v);
-    double wmax = maxElement(s, s->w);
+    double wmax = maxElement(d, d->w);
 
     if (umax > 0) {
         dt = (dt > dx / umax) ? dx / umax : dt;
@@ -282,20 +186,20 @@ void computeTimestep(Solver* s)
         dt = (dt > dz / wmax) ? dz / wmax : dt;
     }
 
-    s->dt = dt * s->tau;
+    d->dt = dt * d->tau;
 }
 
-void setBoundaryConditions(Solver* s)
+void setBoundaryConditions(Discretization* d)
 {
-    int imax = s->grid.imax;
+    int imax = d->grid.imax;
-    int jmax = s->grid.jmax;
+    int jmax = d->grid.jmax;
-    int kmax = s->grid.kmax;
+    int kmax = d->grid.kmax;
 
-    double* u = s->u;
+    double* u = d->u;
-    double* v = s->v;
+    double* v = d->v;
-    double* w = s->w;
+    double* w = d->w;
 
-    switch (s->bcTop) {
+    switch (d->bcTop) {
     case NOSLIP:
         for (int k = 1; k < kmax + 1; k++) {
             for (int i = 1; i < imax + 1; i++) {
@@ -327,7 +231,7 @@ void setBoundaryConditions(Solver* s)
         break;
     }
 
-    switch (s->bcBottom) {
+    switch (d->bcBottom) {
     case NOSLIP:
         for (int k = 1; k < kmax + 1; k++) {
             for (int i = 1; i < imax + 1; i++) {
@@ -359,7 +263,7 @@ void setBoundaryConditions(Solver* s)
         break;
     }
 
-    switch (s->bcLeft) {
+    switch (d->bcLeft) {
     case NOSLIP:
         for (int k = 1; k < kmax + 1; k++) {
             for (int j = 1; j < jmax + 1; j++) {
@@ -391,7 +295,7 @@ void setBoundaryConditions(Solver* s)
         break;
     }
 
-    switch (s->bcRight) {
+    switch (d->bcRight) {
     case NOSLIP:
         for (int k = 1; k < kmax + 1; k++) {
             for (int j = 1; j < jmax + 1; j++) {
@@ -423,7 +327,7 @@ void setBoundaryConditions(Solver* s)
         break;
     }
 
-    switch (s->bcFront) {
+    switch (d->bcFront) {
     case NOSLIP:
         for (int j = 1; j < jmax + 1; j++) {
             for (int i = 1; i < imax + 1; i++) {
@@ -455,7 +359,7 @@ void setBoundaryConditions(Solver* s)
         break;
     }
 
-    switch (s->bcBack) {
+    switch (d->bcBack) {
     case NOSLIP:
         for (int j = 1; j < jmax + 1; j++) {
             for (int i = 1; i < imax + 1; i++) {
@@ -488,23 +392,23 @@ void setBoundaryConditions(Solver* s)
     }
 }
 
-void setSpecialBoundaryCondition(Solver* s)
+void setSpecialBoundaryCondition(Discretization* d)
 {
-    int imax = s->grid.imax;
+    int imax = d->grid.imax;
-    int jmax = s->grid.jmax;
+    int jmax = d->grid.jmax;
-    int kmax = s->grid.kmax;
+    int kmax = d->grid.kmax;
 
-    double mDy = s->grid.dy;
+    double mDy = d->grid.dy;
-    double* u  = s->u;
+    double* u  = d->u;
 
-    if (strcmp(s->problem, "dcavity") == 0) {
+    if (strcmp(d->problem, "dcavity") == 0) {
         for (int k = 1; k < kmax; k++) {
             for (int i = 1; i < imax; i++) {
                 U(i, jmax + 1, k) = 2.0 - U(i, jmax, k);
             }
         }
-    } else if (strcmp(s->problem, "canal") == 0) {
+    } else if (strcmp(d->problem, "canal") == 0) {
-        double ylength = s->grid.ylength;
+        double ylength = d->grid.ylength;
         double y;
 
         for (int k = 1; k < kmax + 1; k++) {
@@ -516,29 +420,29 @@ void setSpecialBoundaryCondition(Solver* s)
     }
 }
 
-void computeFG(Solver* s)
+void computeFG(Discretization* d)
 {
-    int imax = s->grid.imax;
+    int imax = d->grid.imax;
-    int jmax = s->grid.jmax;
+    int jmax = d->grid.jmax;
-    int kmax = s->grid.kmax;
+    int kmax = d->grid.kmax;
 
-    double* u = s->u;
+    double* u = d->u;
-    double* v = s->v;
+    double* v = d->v;
-    double* w = s->w;
+    double* w = d->w;
-    double* f = s->f;
+    double* f = d->f;
-    double* g = s->g;
+    double* g = d->g;
-    double* h = s->h;
+    double* h = d->h;
 
-    double gx = s->gx;
+    double gx = d->gx;
-    double gy = s->gy;
+    double gy = d->gy;
-    double gz = s->gz;
+    double gz = d->gz;
-    double dt = s->dt;
+    double dt = d->dt;
 
-    double gamma     = s->gamma;
+    double gamma     = d->gamma;
-    double inverseRe = 1.0 / s->re;
+    double inverseRe = 1.0 / d->re;
-    double inverseDx = 1.0 / s->grid.dx;
+    double inverseDx = 1.0 / d->grid.dx;
-    double inverseDy = 1.0 / s->grid.dy;
+    double inverseDy = 1.0 / d->grid.dy;
-    double inverseDz = 1.0 / s->grid.dz;
+    double inverseDz = 1.0 / d->grid.dz;
     double du2dx, dv2dy, dw2dz;
     double duvdx, duwdx, duvdy, dvwdy, duwdz, dvwdz;
     double du2dx2, du2dy2, du2dz2;
@@ -705,23 +609,23 @@ void computeFG(Solver* s)
     }
 }
 
-void adaptUV(Solver* s)
+void adaptUV(Discretization* d)
 {
-    int imax = s->grid.imax;
+    int imax = d->grid.imax;
-    int jmax = s->grid.jmax;
+    int jmax = d->grid.jmax;
-    int kmax = s->grid.kmax;
+    int kmax = d->grid.kmax;
 
-    double* p = s->p;
+    double* p = d->p;
-    double* u = s->u;
+    double* u = d->u;
-    double* v = s->v;
+    double* v = d->v;
-    double* w = s->w;
+    double* w = d->w;
-    double* f = s->f;
+    double* f = d->f;
-    double* g = s->g;
+    double* g = d->g;
-    double* h = s->h;
+    double* h = d->h;
 
-    double factorX = s->dt / s->grid.dx;
+    double factorX = d->dt / d->grid.dx;
-    double factorY = s->dt / s->grid.dy;
+    double factorY = d->dt / d->grid.dy;
-    double factorZ = s->dt / s->grid.dz;
+    double factorZ = d->dt / d->grid.dz;
 
     for (int k = 1; k < kmax + 1; k++) {
         for (int j = 1; j < jmax + 1; j++) {

BasicSolver/3D-seq/src/discretization.h

@@ -0,0 +1,41 @@
+/*
+ * Copyright (C)  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 __DISCRETIZATION_H_
+#define __DISCRETIZATION_H_
+
+#include "grid.h"
+#include "parameter.h"
+
+enum BC { NOSLIP = 1, SLIP, OUTFLOW, PERIODIC };
+
+typedef struct {
+    /* geometry and grid information */
+    Grid grid;
+    /* arrays */
+    double *p, *rhs;
+    double *f, *g, *h;
+    double *u, *v, *w;
+    /* parameters */
+    double eps, omega;
+    double re, tau, gamma;
+    double gx, gy, gz;
+    /* time stepping */
+    double dt, te;
+    double dtBound;
+    char* problem;
+    int bcLeft, bcRight, bcBottom, bcTop, bcFront, bcBack;
+} Discretization;
+
+extern void initDiscretization(Discretization*, Parameter*);
+extern void computeRHS(Discretization*);
+extern void normalizePressure(Discretization*);
+extern void computeTimestep(Discretization*);
+extern void setBoundaryConditions(Discretization*);
+extern void setSpecialBoundaryCondition(Discretization*);
+extern void computeFG(Discretization*);
+extern void adaptUV(Discretization*);
+#endif

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

@@ -9,6 +9,7 @@
 #include <unistd.h>
 
 #include "allocate.h"
+#include "discretization.h"
 #include "parameter.h"
 #include "progress.h"
 #include "solver.h"
@@ -17,7 +18,8 @@
 
 #define G(v, i, j, k) v[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)]
 
-static void createBulkArrays(Solver* s, double* pg, double* ug, double* vg, double* wg)
+static void createBulkArrays(
+    Discretization* s, double* pg, double* ug, double* vg, double* wg)
 {
     int imax = s->grid.imax;
     int jmax = s->grid.jmax;
@@ -67,6 +69,7 @@ int main(int argc, char** argv)
 {
     double timeStart, timeStop;
     Parameter p;
+    Discretization d;
     Solver s;
     initParameter(&p);
 
@@ -77,51 +80,53 @@ int main(int argc, char** argv)
 
     readParameter(&p, argv[1]);
     printParameter(&p);
-    initSolver(&s, &p);
+    initDiscretization(&d, &p);
+    initSolver(&s, &d, &p);
 #ifndef VERBOSE
-    initProgress(s.te);
+    initProgress(d.te);
 #endif
 
-    double tau = s.tau;
+    double tau = d.tau;
-    double te  = s.te;
+    double te  = d.te;
     double t   = 0.0;
+    int nt     = 0;
 
     timeStart = getTimeStamp();
     while (t <= te) {
-        if (tau > 0.0) computeTimestep(&s);
+        if (tau > 0.0) computeTimestep(&d);
-        setBoundaryConditions(&s);
+        setBoundaryConditions(&d);
-        setSpecialBoundaryCondition(&s);
+        setSpecialBoundaryCondition(&d);
-        computeFG(&s);
+        computeFG(&d);
-        computeRHS(&s);
+        computeRHS(&d);
-        solve(&s);
+        if (nt % 100 == 0) normalizePressure(&d);
-        adaptUV(&s);
+        solve(&s, d.p, d.rhs);
-        t += s.dt;
+        adaptUV(&d);
+        t += d.dt;
+        nt++;
 
 #ifdef VERBOSE
-        printf("TIME %f , TIMESTEP %f\n", t, s.dt);
+        printf("TIME %f , TIMESTEP %f\n", t, solver.dt);
 #else
         printProgress(t);
 #endif
     }
     timeStop = getTimeStamp();
-#ifndef VERBOSE
     stopProgress();
-#endif
     printf("Solution took %.2fs\n", timeStop - timeStart);
 
     timeStart = getTimeStamp();
     double *pg, *ug, *vg, *wg;
 
-    size_t bytesize = (size_t)(s.grid.imax * s.grid.jmax * s.grid.kmax) * sizeof(double);
+    size_t bytesize = (size_t)(d.grid.imax * d.grid.jmax * d.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);
+    createBulkArrays(&d, pg, ug, vg, wg);
-    VtkOptions opts = { .grid = s.grid };
+    VtkOptions opts = { .grid = d.grid };
-    vtkOpen(&opts, s.problem);
+    vtkOpen(&opts, d.problem);
     vtkScalar(&opts, "pressure", pg);
     vtkVector(&opts, "velocity", (VtkVector) { ug, vg, wg });
     vtkClose(&opts);

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

@@ -26,6 +26,7 @@ void initParameter(Parameter* param)
     param->re      = 100.0;
     param->gamma   = 0.9;
     param->tau     = 0.5;
+    param->levels  = 5;
 }
 
 void readParameter(Parameter* param, const char* filename)
@@ -65,6 +66,7 @@ void readParameter(Parameter* param, const char* filename)
             PARSE_INT(jmax);
             PARSE_INT(kmax);
             PARSE_INT(itermax);
+            PARSE_INT(levels);
             PARSE_REAL(eps);
             PARSE_REAL(omg);
             PARSE_REAL(re);
@@ -123,4 +125,5 @@ void printParameter(Parameter* param)
     printf("\tepsilon (stopping tolerance) : %f\n", param->eps);
     printf("\tgamma (stopping tolerance) : %f\n", param->gamma);
     printf("\tomega (SOR relaxation): %f\n", param->omg);
+    printf("\tMultiGrid levels : %d\n", param->levels);
 }

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

@@ -10,8 +10,8 @@
 typedef struct {
     int imax, jmax, kmax;
     double xlength, ylength, zlength;
-    int itermax;
+    int itermax, levels;
-    double eps, omg;
+    double eps, omg, rho;
     double re, tau, gamma;
     double te, dt;
     double gx, gy, gz;

BasicSolver/3D-seq/src/solver-mg.c

@@ -0,0 +1,250 @@
+/*
+ * Copyright (C) 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 "allocate.h"
+#include "solver.h"
+#include "util.h"
+
+#define FINEST_LEVEL   0
+#define COARSEST_LEVEL (s->levels - 1)
+#define S(i, j, k)     s[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)]
+#define E(i, j, k)     e[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)]
+#define R(i, j, k)     r[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)]
+#define OLD(i, j, k)   old[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)]
+
+static void restrictMG(Solver* s, int level, int imax, int jmax, int kmax)
+{
+    double* r   = s->r[level + 1];
+    double* old = s->r[level];
+
+    for (int k = 1; k < kmax + 1; k++) {
+        for (int j = 1; j < jmax + 1; j++) {
+            for (int i = 1; i < imax + 1; ++i) {
+                R(i, j, k) = (OLD(2 * i - 1, 2 * j - 1, 2 * k) +
+                                 OLD(2 * i, 2 * j - 1, 2 * k) * 2 +
+                                 OLD(2 * i + 1, 2 * j - 1, 2 * k) +
+                                 OLD(2 * i - 1, 2 * j, 2 * k) * 2 +
+                                 OLD(2 * i, 2 * j, 2 * k) * 8 +
+                                 OLD(2 * i + 1, 2 * j, 2 * k) * 2 +
+                                 OLD(2 * i - 1, 2 * j + 1, 2 * k) +
+                                 OLD(2 * i, 2 * j + 1, 2 * k) * 2 +
+                                 OLD(2 * i + 1, 2 * j + 1, 2 * k) +
+
+                                 OLD(2 * i - 1, 2 * j - 1, 2 * k - 1) +
+                                 OLD(2 * i, 2 * j - 1, 2 * k - 1) * 2 +
+                                 OLD(2 * i + 1, 2 * j - 1, 2 * k - 1) +
+                                 OLD(2 * i - 1, 2 * j, 2 * k - 1) * 2 +
+                                 OLD(2 * i, 2 * j, 2 * k - 1) * 4 +
+                                 OLD(2 * i + 1, 2 * j, 2 * k - 1) * 2 +
+                                 OLD(2 * i - 1, 2 * j + 1, 2 * k - 1) +
+                                 OLD(2 * i, 2 * j + 1, 2 * k - 1) * 2 +
+                                 OLD(2 * i + 1, 2 * j + 1, 2 * k - 1) +
+
+                                 OLD(2 * i - 1, 2 * j - 1, 2 * k + 1) +
+                                 OLD(2 * i, 2 * j - 1, 2 * k + 1) * 2 +
+                                 OLD(2 * i + 1, 2 * j - 1, 2 * k + 1) +
+                                 OLD(2 * i - 1, 2 * j, 2 * k + 1) * 2 +
+                                 OLD(2 * i, 2 * j, 2 * k + 1) * 4 +
+                                 OLD(2 * i + 1, 2 * j, 2 * k + 1) * 2 +
+                                 OLD(2 * i - 1, 2 * j + 1, 2 * k + 1) +
+                                 OLD(2 * i, 2 * j + 1, 2 * k + 1) * 2 +
+                                 OLD(2 * i + 1, 2 * j + 1, 2 * k + 1)) /
+                             64.0;
+            }
+        }
+    }
+}
+
+static void prolongate(Solver* s, int level, int imax, int jmax, int kmax)
+{
+    double* old = s->r[level + 1];
+    double* e   = s->r[level];
+
+    for (int k = 2; k < kmax + 1; k += 2) {
+        for (int j = 2; j < jmax + 1; j += 2) {
+            for (int i = 2; i < imax + 1; i += 2) {
+                E(i, j, k) = OLD(i / 2, j / 2, k / 2);
+            }
+        }
+    }
+}
+
+static void correct(Solver* s, double* p, int level, int imax, int jmax, int kmax)
+{
+    double* e = s->e[level];
+
+    for (int k = 1; k < kmax + 1; ++k) {
+        for (int j = 1; j < jmax + 1; ++j) {
+            for (int i = 1; i < imax + 1; ++i) {
+                P(i, j, k) += E(i, j, k);
+            }
+        }
+    }
+}
+
+static void setBoundaryCondition(double* p, int imax, int jmax, int kmax)
+{
+    for (int j = 1; j < jmax + 1; j++) {
+        for (int i = 1; i < imax + 1; i++) {
+            P(i, j, 0)        = P(i, j, 1);
+            P(i, j, kmax + 1) = P(i, j, kmax);
+        }
+    }
+
+    for (int k = 1; k < kmax + 1; k++) {
+        for (int i = 1; i < imax + 1; i++) {
+            P(i, 0, k)        = P(i, 1, k);
+            P(i, jmax + 1, k) = P(i, jmax, k);
+        }
+    }
+
+    for (int k = 1; k < kmax + 1; k++) {
+        for (int j = 1; j < jmax + 1; j++) {
+            P(0, j, k)        = P(1, j, k);
+            P(imax + 1, j, k) = P(imax, j, k);
+        }
+    }
+}
+
+static double smooth(
+    Solver* s, double* p, double* rhs, int level, int imax, int jmax, int kmax)
+{
+    double eps    = s->eps;
+    int itermax   = s->itermax;
+    double dx2    = s->grid->dx * s->grid->dx;
+    double dy2    = s->grid->dy * s->grid->dy;
+    double dz2    = s->grid->dz * s->grid->dz;
+    double idx2   = 1.0 / dx2;
+    double idy2   = 1.0 / dy2;
+    double idz2   = 1.0 / dz2;
+    double factor = s->omega * 0.5 * (dx2 * dy2 * dz2) /
+                    (dy2 * dz2 + dx2 * dz2 + dx2 * dy2);
+    double* r    = s->r[level];
+    double epssq = eps * eps;
+    int it       = 0;
+    int pass, ksw, jsw, isw;
+    double res = 1.0;
+
+    ksw = 1;
+
+    for (pass = 0; pass < 2; pass++) {
+        jsw = ksw;
+
+        for (int k = 1; k < kmax + 1; k++) {
+            isw = jsw;
+            for (int j = 1; j < jmax + 1; j++) {
+                for (int i = isw; i < imax + 1; i += 2) {
+
+                    R(i, j, k) =
+                        RHS(i, j, k) -
+                        ((P(i + 1, j, k) - 2.0 * P(i, j, k) + P(i - 1, j, k)) * idx2 +
+                            (P(i, j + 1, k) - 2.0 * P(i, j, k) + P(i, j - 1, k)) * idy2 +
+                            (P(i, j, k + 1) - 2.0 * P(i, j, k) + P(i, j, k - 1)) * idz2);
+
+                    P(i, j, k) -= (factor * R(i, j, k));
+                    res += (R(i, j, k) * R(i, j, k));
+                }
+                isw = 3 - isw;
+            }
+            jsw = 3 - jsw;
+        }
+        ksw = 3 - ksw;
+    }
+
+    res = res / (double)(imax * jmax * kmax);
+
+    return res;
+}
+
+static double multiGrid(
+    Solver* s, double* p, double* rhs, int level, int imax, int jmax, int kmax)
+{
+    double res = 0.0;
+
+    // coarsest level TODO: Use direct solver?
+    if (level == COARSEST_LEVEL) {
+        for (int i = 0; i < 5; i++) {
+            smooth(s, p, rhs, level, imax, jmax, kmax);
+        }
+        return res;
+    }
+
+    // pre-smoothing TODO: Make smoothing steps configurable?
+    for (int i = 0; i < 5; i++) {
+        smooth(s, p, rhs, level, imax, jmax, kmax);
+        if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax, kmax);
+    }
+
+    // restrict
+    restrictMG(s, level, imax, jmax, kmax);
+
+    // MGSolver on residual and error.
+    // TODO: What if there is a rest?
+    multiGrid(s,
+        s->e[level + 1],
+        s->r[level + 1],
+        level + 1,
+        imax / 2,
+        jmax / 2,
+        kmax / 2);
+
+    // prolongate
+    prolongate(s, level, imax, jmax, kmax);
+
+    // correct p on finer level using residual
+    correct(s, p, level, imax, jmax, kmax);
+    if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax, kmax);
+
+    // post-smoothing
+    for (int i = 0; i < 5; i++) {
+        res = smooth(s, p, rhs, level, imax, jmax, kmax);
+        if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax, kmax);
+    }
+
+    return res;
+}
+
+void initSolver(Solver* s, Discretization* d, Parameter* p)
+{
+    s->eps     = p->eps;
+    s->omega   = p->omg;
+    s->itermax = p->itermax;
+    s->levels  = p->levels;
+    s->grid    = &d->grid;
+
+    int imax   = s->grid->imax;
+    int jmax   = s->grid->jmax;
+    int kmax   = s->grid->kmax;
+    int levels = s->levels;
+    printf("Using Multigrid solver with %d levels\n", levels);
+
+    s->r = malloc(levels * sizeof(double*));
+    s->e = malloc(levels * sizeof(double*));
+
+    size_t size = (imax + 2) * (jmax + 2) * (kmax + 2);
+
+    for (int j = 0; j < levels; j++) {
+        s->r[j] = allocate(64, size * sizeof(double));
+        s->e[j] = allocate(64, size * sizeof(double));
+
+        for (size_t i = 0; i < size; i++) {
+            s->r[j][i] = 0.0;
+            s->e[j][i] = 0.0;
+        }
+    }
+}
+
+void solve(Solver* s, double* p, double* rhs)
+{
+    double res = multiGrid(s, p, rhs, 0, s->grid->imax, s->grid->jmax, s->grid->kmax);
+
+#ifdef VERBOSE
+    printf("Residuum: %.6f\n", res);
+#endif
+}

BasicSolver/3D-seq/src/solver-sor.c

@@ -0,0 +1,99 @@
+/*
+ * Copyright (C) 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 "solver.h"
+#include "util.h"
+
+void initSolver(Solver* s, Discretization* d, Parameter* p)
+{
+    s->grid    = &d->grid;
+    s->itermax = p->itermax;
+    s->eps     = p->eps;
+    s->omega   = p->omg;
+}
+
+void solve(Solver* s, double* p, double* rhs)
+{
+    int imax      = s->grid->imax;
+    int jmax      = s->grid->jmax;
+    int kmax      = s->grid->kmax;
+    double eps    = s->eps;
+    int itermax   = s->itermax;
+    double dx2    = s->grid->dx * s->grid->dx;
+    double dy2    = s->grid->dy * s->grid->dy;
+    double dz2    = s->grid->dz * s->grid->dz;
+    double idx2   = 1.0 / dx2;
+    double idy2   = 1.0 / dy2;
+    double idz2   = 1.0 / dz2;
+    double factor = s->omega * 0.5 * (dx2 * dy2 * dz2) /
+                    (dy2 * dz2 + dx2 * dz2 + dx2 * dy2);
+    double epssq = eps * eps;
+    int it       = 0;
+    double res   = 1.0;
+    int pass, ksw, jsw, isw;
+
+    while ((res >= epssq) && (it < itermax)) {
+        res = 0.0;
+        ksw = 1;
+
+        for (pass = 0; pass < 2; pass++) {
+            jsw = ksw;
+
+            for (int k = 1; k < kmax + 1; k++) {
+                isw = jsw;
+                for (int j = 1; j < jmax + 1; j++) {
+                    for (int i = isw; i < imax + 1; i += 2) {
+
+                        double r =
+                            RHS(i, j, k) -
+                            ((P(i + 1, j, k) - 2.0 * P(i, j, k) + P(i - 1, j, k)) * idx2 +
+                                (P(i, j + 1, k) - 2.0 * P(i, j, k) + P(i, j - 1, k)) *
+                                    idy2 +
+                                (P(i, j, k + 1) - 2.0 * P(i, j, k) + P(i, j, k - 1)) *
+                                    idz2);
+
+                        P(i, j, k) -= (factor * r);
+                        res += (r * r);
+                    }
+                    isw = 3 - isw;
+                }
+                jsw = 3 - jsw;
+            }
+            ksw = 3 - ksw;
+        }
+
+        for (int j = 1; j < jmax + 1; j++) {
+            for (int i = 1; i < imax + 1; i++) {
+                P(i, j, 0)        = P(i, j, 1);
+                P(i, j, kmax + 1) = P(i, j, kmax);
+            }
+        }
+
+        for (int k = 1; k < kmax + 1; k++) {
+            for (int i = 1; i < imax + 1; i++) {
+                P(i, 0, k)        = P(i, 1, k);
+                P(i, jmax + 1, k) = P(i, jmax, k);
+            }
+        }
+
+        for (int k = 1; k < kmax + 1; k++) {
+            for (int j = 1; j < jmax + 1; j++) {
+                P(0, j, k)        = P(1, j, k);
+                P(imax + 1, j, k) = P(imax, j, k);
+            }
+        }
+
+        res = res / (double)(imax * jmax * kmax);
+#ifdef DEBUG
+        printf("%d Residuum: %e\n", it, res);
+#endif
+        it++;
+    }
+
+#ifdef VERBOSE
+    printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
+#endif
+}

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

@@ -6,38 +6,21 @@
  */
 #ifndef __SOLVER_H_
 #define __SOLVER_H_
-
+#include "discretization.h"
 #include "grid.h"
 #include "parameter.h"
 
-enum BC { NOSLIP = 1, SLIP, OUTFLOW, PERIODIC };
-
 typedef struct {
     /* geometry and grid information */
-    Grid grid;
+    Grid* grid;
-    /* arrays */
-    double *p, *rhs;
-    double *f, *g, *h;
-    double *u, *v, *w;
     /* parameters */
-    double eps, omega;
+    double eps, omega, rho;
-    double re, tau, gamma;
-    double gx, gy, gz;
-    /* time stepping */
     int itermax;
-    double dt, te;
+    int levels;
-    double dtBound;
+    double **r, **e;
-    char* problem;
-    int bcLeft, bcRight, bcBottom, bcTop, bcFront, bcBack;
 } Solver;
 
-extern void initSolver(Solver*, Parameter*);
+extern void initSolver(Solver*, Discretization*, Parameter*);
-extern void computeRHS(Solver*);
+extern void solve(Solver*, double*, double*);
-extern void solve(Solver*);
+
-extern void normalizePressure(Solver*);
-extern void computeTimestep(Solver*);
-extern void setBoundaryConditions(Solver*);
-extern void setSpecialBoundaryCondition(Solver*);
-extern void computeFG(Solver*);
-extern void adaptUV(Solver*);
 #endif

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

@@ -19,4 +19,13 @@
 #define ABS(a) ((a) >= 0 ? (a) : -(a))
 #endif
 
+#define P(i, j, k)   p[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)]
+#define F(i, j, k)   f[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)]
+#define G(i, j, k)   g[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)]
+#define H(i, j, k)   h[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)]
+#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 RHS(i, j, k) rhs[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)]
+
 #endif // __UTIL_H_

EnhancedSolver/2D-seq/Makefile

@@ -18,9 +18,10 @@ include $(MAKE_DIR)/include_$(TAG).mk
 INCLUDES  += -I$(SRC_DIR) -I$(BUILD_DIR)
 
 VPATH     = $(SRC_DIR)
-SRC       = $(wildcard $(SRC_DIR)/*.c)
+SRC       = $(filter-out $(wildcard $(SRC_DIR)/*-*.c),$(wildcard $(SRC_DIR)/*.c))
 ASM       = $(patsubst $(SRC_DIR)/%.c, $(BUILD_DIR)/%.s, $(SRC))
 OBJ       = $(patsubst $(SRC_DIR)/%.c, $(BUILD_DIR)/%.o, $(SRC))
+OBJ      += $(BUILD_DIR)/solver-$(SOLVER).o
 SOURCES   = $(SRC) $(wildcard $(SRC_DIR)/*.h)
 CPPFLAGS := $(CPPFLAGS) $(DEFINES) $(OPTIONS) $(INCLUDES)
 
@@ -49,6 +50,8 @@ distclean: clean
 	@rm -f $(TARGET)
 	@rm -f *.dat
 	@rm -f *.png
+	@rm -f ./vis_files/*.dat
+	@rm -f ./vis_files/*.gif
 
 info:
 	$(info $(CFLAGS))

EnhancedSolver/2D-seq/canal.par

@@ -26,13 +26,13 @@ p_init    1.0      # initial value for pressure
 
 xlength    30.0     # domain size in x-direction
 ylength    4.0	   # domain size in y-direction
-imax       200      # number of interior cells in x-direction
+imax       256      # number of interior cells in x-direction
-jmax       50	   # number of interior cells in y-direction
+jmax       64	   # number of interior cells in y-direction
 
 # Time Data:
 # ---------
 
-te      60.0   # final time
+te      80.0   # final time
 dt      0.02    # time stepsize
 tau     0.5     # safety factor for time stepsize control (<0 constant delt)
 
@@ -50,9 +50,9 @@ levels   5         # Multigrid levels
 # -----------------------
 
 numberOfParticles   60
-startTime           100
+startTime           5.0
-injectTimePeriod    2.0
+injectTimePeriod    4.0
-writeTimePeriod     0.5
+writeTimePeriod     1.0
 
 x1                  1.0
 y1                  0.0

EnhancedSolver/2D-seq/config.mk

@@ -1,12 +1,12 @@
 # Supported: GCC, CLANG, ICC
 TAG ?= CLANG
 ENABLE_OPENMP ?= false
+# Supported: sor, mg
+SOLVER ?= mg
+# Run in debug settings
+DEBUG ?= false
 
 #Feature options
 OPTIONS +=  -DARRAY_ALIGNMENT=64
-# OPTIONS +=  -DVERBOSE
+#OPTIONS +=  -DVERBOSE
 #OPTIONS +=  -DDEBUG
-#OPTIONS +=  -DBOUNDCHECK
-#OPTIONS +=  -DVERBOSE_AFFINITY
-#OPTIONS +=  -DVERBOSE_DATASIZE
-#OPTIONS +=  -DVERBOSE_TIMER

EnhancedSolver/2D-seq/dcavity.par

@@ -26,8 +26,8 @@ p_init    0.0		# initial value for pressure
 
 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
+imax       128		# number of interior cells in x-direction
-jmax       100		# number of interior cells in y-direction
+jmax       128		# number of interior cells in y-direction
 
 # Time Data:
 # ---------
@@ -50,7 +50,7 @@ levels   5         # Multigrid levels
 # -----------------------
 
 numberOfParticles   200
-startTime           100
+startTime           2.0
 injectTimePeriod    0.5
 writeTimePeriod     0.2
 

EnhancedSolver/2D-seq/include_CLANG.mk

@@ -2,16 +2,18 @@ CC   = clang
 GCC  = cc
 LINKER = $(CC)
 
-ifeq ($(ENABLE_OPENMP),true)
+ifeq ($(strip $(ENABLE_OPENMP)),true)
 OPENMP   = -fopenmp
 #OPENMP   = -Xpreprocessor -fopenmp #required on Macos with homebrew libomp
 LIBS     = # -lomp
 endif
+ifeq ($(strip $(DEBUG)),true)
+CFLAGS   = -O0 -g -std=c17
+else
+CFLAGS   = -O3 -std=c17 $(OPENMP)
+endif
 
 VERSION  = --version
-# CFLAGS   = -O3 -std=c17 $(OPENMP)
-CFLAGS   = -Ofast -std=c17
-#CFLAGS   = -Ofast -fnt-store=aggressive  -std=c99 $(OPENMP) #AMD CLANG
 LFLAGS   = $(OPENMP) -lm
-DEFINES  = -D_GNU_SOURCE# -DDEBUG
+DEFINES  = -D_GNU_SOURCE
 INCLUDES =

EnhancedSolver/2D-seq/src/affinity.c

@@ -1,61 +0,0 @@
-/*
- * Copyright (C)  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__*/

EnhancedSolver/2D-seq/src/affinity.h

@@ -1,14 +0,0 @@
-/*
- * Copyright (C)  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*/

EnhancedSolver/2D-seq/src/allocate.c

@@ -5,10 +5,13 @@
  * license that can be found in the LICENSE file.
  */
 #include <errno.h>
+#include <stddef.h>
 #include <stdio.h>
 #include <stdlib.h>
 
-void* allocate(int alignment, size_t bytesize)
+#include "allocate.h"
+
+void* allocate(size_t alignment, size_t bytesize)
 {
     int errorCode;
     void* ptr;

EnhancedSolver/2D-seq/src/allocate.h

@@ -8,6 +8,6 @@
 #define __ALLOCATE_H_
 #include <stdlib.h>
 
-extern void* allocate(int alignment, size_t bytesize);
+extern void* allocate(size_t alignment, size_t bytesize);
 
 #endif

EnhancedSolver/2D-seq/src/discretization.c

@@ -0,0 +1,649 @@
+/*
+ * Copyright (C) NHR@FAU, University Erlangen-Nuremberg.
+ * All rights reserved. This file is part of nusif-solver.
+ * Use of this source code is governed by a MIT style
+ * license that can be found in the LICENSE file.
+ */
+#include <float.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "allocate.h"
+#include "discretization.h"
+#include "grid.h"
+#include "parameter.h"
+#include "util.h"
+
+#define S(i, j) s[(j) * (imax + 2) + (i)]
+
+static double distance(double i, double j, double iCenter, double jCenter)
+{
+    return sqrt(pow(iCenter - i, 2) + pow(jCenter - j, 2) * 1.0);
+}
+
+void print(Discretization* d, double* grid)
+{
+    int imax = d->grid.imax;
+    int jmax = d->grid.jmax;
+
+    for (int j = 0; j < jmax + 2; j++) {
+        printf("%02d: ", j);
+        for (int i = 0; i < imax + 2; i++) {
+            printf("%3.2f  ", grid[j * (imax + 2) + i]);
+        }
+        printf("\n");
+    }
+    fflush(stdout);
+}
+
+void printGrid(Discretization* d, int* grid)
+{
+    int imax = d->grid.imax;
+    int jmax = d->grid.jmax;
+
+    for (int j = 0; j < jmax + 2; j++) {
+        printf("%02d: ", j);
+        for (int i = 0; i < imax + 2; i++) {
+            printf("%2d  ", grid[j * (imax + 2) + i]);
+        }
+        printf("\n");
+    }
+    fflush(stdout);
+}
+
+static void printConfig(Discretization* d)
+{
+    printf("Parameters for #%s#\n", d->problem);
+    printf("Boundary conditions Left:%d Right:%d Bottom:%d Top:%d\n",
+        d->bcLeft,
+        d->bcRight,
+        d->bcBottom,
+        d->bcTop);
+    printf("\tReynolds number: %.2f\n", d->re);
+    printf("\tGx Gy: %.2f %.2f\n", d->gx, d->gy);
+    printf("Geometry data:\n");
+    printf("\tDomain box size (x, y): %.2f, %.2f\n", d->grid.xlength, d->grid.ylength);
+    printf("\tCells (x, y): %d, %d\n", d->grid.imax, d->grid.jmax);
+    printf("Timestep parameters:\n");
+    printf("\tDefault stepsize: %.2f, Final time %.2f\n", d->dt, d->te);
+    printf("\tdt bound: %.6f\n", d->dtBound);
+    printf("\tTau factor: %.2f\n", d->tau);
+    printf("Iterative d parameters:\n");
+    printf("\tgamma factor: %f\n", d->gamma);
+}
+
+void initDiscretization(Discretization* d, Parameter* p)
+{
+    d->problem      = p->name;
+    d->bcLeft       = p->bcLeft;
+    d->bcRight      = p->bcRight;
+    d->bcBottom     = p->bcBottom;
+    d->bcTop        = p->bcTop;
+    d->grid.imax    = p->imax;
+    d->grid.jmax    = p->jmax;
+    d->grid.xlength = p->xlength;
+    d->grid.ylength = p->ylength;
+    d->grid.dx      = p->xlength / p->imax;
+    d->grid.dy      = p->ylength / p->jmax;
+    d->re           = p->re;
+    d->gx           = p->gx;
+    d->gy           = p->gy;
+    d->dt           = p->dt;
+    d->te           = p->te;
+    d->tau          = p->tau;
+    d->gamma        = p->gamma;
+
+    int imax    = d->grid.imax;
+    int jmax    = d->grid.jmax;
+    size_t size = (imax + 2) * (jmax + 2) * sizeof(double);
+    d->u        = allocate(64, size);
+    d->v        = allocate(64, size);
+    d->grid.s   = allocate(64, size);
+    d->p        = allocate(64, size);
+    d->rhs      = allocate(64, size);
+    d->f        = allocate(64, size);
+    d->g        = allocate(64, size);
+
+    for (int i = 0; i < (imax + 2) * (jmax + 2); i++) {
+
+        d->u[i]      = p->u_init;
+        d->v[i]      = p->v_init;
+        d->p[i]      = p->p_init;
+        d->rhs[i]    = 0.0;
+        d->f[i]      = 0.0;
+        d->g[i]      = 0.0;
+        d->grid.s[i] = FLUID;
+    }
+
+    double dx        = d->grid.dx;
+    double dy        = d->grid.dy;
+    double invSqrSum = 1.0 / (dx * dx) + 1.0 / (dy * dy);
+    d->dtBound       = 0.5 * d->re * 1.0 / invSqrSum;
+
+    double xCenter = 0, yCenter = 0, radius = 0;
+    double x1 = 0, x2 = 0, y1 = 0, y2 = 0;
+
+    int* s = d->grid.s;
+
+    switch (p->shape) {
+    case NOSHAPE:
+        break;
+    case RECT:
+        x1 = p->xCenter - p->xRectLength / 2;
+        x2 = p->xCenter + p->xRectLength / 2;
+        y1 = p->yCenter - p->yRectLength / 2;
+        y2 = p->yCenter + p->yRectLength / 2;
+
+        for (int j = 1; j < jmax + 1; j++) {
+            for (int i = 1; i < imax + 1; i++) {
+                if ((x1 <= (i * dx)) && ((i * dx) <= x2) && (y1 <= (j * dy)) &&
+                    ((j * dy) <= y2)) {
+                    S(i, j) = OBSTACLE;
+                }
+            }
+        }
+        break;
+    case CIRCLE:
+        xCenter = p->xCenter;
+        yCenter = p->yCenter;
+        radius  = p->circleRadius;
+
+        for (int j = 1; j < jmax + 1; j++) {
+            for (int i = 1; i < imax + 1; i++) {
+                if (distance((i * dx), (j * dy), xCenter, yCenter) <= radius) {
+                    S(i, j) = OBSTACLE;
+                }
+            }
+        }
+        break;
+    }
+
+    if (p->shape != NOSHAPE) {
+        for (int j = 1; j < jmax + 1; j++) {
+            for (int i = 1; i < imax + 1; i++) {
+
+                if (S(i, j - 1) == FLUID && S(i, j + 1) == OBSTACLE &&
+                    S(i, j) == OBSTACLE)
+                    S(i, j) = BOTTOM; // TOP
+                if (S(i - 1, j) == FLUID && S(i + 1, j) == OBSTACLE &&
+                    S(i, j) == OBSTACLE)
+                    S(i, j) = LEFT; // LEFT
+                if (S(i + 1, j) == FLUID && S(i - 1, j) == OBSTACLE &&
+                    S(i, j) == OBSTACLE)
+                    S(i, j) = RIGHT; // RIGHT
+                if (S(i, j + 1) == FLUID && S(i, j - 1) == OBSTACLE &&
+                    S(i, j) == OBSTACLE)
+                    S(i, j) = TOP; // BOTTOM
+                if (S(i - 1, j - 1) == FLUID && S(i, j - 1) == FLUID &&
+                    S(i - 1, j) == FLUID && S(i + 1, j + 1) == OBSTACLE &&
+                    (S(i, j) == OBSTACLE || S(i, j) == LEFT || S(i, j) == BOTTOM))
+                    S(i, j) = BOTTOMLEFT; // TOPLEFT
+                if (S(i + 1, j - 1) == FLUID && S(i, j - 1) == FLUID &&
+                    S(i + 1, j) == FLUID && S(i - 1, j + 1) == OBSTACLE &&
+                    (S(i, j) == OBSTACLE || S(i, j) == RIGHT || S(i, j) == BOTTOM))
+                    S(i, j) = BOTTOMRIGHT; // TOPRIGHT
+                if (S(i - 1, j + 1) == FLUID && S(i - 1, j) == FLUID &&
+                    S(i, j + 1) == FLUID && S(i + 1, j - 1) == OBSTACLE &&
+                    (S(i, j) == OBSTACLE || S(i, j) == LEFT || S(i, j) == TOP))
+                    S(i, j) = TOPLEFT; // BOTTOMLEFT
+                if (S(i + 1, j + 1) == FLUID && S(i + 1, j) == FLUID &&
+                    S(i, j + 1) == FLUID && S(i - 1, j - 1) == OBSTACLE &&
+                    (S(i, j) == OBSTACLE || S(i, j) == RIGHT || S(i, j) == TOP))
+                    S(i, j) = TOPRIGHT; // BOTTOMRIGHT
+            }
+        }
+    }
+
+#ifdef VERBOSE
+    printConfig(solver);
+#endif
+}
+
+static double maxElement(Discretization* d, double* m)
+{
+    int size      = (d->grid.imax + 2) * (d->grid.jmax + 2);
+    double maxval = DBL_MIN;
+
+    for (int i = 0; i < size; i++) {
+        maxval = MAX(maxval, fabs(m[i]));
+    }
+
+    return maxval;
+}
+
+void computeRHS(Discretization* d)
+{
+    int imax    = d->grid.imax;
+    int jmax    = d->grid.jmax;
+    double idx  = 1.0 / d->grid.dx;
+    double idy  = 1.0 / d->grid.dy;
+    double idt  = 1.0 / d->dt;
+    double* rhs = d->rhs;
+    double* f   = d->f;
+    double* g   = d->g;
+    int* s      = d->grid.s;
+
+    for (int j = 1; j < jmax + 1; j++) {
+        for (int i = 1; i < imax + 1; i++) {
+            RHS(i, j) = idt *
+                        ((F(i, j) - F(i - 1, j)) * idx + (G(i, j) - G(i, j - 1)) * idy);
+        }
+    }
+}
+
+void normalizePressure(Discretization* d)
+{
+    int size    = (d->grid.imax + 2) * (d->grid.jmax + 2);
+    double* p   = d->p;
+    double avgP = 0.0;
+
+    for (int i = 0; i < size; i++) {
+        avgP += p[i];
+    }
+    avgP /= size;
+
+    for (int i = 0; i < size; i++) {
+        p[i] = p[i] - avgP;
+    }
+}
+
+void computeTimestep(Discretization* d)
+{
+    double dt   = d->dtBound;
+    double dx   = d->grid.dx;
+    double dy   = d->grid.dy;
+    double umax = maxElement(d, d->u);
+    double vmax = maxElement(d, d->v);
+
+    if (umax > 0) {
+        dt = (dt > dx / umax) ? dx / umax : dt;
+    }
+    if (vmax > 0) {
+        dt = (dt > dy / vmax) ? dy / vmax : dt;
+    }
+
+    d->dt = dt * d->tau;
+}
+
+void setBoundaryConditions(Discretization* d)
+{
+    int imax  = d->grid.imax;
+    int jmax  = d->grid.jmax;
+    double* u = d->u;
+    double* v = d->v;
+
+    // Left boundary
+    switch (d->bcLeft) {
+    case NOSLIP:
+        for (int j = 1; j < jmax + 1; j++) {
+            U(0, j) = 0.0;
+            V(0, j) = -V(1, j);
+        }
+        break;
+    case SLIP:
+        for (int j = 1; j < jmax + 1; j++) {
+            U(0, j) = 0.0;
+            V(0, j) = V(1, j);
+        }
+        break;
+    case OUTFLOW:
+        for (int j = 1; j < jmax + 1; j++) {
+            U(0, j) = U(1, j);
+            V(0, j) = V(1, j);
+        }
+        break;
+    case PERIODIC:
+        break;
+    }
+
+    // Right boundary
+    switch (d->bcRight) {
+    case NOSLIP:
+        for (int j = 1; j < jmax + 1; j++) {
+            U(imax, j)     = 0.0;
+            V(imax + 1, j) = -V(imax, j);
+        }
+        break;
+    case SLIP:
+        for (int j = 1; j < jmax + 1; j++) {
+            U(imax, j)     = 0.0;
+            V(imax + 1, j) = V(imax, j);
+        }
+        break;
+    case OUTFLOW:
+        for (int j = 1; j < jmax + 1; j++) {
+            U(imax, j)     = U(imax - 1, j);
+            V(imax + 1, j) = V(imax, j);
+        }
+        break;
+    case PERIODIC:
+        break;
+    }
+
+    // Bottom boundary
+    switch (d->bcBottom) {
+    case NOSLIP:
+        for (int i = 1; i < imax + 1; i++) {
+            V(i, 0) = 0.0;
+            U(i, 0) = -U(i, 1);
+        }
+        break;
+    case SLIP:
+        for (int i = 1; i < imax + 1; i++) {
+            V(i, 0) = 0.0;
+            U(i, 0) = U(i, 1);
+        }
+        break;
+    case OUTFLOW:
+        for (int i = 1; i < imax + 1; i++) {
+            U(i, 0) = U(i, 1);
+            V(i, 0) = V(i, 1);
+        }
+        break;
+    case PERIODIC:
+        break;
+    }
+
+    // Top boundary
+    switch (d->bcTop) {
+    case NOSLIP:
+        for (int i = 1; i < imax + 1; i++) {
+            V(i, jmax)     = 0.0;
+            U(i, jmax + 1) = -U(i, jmax);
+        }
+        break;
+    case SLIP:
+        for (int i = 1; i < imax + 1; i++) {
+            V(i, jmax)     = 0.0;
+            U(i, jmax + 1) = U(i, jmax);
+        }
+        break;
+    case OUTFLOW:
+        for (int i = 1; i < imax + 1; i++) {
+            U(i, jmax + 1) = U(i, jmax);
+            V(i, jmax)     = V(i, jmax - 1);
+        }
+        break;
+    case PERIODIC:
+        break;
+    }
+}
+
+void setSpecialBoundaryCondition(Discretization* d)
+{
+    int imax   = d->grid.imax;
+    int jmax   = d->grid.jmax;
+    double mDy = d->grid.dy;
+    double* u  = d->u;
+    int* s     = d->grid.s;
+
+    if (strcmp(d->problem, "dcavity") == 0) {
+        for (int i = 1; i < imax; i++) {
+            U(i, jmax + 1) = 2.0 - U(i, jmax);
+        }
+    } else if (strcmp(d->problem, "canal") == 0) {
+        double ylength = d->grid.ylength;
+        double y;
+
+        for (int j = 1; j < jmax + 1; j++) {
+            y       = mDy * (j - 0.5);
+            U(0, j) = y * (ylength - y) * 4.0 / (ylength * ylength);
+        }
+    } else if (strcmp(d->problem, "backstep") == 0) {
+        for (int j = 1; j < jmax + 1; j++) {
+            if (S(0, j) == FLUID) U(0, j) = 1.0;
+        }
+    } else if (strcmp(d->problem, "karman") == 0) {
+        for (int j = 1; j < jmax + 1; j++) {
+            U(0, j) = 1.0;
+        }
+    }
+}
+
+void setObjectBoundaryCondition(Discretization* d)
+{
+    int imax  = d->grid.imax;
+    int jmax  = d->grid.jmax;
+    double* u = d->u;
+    double* v = d->v;
+    int* s    = d->grid.s;
+
+    for (int j = 1; j < jmax + 1; j++) {
+        for (int i = 1; i < imax + 1; i++) {
+            switch (S(i, j)) {
+            case TOP:
+                U(i, j)     = -U(i, j + 1);
+                U(i - 1, j) = -U(i - 1, j + 1);
+                V(i, j)     = 0.0;
+                break;
+            case BOTTOM:
+                U(i, j)     = -U(i, j - 1);
+                U(i - 1, j) = -U(i - 1, j - 1);
+                V(i, j)     = 0.0;
+                break;
+            case LEFT:
+                U(i - 1, j) = 0.0;
+                V(i, j)     = -V(i - 1, j);
+                V(i, j - 1) = -V(i - 1, j - 1);
+                break;
+            case RIGHT:
+                U(i, j)     = 0.0;
+                V(i, j)     = -V(i + 1, j);
+                V(i, j - 1) = -V(i + 1, j - 1);
+                break;
+            case TOPLEFT:
+                U(i, j)     = -U(i, j + 1);
+                U(i - 1, j) = 0.0;
+                V(i, j)     = 0.0;
+                V(i, j - 1) = -V(i - 1, j - 1);
+                break;
+            case TOPRIGHT:
+                U(i, j)     = 0.0;
+                U(i - 1, j) = -U(i - 1, j + 1);
+                V(i, j)     = 0.0;
+                V(i, j - 1) = -V(i + 1, j - 1);
+                break;
+            case BOTTOMLEFT:
+                U(i, j)     = -U(i, j - 1);
+                U(i - 1, j) = 0.0;
+                V(i, j)     = -V(i - 1, j);
+                V(i, j - 1) = 0.0;
+                break;
+            case BOTTOMRIGHT:
+                U(i, j)     = 0.0;
+                U(i - 1, j) = -U(i - 1, j - 1);
+                V(i, j)     = -V(i, j + 1);
+                V(i, j - 1) = 0.0;
+                break;
+            }
+        }
+    }
+}
+
+void computeFG(Discretization* d)
+{
+    double* u        = d->u;
+    double* v        = d->v;
+    double* f        = d->f;
+    double* g        = d->g;
+    int* s           = d->grid.s;
+    int imax         = d->grid.imax;
+    int jmax         = d->grid.jmax;
+    double gx        = d->gx;
+    double gy        = d->gy;
+    double gamma     = d->gamma;
+    double dt        = d->dt;
+    double inverseRe = 1.0 / d->re;
+    double inverseDx = 1.0 / d->grid.dx;
+    double inverseDy = 1.0 / d->grid.dy;
+    double du2dx, dv2dy, duvdx, duvdy;
+    double du2dx2, du2dy2, dv2dx2, dv2dy2;
+
+    for (int j = 1; j < jmax + 1; j++) {
+        for (int i = 1; i < imax + 1; i++) {
+            if (S(i, j) == FLUID) {
+                du2dx = inverseDx * 0.25 *
+                            ((U(i, j) + U(i + 1, j)) * (U(i, j) + U(i + 1, j)) -
+                                (U(i, j) + U(i - 1, j)) * (U(i, j) + U(i - 1, j))) +
+                        gamma * inverseDx * 0.25 *
+                            (fabs(U(i, j) + U(i + 1, j)) * (U(i, j) - U(i + 1, j)) +
+                                fabs(U(i, j) + U(i - 1, j)) * (U(i, j) - U(i - 1, j)));
+
+                duvdy = inverseDy * 0.25 *
+                            ((V(i, j) + V(i + 1, j)) * (U(i, j) + U(i, j + 1)) -
+                                (V(i, j - 1) + V(i + 1, j - 1)) *
+                                    (U(i, j) + U(i, j - 1))) +
+                        gamma * inverseDy * 0.25 *
+                            (fabs(V(i, j) + V(i + 1, j)) * (U(i, j) - U(i, j + 1)) +
+                                fabs(V(i, j - 1) + V(i + 1, j - 1)) *
+                                    (U(i, j) - U(i, j - 1)));
+
+                du2dx2 = inverseDx * inverseDx *
+                         (U(i + 1, j) - 2.0 * U(i, j) + U(i - 1, j));
+                du2dy2 = inverseDy * inverseDy *
+                         (U(i, j + 1) - 2.0 * U(i, j) + U(i, j - 1));
+                F(i, j) = U(i, j) +
+                          dt * (inverseRe * (du2dx2 + du2dy2) - du2dx - duvdy + gx);
+
+                duvdx = inverseDx * 0.25 *
+                            ((U(i, j) + U(i, j + 1)) * (V(i, j) + V(i + 1, j)) -
+                                (U(i - 1, j) + U(i - 1, j + 1)) *
+                                    (V(i, j) + V(i - 1, j))) +
+                        gamma * inverseDx * 0.25 *
+                            (fabs(U(i, j) + U(i, j + 1)) * (V(i, j) - V(i + 1, j)) +
+                                fabs(U(i - 1, j) + U(i - 1, j + 1)) *
+                                    (V(i, j) - V(i - 1, j)));
+
+                dv2dy = inverseDy * 0.25 *
+                            ((V(i, j) + V(i, j + 1)) * (V(i, j) + V(i, j + 1)) -
+                                (V(i, j) + V(i, j - 1)) * (V(i, j) + V(i, j - 1))) +
+                        gamma * inverseDy * 0.25 *
+                            (fabs(V(i, j) + V(i, j + 1)) * (V(i, j) - V(i, j + 1)) +
+                                fabs(V(i, j) + V(i, j - 1)) * (V(i, j) - V(i, j - 1)));
+
+                dv2dx2 = inverseDx * inverseDx *
+                         (V(i + 1, j) - 2.0 * V(i, j) + V(i - 1, j));
+                dv2dy2 = inverseDy * inverseDy *
+                         (V(i, j + 1) - 2.0 * V(i, j) + V(i, j - 1));
+                G(i, j) = V(i, j) +
+                          dt * (inverseRe * (dv2dx2 + dv2dy2) - duvdx - dv2dy + gy);
+            } else {
+                switch (S(i, j)) {
+                case TOP:
+                    G(i, j) = V(i, j);
+                    break;
+                case BOTTOM:
+                    G(i, j - 1) = V(i, j - 1);
+                    break;
+                case LEFT:
+                    F(i - 1, j) = U(i - 1, j);
+                    break;
+                case RIGHT:
+                    F(i, j) = U(i, j);
+                    break;
+                case TOPLEFT:
+                    F(i - 1, j) = U(i - 1, j);
+                    G(i, j)     = V(i, j);
+                    break;
+                case TOPRIGHT:
+                    F(i, j) = U(i, j);
+                    G(i, j) = V(i, j);
+                    break;
+                case BOTTOMLEFT:
+                    F(i - 1, j) = U(i - 1, j);
+                    G(i, j - 1) = V(i, j - 1);
+                    break;
+                case BOTTOMRIGHT:
+                    F(i, j)     = U(i, j);
+                    G(i, j - 1) = V(i, j - 1);
+                    break;
+                }
+            }
+        }
+    }
+
+    /* ---------------------- boundary of F --------------------------- */
+    for (int j = 1; j < jmax + 1; j++) {
+        F(0, j)    = U(0, j);
+        F(imax, j) = U(imax, j);
+    }
+
+    /* ---------------------- boundary of G --------------------------- */
+    for (int i = 1; i < imax + 1; i++) {
+        G(i, 0)    = V(i, 0);
+        G(i, jmax) = V(i, jmax);
+    }
+}
+
+void adaptUV(Discretization* d)
+{
+    int imax       = d->grid.imax;
+    int jmax       = d->grid.jmax;
+    double* p      = d->p;
+    double* u      = d->u;
+    double* v      = d->v;
+    double* f      = d->f;
+    double* g      = d->g;
+    double factorX = d->dt / d->grid.dx;
+    double factorY = d->dt / d->grid.dy;
+
+    for (int j = 1; j < jmax + 1; j++) {
+        for (int i = 1; i < imax + 1; i++) {
+            U(i, j) = F(i, j) - (P(i + 1, j) - P(i, j)) * factorX;
+            V(i, j) = G(i, j) - (P(i, j + 1) - P(i, j)) * factorY;
+        }
+    }
+}
+
+void writeResult(Discretization* d)
+{
+    int imax  = d->grid.imax;
+    int jmax  = d->grid.jmax;
+    double dx = d->grid.dx;
+    double dy = d->grid.dy;
+    double* p = d->p;
+    double* u = d->u;
+    double* v = d->v;
+    double x = 0.0, y = 0.0;
+
+    FILE* fp;
+    fp = fopen("pressure.dat", "w");
+
+    if (fp == NULL) {
+        printf("Error!\n");
+        exit(EXIT_FAILURE);
+    }
+
+    for (int j = 1; j < jmax + 1; j++) {
+        y = (double)(j - 0.5) * dy;
+        for (int i = 1; i < imax + 1; i++) {
+            x = (double)(i - 0.5) * dx;
+            fprintf(fp, "%.2f %.2f %f\n", x, y, P(i, j));
+        }
+        fprintf(fp, "\n");
+    }
+
+    fclose(fp);
+
+    fp = fopen("velocity.dat", "w");
+
+    if (fp == NULL) {
+        printf("Error!\n");
+        exit(EXIT_FAILURE);
+    }
+
+    for (int j = 1; j < jmax + 1; j++) {
+        y = dy * (j - 0.5);
+        for (int i = 1; i < imax + 1; i++) {
+            x           = dx * (i - 0.5);
+            double velU = (U(i, j) + U(i - 1, j)) / 2.0;
+            double velV = (V(i, j) + V(i, j - 1)) / 2.0;
+            double len  = sqrt((velU * velU) + (velV * velV));
+            fprintf(fp, "%.2f %.2f %f %f %f\n", x, y, velU, velV, len);
+        }
+    }
+
+    fclose(fp);
+}

EnhancedSolver/2D-seq/src/discretization.h

@@ -0,0 +1,44 @@
+/*
+ * Copyright (C) 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 __DISCRETIZATION_H_
+#define __DISCRETIZATION_H_
+#include "grid.h"
+#include "parameter.h"
+
+enum BC { NOSLIP = 1, SLIP, OUTFLOW, PERIODIC };
+
+typedef struct {
+    /* geometry and grid information */
+    Grid grid;
+    /* arrays */
+    double *p, *rhs;
+    double *f, *g;
+    double *u, *v;
+    /* parameters */
+    double rho;
+    double re, tau, gamma;
+    double gx, gy;
+    /* time stepping */
+    double dt, te;
+    double dtBound;
+    char* problem;
+    int bcLeft, bcRight, bcBottom, bcTop;
+} Discretization;
+
+extern void initDiscretization(Discretization*, Parameter*);
+extern void computeRHS(Discretization*);
+extern void normalizePressure(Discretization*);
+extern void computeTimestep(Discretization*);
+extern void setBoundaryConditions(Discretization*);
+extern void setSpecialBoundaryCondition(Discretization*);
+extern void setObjectBoundaryCondition(Discretization*);
+extern void computeFG(Discretization*);
+extern void adaptUV(Discretization*);
+extern void writeResult(Discretization*);
+extern void print(Discretization*, double*);
+extern void printGrid(Discretization*, int*);
+#endif

EnhancedSolver/2D-seq/src/grid.h

@@ -7,10 +7,32 @@
 #ifndef __GRID_H_
 #define __GRID_H_
 
+#define S(i, j) s[(j) * (imax + 2) + (i)]
+
+enum OBJECTBOUNDARY {
+    FLUID = 0,
+    TOP,
+    BOTTOM,
+    LEFT,
+    RIGHT,
+    TOPLEFT,
+    BOTTOMLEFT,
+    TOPRIGHT,
+    BOTTOMRIGHT,
+    OBSTACLE
+};
+
+enum SHAPE { NOSHAPE = 0, RECT, CIRCLE };
+
 typedef struct {
     double dx, dy;
     int imax, jmax;
     double xlength, ylength;
+    int* s;
 } Grid;
 
+static inline int gridIsFluid(Grid* g, int i, int j)
+{
+    return g->s[j * (g->imax + 2) + i] == FLUID;
+}
 #endif // __GRID_H_

EnhancedSolver/2D-seq/src/main.c

@@ -8,6 +8,7 @@
 #include <stdlib.h>
 #include <unistd.h>
 
+#include "discretization.h"
 #include "parameter.h"
 #include "particletracing.h"
 #include "progress.h"
@@ -16,52 +17,51 @@
 
 int main(int argc, char** argv)
 {
-    double S, E;
+    double timeStart, timeStop;
-    Parameter params;
+    Parameter p;
-    Solver solver;
+    Discretization d;
+    Solver s;
     ParticleTracer particletracer;
-    initParameter(&params);
+
+    initParameter(&p);
 
     if (argc != 2) {
         printf("Usage: %s <configFile>\n", argv[0]);
         exit(EXIT_SUCCESS);
     }
 
-    readParameter(&params, argv[1]);
+    readParameter(&p, argv[1]);
-    printParameter(&params);
+    printParameter(&p);
-    initSolver(&solver, &params);
+    initDiscretization(&d, &p);
-    printf("initsolver done\n");
+    initSolver(&s, &d, &p);
-
+    initParticleTracer(&particletracer, &d.grid, &p);
-    initParticleTracer(&particletracer, &solver.grid, &params);
     printParticleTracerParameters(&particletracer);
 
 #ifndef VERBOSE
-    initProgress(solver.te);
+    initProgress(d.te);
 #endif
 
-    double tau = solver.tau;
+    double tau = d.tau;
-    double te  = solver.te;
+    double te  = d.te;
     double t   = 0.0;
     int nt     = 0;
 
-    S = getTimeStamp();
+    timeStart = getTimeStamp();
+
     while (t <= te) {
-        if (tau > 0.0) computeTimestep(&solver);
+        if (tau > 0.0) computeTimestep(&d);
-        setBoundaryConditions(&solver);
+        setBoundaryConditions(&d);
-        setSpecialBoundaryCondition(&solver);
+        setSpecialBoundaryCondition(&d);
-        setObjectBoundaryCondition(&solver);
+        setObjectBoundaryCondition(&d);
-        computeFG(&solver);
-        computeRHS(&solver);
-        if (nt % 100 == 0) normalizePressure(&solver);
-        multiGrid(&solver);
 
-        adaptUV(&solver);
+        computeFG(&d);
+        computeRHS(&d);
+        if (nt % 100 == 0) normalizePressure(&d);
+        solve(&s, d.p, d.rhs);
+        adaptUV(&d);
+        trace(&particletracer, d.u, d.v, d.dt, t);
 
-        /* Added function for particle tracing. Will inject and advance particles as per
+        t += d.dt;
-         * timePeriod */
-        trace(&particletracer, solver.u, solver.v, solver.s, t);
-
-        t += solver.dt;
         nt++;
 
 #ifdef VERBOSE
@@ -70,14 +70,12 @@ int main(int argc, char** argv)
         printProgress(t);
 #endif
     }
-    printf("Total particles : %d\n", particletracer.totalParticles);
 
-    E = getTimeStamp();
+    timeStop = getTimeStamp();
+
     stopProgress();
-
     freeParticles(&particletracer);
-
+    printf("Solution took %.2fs\n", timeStop - timeStart);
-    printf("Solution took %.2fs\n", E - S);
+    writeResult(&d);
-    writeResult(&solver);
     return EXIT_SUCCESS;
 }

EnhancedSolver/2D-seq/src/particletracing.c

@@ -8,14 +8,13 @@
 #include <stdlib.h>
 #include <string.h>
 
+#include "grid.h"
 #include "vtkWriter.h"
 
 #define U(i, j) u[(j) * (imax + 2) + (i)]
 #define V(i, j) v[(j) * (imax + 2) + (i)]
 #define S(i, j) s[(j) * (imax + 2) + (i)]
 
-static int ts = 0;
-
 void printParticles(ParticleTracer* p)
 {
     for (int i = 0; i < p->totalParticles; ++i) {
@@ -25,29 +24,28 @@ void printParticles(ParticleTracer* p)
             p->particlePool[i].flag);
     }
 }
-void injectParticles(ParticleTracer* p)
+
+static void injectParticles(ParticleTracer* p)
 {
-    for (int i = 0; i < p->numberOfParticles; ++i) {
+    if (p->totalParticles + p->numParticlesInLine > p->numAllocatedParticles) {
+        return;
+    }
+    for (int i = 0; i < p->numParticlesInLine; ++i) {
         p->particlePool[p->pointer].x    = p->linSpaceLine[i].x;
         p->particlePool[p->pointer].y    = p->linSpaceLine[i].y;
         p->particlePool[p->pointer].flag = true;
-        ++(p->pointer);
+        p->pointer++;
-        ++(p->totalParticles);
+        p->totalParticles++;
     }
 }
 
-void advanceParticles(ParticleTracer* p,
+static void advanceParticles(
-    double* restrict u,
+    ParticleTracer* p, double* restrict u, double* restrict v, double dt)
-    double* restrict v,
-    int* restrict s,
-    double time)
 {
     int imax       = p->grid->imax;
     int jmax       = p->grid->jmax;
-
     double dx      = p->grid->dx;
     double dy      = p->grid->dy;
-
     double xlength = p->grid->xlength;
     double ylength = p->grid->ylength;
 
@@ -64,14 +62,14 @@ void advanceParticles(ParticleTracer* p,
             double x2 = (double)iCoord * dx;
             double y2 = ((double)jCoord - 0.5) * dy;
 
-            double u_n = (1.0 / (dx * dy)) *
+            double intU = (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 + p->dt * u_n;
+            double newX          = x + dt * intU;
-            p->particlePool[i].x = new_x;
+            p->particlePool[i].x = newX;
 
             iCoord = (int)((x + 0.5 * dx) / dx) + 1;
             jCoord = (int)(y / dy) + 1;
@@ -81,50 +79,58 @@ void advanceParticles(ParticleTracer* p,
             x2 = ((double)iCoord - 0.5) * dx;
             y2 = (double)jCoord * dy;
 
-            double v_n = (1.0 / (dx * dy)) *
+            double intV = (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 + p->dt * v_n;
+            double newY          = y + dt * intV;
-            p->particlePool[i].y = new_y;
+            p->particlePool[i].y = newY;
 
-            // printf("\tOld X : %.2f, New X : %.2f, iCoord : %d\n\tOld Y : %.2f, New Y :
+            if (((newX < 0.0) || (newX > xlength) || (newY < 0.0) || (newY > ylength))) {
-            // %.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))) {
                 p->particlePool[i].flag = false;
+                p->removedParticles++;
             }
-            int i_new = new_x / dx, j_new = new_y / dy;
+
-            if (S(i_new, j_new) != NONE) {
+            int newI = newX / dx, newJ = newY / dy;
+
+            if (!gridIsFluid(p->grid, newI, newJ)) {
                 p->particlePool[i].flag = false;
+                p->removedParticles++;
+                printf("Forbidden movement of particle into obstacle!\n");
             }
         }
     }
 }
 
-void freeParticles(ParticleTracer* p)
+static void compress(ParticleTracer* p)
 {
-    if (p->particlePool != NULL) {
+    Particle* memPool = p->particlePool;
-        free(p->particlePool);
+    Particle tempPool[p->totalParticles];
-        free(p->linSpaceLine);
+    int totalParticles = 0;
+
+    printf("Performing compression ...");
+
+    for (int i = 0; i < p->totalParticles; i++) {
+        if (memPool[i].flag == 1) {
+            tempPool[totalParticles].x    = memPool[i].x;
+            tempPool[totalParticles].y    = memPool[i].y;
+            tempPool[totalParticles].flag = memPool[i].flag;
+            totalParticles++;
         }
+    }
+
+    printf(" remove %d particles\n", p->totalParticles - totalParticles);
+    p->totalParticles   = totalParticles;
+    p->removedParticles = 0;
+    p->pointer          = totalParticles + 1;
+    memcpy(p->particlePool, tempPool, totalParticles * sizeof(Particle));
 }
 
 void writeParticles(ParticleTracer* p)
 {
+    static int ts   = 0;
     VtkOptions opts = { .particletracer = p };
 
     char filename[50];
@@ -154,44 +160,43 @@ void writeParticles(ParticleTracer* p)
     ++ts;
 }
 
-void initParticleTracer(ParticleTracer* p, Grid* grid, Parameter* params)
+void initParticleTracer(ParticleTracer* pt, Grid* g, Parameter* p)
 {
-    p->numberOfParticles = params->numberOfParticles;
+    pt->numParticlesInLine = p->numberOfParticles;
-    p->startTime         = params->startTime;
+    pt->startTime          = p->startTime;
-    p->injectTimePeriod  = params->injectTimePeriod;
+    pt->injectTimePeriod   = p->injectTimePeriod;
-    p->writeTimePeriod   = params->writeTimePeriod;
+    pt->writeTimePeriod    = p->writeTimePeriod;
+    pt->grid               = g;
 
-    p->dt   = params->dt;
+    pt->x1 = p->x1;
-    p->grid = grid;
+    pt->y1 = p->y1;
+    pt->x2 = p->x2;
+    pt->y2 = p->y2;
 
-    p->x1 = params->x1;
+    pt->lastInjectTime = p->startTime;
-    p->y1 = params->y1;
+    pt->lastWriteTime  = p->startTime;
-    p->x2 = params->x2;
-    p->y2 = params->y2;
 
-    p->lastInjectTime = params->startTime;
+    pt->pointer          = 0;
-    p->lastUpdateTime = params->startTime;
+    pt->removedParticles = 0;
-    p->lastWriteTime  = params->startTime;
+    pt->totalParticles   = 0;
 
-    p->pointer        = 0;
+    if (p->te > p->startTime) {
-    p->totalParticles = 0;
+        pt->numAllocatedParticles = ((p->te - p->startTime) / p->injectTimePeriod) *
+                                    p->numberOfParticles;
+        pt->numAllocatedParticles += (2 * p->numberOfParticles);
 
-    if (params->te > params->startTime) {
+        pt->particlePool = malloc(sizeof(Particle) * pt->numAllocatedParticles);
-        p->estimatedNumParticles = ((params->te - params->startTime) + 2) *
+        pt->linSpaceLine = malloc(sizeof(Particle) * pt->numParticlesInLine);
-                                   params->numberOfParticles;
 
-        p->particlePool = malloc(sizeof(Particle) * p->estimatedNumParticles);
+        for (int i = 0; i < pt->numParticlesInLine; ++i) {
-        p->linSpaceLine = malloc(sizeof(Particle) * p->numberOfParticles);
+            double spacing           = (double)i / (double)(pt->numParticlesInLine - 1);
-
+            pt->linSpaceLine[i].x    = spacing * pt->x1 + (1.0 - spacing) * pt->x2;
-        for (int i = 0; i < p->numberOfParticles; ++i) {
+            pt->linSpaceLine[i].y    = spacing * pt->y1 + (1.0 - spacing) * pt->y2;
-            double spacing          = (double)i / (double)(p->numberOfParticles - 1);
+            pt->linSpaceLine[i].flag = true;
-            p->linSpaceLine[i].x    = spacing * p->x1 + (1.0 - spacing) * p->x2;
-            p->linSpaceLine[i].y    = spacing * p->y1 + (1.0 - spacing) * p->y2;
-            p->linSpaceLine[i].flag = true;
         }
     } else {
-        p->particlePool = NULL;
+        pt->particlePool = NULL;
-        p->linSpaceLine = NULL;
+        pt->linSpaceLine = NULL;
     }
 }
 
@@ -199,7 +204,7 @@ void printParticleTracerParameters(ParticleTracer* p)
 {
     printf("Particle Tracing data:\n");
     printf("\tNumber of particles : %d being injected for every period of %.2f\n",
-        p->numberOfParticles,
+        p->numParticlesInLine,
         p->injectTimePeriod);
     printf("\tstartTime : %.2f\n", p->startTime);
     printf("\t(Line along which the particles are to be injected) \n\tx1 : %.2f, y1 : "
@@ -209,43 +214,33 @@ void printParticleTracerParameters(ParticleTracer* p)
         p->x2,
         p->y2);
     printf("\tPointer : %d, TotalParticles : %d\n", p->pointer, p->totalParticles);
-    printf("\tdt : %.2f, dx : %.2f, dy : %.2f\n", p->dt, p->grid->dx, p->grid->dy);
 }
 
-void trace(ParticleTracer* p, double* u, double* v, int* s, double time)
+void trace(ParticleTracer* p, double* u, double* v, double dt, double time)
 {
     if (time >= p->startTime) {
-        // printParticles(particletracer);
         if ((time - p->lastInjectTime) >= p->injectTimePeriod) {
             injectParticles(p);
             p->lastInjectTime = time;
         }
+
         if ((time - p->lastWriteTime) >= p->writeTimePeriod) {
             writeParticles(p);
             p->lastWriteTime = time;
         }
-        advanceParticles(p, u, v, s, time);
+
+        advanceParticles(p, u, v, dt);
+
+        if (p->removedParticles > (p->totalParticles * 0.2)) {
             compress(p);
-        p->lastUpdateTime = time;
+        }
     }
 }
 
-void compress(ParticleTracer* p)
+void freeParticles(ParticleTracer* p)
 {
-    Particle* memPool = p->particlePool;
+    if (p->particlePool != NULL) {
-    Particle tempPool[p->totalParticles];
+        free(p->particlePool);
-    int totalParticles = 0;
+        free(p->linSpaceLine);
-
-    for (int i = 0; i < p->totalParticles; ++i) {
-        if (memPool[i].flag == 1) {
-            tempPool[totalParticles].x    = memPool[i].x;
-            tempPool[totalParticles].y    = memPool[i].y;
-            tempPool[totalParticles].flag = memPool[i].flag;
-            ++totalParticles;
     }
-    }
-
-    p->totalParticles = totalParticles;
-    p->pointer        = totalParticles + 1;
-    memcpy(p->particlePool, tempPool, totalParticles * sizeof(Particle));
 }

EnhancedSolver/2D-seq/src/particletracing.h

@@ -19,13 +19,12 @@ typedef struct {
 } Particle;
 
 typedef struct {
-    int numberOfParticles, totalParticles;
+    int numParticlesInLine, removedParticles, totalParticles;
     double startTime, injectTimePeriod, writeTimePeriod;
     double lastInjectTime, lastUpdateTime, lastWriteTime;
 
-    int estimatedNumParticles;
+    int numAllocatedParticles;
 
-    double dt;
     Particle* linSpaceLine;
     Particle* particlePool;
 
@@ -35,12 +34,8 @@ typedef struct {
 } ParticleTracer;
 
 extern void initParticleTracer(ParticleTracer*, Grid*, Parameter*);
-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*, double, double);
-extern void trace(ParticleTracer*, double*, double*, int*, double);
-extern void compress(ParticleTracer*);
 #endif

EnhancedSolver/2D-seq/src/progress.h

@@ -9,6 +9,6 @@
 
 extern void initProgress(double);
 extern void printProgress(double);
-extern void stopProgress();
+extern void stopProgress(void);
 
 #endif

EnhancedSolver/2D-seq/src/solver-mg.c

@@ -0,0 +1,186 @@
+/*
+ * Copyright (C) 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 "allocate.h"
+#include "solver.h"
+#include "util.h"
+
+#define FINEST_LEVEL   0
+#define COARSEST_LEVEL (s->levels - 1)
+#define S(i, j)        s[(j) * (imax + 2) + (i)]
+#define E(i, j)        e[(j) * (imax + 2) + (i)]
+#define R(i, j)        r[(j) * (imax + 2) + (i)]
+#define OLD(i, j)      old[(j) * (imax + 2) + (i)]
+
+static void restrictMG(Solver* s, int level, int imax, int jmax)
+{
+    double* r   = s->r[level + 1];
+    double* old = s->r[level];
+
+    for (int j = 1; j < jmax + 1; j++) {
+        for (int i = 1; i < imax + 1; i++) {
+            R(i, j) = (OLD(2 * i - 1, 2 * j - 1) + OLD(2 * i, 2 * j - 1) * 2 +
+                          OLD(2 * i + 1, 2 * j - 1) + OLD(2 * i - 1, 2 * j) * 2 +
+                          OLD(2 * i, 2 * j) * 4 + OLD(2 * i + 1, 2 * j) * 2 +
+                          OLD(2 * i - 1, 2 * j + 1) + OLD(2 * i, 2 * j + 1) * 2 +
+                          OLD(2 * i + 1, 2 * j + 1)) /
+                      16.0;
+        }
+    }
+}
+
+static void prolongate(Solver* s, int level, int imax, int jmax)
+{
+    double* old = s->r[level + 1];
+    double* e   = s->r[level];
+
+    for (int j = 2; j < jmax + 1; j += 2) {
+        for (int i = 2; i < imax + 1; i += 2) {
+            E(i, j) = OLD(i / 2, j / 2);
+        }
+    }
+}
+
+static void correct(Solver* s, double* p, int level, int imax, int jmax)
+{
+    double* e = s->e[level];
+
+    for (int j = 1; j < jmax + 1; ++j) {
+        for (int i = 1; i < imax + 1; ++i) {
+            P(i, j) += E(i, j);
+        }
+    }
+}
+
+static void setBoundaryCondition(double* p, int imax, int 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);
+    }
+}
+
+static double smooth(Solver* s, double* p, double* rhs, int level, int imax, int jmax)
+{
+    double dx2    = s->grid->dx * s->grid->dx;
+    double dy2    = s->grid->dy * s->grid->dy;
+    double idx2   = 1.0 / dx2;
+    double idy2   = 1.0 / dy2;
+    double factor = s->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2);
+    double* r     = s->r[level];
+    double res    = 1.0;
+    int pass, jsw, isw;
+
+    jsw = 1;
+
+    for (pass = 0; pass < 2; pass++) {
+        isw = jsw;
+
+        for (int j = 1; j < jmax + 1; j++) {
+            for (int i = isw; i < imax + 1; i += 2) {
+
+                R(i, j) = RHS(i, j) -
+                          ((P(i + 1, j) - 2.0 * P(i, j) + P(i - 1, j)) * idx2 +
+                              (P(i, j + 1) - 2.0 * P(i, j) + P(i, j - 1)) * idy2);
+
+                P(i, j) -= (factor * R(i, j));
+                res += (R(i, j) * R(i, j));
+            }
+            isw = 3 - isw;
+        }
+        jsw = 3 - jsw;
+    }
+
+    res = res / (double)(imax * jmax);
+    return res;
+}
+
+static double multiGrid(Solver* s, double* p, double* rhs, int level, int imax, int jmax)
+{
+    double res = 0.0;
+
+    // coarsest level TODO: Use direct solver?
+    if (level == COARSEST_LEVEL) {
+        for (int i = 0; i < 5; i++) {
+            smooth(s, p, rhs, level, imax, jmax);
+        }
+        return res;
+    }
+
+    // pre-smoothing TODO: Make smoothing steps configurable?
+    for (int i = 0; i < 5; i++) {
+        smooth(s, p, rhs, level, imax, jmax);
+        if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax);
+    }
+
+    // restrict
+    restrictMG(s, level, imax, jmax);
+
+    // MGSolver on residual and error.
+    // TODO: What if there is a rest?
+    multiGrid(s, s->e[level + 1], s->r[level + 1], level + 1, imax / 2, jmax / 2);
+
+    // prolongate
+    prolongate(s, level, imax, jmax);
+
+    // correct p on finer level using residual
+    correct(s, p, level, imax, jmax);
+    if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax);
+
+    // post-smoothing
+    for (int i = 0; i < 5; i++) {
+        res = smooth(s, p, rhs, level, imax, jmax);
+        if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax);
+    }
+
+    return res;
+}
+
+void initSolver(Solver* s, Discretization* d, Parameter* p)
+{
+    s->eps     = p->eps;
+    s->omega   = p->omg;
+    s->itermax = p->itermax;
+    s->levels  = p->levels;
+    s->grid    = &d->grid;
+
+    int imax   = s->grid->imax;
+    int jmax   = s->grid->jmax;
+    int levels = s->levels;
+    printf("Using Multigrid solver with %d levels\n", levels);
+
+    s->r = malloc(levels * sizeof(double*));
+    s->e = malloc(levels * sizeof(double*));
+
+    size_t size = (imax + 2) * (jmax + 2) * sizeof(double);
+
+    for (int j = 0; j < levels; j++) {
+        s->r[j] = allocate(64, size);
+        s->e[j] = allocate(64, size);
+
+        for (int i = 0; i < (imax + 2) * (jmax + 2); i++) {
+            s->r[j][i] = 0.0;
+            s->e[j][i] = 0.0;
+        }
+    }
+}
+
+void solve(Solver* s, double* p, double* rhs)
+{
+    double res = multiGrid(s, p, rhs, 0, s->grid->imax, s->grid->jmax);
+
+#ifdef VERBOSE
+    printf("Residuum: %.6f\n", res);
+#endif
+}

EnhancedSolver/2D-seq/src/solver-sor.c

@@ -0,0 +1,128 @@
+/*
+ * Copyright (C)  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 "solver.h"
+#include "util.h"
+
+void initSolver(Solver* s, Discretization* d, Parameter* p)
+{
+    s->grid    = &d->grid;
+    s->itermax = p->itermax;
+    s->eps     = p->eps;
+    s->omega   = p->omg;
+}
+
+void solveSOR(Solver* solver, double* p, double* rhs)
+{
+    int imax      = solver->grid->imax;
+    int jmax      = solver->grid->jmax;
+    double eps    = solver->eps;
+    int itermax   = solver->itermax;
+    double dx2    = solver->grid->dx * solver->grid->dx;
+    double dy2    = solver->grid->dy * solver->grid->dy;
+    double idx2   = 1.0 / dx2;
+    double idy2   = 1.0 / dy2;
+    double factor = solver->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2);
+    double epssq  = eps * eps;
+    int it        = 0;
+    double res    = 1.0;
+
+    while ((res >= epssq) && (it < itermax)) {
+        res = 0.0;
+
+        for (int j = 1; j < jmax + 1; j++) {
+            for (int i = 1; i < imax + 1; i++) {
+
+                double r = RHS(i, j) -
+                           ((P(i + 1, j) - 2.0 * P(i, j) + P(i - 1, j)) * idx2 +
+                               (P(i, j + 1) - 2.0 * P(i, j) + P(i, j - 1)) * idy2);
+
+                P(i, j) -= (factor * r);
+                res += (r * r);
+            }
+        }
+
+        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);
+        }
+
+        res = res / (double)(imax * jmax);
+#ifdef DEBUG
+        printf("%d Residuum: %e\n", it, res);
+#endif
+        it++;
+    }
+
+#ifdef VERBOSE
+    printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
+#endif
+}
+
+void solve(Solver* solver, double* p, double* rhs)
+{
+    int imax      = solver->grid->imax;
+    int jmax      = solver->grid->jmax;
+    double eps    = solver->eps;
+    int itermax   = solver->itermax;
+    double dx2    = solver->grid->dx * solver->grid->dx;
+    double dy2    = solver->grid->dy * solver->grid->dy;
+    double idx2   = 1.0 / dx2;
+    double idy2   = 1.0 / dy2;
+    double factor = solver->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2);
+    double epssq  = eps * eps;
+    int it        = 0;
+    double res    = 1.0;
+    int pass, jsw, isw;
+
+    while ((res >= epssq) && (it < itermax)) {
+        res = 0.0;
+        jsw = 1;
+
+        for (pass = 0; pass < 2; pass++) {
+            isw = jsw;
+
+            for (int j = 1; j < jmax + 1; j++) {
+                for (int i = isw; i < imax + 1; i += 2) {
+
+                    double r = RHS(i, j) -
+                               ((P(i + 1, j) - 2.0 * P(i, j) + P(i - 1, j)) * idx2 +
+                                   (P(i, j + 1) - 2.0 * P(i, j) + P(i, j - 1)) * idy2);
+
+                    P(i, j) -= (factor * r);
+                    res += (r * r);
+                }
+                isw = 3 - isw;
+            }
+            jsw = 3 - jsw;
+        }
+
+        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);
+        }
+
+        res = res / (double)(imax * jmax);
+#ifdef DEBUG
+        printf("%d Residuum: %e\n", it, res);
+#endif
+        it++;
+    }
+
+#ifdef VERBOSE
+    printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
+#endif
+}

EnhancedSolver/2D-seq/src/solver.c

@@ -1,900 +0,0 @@
-/*
- * Copyright (C)  NHR@FAU, University Erlangen-Nuremberg.
- * All rights reserved. This file is part of nusif-solver.
- * Use of this source code is governed by a MIT style
- * license that can be found in the LICENSE file.
- */
-#include <float.h>
-#include <math.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-
-#include "allocate.h"
-#include "parameter.h"
-#include "solver.h"
-#include "util.h"
-
-#define P(i, j)    p[(j) * (imax + 2) + (i)]
-#define F(i, j)    f[(j) * (imax + 2) + (i)]
-#define G(i, j)    g[(j) * (imax + 2) + (i)]
-#define U(i, j)    u[(j) * (imax + 2) + (i)]
-#define V(i, j)    v[(j) * (imax + 2) + (i)]
-#define S(i, j)    s[(j) * (imax + 2) + (i)]
-#define E(i, j)    e[(j) * (imax + 2) + (i)]
-#define R(i, j)    r[(j) * (imax + 2) + (i)]
-#define oldR(i, j) oldr[(j) * (imax + 2) + (i)]
-#define oldE(i, j) olde[(j) * (imax + 2) + (i)]
-#define RHS(i, j)  rhs[(j) * (imax + 2) + (i)]
-
-static double distance(double i, double j, double iCenter, double jCenter)
-{
-    return sqrt(pow(iCenter - i, 2) + pow(jCenter - j, 2) * 1.0);
-}
-
-void print(Solver* solver, double* grid)
-{
-    int imax = solver->grid.imax;
-    int jmax = solver->grid.jmax;
-
-    for (int j = 0; j < jmax + 2; j++) {
-        printf("%02d: ", j);
-        for (int i = 0; i < imax + 2; i++) {
-            printf("%3.2f  ", grid[j * (imax + 2) + i]);
-        }
-        printf("\n");
-    }
-    fflush(stdout);
-}
-
-void printGrid(Solver* solver, int* grid)
-{
-    int imax = solver->grid.imax;
-    int jmax = solver->grid.jmax;
-
-    for (int j = 0; j < jmax + 2; j++) {
-        printf("%02d: ", j);
-        for (int i = 0; i < imax + 2; i++) {
-            printf("%2d  ", grid[j * (imax + 2) + i]);
-        }
-        printf("\n");
-    }
-    fflush(stdout);
-}
-
-static void printConfig(Solver* solver)
-{
-    printf("Parameters for #%s#\n", solver->problem);
-    printf("Boundary conditions Left:%d Right:%d Bottom:%d Top:%d\n",
-        solver->bcLeft,
-        solver->bcRight,
-        solver->bcBottom,
-        solver->bcTop);
-    printf("\tReynolds number: %.2f\n", solver->re);
-    printf("\tGx Gy: %.2f %.2f\n", solver->gx, solver->gy);
-    printf("Geometry data:\n");
-    printf("\tDomain box size (x, y): %.2f, %.2f\n",
-        solver->grid.xlength,
-        solver->grid.ylength);
-    printf("\tCells (x, y): %d, %d\n", solver->grid.imax, solver->grid.jmax);
-    printf("Timestep parameters:\n");
-    printf("\tDefault stepsize: %.2f, Final time %.2f\n", solver->dt, solver->te);
-    printf("\tdt bound: %.6f\n", solver->dtBound);
-    printf("\tTau factor: %.2f\n", solver->tau);
-    printf("Iterative solver parameters:\n");
-    printf("\tMax iterations: %d\n", solver->itermax);
-    printf("\tepsilon (stopping tolerance) : %f\n", solver->eps);
-    printf("\tgamma factor: %f\n", solver->gamma);
-    printf("\tomega (SOR relaxation): %f\n", solver->omega);
-}
-
-void initSolver(Solver* solver, Parameter* params)
-{
-    solver->problem      = params->name;
-    solver->bcLeft       = params->bcLeft;
-    solver->bcRight      = params->bcRight;
-    solver->bcBottom     = params->bcBottom;
-    solver->bcTop        = params->bcTop;
-    solver->grid.imax    = params->imax;
-    solver->grid.jmax    = params->jmax;
-    solver->grid.xlength = params->xlength;
-    solver->grid.ylength = params->ylength;
-    solver->grid.dx      = params->xlength / params->imax;
-    solver->grid.dy      = params->ylength / params->jmax;
-    solver->eps          = params->eps;
-    solver->omega        = params->omg;
-    solver->itermax      = params->itermax;
-    solver->re           = params->re;
-    solver->gx           = params->gx;
-    solver->gy           = params->gy;
-    solver->dt           = params->dt;
-    solver->te           = params->te;
-    solver->tau          = params->tau;
-    solver->gamma        = params->gamma;
-    solver->rho          = params->rho;
-    solver->levels       = params->levels;
-    solver->currentlevel = 0;
-
-    int imax   = solver->grid.imax;
-    int jmax   = solver->grid.jmax;
-    int levels = solver->levels;
-
-    size_t size_level = levels * (imax + 2) * (jmax + 2) * sizeof(double);
-
-    size_t size = (imax + 2) * (jmax + 2) * sizeof(double);
-    solver->u   = allocate(64, size);
-    solver->v   = allocate(64, size);
-    solver->s   = allocate(64, size);
-    solver->p   = allocate(64, size);
-    solver->rhs = allocate(64, size);
-    solver->f   = allocate(64, size);
-    solver->g   = allocate(64, size);
-
-    solver->r = malloc(levels * sizeof(double*));
-    solver->e = malloc(levels * sizeof(double*));
-
-    for (int j = 0; j < levels; ++j) {
-
-        solver->r[j] = allocate(64, size);
-        solver->e[j] = allocate(64, size);
-    }
-
-    for (int i = 0; i < (imax + 2) * (jmax + 2); i++) {
-
-        solver->u[i]   = params->u_init;
-        solver->v[i]   = params->v_init;
-        solver->p[i]   = params->p_init;
-        solver->rhs[i] = 0.0;
-        solver->f[i]   = 0.0;
-        solver->g[i]   = 0.0;
-        solver->s[i]   = NONE;
-        for (int j = 0; j < levels; ++j) {
-
-            solver->r[j][i] = 0.0;
-            solver->e[j][i] = 0.0;
-        }
-    }
-
-    double dx        = solver->grid.dx;
-    double dy        = solver->grid.dy;
-    double invSqrSum = 1.0 / (dx * dx) + 1.0 / (dy * dy);
-    solver->dtBound  = 0.5 * solver->re * 1.0 / invSqrSum;
-
-    double xCenter = 0, yCenter = 0, radius = 0;
-    double x1 = 0, x2 = 0, y1 = 0, y2 = 0;
-
-    int* s = solver->s;
-
-    switch (params->shape) {
-    case NOSHAPE:
-        break;
-    case RECT:
-        x1 = params->xCenter - params->xRectLength / 2;
-        x2 = params->xCenter + params->xRectLength / 2;
-        y1 = params->yCenter - params->yRectLength / 2;
-        y2 = params->yCenter + params->yRectLength / 2;
-
-        for (int j = 1; j < jmax + 1; j++) {
-            for (int i = 1; i < imax + 1; i++) {
-                if ((x1 <= (i * dx)) && ((i * dx) <= x2) && (y1 <= (j * dy)) &&
-                    ((j * dy) <= y2)) {
-                    S(i, j) = LOCAL;
-                }
-            }
-        }
-
-        break;
-    case CIRCLE:
-        xCenter = params->xCenter;
-        yCenter = params->yCenter;
-        radius  = params->circleRadius;
-
-        for (int j = 1; j < jmax + 1; j++) {
-            for (int i = 1; i < imax + 1; i++) {
-                if (distance((i * dx), (j * dy), xCenter, yCenter) <= radius) {
-                    S(i, j) = LOCAL;
-                }
-            }
-        }
-
-        break;
-    default:
-        break;
-    }
-
-    if (params->shape != NOSHAPE) {
-        for (int j = 1; j < jmax + 1; j++) {
-            for (int i = 1; i < imax + 1; i++) {
-
-                if (S(i, j - 1) == NONE && S(i, j + 1) == LOCAL && S(i, j) == LOCAL)
-                    S(i, j) = BOTTOM; // TOP
-                if (S(i - 1, j) == NONE && S(i + 1, j) == LOCAL && S(i, j) == LOCAL)
-                    S(i, j) = LEFT; // LEFT
-                if (S(i + 1, j) == NONE && S(i - 1, j) == LOCAL && S(i, j) == LOCAL)
-                    S(i, j) = RIGHT; // RIGHT
-                if (S(i, j + 1) == NONE && S(i, j - 1) == LOCAL && S(i, j) == LOCAL)
-                    S(i, j) = TOP; // BOTTOM
-                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) == LEFT || S(i, j) == BOTTOM))
-                    S(i, j) = BOTTOMLEFT; // TOPLEFT
-                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) == RIGHT || S(i, j) == BOTTOM))
-                    S(i, j) = BOTTOMRIGHT; // TOPRIGHT
-                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) == LEFT || S(i, j) == TOP))
-                    S(i, j) = TOPLEFT; // BOTTOMLEFT
-                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) == RIGHT || S(i, j) == TOP))
-                    S(i, j) = TOPRIGHT; // BOTTOMRIGHT
-            }
-        }
-    }
-
-#ifdef VERBOSE
-    printConfig(solver);
-#endif
-}
-
-static double maxElement(Solver* solver, double* m)
-{
-    int size      = (solver->grid.imax + 2) * (solver->grid.jmax + 2);
-    double maxval = DBL_MIN;
-
-    for (int i = 0; i < size; i++) {
-        maxval = MAX(maxval, fabs(m[i]));
-    }
-
-    return maxval;
-}
-
-void computeRHS(Solver* solver)
-{
-    int imax    = solver->grid.imax;
-    int jmax    = solver->grid.jmax;
-    double idx  = 1.0 / solver->grid.dx;
-    double idy  = 1.0 / solver->grid.dy;
-    double idt  = 1.0 / solver->dt;
-    double* rhs = solver->rhs;
-    double* f   = solver->f;
-    double* g   = solver->g;
-    int* s      = solver->s;
-
-    for (int j = 1; j < jmax + 1; j++) {
-        for (int i = 1; i < imax + 1; i++) {
-            RHS(i, j) = idt *
-                        ((F(i, j) - F(i - 1, j)) * idx + (G(i, j) - G(i, j - 1)) * idy);
-        }
-    }
-}
-
-void normalizePressure(Solver* solver)
-{
-    int size    = (solver->grid.imax + 2) * (solver->grid.jmax + 2);
-    double* p   = solver->p;
-    double avgP = 0.0;
-
-    for (int i = 0; i < size; i++) {
-        avgP += p[i];
-    }
-    avgP /= size;
-
-    for (int i = 0; i < size; i++) {
-        p[i] = p[i] - avgP;
-    }
-}
-
-void computeTimestep(Solver* solver)
-{
-    double dt   = solver->dtBound;
-    double dx   = solver->grid.dx;
-    double dy   = solver->grid.dy;
-    double umax = maxElement(solver, solver->u);
-    double vmax = maxElement(solver, solver->v);
-
-    if (umax > 0) {
-        dt = (dt > dx / umax) ? dx / umax : dt;
-    }
-    if (vmax > 0) {
-        dt = (dt > dy / vmax) ? dy / vmax : dt;
-    }
-
-    solver->dt = dt * solver->tau;
-}
-
-void setBoundaryConditions(Solver* solver)
-{
-    int imax  = solver->grid.imax;
-    int jmax  = solver->grid.jmax;
-    double* u = solver->u;
-    double* v = solver->v;
-
-    // Left boundary
-    switch (solver->bcLeft) {
-    case NOSLIP:
-        for (int j = 1; j < jmax + 1; j++) {
-            U(0, j) = 0.0;
-            V(0, j) = -V(1, j);
-        }
-        break;
-    case SLIP:
-        for (int j = 1; j < jmax + 1; j++) {
-            U(0, j) = 0.0;
-            V(0, j) = V(1, j);
-        }
-        break;
-    case OUTFLOW:
-        for (int j = 1; j < jmax + 1; j++) {
-            U(0, j) = U(1, j);
-            V(0, j) = V(1, j);
-        }
-        break;
-    case PERIODIC:
-        break;
-    }
-
-    // Right boundary
-    switch (solver->bcRight) {
-    case NOSLIP:
-        for (int j = 1; j < jmax + 1; j++) {
-            U(imax, j)     = 0.0;
-            V(imax + 1, j) = -V(imax, j);
-        }
-        break;
-    case SLIP:
-        for (int j = 1; j < jmax + 1; j++) {
-            U(imax, j)     = 0.0;
-            V(imax + 1, j) = V(imax, j);
-        }
-        break;
-    case OUTFLOW:
-        for (int j = 1; j < jmax + 1; j++) {
-            U(imax, j)     = U(imax - 1, j);
-            V(imax + 1, j) = V(imax, j);
-        }
-        break;
-    case PERIODIC:
-        break;
-    }
-
-    // Bottom boundary
-    switch (solver->bcBottom) {
-    case NOSLIP:
-        for (int i = 1; i < imax + 1; i++) {
-            V(i, 0) = 0.0;
-            U(i, 0) = -U(i, 1);
-        }
-        break;
-    case SLIP:
-        for (int i = 1; i < imax + 1; i++) {
-            V(i, 0) = 0.0;
-            U(i, 0) = U(i, 1);
-        }
-        break;
-    case OUTFLOW:
-        for (int i = 1; i < imax + 1; i++) {
-            U(i, 0) = U(i, 1);
-            V(i, 0) = V(i, 1);
-        }
-        break;
-    case PERIODIC:
-        break;
-    }
-
-    // Top boundary
-    switch (solver->bcTop) {
-    case NOSLIP:
-        for (int i = 1; i < imax + 1; i++) {
-            V(i, jmax)     = 0.0;
-            U(i, jmax + 1) = -U(i, jmax);
-        }
-        break;
-    case SLIP:
-        for (int i = 1; i < imax + 1; i++) {
-            V(i, jmax)     = 0.0;
-            U(i, jmax + 1) = U(i, jmax);
-        }
-        break;
-    case OUTFLOW:
-        for (int i = 1; i < imax + 1; i++) {
-            U(i, jmax + 1) = U(i, jmax);
-            V(i, jmax)     = V(i, jmax - 1);
-        }
-        break;
-    case PERIODIC:
-        break;
-    }
-}
-
-void setSpecialBoundaryCondition(Solver* solver)
-{
-    int imax   = solver->grid.imax;
-    int jmax   = solver->grid.jmax;
-    double mDy = solver->grid.dy;
-    double* u  = solver->u;
-    int* s     = solver->s;
-
-    if (strcmp(solver->problem, "dcavity") == 0) {
-        for (int i = 1; i < imax; i++) {
-            U(i, jmax + 1) = 2.0 - U(i, jmax);
-        }
-    } else if (strcmp(solver->problem, "canal") == 0) {
-        double ylength = solver->grid.ylength;
-        double y;
-
-        for (int j = 1; j < jmax + 1; j++) {
-            y       = mDy * (j - 0.5);
-            U(0, j) = y * (ylength - y) * 4.0 / (ylength * ylength);
-        }
-    } else if (strcmp(solver->problem, "backstep") == 0) {
-        for (int j = 1; j < jmax + 1; j++) {
-            if (S(0, j) == NONE) U(0, j) = 1.0;
-        }
-    } else if (strcmp(solver->problem, "karman") == 0) {
-        for (int j = 1; j < jmax + 1; j++) {
-            U(0, j) = 1.0;
-        }
-    }
-}
-
-void setObjectBoundaryCondition(Solver* solver)
-{
-    int imax  = solver->grid.imax;
-    int jmax  = solver->grid.jmax;
-    double* u = solver->u;
-    double* v = solver->v;
-    int* s    = solver->s;
-
-    for (int j = 1; j < jmax + 1; j++) {
-        for (int i = 1; i < imax + 1; i++) {
-            switch (S(i, j)) {
-            case TOP:
-                U(i, j)     = -U(i, j + 1);
-                U(i - 1, j) = -U(i - 1, j + 1);
-                V(i, j)     = 0.0;
-                break;
-            case BOTTOM:
-                U(i, j)     = -U(i, j - 1);
-                U(i - 1, j) = -U(i - 1, j - 1);
-                V(i, j)     = 0.0;
-                break;
-            case LEFT:
-                U(i - 1, j) = 0.0;
-                V(i, j)     = -V(i - 1, j);
-                V(i, j - 1) = -V(i - 1, j - 1);
-                break;
-            case RIGHT:
-                U(i, j)     = 0.0;
-                V(i, j)     = -V(i + 1, j);
-                V(i, j - 1) = -V(i + 1, j - 1);
-                break;
-            case TOPLEFT:
-                U(i, j)     = -U(i, j + 1);
-                U(i - 1, j) = 0.0;
-                V(i, j)     = 0.0;
-                V(i, j - 1) = -V(i - 1, j - 1);
-                break;
-            case TOPRIGHT:
-                U(i, j)     = 0.0;
-                U(i - 1, j) = -U(i - 1, j + 1);
-                V(i, j)     = 0.0;
-                V(i, j - 1) = -V(i + 1, j - 1);
-                break;
-            case BOTTOMLEFT:
-                U(i, j)     = -U(i, j - 1);
-                U(i - 1, j) = 0.0;
-                V(i, j)     = -V(i - 1, j);
-                V(i, j - 1) = 0.0;
-                break;
-            case BOTTOMRIGHT:
-                U(i, j)     = 0.0;
-                U(i - 1, j) = -U(i - 1, j - 1);
-                V(i, j)     = -V(i, j + 1);
-                V(i, j - 1) = 0.0;
-                break;
-            }
-        }
-    }
-}
-
-void computeFG(Solver* solver)
-{
-    double* u        = solver->u;
-    double* v        = solver->v;
-    double* f        = solver->f;
-    double* g        = solver->g;
-    int* s           = solver->s;
-    int imax         = solver->grid.imax;
-    int jmax         = solver->grid.jmax;
-    double gx        = solver->gx;
-    double gy        = solver->gy;
-    double gamma     = solver->gamma;
-    double dt        = solver->dt;
-    double inverseRe = 1.0 / solver->re;
-    double inverseDx = 1.0 / solver->grid.dx;
-    double inverseDy = 1.0 / solver->grid.dy;
-    double du2dx, dv2dy, duvdx, duvdy;
-    double du2dx2, du2dy2, dv2dx2, dv2dy2;
-
-    for (int j = 1; j < jmax + 1; j++) {
-        for (int i = 1; i < imax + 1; i++) {
-            if (S(i, j) == NONE) {
-                du2dx = inverseDx * 0.25 *
-                            ((U(i, j) + U(i + 1, j)) * (U(i, j) + U(i + 1, j)) -
-                                (U(i, j) + U(i - 1, j)) * (U(i, j) + U(i - 1, j))) +
-                        gamma * inverseDx * 0.25 *
-                            (fabs(U(i, j) + U(i + 1, j)) * (U(i, j) - U(i + 1, j)) +
-                                fabs(U(i, j) + U(i - 1, j)) * (U(i, j) - U(i - 1, j)));
-
-                duvdy = inverseDy * 0.25 *
-                            ((V(i, j) + V(i + 1, j)) * (U(i, j) + U(i, j + 1)) -
-                                (V(i, j - 1) + V(i + 1, j - 1)) *
-                                    (U(i, j) + U(i, j - 1))) +
-                        gamma * inverseDy * 0.25 *
-                            (fabs(V(i, j) + V(i + 1, j)) * (U(i, j) - U(i, j + 1)) +
-                                fabs(V(i, j - 1) + V(i + 1, j - 1)) *
-                                    (U(i, j) - U(i, j - 1)));
-
-                du2dx2 = inverseDx * inverseDx *
-                         (U(i + 1, j) - 2.0 * U(i, j) + U(i - 1, j));
-                du2dy2 = inverseDy * inverseDy *
-                         (U(i, j + 1) - 2.0 * U(i, j) + U(i, j - 1));
-                F(i, j) = U(i, j) +
-                          dt * (inverseRe * (du2dx2 + du2dy2) - du2dx - duvdy + gx);
-
-                duvdx = inverseDx * 0.25 *
-                            ((U(i, j) + U(i, j + 1)) * (V(i, j) + V(i + 1, j)) -
-                                (U(i - 1, j) + U(i - 1, j + 1)) *
-                                    (V(i, j) + V(i - 1, j))) +
-                        gamma * inverseDx * 0.25 *
-                            (fabs(U(i, j) + U(i, j + 1)) * (V(i, j) - V(i + 1, j)) +
-                                fabs(U(i - 1, j) + U(i - 1, j + 1)) *
-                                    (V(i, j) - V(i - 1, j)));
-
-                dv2dy = inverseDy * 0.25 *
-                            ((V(i, j) + V(i, j + 1)) * (V(i, j) + V(i, j + 1)) -
-                                (V(i, j) + V(i, j - 1)) * (V(i, j) + V(i, j - 1))) +
-                        gamma * inverseDy * 0.25 *
-                            (fabs(V(i, j) + V(i, j + 1)) * (V(i, j) - V(i, j + 1)) +
-                                fabs(V(i, j) + V(i, j - 1)) * (V(i, j) - V(i, j - 1)));
-
-                dv2dx2 = inverseDx * inverseDx *
-                         (V(i + 1, j) - 2.0 * V(i, j) + V(i - 1, j));
-                dv2dy2 = inverseDy * inverseDy *
-                         (V(i, j + 1) - 2.0 * V(i, j) + V(i, j - 1));
-                G(i, j) = V(i, j) +
-                          dt * (inverseRe * (dv2dx2 + dv2dy2) - duvdx - dv2dy + gy);
-            } else {
-                switch (S(i, j)) {
-                case TOP:
-                    G(i, j) = V(i, j);
-                    break;
-                case BOTTOM:
-                    G(i, j - 1) = V(i, j - 1);
-                    break;
-                case LEFT:
-                    F(i - 1, j) = U(i - 1, j);
-                    break;
-                case RIGHT:
-                    F(i, j) = U(i, j);
-                    break;
-                case TOPLEFT:
-                    F(i - 1, j) = U(i - 1, j);
-                    G(i, j)     = V(i, j);
-                    break;
-                case TOPRIGHT:
-                    F(i, j) = U(i, j);
-                    G(i, j) = V(i, j);
-                    break;
-                case BOTTOMLEFT:
-                    F(i - 1, j) = U(i - 1, j);
-                    G(i, j - 1) = V(i, j - 1);
-                    break;
-                case BOTTOMRIGHT:
-                    F(i, j)     = U(i, j);
-                    G(i, j - 1) = V(i, j - 1);
-                    break;
-                }
-            }
-        }
-    }
-
-    /* ---------------------- boundary of F --------------------------- */
-    for (int j = 1; j < jmax + 1; j++) {
-        F(0, j)    = U(0, j);
-        F(imax, j) = U(imax, j);
-    }
-
-    /* ---------------------- boundary of G --------------------------- */
-    for (int i = 1; i < imax + 1; i++) {
-        G(i, 0)    = V(i, 0);
-        G(i, jmax) = V(i, jmax);
-    }
-}
-
-void adaptUV(Solver* solver)
-{
-    int imax       = solver->grid.imax;
-    int jmax       = solver->grid.jmax;
-    double* p      = solver->p;
-    double* u      = solver->u;
-    double* v      = solver->v;
-    int* s         = solver->s;
-    double* f      = solver->f;
-    double* g      = solver->g;
-    double factorX = solver->dt / solver->grid.dx;
-    double factorY = solver->dt / solver->grid.dy;
-
-    double val = 0;
-
-    for (int j = 1; j < jmax + 1; j++) {
-        for (int i = 1; i < imax + 1; i++) {
-            U(i, j) = F(i, j) - (P(i + 1, j) - P(i, j)) * factorX;
-            V(i, j) = G(i, j) - (P(i, j + 1) - P(i, j)) * factorY;
-        }
-    }
-}
-
-double smoothRB(Solver* solver)
-{
-    int imax      = solver->grid.imax;
-    int jmax      = solver->grid.jmax;
-    double eps    = solver->eps;
-    int itermax   = solver->itermax;
-    double dx2    = solver->grid.dx * solver->grid.dx;
-    double dy2    = solver->grid.dy * solver->grid.dy;
-    double idx2   = 1.0 / dx2;
-    double idy2   = 1.0 / dy2;
-    double factor = solver->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2);
-    double* p     = solver->p;
-    double* r     = solver->r[solver->currentlevel];
-    double* rhs   = solver->rhs;
-    double epssq  = eps * eps;
-    int it        = 0;
-    double res    = 1.0;
-    int pass, jsw, isw;
-
-    jsw = 1;
-
-    for (pass = 0; pass < 2; pass++) {
-        isw = jsw;
-
-        for (int j = 1; j < jmax + 1; j++) {
-            for (int i = isw; i < imax + 1; i += 2) {
-
-                R(i, j) = RHS(i, j) -
-                          ((P(i + 1, j) - 2.0 * P(i, j) + P(i - 1, j)) * idx2 +
-                              (P(i, j + 1) - 2.0 * P(i, j) + P(i, j - 1)) * idy2);
-
-                P(i, j) -= (factor * R(i, j));
-                res += (R(i, j) * R(i, j));
-            }
-            isw = 3 - isw;
-        }
-        jsw = 3 - jsw;
-    }
-
-    res = res / (double)(imax * jmax);
-    return res;
-}
-
-void multiGrid(Solver* solver)
-{
-    double res = 0.0;
-    int imax   = solver->grid.imax;
-    int jmax   = solver->grid.jmax;
-    if (solver->currentlevel == (solver->levels - 1)) {
-        for (int i = 0; i < 5; i++) {
-            smoothRB(solver);
-        }
-        return;
-    }
-
-    for (int i = 0; i < 5; i++) {
-        smoothRB(solver);
-        if (solver->currentlevel == 0) {
-
-            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);
-
-    // restrict
-    restrictMG(solver);
-
-    coarseSolver.p   = solver->e[coarseSolver.currentlevel];
-    coarseSolver.rhs = solver->r[coarseSolver.currentlevel];
-    coarseSolver.grid.imax /= 2;
-    coarseSolver.grid.jmax /= 2;
-
-    // MGSolver on residual and error.
-    multiGrid(&coarseSolver);
-
-    // prolongate
-    prolongate(solver);
-
-    // correct p on finest level using residual
-    correct(solver);
-
-    if (solver->currentlevel == 0) {
-
-        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++) {
-        res = smoothRB(solver);
-        if (solver->currentlevel == 0) {
-
-            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
-    if (solver->currentlevel == 0) {
-        printf("Residuum: %.6f\n", res);
-    }
-#endif
-}
-
-void restrictMG(Solver* solver)
-{
-    int imax     = solver->grid.imax;
-    int jmax     = solver->grid.jmax;
-    double* r    = solver->r[solver->currentlevel + 1];
-    double* oldr = solver->r[solver->currentlevel];
-
-    for (int j = 1; j < jmax + 1; j++) {
-        for (int i = 1; i < imax + 1; ++i) {
-            R(i, j) = (oldR(2 * i - 1, 2 * j - 1) + oldR(2 * i, 2 * j - 1) * 2 +
-                          oldR(2 * i + 1, 2 * j - 1) + oldR(2 * i - 1, 2 * j) * 2 +
-                          oldR(2 * i, 2 * j) * 4 + oldR(2 * i + 1, 2 * j) * 2 +
-                          oldR(2 * i - 1, 2 * j + 1) + oldR(2 * i, 2 * j + 1) * 2 +
-                          oldR(2 * i + 1, 2 * j + 1)) /
-                      16.0;
-        }
-    }
-}
-
-void prolongate(Solver* solver)
-{
-    int imax     = solver->grid.imax;
-    int jmax     = solver->grid.jmax;
-    double* olde = solver->r[solver->currentlevel + 1];
-    double* e    = solver->r[solver->currentlevel];
-
-    for (int j = 2; j < jmax + 1; j += 2) {
-        for (int i = 2; i < imax + 1; i += 2) {
-            E(i, j) = oldE(i / 2, j / 2);
-        }
-    }
-}
-
-void correct(Solver* solver)
-{
-    int imax  = solver->grid.imax;
-    int jmax  = solver->grid.jmax;
-    double* p = solver->p;
-    double* e = solver->e[solver->currentlevel];
-    for (int j = 1; j < jmax + 1; ++j) {
-        for (int i = 1; i < imax + 1; ++i) {
-            P(i, j) += E(i, j);
-        }
-    }
-}
-
-Solver copySolver(Solver* solver)
-{
-    Solver newSolver;
-    newSolver.problem      = solver->problem;
-    newSolver.bcLeft       = solver->bcLeft;
-    newSolver.bcRight      = solver->bcRight;
-    newSolver.bcBottom     = solver->bcBottom;
-    newSolver.bcTop        = solver->bcTop;
-    newSolver.grid.imax    = solver->grid.imax;
-    newSolver.grid.jmax    = solver->grid.jmax;
-    newSolver.grid.xlength = solver->grid.xlength;
-    newSolver.grid.ylength = solver->grid.ylength;
-    newSolver.grid.dx      = solver->grid.xlength / solver->grid.imax;
-    newSolver.grid.dy      = solver->grid.ylength / solver->grid.jmax;
-    newSolver.eps          = solver->eps;
-    newSolver.omega        = solver->omega;
-    newSolver.itermax      = solver->itermax;
-    newSolver.re           = solver->re;
-    newSolver.gx           = solver->gx;
-    newSolver.gy           = solver->gy;
-    newSolver.dt           = solver->dt;
-    newSolver.te           = solver->te;
-    newSolver.tau          = solver->tau;
-    newSolver.gamma        = solver->gamma;
-    newSolver.rho          = solver->rho;
-    newSolver.levels       = solver->levels;
-    newSolver.currentlevel = solver->currentlevel + 1;
-
-    newSolver.r = solver->r;
-    newSolver.e = solver->e;
-
-    return newSolver;
-}
-
-void writeResult(Solver* solver)
-{
-    int imax  = solver->grid.imax;
-    int jmax  = solver->grid.jmax;
-    double dx = solver->grid.dx;
-    double dy = solver->grid.dy;
-    double* p = solver->p;
-    double* u = solver->u;
-    double* v = solver->v;
-    double x = 0.0, y = 0.0;
-
-    FILE* fp;
-    fp = fopen("pressure.dat", "w");
-
-    if (fp == NULL) {
-        printf("Error!\n");
-        exit(EXIT_FAILURE);
-    }
-
-    for (int j = 1; j < jmax + 1; j++) {
-        y = (double)(j - 0.5) * dy;
-        for (int i = 1; i < imax + 1; i++) {
-            x = (double)(i - 0.5) * dx;
-            fprintf(fp, "%.2f %.2f %f\n", x, y, P(i, j));
-        }
-        fprintf(fp, "\n");
-    }
-
-    fclose(fp);
-
-    fp = fopen("velocity.dat", "w");
-
-    if (fp == NULL) {
-        printf("Error!\n");
-        exit(EXIT_FAILURE);
-    }
-
-    for (int j = 1; j < jmax + 1; j++) {
-        y = dy * (j - 0.5);
-        for (int i = 1; i < imax + 1; i++) {
-            x            = dx * (i - 0.5);
-            double vel_u = (U(i, j) + U(i - 1, j)) / 2.0;
-            double vel_v = (V(i, j) + V(i, j - 1)) / 2.0;
-            double len   = sqrt((vel_u * vel_u) + (vel_v * vel_v));
-            fprintf(fp, "%.2f %.2f %f %f %f\n", x, y, vel_u, vel_v, len);
-        }
-    }
-
-    fclose(fp);
-}

EnhancedSolver/2D-seq/src/solver.h

@@ -6,62 +6,21 @@
  */
 #ifndef __SOLVER_H_
 #define __SOLVER_H_
+#include "discretization.h"
 #include "grid.h"
 #include "parameter.h"
 
-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 */
-    Grid grid;
+    Grid* grid;
-    /* arrays */
-    double *p, *rhs, **r, **e;
-    double *f, *g;
-    double *u, *v;
-    int* s;
     /* parameters */
     double eps, omega, rho;
-    double re, tau, gamma;
+    int itermax;
-    double gx, gy;
+    int levels;
-    /* time stepping */
+    double **r, **e;
-    int itermax, levels, currentlevel;
-    double dt, te;
-    double dtBound;
-    char* problem;
-    int bcLeft, bcRight, bcBottom, bcTop;
 } Solver;
 
-extern void initSolver(Solver*, Parameter*);
+extern void initSolver(Solver*, Discretization*, Parameter*);
-extern void computeRHS(Solver*);
+extern void solve(Solver*, double*, double*);
-extern double smoothRB(Solver*);
-extern void restrictMG(Solver*);
-extern void prolongate(Solver*);
-extern void correct(Solver*);
-extern Solver copySolver(Solver*);
-extern void multiGrid(Solver*);
-extern void normalizePressure(Solver*);
-extern void computeTimestep(Solver*);
-extern void setBoundaryConditions(Solver*);
-extern void setSpecialBoundaryCondition(Solver*);
-extern void setObjectBoundaryCondition(Solver*);
-extern void computeFG(Solver*);
-extern void adaptUV(Solver*);
-extern void writeResult(Solver*);
-extern void print(Solver*, double*);
-extern void printGrid(Solver*, int*);
 
 #endif

EnhancedSolver/2D-seq/src/timing.c

@@ -7,18 +7,16 @@
 #include <stdlib.h>
 #include <time.h>
 
-double getTimeStamp()
+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()
+double getTimeResolution(void)
 {
     struct timespec ts;
     clock_getres(CLOCK_MONOTONIC, &ts);
     return (double)ts.tv_sec + (double)ts.tv_nsec * 1.e-9;
 }
-
-double getTimeStamp_() { return getTimeStamp(); }

EnhancedSolver/2D-seq/src/timing.h

@@ -7,8 +7,7 @@
 #ifndef __TIMING_H_
 #define __TIMING_H_
 
-extern double getTimeStamp();
+extern double getTimeStamp(void);
-extern double getTimeResolution();
+extern double getTimeResolution(void);
-extern double getTimeStamp_();
 
 #endif // __TIMING_H_

EnhancedSolver/2D-seq/src/util.h

@@ -7,8 +7,7 @@
 #ifndef __UTIL_H_
 #define __UTIL_H_
 #define HLINE                                                                            \
-    "------------------------------------------------------------------------"           \
+    "----------------------------------------------------------------------------\n"
-    "----\n"
 
 #ifndef MIN
 #define MIN(x, y) ((x) < (y) ? (x) : (y))
@@ -20,4 +19,11 @@
 #define ABS(a) ((a) >= 0 ? (a) : -(a))
 #endif
 
+#define P(i, j)   p[(j) * (imax + 2) + (i)]
+#define F(i, j)   f[(j) * (imax + 2) + (i)]
+#define G(i, j)   g[(j) * (imax + 2) + (i)]
+#define U(i, j)   u[(j) * (imax + 2) + (i)]
+#define V(i, j)   v[(j) * (imax + 2) + (i)]
+#define RHS(i, j) rhs[(j) * (imax + 2) + (i)]
+
 #endif // __UTIL_H_

EnhancedSolver/2D-seq/vis_files/canal.plot

@@ -0,0 +1,10 @@
+unset border; unset tics; unset key;
+set term gif animate delay 30
+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