EnhancedSolver completed

EnhancedSolver/2D-mpi/Makefile

@@ -45,6 +45,7 @@ vis:
 	$(info ===>  GENERATE VISUALIZATION)
 	@gnuplot -e "filename='pressure.dat'" ./surface.plot
 	@gnuplot -e "filename='velocity.dat'" ./vector.plot
+	@gnuplot -e "filename='residual.dat'" ./residual.plot
 
 vis_clean:
 	$(info ===>  CLEAN VISUALIZATION)

EnhancedSolver/2D-mpi/backstep.par

@@ -26,8 +26,8 @@ p_init    1.0      # initial value for pressure
 
 xlength    7.0     # domain size in x-direction
 ylength    1.5	   # domain size in y-direction
-imax       210     # number of interior cells in x-direction
+imax       200     # number of interior cells in x-direction
-jmax       45	   # number of interior cells in y-direction
+jmax       40	   # number of interior cells in y-direction
 
 # Time Data:
 # ---------
@@ -42,20 +42,20 @@ tau     0.5     # safety factor for time stepsize control (<0 constant delt)
 itermax  500       # maximal number of pressure iteration in one time step
 eps      0.0001    # stopping tolerance for pressure iteration
 rho      0.52
-omg      1.8       # relaxation parameter for SOR iteration
+omg      1.7       # relaxation parameter for SOR iteration
 gamma    0.9       # upwind differencing factor gamma
 
 # Multigrid data:
 # ---------
 
 levels        3         # Multigrid levels
-presmooth     15        # Pre-smoothning iterations
+presmooth     55        # Pre-smoothning iterations
 postsmooth    5         # Post-smoothning iterations
 
 # Particle Tracing Data:
 # -----------------------
 
-numberOfParticles   200
+numberOfParticles   500
 startTime           0     #if you want to see particles trapped in recirculation zone, startTime should be set to 0 
 injectTimePeriod    1.0
 writeTimePeriod     0.5

EnhancedSolver/2D-mpi/canal.par

@@ -39,7 +39,7 @@ tau      0.5     # safety factor for time stepsize control (<0 constant delt)
 # Multigrid data:
 # ---------
 
-levels        2         # Multigrid levels
+levels        3         # Multigrid levels
 presmooth     5         # Pre-smoothning iterations
 postsmooth    5         # Post-smoothning iterations
 
@@ -48,7 +48,7 @@ postsmooth    5         # Post-smoothning iterations
 
 itermax       500       # maximal number of pressure iteration in one time step
 eps           0.00001   # stopping tolerance for pressure iteration
-omg           1.8       # relaxation parameter for SOR iteration
+omg           1.7       # relaxation parameter for SOR iteration
 gamma         0.9       # upwind differencing factor gamma
 
 # Particle Tracing Data:

EnhancedSolver/2D-mpi/dcavity.par

@@ -42,14 +42,14 @@ 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.8		# relaxation parameter for SOR iteration
+omg      1.7		# relaxation parameter for SOR iteration
 gamma    0.9		# upwind differencing factor gamma
 
 # Multigrid data:
 # ---------
 
 levels        3         # Multigrid levels
-presmooth     10        # Pre-smoothning iterations
+presmooth     20        # Pre-smoothning iterations
 postsmooth    5         # Post-smoothning iterations
 
 # Particle Tracing Data:

EnhancedSolver/2D-mpi/residual.plot

@@ -0,0 +1,9 @@
+set terminal png size 1800,768 enhanced font ,12
+set output 'residual.png'
+set datafile separator whitespace
+set xlabel "Timestep"
+set ylabel "Residual"
+
+set logscale y 2
+
+plot 'residual.dat' using 1:2 title "Residual"

EnhancedSolver/2D-mpi/src/comm-v1.c

@@ -14,7 +14,7 @@ static int sum(int* sizes, int position)
 {
     int sum = 0;
 
-    for (int i = 0; i < position; i++) {
+    for (int i = 0; i < position; i += position) {
         sum += sizes[i];
     }
 
@@ -174,6 +174,17 @@ void commPartition(Comm* c, int jmax, int imax)
 #ifdef _MPI
     c->imaxLocal = imax;
     c->jmaxLocal = sizeOfRank(c->coords[JDIM], c->size, jmax);
+
+    c->neighbours[BOTTOM] = c->rank == 0 ? -1 : c->rank - 1;
+    c->neighbours[TOP]    = c->rank == (c->size - 1) ? -1 : c->rank + 1;
+    c->neighbours[LEFT]   = -1;
+    c->neighbours[RIGHT]  = -1;
+
+    c->coords[IDIM] = 0;
+    c->coords[JDIM] = c->rank;
+
+    c->dims[IDIM] = 1;
+    c->dims[JDIM] = c->size;
 #else
     c->imaxLocal = imax;
     c->jmaxLocal = jmax;
@@ -182,15 +193,33 @@ void commPartition(Comm* c, int jmax, int imax)
 
 void commUpdateDatatypes(Comm* oldcomm, Comm* newcomm, int imaxLocal, int jmaxLocal)
 {
-    Comm* newcomm;
 
 #if defined _MPI
     newcomm->comm = MPI_COMM_NULL;
-    int result    = MPI_Comm_dup(oldcomm->comm, &newcomm->comm);
+    int result    = MPI_Comm_dup(MPI_COMM_WORLD, &newcomm->comm);
 
     if (result == MPI_ERR_COMM) {
         printf("\nNull communicator. Duplication failed !!\n");
     }
+
+    newcomm->rank = oldcomm->rank;
+    newcomm->size = oldcomm->size;
+
+    newcomm->imaxLocal = imaxLocal / 2;
+    newcomm->jmaxLocal = jmaxLocal / 2;
+
+    newcomm->neighbours[BOTTOM] = newcomm->rank == 0 ? -1 : newcomm->rank - 1;
+    newcomm->neighbours[TOP]    = newcomm->rank == (newcomm->size - 1) ? -1 : newcomm->rank + 1;
+    newcomm->neighbours[LEFT]   = -1;
+    newcomm->neighbours[RIGHT]  = -1;
+
+    newcomm->coords[IDIM] = 0;
+    newcomm->coords[JDIM] = newcomm->rank;
+
+    newcomm->dims[IDIM] = 1;
+    newcomm->dims[JDIM] = newcomm->size;
+
+
 #endif
     newcomm->imaxLocal = imaxLocal;
     newcomm->jmaxLocal = jmaxLocal;

EnhancedSolver/2D-mpi/src/comm-v2.c

@@ -4,18 +4,18 @@
  * Use of this source code is governed by a MIT style
  * license that can be found in the LICENSE file.
  */
+#include "comm.h"
 #include <stdio.h>
 #include <stdlib.h>
-
+#include <string.h>
-#include "comm.h"
 
 #ifdef _MPI
 // subroutines local to this module
-static int sum(int* sizes, int position)
+static int sum(int* sizes, int init, int offset, int coord)
 {
     int sum = 0;
 
-    for (int i = 0; i < position; i++) {
+    for (int i = init - offset; coord > 0; i -= offset, --coord) {
         sum += sizes[i];
     }
 
@@ -79,8 +79,8 @@ static void assembleResult(Comm* c, double* src, double* dst, int imax, int jmax
             int newSizes[NDIMS] = { newSizesJ[i], newSizesI[i] };
             int coords[NDIMS];
             MPI_Cart_coords(c->comm, i, NDIMS, coords);
-            int starts[NDIMS] = { sum(newSizesJ, coords[JDIM]),
+            int starts[NDIMS] = { sum(newSizesJ, i, 1, coords[JDIM]),
-                sum(newSizesI, coords[IDIM]) };
+                sum(newSizesI, i, c->dims[JDIM], coords[IDIM]) };
             printf(
                 "Rank: %d, Coords(i,j): %d %d, Size(i,j): %d %d, Target Size(i,j): %d %d "
                 "Starts(i,j): %d %d\n",
@@ -296,8 +296,15 @@ void commUpdateDatatypes(Comm* oldcomm, Comm* newcomm, int imaxLocal, int jmaxLo
         printf("\nNull communicator. Duplication failed !!\n");
     }
 
-    newcomm->imaxLocal = imaxLocal;
+    newcomm->rank = oldcomm->rank;
-    newcomm->jmaxLocal = jmaxLocal;
+    newcomm->size = oldcomm->size;
+
+    memcpy(&newcomm->neighbours, &oldcomm->neighbours, sizeof(oldcomm->neighbours));
+    memcpy(&newcomm->coords, &oldcomm->coords, sizeof(oldcomm->coords));
+    memcpy(&newcomm->dims, &oldcomm->dims, sizeof(oldcomm->dims));
+
+    newcomm->imaxLocal = imaxLocal/2;
+    newcomm->jmaxLocal = jmaxLocal/2;
 
     MPI_Datatype jBufferType;
     MPI_Type_contiguous(imaxLocal, MPI_DOUBLE, &jBufferType);

EnhancedSolver/2D-mpi/src/comm-v3.c

@@ -11,11 +11,11 @@
 
 #ifdef _MPI
 // subroutines local to this module
-static int sum(int* sizes, int position)
+static int sum(int* sizes, int init, int offset, int coord)
 {
     int sum = 0;
 
-    for (int i = 0; i < position; i++) {
+    for (int i = init - offset; coord > 0; i -= offset, --coord) {
         sum += sizes[i];
     }
 
@@ -79,8 +79,8 @@ static void assembleResult(Comm* c, double* src, double* dst, int imax, int jmax
             int newSizes[NDIMS] = { newSizesJ[i], newSizesI[i] };
             int coords[NDIMS];
             MPI_Cart_coords(c->comm, i, NDIMS, coords);
-            int starts[NDIMS] = { sum(newSizesJ, coords[JDIM]),
+            int starts[NDIMS] = { sum(newSizesJ, i, 1, coords[JDIM]),
-                sum(newSizesI, coords[IDIM]) };
+                sum(newSizesI, i, c->dims[JDIM], coords[IDIM]) };
             printf(
                 "Rank: %d, Coords(i,j): %d %d, Size(i,j): %d %d, Target Size(i,j): %d %d "
                 "Starts(i,j): %d %d\n",
@@ -277,8 +277,12 @@ void commUpdateDatatypes(Comm* oldcomm, Comm* newcomm, int imaxLocal, int jmaxLo
         printf("\nNull communicator. Duplication failed !!\n");
     }
 
-    newcomm->imaxLocal = imaxLocal;
+    newcomm->rank = oldcomm->rank;
-    newcomm->jmaxLocal = jmaxLocal;
+    newcomm->size = oldcomm->size;
+    
+
+    newcomm->imaxLocal = imaxLocal / 2;
+    newcomm->jmaxLocal = jmaxLocal / 2;
 
     MPI_Datatype jBufferType;
     MPI_Type_contiguous(imaxLocal, MPI_DOUBLE, &jBufferType);

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

@@ -140,7 +140,6 @@ void initDiscretiztion(Discretization* d, Parameter* params)
 #ifdef _MPI
     MPI_Allgather(&d->xLocal, 1, MPI_DOUBLE, xLocal, 1, MPI_DOUBLE, d->comm.comm);
     MPI_Allgather(&d->yLocal, 1, MPI_DOUBLE, yLocal, 1, MPI_DOUBLE, d->comm.comm);
-#endif
 
     d->xOffset    = sumOffset(xLocal,
         d->comm.rank,
@@ -149,6 +148,13 @@ void initDiscretiztion(Discretization* d, Parameter* params)
     d->yOffset    = sumOffset(yLocal, d->comm.rank, 1, d->comm.coords[JDIM]);
     d->xOffsetEnd = d->xOffset + d->xLocal;
     d->yOffsetEnd = d->yOffset + d->yLocal;
+#else
+
+    d->xOffset    = 0;
+    d->yOffset    = 0;
+    d->xOffsetEnd = d->xOffset + d->xLocal;
+    d->yOffsetEnd = d->yOffset + d->yLocal;
+#endif
 
     printf("Rank : %d, xOffset : %.2f, yOffset : %.2f, xOffsetEnd : %.2f, yOffsetEnd : "
            "%.2f\n",

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

@@ -50,6 +50,8 @@ int main(int argc, char** argv)
 
     commInit(&d.comm, argc, argv);
     initParameter(&p);
+    FILE* fp;
+    fp = initResidualWriter();
 
     if (argc != 2) {
         printf("Usage: %s <configFile>\n", argv[0]);
@@ -83,7 +85,7 @@ int main(int argc, char** argv)
     double tau = d.tau;
     double te  = d.te;
     double t   = 0.0;
-
+    double res = 0.0;
     timeStart = getTimeStamp();
     while (t <= te) {
         if (tau > 0.0) computeTimestep(&d);
@@ -92,9 +94,12 @@ int main(int argc, char** argv)
         setObjectBoundaryCondition(&d);
         computeFG(&d);
         computeRHS(&d);
-        solve(&s, d.p, d.rhs);
+        res = solve(&s, d.p, d.rhs);
         adaptUV(&d);
         trace(&particletracer, &d, t);
+
+        writeResidual(fp, t, res);
+
         t += d.dt;
 
 #ifdef VERBOSE
@@ -112,7 +117,7 @@ int main(int argc, char** argv)
     if (commIsMaster(s.comm)) {
         printf("Solution took %.2fs\n", timeStop - timeStart);
     }
-
+    fclose(fp);
     freeParticles(&particletracer);
     writeResults(&d);
     commFinalize(s.comm);

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

@@ -105,6 +105,7 @@ void initParticleTracer(
 
     particletracer->pointer          = 0;
     particletracer->totalParticles   = 0;
+    particletracer->removedParticles = 0;
 
     particletracer->imax = params->imax;
     particletracer->jmax = params->jmax;
@@ -199,14 +200,7 @@ void initParticleTracer(
 
     memcpy(particletracer->offset, offset, sizeof(offset));
 
-    for (int i = 0; i < particletracer->numberOfParticles; ++i) {
+    particleRandomizer(particletracer);
-        double spacing = (double)i / (double)(particletracer->numberOfParticles - 1);
-        particletracer->linSpaceLine[i].x = spacing * particletracer->x1 +
-                                            (1.0 - spacing) * particletracer->x2;
-        particletracer->linSpaceLine[i].y = spacing * particletracer->y1 +
-                                            (1.0 - spacing) * particletracer->y2;
-        particletracer->linSpaceLine[i].flag = true;
-    }
 
 #ifdef _MPI
     // Create the mpi_particle datatype
@@ -252,9 +246,15 @@ void printParticles(ParticleTracer* particletracer)
             particletracer->yOffsetEnd);
     }
 }
-void injectParticles(ParticleTracer* particletracer)
+void injectParticles(ParticleTracer* particletracer, double* s)
 {
     double x, y;
+    compress(particletracer);
+    particleRandomizer(particletracer);
+
+    int imaxLocal = particletracer->imaxLocal;
+    int jmaxLocal = particletracer->jmaxLocal;
+
     for (int i = 0; i < particletracer->numberOfParticles; ++i) {
         x = particletracer->linSpaceLine[i].x;
         y = particletracer->linSpaceLine[i].y;
@@ -263,12 +263,23 @@ void injectParticles(ParticleTracer* particletracer)
 
             particletracer->particlePool[particletracer->pointer].x = x;
             particletracer->particlePool[particletracer->pointer].y = y;
+
+            int i = particletracer->particlePool[particletracer->pointer].x /
+                    particletracer->dx;
+            int j = particletracer->particlePool[particletracer->pointer].y /
+                    particletracer->dy;
+
+            int iOffset = particletracer->xOffset / particletracer->dx,
+                jOffset = particletracer->yOffset / particletracer->dy;
+
+            if (S(i - iOffset, j - jOffset) == FLUID) {
                 particletracer->particlePool[particletracer->pointer].flag = true;
                 ++(particletracer->pointer);
                 ++(particletracer->totalParticles);
             }
         }
     }
+}
 
 void advanceParticles(ParticleTracer* particletracer,
     double* restrict u,
@@ -289,18 +300,16 @@ void advanceParticles(ParticleTracer* particletracer,
     double xlength = particletracer->xlength;
     double ylength = particletracer->ylength;
 
-    Particle buff[particletracer->size][100];
+    Particle buff[particletracer->size][(particletracer->estimatedNumParticles)];
+    memset(buff, 0, sizeof(buff));
+    Particle recvbuff[particletracer->size][(particletracer->estimatedNumParticles)];
+    memset(buff, 0, sizeof(recvbuff));
 
-    for (int i = 0; i < particletracer->size; ++i) {
+    int particleBufIndex[particletracer->size],
-        for (int j = 0; j < 100; ++j) {
+        recvparticleBufIndex[particletracer->size];
-            buff[i][j].x    = 0.0;
-            buff[i][j].y    = 0.0;
-            buff[i][j].flag = false;
-        }
-    }
-    int particleBufIndex[particletracer->size];
 
     memset(particleBufIndex, 0, sizeof(particleBufIndex));
+    memset(recvparticleBufIndex, 0, sizeof(recvparticleBufIndex));
 
     for (int i = 0; i < particletracer->totalParticles; ++i) {
         if (particletracer->particlePool[i].flag == true) {
@@ -345,12 +354,12 @@ void advanceParticles(ParticleTracer* particletracer,
             particletracer->particlePool[i].y = new_y;
 
             if (((new_x < particletracer->xOffset) ||
-                    (new_x >= particletracer->xOffsetEnd) ||
+                    (new_x > particletracer->xOffsetEnd) ||
                     (new_y < particletracer->yOffset) ||
-                    (new_y >= particletracer->yOffsetEnd))) {
+                    (new_y > particletracer->yOffsetEnd))) {
                 // New logic to transfer particles to neighbouring ranks or discard the
                 // particle.
-
+#ifdef _MPI
                 for (int i = 0; i < particletracer->size; ++i) {
                     if ((new_x >=
                             particletracer->offset[i + particletracer->size * XOFFSET]) &&
@@ -367,6 +376,7 @@ void advanceParticles(ParticleTracer* particletracer,
                         ++particleBufIndex[i];
                     }
                 }
+#endif
                 particletracer->particlePool[i].flag = false;
                 particletracer->removedParticles++;
             }
@@ -384,14 +394,43 @@ void advanceParticles(ParticleTracer* particletracer,
 
 #ifdef _MPI
     for (int i = 0; i < particletracer->size; ++i) {
-        if (i != particletracer->rank) {
+        if (i != comm->rank) {
-            MPI_Send(buff[i], 100, particletracer->mpi_particle, i, 0, comm->comm);
+            MPI_Send(&particleBufIndex[i], 1, MPI_INT, i, 0, comm->comm);
         }
     }
     for (int i = 0; i < particletracer->size; ++i) {
         if (i != particletracer->rank) {
-            MPI_Recv(buff[i],
+            MPI_Recv(&recvparticleBufIndex[i],
-                100,
+                1,
+                MPI_INT,
+                i,
+                0,
+                comm->comm,
+                MPI_STATUS_IGNORE);
+
+            // if (0 !=recvparticleBufIndex[i]) {
+            //     printf("Rank %d will receive %d particles from rank %d\n",
+            //         particletracer->rank,
+            //         recvparticleBufIndex[i],
+            //         i);
+            // }
+        }
+    }
+
+    for (int i = 0; i < particletracer->size; ++i) {
+        if (i != particletracer->rank) {
+            MPI_Send(buff[i],
+                particleBufIndex[i],
+                particletracer->mpi_particle,
+                i,
+                0,
+                comm->comm);
+        }
+    }
+    for (int i = 0; i < particletracer->size; ++i) {
+        if (i != particletracer->rank) {
+            MPI_Recv(recvbuff[i],
+                recvparticleBufIndex[i],
                 particletracer->mpi_particle,
                 i,
                 0,
@@ -401,11 +440,10 @@ void advanceParticles(ParticleTracer* particletracer,
     }
     for (int i = 0; i < particletracer->size; ++i) {
         if (i != particletracer->rank) {
-            for (int j = 0; j < 100; ++j) {
+            for (int j = 0; j < recvparticleBufIndex[i]; ++j) {
-                if (buff[i][j].flag == true) {
+                particletracer->particlePool[particletracer->pointer].x = recvbuff[i][j]
-                    particletracer->particlePool[particletracer->pointer].x = buff[i][j]
                                                                               .x;
-                    particletracer->particlePool[particletracer->pointer].y = buff[i][j]
+                particletracer->particlePool[particletracer->pointer].y = recvbuff[i][j]
                                                                               .y;
                 particletracer->particlePool[particletracer->pointer].flag = true;
                 ++(particletracer->pointer);
@@ -413,7 +451,6 @@ void advanceParticles(ParticleTracer* particletracer,
             }
         }
     }
-    }
 #endif
 }
 
@@ -427,7 +464,7 @@ void freeParticles(ParticleTracer* particletracer)
 void writeParticles(ParticleTracer* particletracer, Comm* comm)
 {
     int collectedBuffIndex[particletracer->size];
-
+    compress(particletracer);
 #ifdef _MPI
 
     MPI_Gather(&particletracer->totalParticles,
@@ -480,13 +517,12 @@ void writeParticles(ParticleTracer* particletracer, Comm* comm)
         // fprintf(fp, "CYCLE 1 1 int\n");
         // fprintf(fp, "1\n");
 
+#ifdef _MPI
         int overallTotalParticles = sum(collectedBuffIndex, particletracer->size);
 
         // fprintf(fp, "POINTS %d float\n", overallTotalParticles);
 
         // printf("Total particles : %d\n", overallTotalParticles);
-
-#ifdef _MPI
         for (int i = 1; i < particletracer->size; ++i) {
             Particle recvBuff[collectedBuffIndex[i]];
             MPI_Recv(&recvBuff,
@@ -501,10 +537,17 @@ void writeParticles(ParticleTracer* particletracer, Comm* comm)
                 double x = recvBuff[j].x;
                 double y = recvBuff[j].y;
                 fprintf(fp, "%f %f\n", x, y);
+
                 // printf("Rank : 0 receiving from rank %d X : %.2f, Y : %.2f with totalpt
                 // : %d\n", i, x, y, particletracer->totalParticles);
             }
         }
+#else
+        int overallTotalParticles = particletracer->totalParticles;
+
+        // fprintf(fp, "POINTS %d float\n", overallTotalParticles);
+
+        // printf("Total particles : %d\n", overallTotalParticles);
 #endif
         for (int i = 0; i < particletracer->totalParticles; ++i) {
             double x = particletracer->particlePool[i].x;
@@ -568,7 +611,7 @@ void trace(ParticleTracer* particletracer, Discretization* d, double time)
 {
     if (time >= particletracer->startTime) {
         if ((time - particletracer->lastInjectTime) >= particletracer->injectTimePeriod) {
-            injectParticles(particletracer);
+            injectParticles(particletracer, d->grid.s);
             particletracer->lastInjectTime = time;
         }
         if ((time - particletracer->lastWriteTime) >= particletracer->writeTimePeriod) {
@@ -591,7 +634,7 @@ void compress(ParticleTracer* particletracer)
     Particle tempPool[particletracer->totalParticles];
     int totalParticles = 0;
 
-    printf("Performing compression ...");
+    // printf("\nPerforming compression ...");
 
     for (int i = 0; i < particletracer->totalParticles; ++i) {
         if (memPool[i].flag == true) {
@@ -602,11 +645,31 @@ void compress(ParticleTracer* particletracer)
         }
     }
 
-    printf(" remove %d particles\n", particletracer->totalParticles - totalParticles);
+    // printf(" remove %d particles\n", particletracer->totalParticles - totalParticles);
 
     particletracer->totalParticles   = totalParticles;
     particletracer->removedParticles = 0;
-    particletracer->pointer          = totalParticles + 1;
+    particletracer->pointer          = totalParticles;
 
     memcpy(particletracer->particlePool, tempPool, totalParticles * sizeof(Particle));
 }
+
+void particleRandomizer(ParticleTracer* particletracer)
+{
+    memset(particletracer->linSpaceLine,
+        0,
+        sizeof(Particle) * particletracer->numberOfParticles);
+
+    for (int i = 0; i < particletracer->numberOfParticles; ++i) {
+
+        particletracer->linSpaceLine[i].x = (((double)rand() / RAND_MAX) *
+                                                (particletracer->x2 -
+                                                    particletracer->x1)) +
+                                            particletracer->x1;
+        particletracer->linSpaceLine[i].y = (((double)rand() / RAND_MAX) *
+                                                (particletracer->y2 -
+                                                    particletracer->y1)) +
+                                            particletracer->y1;
+        particletracer->linSpaceLine[i].flag = true;
+    }
+}

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

@@ -52,7 +52,7 @@ typedef struct {
 } ParticleTracer;
 
 extern void initParticleTracer(ParticleTracer*, Discretization*, Parameter*);
-extern void injectParticles(ParticleTracer*);
+extern void injectParticles(ParticleTracer*, double*);
 extern void advanceParticles(ParticleTracer*, double*, double*, double*, Comm*, double);
 extern void freeParticles(ParticleTracer*);
 extern void writeParticles(ParticleTracer*, Comm*);
@@ -60,5 +60,6 @@ extern void printParticleTracerParameters(ParticleTracer*);
 extern void printParticles(ParticleTracer*);
 extern void trace(ParticleTracer*, Discretization*, double);
 extern void compress(ParticleTracer*);
+extern void particleRandomizer(ParticleTracer*);
 
 #endif

EnhancedSolver/2D-mpi/src/progress.c

@@ -6,8 +6,8 @@
  */
 #include <math.h>
 #include <mpi.h>
-#include <stdio.h>
 #include <string.h>
+#include <stdlib.h>
 
 #include "progress.h"
 
@@ -49,3 +49,22 @@ void stopProgress()
     printf("\n");
     fflush(stdout);
 }
+
+FILE* initResidualWriter()
+{
+    FILE* fp;
+    fp = fopen("residual.dat", "w");
+
+    if (fp == NULL) {
+        printf("Error!\n");
+        exit(-1);
+    }
+
+    return fp;
+
+}
+
+void writeResidual(FILE* fp, double ts, double res)
+{
+    fprintf(fp, "%f, %f\n", ts, res);
+}

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

@@ -4,11 +4,14 @@
  * Use of this source code is governed by a MIT-style
  * license that can be found in the LICENSE file.
  */
+#include <stdio.h>
+
 #ifndef __PROGRESS_H_
 #define __PROGRESS_H_
 
 extern void initProgress(double);
 extern void printProgress(double);
 extern void stopProgress();
-
+extern FILE* initResidualWriter(void);
+extern void writeResidual(FILE*, double, double);
 #endif

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

@@ -312,7 +312,7 @@ void initSolver(Solver* s, Discretization* d, Parameter* p)
     }
 }
 
-void solve(Solver* s, double* p, double* rhs)
+double solve(Solver* s, double* p, double* rhs)
 {
     double res = multiGrid(s, p, rhs, 0, s->comm);
 
@@ -321,4 +321,5 @@ void solve(Solver* s, double* p, double* rhs)
         printf("Residuum: %.6f\n", res);
     }
 #endif
+return res;
 }

EnhancedSolver/2D-mpi/src/solver-rb.c

@@ -48,7 +48,7 @@ void initSolver(Solver* s, Discretization* d, Parameter* p)
     s->comm    = &d->comm;
 }
 
-void solve(Solver* s, double* p, double* rhs)
+double solve(Solver* s, double* p, double* rhs)
 {
     int imax      = s->grid->imax;
     int jmax      = s->grid->jmax;
@@ -126,4 +126,6 @@ void solve(Solver* s, double* p, double* rhs)
         printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
     }
 #endif
+
+return res;
 }

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

@@ -24,7 +24,7 @@ typedef struct {
     Comm* comm;
 } Solver;
 
-void initSolver(Solver*, Discretization*, Parameter*);
+extern void initSolver(Solver*, Discretization*, Parameter*);
-void solve(Solver*, double*, double*);
+extern double solve(Solver*, double*, double*);
-void printSolver(Solver* , double*, char*);
+extern void printSolver(Solver* , double*, char*);
 #endif

EnhancedSolver/2D-mpi/vis_files/backstep_animate.plot

@@ -8,7 +8,7 @@ set size ratio -1
 set object 1 rect from 0.0,0.0 to 1.0,0.5 lw 5
 
 
-do for [ts=0:300] {
+do for [ts=0:550] {
     plot "particles_".ts.".dat" with points pointtype 7
 }
 unset output

EnhancedSolver/2D-seq/Makefile

@@ -38,9 +38,9 @@ $(BUILD_DIR)/%.s:  %.c
 	$(info ===>  GENERATE ASM  $@)
 	$(CC) -S $(CPPFLAGS) $(CFLAGS) $< -o $@
 
-.PHONY: clean distclean tags info asm format
+.PHONY: clean distclean vis vis_clean tags info asm format
 
-clean:
+clean: vis_clean
 	$(info ===>  CLEAN)
 	@rm -rf $(BUILD_DIR)
 	@rm -f tags
@@ -74,6 +74,7 @@ vis:
 	$(info ===>  GENERATE VISUALIZATION)
 	@gnuplot -e "filename='pressure.dat'" ./surface.plot
 	@gnuplot -e "filename='velocity.dat'" ./vector.plot
+	@gnuplot -e "filename='residual.dat'" ./residual.plot
 
 vis_clean:
 	$(info ===>  CLEAN VISUALIZATION)

EnhancedSolver/2D-seq/backstep.par

@@ -15,7 +15,7 @@ bcRight    3			#
 gx    0.0      # Body forces (e.g. gravity)
 gy    0.0      #
 
-re    20000.0	   # Reynolds number
+re    36000.0	   # Reynolds number
 
 u_init    1.0      # initial value for velocity in x-direction
 v_init    0.0      # initial value for velocity in y-direction
@@ -56,7 +56,7 @@ postsmooth    5         # Post-smoothning iterations
 # -----------------------
 
 numberOfParticles   200
-startTime           10
+startTime           0
 injectTimePeriod    1.0
 writeTimePeriod     0.5
 

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       256      # number of interior cells in x-direction
+imax       200      # number of interior cells in x-direction
-jmax       64	   # number of interior cells in y-direction
+jmax       40	    # number of interior cells in y-direction
 
 # Time Data:
 # ---------
 
-te      80.0   # final time
+te      60.0   # final time
 dt      0.02    # time stepsize
 tau     0.5     # safety factor for time stepsize control (<0 constant delt)
 

EnhancedSolver/2D-seq/dcavity.par

@@ -49,7 +49,7 @@ gamma    0.9		# upwind differencing factor gamma
 # ---------
 
 levels        3         # Multigrid levels
-presmooth     10        # Pre-smoothning iterations
+presmooth     20        # Pre-smoothning iterations
 postsmooth    5         # Post-smoothning iterations
 
 # Particle Tracing Data:

EnhancedSolver/2D-seq/residual.plot

@@ -0,0 +1,9 @@
+set terminal png size 1800,768 enhanced font ,12
+set output 'residual.png'
+set datafile separator whitespace
+set xlabel "Timestep"
+set ylabel "Residual"
+
+set logscale y 2
+
+plot 'residual.dat' using 1:2 title "Residual"

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

@@ -24,6 +24,8 @@ int main(int argc, char** argv)
     ParticleTracer particletracer;
 
     initParameter(&p);
+    FILE* fp;
+    fp = initResidualWriter();
 
     if (argc != 2) {
         printf("Usage: %s <configFile>\n", argv[0]);
@@ -45,6 +47,7 @@ int main(int argc, char** argv)
     double te  = d.te;
     double t   = 0.0;
     int nt     = 0;
+    double res = 0.0;
 
     timeStart  = getTimeStamp();
 
@@ -57,10 +60,12 @@ int main(int argc, char** argv)
         computeFG(&d);
         computeRHS(&d);
         if (nt % 100 == 0) normalizePressure(&d);
-        solve(&s, d.p, d.rhs);
+        res = solve(&s, d.p, d.rhs);
         adaptUV(&d);
         trace(&particletracer, d.u, d.v, d.dt, t);
 
+        writeResidual(fp, t, res);
+
         t += d.dt;
         nt++;
 
@@ -73,6 +78,7 @@ int main(int argc, char** argv)
 
     timeStop = getTimeStamp();
 
+    fclose(fp);
     stopProgress();
     freeParticles(&particletracer);
     printf("Solution took %.2fs\n", timeStop - timeStart);

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

@@ -26,6 +26,8 @@ void initParameter(Parameter* param)
     param->tau        = 0.5;
     param->rho        = 0.99;
     param->levels     = 5;
+    param->presmooth  = 5;
+    param->postsmooth = 5;
 }
 
 void readParameter(Parameter* param, const char* filename)
@@ -64,6 +66,8 @@ void readParameter(Parameter* param, const char* filename)
             PARSE_INT(jmax);
             PARSE_INT(itermax);
             PARSE_INT(levels);
+            PARSE_INT(presmooth);
+            PARSE_INT(postsmooth);
             PARSE_REAL(eps);
             PARSE_REAL(omg);
             PARSE_REAL(re);

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

@@ -10,7 +10,7 @@
 typedef struct {
     double xlength, ylength;
     int imax, jmax;
-    int itermax, levels;
+    int itermax, levels, presmooth, postsmooth;
     double eps, omg, rho;
     double re, tau, gamma;
     double te, dt;

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

@@ -110,7 +110,7 @@ static void compress(ParticleTracer* p)
     Particle tempPool[p->totalParticles];
     int totalParticles = 0;
 
-    printf("Performing compression ...");
+    // printf("Performing compression ...");
 
     for (int i = 0; i < p->totalParticles; i++) {
         if (memPool[i].flag == 1) {
@@ -121,7 +121,7 @@ static void compress(ParticleTracer* p)
         }
     }
 
-    printf(" remove %d particles\n", p->totalParticles - totalParticles);
+    // printf(" remove %d particles\n", p->totalParticles - totalParticles);
     p->totalParticles   = totalParticles;
     p->removedParticles = 0;
     p->pointer          = totalParticles + 1;
@@ -131,6 +131,8 @@ static void compress(ParticleTracer* p)
 void writeParticles(ParticleTracer* p)
 {
     static int ts = 0;
+    compress(p);
+    
     VtkOptions opts = { .particletracer = p };
 
     char filename[50];

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

@@ -49,3 +49,21 @@ void stopProgress()
     printf("\n");
     fflush(stdout);
 }
+FILE* initResidualWriter()
+{
+    FILE* fp;
+    fp = fopen("residual.dat", "w");
+
+    if (fp == NULL) {
+        printf("Error!\n");
+        exit(EXIT_FAILURE);
+    }
+
+    return fp;
+
+}
+
+void writeResidual(FILE* fp, double ts, double res)
+{
+    fprintf(fp, "%f, %f\n", ts, res);
+}

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

@@ -10,5 +10,6 @@
 extern void initProgress(double);
 extern void printProgress(double);
 extern void stopProgress(void);
-
+extern FILE* initResidualWriter(void);
+extern void writeResidual(FILE*, double, double);
 #endif

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

@@ -71,7 +71,38 @@ static void setBoundaryCondition(double* p, int imax, int jmax)
     }
 }
 
-static double smooth(Solver* s, double* p, double* rhs, int level, int imax, int jmax)
+static void 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) {
+
+                P(i, j) -= factor * (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));
+
+            }
+            isw = 3 - isw;
+        }
+        jsw = 3 - jsw;
+    }
+
+}
+
+static double calculateResidual(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;
@@ -94,7 +125,6 @@ static double smooth(Solver* s, double* p, double* rhs, int level, int imax, int
                           ((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;
@@ -110,7 +140,7 @@ static double multiGrid(Solver* s, double* p, double* rhs, int level, int imax,
 {
     double res = 0.0;
 
-    // coarsest level TODO: Use direct solver?
+    // coarsest level
     if (level == COARSEST_LEVEL) {
         for (int i = 0; i < 5; i++) {
             smooth(s, p, rhs, level, imax, jmax);
@@ -118,17 +148,18 @@ static double multiGrid(Solver* s, double* p, double* rhs, int level, int imax,
         return res;
     }
 
-    // pre-smoothing TODO: Make smoothing steps configurable?
+    // pre-smoothing
-    for (int i = 0; i < 5; i++) {
+    for (int i = 0; i < s->presmooth; i++) {
         smooth(s, p, rhs, level, imax, jmax);
         if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax);
     }
 
+    res = calculateResidual(s, p, rhs, level, 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
@@ -139,8 +170,8 @@ static double multiGrid(Solver* s, double* p, double* rhs, int level, int imax,
     if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax);
 
     // post-smoothing
-    for (int i = 0; i < 5; i++) {
+    for (int i = 0; i < s->postsmooth; i++) {
-        res = smooth(s, p, rhs, level, imax, jmax);
+        smooth(s, p, rhs, level, imax, jmax);
         if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax);
     }
 
@@ -154,6 +185,8 @@ void initSolver(Solver* s, Discretization* d, Parameter* p)
     s->itermax = p->itermax;
     s->levels  = p->levels;
     s->grid    = &d->grid;
+    s->presmooth = p->presmooth;
+    s->postsmooth = p->postsmooth;
 
     int imax   = s->grid->imax;
     int jmax   = s->grid->jmax;
@@ -176,11 +209,13 @@ void initSolver(Solver* s, Discretization* d, Parameter* p)
     }
 }
 
-void solve(Solver* s, double* p, double* rhs)
+double 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
+
+return res;
 }

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

@@ -18,18 +18,9 @@ void initSolver(Solver* s, Discretization* d, Parameter* p)
     Grid* g  = s->grid;
     int imax = s->grid->imax;
     int jmax = s->grid->jmax;
-
-    s->totalFluidCells = 0;
-    for (int j = 0; j < jmax + 2; j++) {
-        for (int i = 0; i < imax + 2; i++) {
-            if (gridIsFluid(g, i, j)) {
-                s->totalFluidCells++;
-            }
-        }
-    }
 }
 
-void solve(Solver* s, double* p, double* rhs)
+double solve(Solver* s, double* p, double* rhs)
 {
     int imax      = s->grid->imax;
     int jmax      = s->grid->jmax;
@@ -55,16 +46,14 @@ void solve(Solver* s, double* p, double* rhs)
 
             for (int j = 1; j < jmax + 1; j++) {
                 for (int i = isw; i < imax + 1; i += 2) {
-                    if (gridIsFluid(g, i, j)) {
+
                     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)) *
+                                   (P(i, j + 1) - 2.0 * P(i, j) + P(i, j - 1)) * idy2);
-                                           idy2);
 
                     P(i, j) -= (factor * r);
                     res += (r * r);
                 }
-                }
                 isw = 3 - isw;
             }
             jsw = 3 - jsw;
@@ -80,7 +69,7 @@ void solve(Solver* s, double* p, double* rhs)
             P(imax + 1, j) = P(imax, j);
         }
 
-        res = res / (double)(s->totalFluidCells);
+        res = res / (double)(imax * jmax);
 #ifdef DEBUG
         printf("%d Residuum: %e\n", it, res);
 #endif
@@ -90,4 +79,6 @@ void solve(Solver* s, double* p, double* rhs)
 #ifdef VERBOSE
     printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
 #endif
+
+    return res;
 }

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

@@ -18,18 +18,9 @@ void initSolver(Solver* s, Discretization* d, Parameter* p)
     Grid* g  = s->grid;
     int imax = s->grid->imax;
     int jmax = s->grid->jmax;
-
-    s->totalFluidCells = 0;
-    for (int j = 0; j < jmax + 2; j++) {
-        for (int i = 0; i < imax + 2; i++) {
-            if (gridIsFluid(g, i, j)) {
-                s->totalFluidCells++;
-            }
-        }
-    }
 }
 
-void solve(Solver* s, double* p, double* rhs)
+double solve(Solver* s, double* p, double* rhs)
 {
     int imax      = s->grid->imax;
     int jmax      = s->grid->jmax;
@@ -50,7 +41,6 @@ void solve(Solver* s, double* p, double* rhs)
 
         for (int j = 1; j < jmax + 1; j++) {
             for (int i = 1; i < imax + 1; i++) {
-                if (gridIsFluid(g, i, j)) {
                 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);
@@ -59,7 +49,6 @@ void solve(Solver* s, double* p, double* rhs)
                 res += (r * r);
             }
         }
-        }
 
         for (int i = 1; i < imax + 1; i++) {
             P(i, 0)        = P(i, 1);
@@ -71,7 +60,7 @@ void solve(Solver* s, double* p, double* rhs)
             P(imax + 1, j) = P(imax, j);
         }
 
-        res = res / (double)(s->totalFluidCells);
+        res = res / (double)(imax * jmax);
 #ifdef DEBUG
         printf("%d Residuum: %e\n", it, res);
 #endif
@@ -81,4 +70,6 @@ void solve(Solver* s, double* p, double* rhs)
 #ifdef VERBOSE
     printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
 #endif
+
+    return res;
 }

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

@@ -16,12 +16,12 @@ typedef struct {
     /* parameters */
     double eps, omega, rho;
     int itermax;
-    int levels;
+    int levels, presmooth, postsmooth;
     int totalFluidCells;
     double **r, **e;
 } Solver;
 
 extern void initSolver(Solver*, Discretization*, Parameter*);
-extern void solve(Solver*, double*, double*);
+extern double solve(Solver*, double*, double*);
 
 #endif

EnhancedSolver/3D-mpi/Makefile

@@ -0,0 +1,96 @@
+#=======================================================================================
+# 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.
+#=======================================================================================
+
+#CONFIGURE BUILD SYSTEM
+TARGET	   = exe-$(TAG)
+BUILD_DIR  = ./$(TAG)
+SRC_DIR    = ./src
+MAKE_DIR   = ./
+Q         ?= @
+
+#DO NOT EDIT BELOW
+include $(MAKE_DIR)/config.mk
+include $(MAKE_DIR)/include_$(TAG).mk
+INCLUDES  += -I$(SRC_DIR) -I$(BUILD_DIR)
+
+VPATH     = $(SRC_DIR)
+SRC       = $(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)/vtkWriter-$(VTK_OUTPUT_FMT).o
+OBJ      += $(BUILD_DIR)/solver-$(SOLVER).o
+SOURCES   = $(SRC) $(wildcard $(SRC_DIR)/*.h)
+ifeq ($(VTK_OUTPUT_FMT),mpi)
+DEFINES  += -D_VTK_WRITER_MPI
+endif
+
+CPPFLAGS := $(CPPFLAGS) $(DEFINES) $(OPTIONS) $(INCLUDES)
+
+${TARGET}: sanity-checks $(BUILD_DIR) $(OBJ)
+	$(info ===>  LINKING  $(TARGET))
+	$(Q)${LINKER} ${LFLAGS} -o $(TARGET) $(OBJ) $(LIBS)
+
+$(BUILD_DIR)/%.o:  %.c $(MAKE_DIR)/include_$(TAG).mk $(MAKE_DIR)/config.mk
+	$(info ===>  COMPILE  $@)
+	$(CC) -c $(CPPFLAGS) $(CFLAGS) $< -o $@
+	$(Q)$(GCC) $(CPPFLAGS) -MT $(@:.d=.o) -MM  $< > $(BUILD_DIR)/$*.d
+
+$(BUILD_DIR)/%.s:  %.c
+	$(info ===>  GENERATE ASM  $@)
+	$(CC) -S $(CPPFLAGS) $(CFLAGS) $< -o $@
+
+.PHONY: clean distclean vis vis_clean tags info asm format
+
+vis:
+	$(info ===>  GENERATE VISUALIZATION)
+	@gnuplot -e "filename='residual.dat'" ./residual.plot
+
+vis_clean:
+	$(info ===>  CLEAN VISUALIZATION)
+	@rm -f *.dat
+	@rm -f *.vtk
+	@rm -f vtk_files/*.vtk
+	@rm -f *.png
+
+clean: vis_clean
+	$(info ===>  CLEAN)
+	@rm -rf $(BUILD_DIR)
+	@rm -f tags
+
+distclean: clean
+	$(info ===>  DIST CLEAN)
+	@rm -f $(TARGET)
+
+info:
+	$(info $(CFLAGS))
+	$(Q)$(CC) $(VERSION)
+
+asm:  $(BUILD_DIR) $(ASM)
+
+tags:
+	$(info ===>  GENERATE TAGS)
+	$(Q)ctags -R
+
+format:
+	@for src in $(SOURCES) ; do \
+		echo "Formatting $$src" ; \
+		clang-format -i $$src ; \
+	done
+	@echo "Done"
+
+sanity-checks:
+ifeq ($(VTK_OUTPUT_FMT),mpi)
+ifeq ($(ENABLE_MPI),false)
+	$(error VTK_OUTPUT_FMT mpi only supported for ENABLE_MPI true!)
+endif
+endif
+
+
+$(BUILD_DIR):
+	@mkdir $(BUILD_DIR)
+
+-include $(OBJ:.o=.d)

EnhancedSolver/3D-mpi/README.md

@@ -0,0 +1,78 @@
+# C source skeleton
+
+## Build
+
+1. Configure the toolchain and additional options in `config.mk`:
+```
+# Supported: GCC, CLANG, ICC
+TAG ?= GCC
+ENABLE_OPENMP ?= false
+
+OPTIONS +=  -DARRAY_ALIGNMENT=64
+#OPTIONS +=  -DVERBOSE
+#OPTIONS +=  -DVERBOSE_AFFINITY
+#OPTIONS +=  -DVERBOSE_DATASIZE
+#OPTIONS +=  -DVERBOSE_TIMER
+```
+
+The verbosity options enable detailed output about solver, affinity settings, allocation sizes and timer resolution.
+For debugging you may want to set the VERBOSE option:
+```
+# Supported: GCC, CLANG, ICC
+TAG ?= GCC
+ENABLE_OPENMP ?= false
+
+OPTIONS +=  -DARRAY_ALIGNMENT=64
+OPTIONS +=  -DVERBOSE
+#OPTIONS +=  -DVERBOSE_AFFINITY
+#OPTIONS +=  -DVERBOSE_DATASIZE
+#OPTIONS +=  -DVERBOSE_TIMER
+`
+
+2. Build with:
+```
+make
+```
+
+You can build multiple toolchains in the same directory, but notice that the Makefile is only acting on the one currently set.
+Intermediate build results are located in the `<TOOLCHAIN>` directory.
+
+To output the executed commands use:
+```
+make Q=
+```
+
+3. Clean up with:
+```
+make clean
+```
+to clean intermediate build results.
+
+```
+make distclean
+```
+to clean intermediate build results and binary.
+
+4. (Optional) Generate assembler:
+```
+make asm
+```
+The assembler files will also be located in the `<TOOLCHAIN>` directory.
+
+## Usage
+
+You have to provide a parameter file describing the problem you want to solve:
+```
+./exe-CLANG  dcavity.par
+```
+
+Examples are given in in dcavity (a lid driven cavity test case) and canal (simulating a empty canal).
+
+You can plot the resulting velocity and pressure fields using gnuplot:
+```
+gnuplot vector.plot
+```
+and for the pressure:
+```
+gnuplot surface.plot
+```

EnhancedSolver/3D-mpi/backstep.par

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

EnhancedSolver/3D-mpi/canal.par

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

EnhancedSolver/3D-mpi/config.mk

@@ -0,0 +1,18 @@
+# Supported: GCC, CLANG, ICC
+TAG ?= ICC
+# Supported: true, false
+ENABLE_MPI ?= true
+ENABLE_OPENMP ?= false
+# Supported: rb, mg
+SOLVER ?= mg
+# Supported: seq, mpi
+VTK_OUTPUT_FMT ?= seq
+
+#Feature options
+OPTIONS +=  -DARRAY_ALIGNMENT=64
+OPTIONS +=  -DVERBOSE
+#OPTIONS +=  -DDEBUG
+#OPTIONS +=  -DBOUNDCHECK
+#OPTIONS +=  -DVERBOSE_AFFINITY
+#OPTIONS +=  -DVERBOSE_DATASIZE
+#OPTIONS +=  -DVERBOSE_TIMER

EnhancedSolver/3D-mpi/dcavity.par

@@ -0,0 +1,88 @@
+#==============================================================================
+#                              Driven Cavity
+#==============================================================================
+
+# Problem specific Data:
+# ---------------------
+
+name dcavity        # name of flow setup
+
+bcLeft    1			#  flags for boundary conditions
+bcRight   1			#  1 = no-slip      3 = outflow
+bcBottom  1			#  2 = free-slip    4 = periodic
+bcTop     1			#
+bcFront   1			#
+bcBack    1			#
+
+gx    0.0			# Body forces (e.g. gravity)
+gy    0.0			#
+gz    0.0			#
+
+re    1000.0		    # Reynolds number
+
+u_init    0.0		# initial value for velocity in x-direction
+v_init    0.0		# initial value for velocity in y-direction
+w_init    0.0		# initial value for velocity in z-direction
+p_init    0.0		# initial value for pressure
+
+# Geometry Data:
+# -------------
+
+xlength    1.0		# domain size in x-direction
+ylength    1.0		# domain size in y-direction
+zlength    1.0		# domain size in z-direction
+imax       80		# number of interior cells in x-direction
+jmax       80		# number of interior cells in y-direction
+kmax       80		# number of interior cells in z-direction
+
+# Time Data:
+# ---------
+
+te       10.0		# final time
+dt       0.02	    # time stepsize
+tau      0.5		# safety factor for time stepsize control (<0 constant delt)
+
+# Pressure Iteration Data:
+# -----------------------
+
+itermax  1000		# maximal number of pressure iteration in one time step
+eps      0.001		# stopping tolerance for pressure iteration
+rho      0.5
+omg      1.7		# relaxation parameter for SOR iteration
+gamma    0.9		# upwind differencing factor gamma
+
+# Multigrid data:
+# ---------
+
+levels        2         # Multigrid levels
+presmooth     20        # Pre-smoothning iterations
+postsmooth    5         # Post-smoothning iterations
+
+# Particle Tracing Data:
+# -----------------------
+
+numberOfParticles   200
+startTime           0
+injectTimePeriod    1.0
+writeTimePeriod     0.1
+
+x1                  0.0
+y1                  0.1
+z1                  0.0
+x2                  1.0
+y2                  0.1
+z2                  1.0
+
+# Obstacle Geometry Data:
+# -----------------------
+# Shape 0 disable, 1 Rectangle/Square, 2 Circle
+
+shape               0 
+xCenter             0.5
+yCenter             0.5
+zCenter             0.5
+xRectLength         0.25
+yRectLength         0.25
+zRectLength         0.25
+circleRadius        0.25
+#===============================================================================

EnhancedSolver/3D-mpi/include_CLANG.mk

@@ -0,0 +1,21 @@
+ifeq ($(ENABLE_MPI),true)
+CC   = mpicc
+DEFINES  = -D_MPI
+else
+CC   = cc
+endif
+
+GCC  = cc
+LINKER = $(CC)
+
+ifeq ($(ENABLE_OPENMP),true)
+OPENMP   = -fopenmp
+#OPENMP   = -Xpreprocessor -fopenmp #required on Macos with homebrew libomp
+LIBS     = # -lomp
+endif
+
+VERSION  = --version
+CFLAGS   = -Ofast -std=c17
+LFLAGS   = $(OPENMP) -lm
+DEFINES  += -D_GNU_SOURCE# -DDEBUG
+INCLUDES = -I/opt/homebrew/include

EnhancedSolver/3D-mpi/include_GCC.mk

@@ -0,0 +1,20 @@
+ifeq ($(ENABLE_MPI),true)
+CC   = mpicc
+DEFINES  = -D_MPI
+else
+CC   = gcc
+endif
+
+GCC  = gcc
+LINKER = $(CC)
+
+ifeq ($(ENABLE_OPENMP),true)
+OPENMP   = -fopenmp
+endif
+
+VERSION  = --version
+CFLAGS   = -Ofast -ffreestanding -std=c99 $(OPENMP)
+LFLAGS   = $(OPENMP)
+DEFINES  += -D_GNU_SOURCE
+INCLUDES =
+LIBS     =

EnhancedSolver/3D-mpi/include_ICC.mk

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

EnhancedSolver/3D-mpi/karman.par

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

EnhancedSolver/3D-mpi/residual.plot

@@ -0,0 +1,9 @@
+set terminal png size 1800,768 enhanced font ,12
+set output 'residual.png'
+set datafile separator whitespace
+set xlabel "Timestep"
+set ylabel "Residual"
+
+set logscale y 2
+
+plot 'residual.dat' using 1:2 title "Residual"

EnhancedSolver/3D-mpi/src/allocate.c

@@ -0,0 +1,38 @@
+/*
+ * 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.
+ */
+#include <errno.h>
+#include <stddef.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "allocate.h"
+
+void* allocate(size_t alignment, size_t bytesize)
+{
+    int errorCode;
+    void* ptr;
+
+    errorCode = posix_memalign(&ptr, alignment, bytesize);
+
+    if (errorCode) {
+        if (errorCode == EINVAL) {
+            fprintf(stderr, "Error: Alignment parameter is not a power of two\n");
+            exit(EXIT_FAILURE);
+        }
+        if (errorCode == ENOMEM) {
+            fprintf(stderr, "Error: Insufficient memory to fulfill the request\n");
+            exit(EXIT_FAILURE);
+        }
+    }
+
+    if (ptr == NULL) {
+        fprintf(stderr, "Error: posix_memalign failed!\n");
+        exit(EXIT_FAILURE);
+    }
+
+    return ptr;
+}

EnhancedSolver/3D-mpi/src/allocate.h

@@ -0,0 +1,13 @@
+/*
+ * 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 __ALLOCATE_H_
+#define __ALLOCATE_H_
+#include <stdlib.h>
+
+extern void* allocate(size_t alignment, size_t bytesize);
+
+#endif

EnhancedSolver/3D-mpi/src/comm.c

@@ -0,0 +1,647 @@
+/*
+ * 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 <stddef.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "allocate.h"
+#include "comm.h"
+
+#if defined(_MPI)
+// subroutines local to this module
+static int sizeOfRank(int rank, int size, int N)
+{
+    return N / size + ((N % size > rank) ? 1 : 0);
+}
+
+static void setupCommunication(Comm* c, Direction direction, int layer)
+{
+    int imaxLocal = c->imaxLocal;
+    int jmaxLocal = c->jmaxLocal;
+    int kmaxLocal = c->kmaxLocal;
+
+    size_t dblsize = sizeof(double);
+    int sizes[NDIMS];
+    int subSizes[NDIMS];
+    int starts[NDIMS];
+    int offset = 0;
+
+    sizes[IDIM] = imaxLocal + 2;
+    sizes[JDIM] = jmaxLocal + 2;
+    sizes[KDIM] = kmaxLocal + 2;
+
+    if (layer == HALO) {
+        offset = 1;
+    }
+
+    switch (direction) {
+    case LEFT:
+        subSizes[IDIM] = 1;
+        subSizes[JDIM] = jmaxLocal;
+        subSizes[KDIM] = kmaxLocal;
+        starts[IDIM]   = 1 - offset;
+        starts[JDIM]   = 1;
+        starts[KDIM]   = 1;
+        break;
+    case RIGHT:
+        subSizes[IDIM] = 1;
+        subSizes[JDIM] = jmaxLocal;
+        subSizes[KDIM] = kmaxLocal;
+        starts[IDIM]   = imaxLocal + offset;
+        starts[JDIM]   = 1;
+        starts[KDIM]   = 1;
+        break;
+    case BOTTOM:
+        subSizes[IDIM] = imaxLocal;
+        subSizes[JDIM] = 1;
+        subSizes[KDIM] = kmaxLocal;
+        starts[IDIM]   = 1;
+        starts[JDIM]   = 1 - offset;
+        starts[KDIM]   = 1;
+        break;
+    case TOP:
+        subSizes[IDIM] = imaxLocal;
+        subSizes[JDIM] = 1;
+        subSizes[KDIM] = kmaxLocal;
+        starts[IDIM]   = 1;
+        starts[JDIM]   = jmaxLocal + offset;
+        starts[KDIM]   = 1;
+        break;
+    case FRONT:
+        subSizes[IDIM] = imaxLocal;
+        subSizes[JDIM] = jmaxLocal;
+        subSizes[KDIM] = 1;
+        starts[IDIM]   = 1;
+        starts[JDIM]   = 1;
+        starts[KDIM]   = 1 - offset;
+        break;
+    case BACK:
+        subSizes[IDIM] = imaxLocal;
+        subSizes[JDIM] = jmaxLocal;
+        subSizes[KDIM] = 1;
+        starts[IDIM]   = 1;
+        starts[JDIM]   = 1;
+        starts[KDIM]   = kmaxLocal + offset;
+        break;
+    case NDIRS:
+        printf("ERROR!\n");
+        break;
+    }
+
+    if (layer == HALO) {
+        MPI_Type_create_subarray(NDIMS,
+            sizes,
+            subSizes,
+            starts,
+            MPI_ORDER_C,
+            MPI_DOUBLE,
+            &c->rbufferTypes[direction]);
+        MPI_Type_commit(&c->rbufferTypes[direction]);
+    } else if (layer == BULK) {
+        MPI_Type_create_subarray(NDIMS,
+            sizes,
+            subSizes,
+            starts,
+            MPI_ORDER_C,
+            MPI_DOUBLE,
+            &c->sbufferTypes[direction]);
+        MPI_Type_commit(&c->sbufferTypes[direction]);
+    }
+}
+
+static void assembleResult(Comm* c,
+    double* src,
+    double* dst,
+    int imaxLocal[],
+    int jmaxLocal[],
+    int kmaxLocal[],
+    int offset[],
+    int kmax,
+    int jmax,
+    int imax)
+{
+    int numRequests = 1;
+
+    if (c->rank == 0) {
+        numRequests = c->size + 1;
+    }
+
+    MPI_Request requests[numRequests];
+
+    /* all ranks send their interpolated bulk array */
+    MPI_Isend(src,
+        c->imaxLocal * c->jmaxLocal * c->kmaxLocal,
+        MPI_DOUBLE,
+        0,
+        0,
+        c->comm,
+        &requests[0]);
+
+    /* rank 0 assembles the subdomains */
+    if (c->rank == 0) {
+        for (int i = 0; i < c->size; i++) {
+            MPI_Datatype domainType;
+            int oldSizes[NDIMS] = { kmax, jmax, imax };
+            int newSizes[NDIMS] = { kmaxLocal[i], jmaxLocal[i], imaxLocal[i] };
+            int starts[NDIMS]   = { offset[i * NDIMS + KDIM],
+                  offset[i * NDIMS + JDIM],
+                  offset[i * NDIMS + IDIM] };
+            MPI_Type_create_subarray(NDIMS,
+                oldSizes,
+                newSizes,
+                starts,
+                MPI_ORDER_C,
+                MPI_DOUBLE,
+                &domainType);
+            MPI_Type_commit(&domainType);
+
+            MPI_Irecv(dst, 1, domainType, i, 0, c->comm, &requests[i + 1]);
+            MPI_Type_free(&domainType);
+        }
+    }
+
+    MPI_Waitall(numRequests, requests, MPI_STATUSES_IGNORE);
+}
+
+// subroutines local to this module
+static int sum(int* sizes, int init, int offset, int coord)
+{
+    int sum = 0;
+
+    for (int i = init - offset; coord > 0; i -= offset, --coord) {
+        sum += sizes[i];
+    }
+
+    return sum;
+}
+#endif // defined _MPI
+
+// exported subroutines
+void commReduction(double* v, int op)
+{
+#if defined(_MPI)
+    if (op == MAX) {
+        MPI_Allreduce(MPI_IN_PLACE, v, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
+    } else if (op == SUM) {
+        MPI_Allreduce(MPI_IN_PLACE, v, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+    }
+#endif
+}
+
+int commIsBoundary(Comm* c, Direction direction)
+{
+#if defined(_MPI)
+    switch (direction) {
+    case LEFT:
+        return c->coords[ICORD] == 0;
+        break;
+    case RIGHT:
+        return c->coords[ICORD] == (c->dims[ICORD] - 1);
+        break;
+    case BOTTOM:
+        return c->coords[JCORD] == 0;
+        break;
+    case TOP:
+        return c->coords[JCORD] == (c->dims[JCORD] - 1);
+        break;
+    case FRONT:
+        return c->coords[KCORD] == 0;
+        break;
+    case BACK:
+        return c->coords[KCORD] == (c->dims[KCORD] - 1);
+        break;
+    case NDIRS:
+        printf("ERROR!\n");
+        break;
+    }
+#endif
+
+    return 1;
+}
+
+void commExchange(Comm* c, double* grid)
+{
+#if defined(_MPI)
+    int counts[6]      = { 1, 1, 1, 1, 1, 1 };
+    MPI_Aint displs[6] = { 0, 0, 0, 0, 0, 0 };
+
+    MPI_Neighbor_alltoallw(grid,
+        counts,
+        displs,
+        c->sbufferTypes,
+        grid,
+        counts,
+        displs,
+        c->rbufferTypes,
+        c->comm);
+#endif
+}
+
+void commShift(Comm* c, double* f, double* g, double* h)
+{
+#if defined(_MPI)
+    MPI_Request requests[6] = { MPI_REQUEST_NULL,
+        MPI_REQUEST_NULL,
+        MPI_REQUEST_NULL,
+        MPI_REQUEST_NULL,
+        MPI_REQUEST_NULL,
+        MPI_REQUEST_NULL };
+
+    /* shift G */
+    /* receive ghost cells from bottom neighbor */
+    MPI_Irecv(g,
+        1,
+        c->rbufferTypes[BOTTOM],
+        c->neighbours[BOTTOM],
+        0,
+        c->comm,
+        &requests[0]);
+
+    /* send ghost cells to top neighbor */
+    MPI_Isend(g, 1, c->sbufferTypes[TOP], c->neighbours[TOP], 0, c->comm, &requests[1]);
+
+    /* shift F */
+    /* receive ghost cells from left neighbor */
+    MPI_Irecv(f, 1, c->rbufferTypes[LEFT], c->neighbours[LEFT], 1, c->comm, &requests[2]);
+
+    /* send ghost cells to right neighbor */
+    MPI_Isend(f,
+        1,
+        c->sbufferTypes[RIGHT],
+        c->neighbours[RIGHT],
+        1,
+        c->comm,
+        &requests[3]);
+
+    /* shift H */
+    /* receive ghost cells from front neighbor */
+    MPI_Irecv(h,
+        1,
+        c->rbufferTypes[FRONT],
+        c->neighbours[FRONT],
+        2,
+        c->comm,
+        &requests[4]);
+
+    /* send ghost cells to back neighbor */
+    MPI_Isend(h, 1, c->sbufferTypes[BACK], c->neighbours[BACK], 2, c->comm, &requests[5]);
+
+    MPI_Waitall(6, requests, MPI_STATUSES_IGNORE);
+#endif
+}
+
+void commGetOffsets(Comm* c, int offsets[], int kmax, int jmax, int imax)
+{
+#if defined(_MPI)
+    int sum = 0;
+
+    for (int i = 0; i < c->coords[ICORD]; i++) {
+        sum += sizeOfRank(i, c->dims[ICORD], imax);
+    }
+    offsets[IDIM] = sum;
+    sum           = 0;
+
+    for (int i = 0; i < c->coords[JCORD]; i++) {
+        sum += sizeOfRank(i, c->dims[JCORD], jmax);
+    }
+    offsets[JDIM] = sum;
+    sum           = 0;
+
+    for (int i = 0; i < c->coords[KCORD]; i++) {
+        sum += sizeOfRank(i, c->dims[KCORD], kmax);
+    }
+    offsets[KDIM] = sum;
+#endif
+}
+
+#define G(v, i, j, k)                                                                    \
+    v[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)]
+
+void commCollectResult(Comm* c,
+    double* ug,
+    double* vg,
+    double* wg,
+    double* pg,
+    double* u,
+    double* v,
+    double* w,
+    double* p,
+    int kmax,
+    int jmax,
+    int imax)
+{
+    int imaxLocal = c->imaxLocal;
+    int jmaxLocal = c->jmaxLocal;
+    int kmaxLocal = c->kmaxLocal;
+#if defined(_MPI)
+    int offset[c->size * NDIMS];
+    int imaxLocalAll[c->size];
+    int jmaxLocalAll[c->size];
+    int kmaxLocalAll[c->size];
+
+    MPI_Gather(&imaxLocal, 1, MPI_INT, imaxLocalAll, 1, MPI_INT, 0, MPI_COMM_WORLD);
+    MPI_Gather(&jmaxLocal, 1, MPI_INT, jmaxLocalAll, 1, MPI_INT, 0, MPI_COMM_WORLD);
+    MPI_Gather(&kmaxLocal, 1, MPI_INT, kmaxLocalAll, 1, MPI_INT, 0, MPI_COMM_WORLD);
+
+    if (c->rank == 0) {
+        for (int i = 0; i < c->size; i++) {
+            int coords[NCORDS];
+            MPI_Cart_coords(c->comm, i, NDIMS, coords);
+            offset[i * NDIMS + IDIM] = sum(imaxLocalAll,
+                i,
+                c->dims[IDIM] * c->dims[JDIM],
+                coords[ICORD]);
+            offset[i * NDIMS + JDIM] = sum(jmaxLocalAll, i, c->dims[IDIM], coords[JCORD]);
+            offset[i * NDIMS + KDIM] = sum(kmaxLocalAll, i, 1, coords[KCORD]);
+
+            printf("Rank: %d, Coords(k,j,i): %d %d %d, Size(k,j,i): %d %d %d, "
+                   "Offset(k,j,i): %d %d %d\n",
+                i,
+                coords[KCORD],
+                coords[JCORD],
+                coords[ICORD],
+                kmaxLocalAll[i],
+                jmaxLocalAll[i],
+                imaxLocalAll[i],
+                offset[i * NDIMS + KDIM],
+                offset[i * NDIMS + JDIM],
+                offset[i * NDIMS + IDIM]);
+        }
+    }
+
+    size_t bytesize = imaxLocal * jmaxLocal * kmaxLocal * sizeof(double);
+    double* tmp     = allocate(64, bytesize);
+    int idx         = 0;
+
+    /* collect P */
+    for (int k = 1; k < kmaxLocal + 1; k++) {
+        for (int j = 1; j < jmaxLocal + 1; j++) {
+            for (int i = 1; i < imaxLocal + 1; i++) {
+                tmp[idx++] = G(p, i, j, k);
+            }
+        }
+    }
+
+    assembleResult(c,
+        tmp,
+        pg,
+        imaxLocalAll,
+        jmaxLocalAll,
+        kmaxLocalAll,
+        offset,
+        kmax,
+        jmax,
+        imax);
+
+    /* collect U */
+    idx = 0;
+
+    for (int k = 1; k < kmaxLocal + 1; k++) {
+        for (int j = 1; j < jmaxLocal + 1; j++) {
+            for (int i = 1; i < imaxLocal + 1; i++) {
+                tmp[idx++] = (G(u, i, j, k) + G(u, i - 1, j, k)) / 2.0;
+            }
+        }
+    }
+
+    assembleResult(c,
+        tmp,
+        ug,
+        imaxLocalAll,
+        jmaxLocalAll,
+        kmaxLocalAll,
+        offset,
+        kmax,
+        jmax,
+        imax);
+
+    /* collect V */
+    idx = 0;
+
+    for (int k = 1; k < kmaxLocal + 1; k++) {
+        for (int j = 1; j < jmaxLocal + 1; j++) {
+            for (int i = 1; i < imaxLocal + 1; i++) {
+                tmp[idx++] = (G(v, i, j, k) + G(v, i, j - 1, k)) / 2.0;
+            }
+        }
+    }
+
+    assembleResult(c,
+        tmp,
+        vg,
+        imaxLocalAll,
+        jmaxLocalAll,
+        kmaxLocalAll,
+        offset,
+        kmax,
+        jmax,
+        imax);
+
+    /* collect W */
+    idx = 0;
+
+    for (int k = 1; k < kmaxLocal + 1; k++) {
+        for (int j = 1; j < jmaxLocal + 1; j++) {
+            for (int i = 1; i < imaxLocal + 1; i++) {
+                tmp[idx++] = (G(w, i, j, k) + G(w, i, j, k - 1)) / 2.0;
+            }
+        }
+    }
+
+    assembleResult(c,
+        tmp,
+        wg,
+        imaxLocalAll,
+        jmaxLocalAll,
+        kmaxLocalAll,
+        offset,
+        kmax,
+        jmax,
+        imax);
+
+    free(tmp);
+#else
+    int idx = 0;
+
+    for (int k = 1; k < kmaxLocal + 1; k++) {
+        for (int j = 1; j < jmaxLocal + 1; j++) {
+            for (int i = 1; i < imaxLocal + 1; i++) {
+                pg[idx++] = G(p, i, j, k);
+            }
+        }
+    }
+
+    idx = 0;
+
+    for (int k = 1; k < kmaxLocal + 1; k++) {
+        for (int j = 1; j < jmaxLocal + 1; j++) {
+            for (int i = 1; i < imaxLocal + 1; i++) {
+                ug[idx++] = (G(u, i, j, k) + G(u, i - 1, j, k)) / 2.0;
+            }
+        }
+    }
+
+    idx = 0;
+
+    for (int k = 1; k < kmaxLocal + 1; k++) {
+        for (int j = 1; j < jmaxLocal + 1; j++) {
+            for (int i = 1; i < imaxLocal + 1; i++) {
+                vg[idx++] = (G(v, i, j, k) + G(v, i, j - 1, k)) / 2.0;
+            }
+        }
+    }
+
+    idx = 0;
+
+    for (int k = 1; k < kmaxLocal + 1; k++) {
+        for (int j = 1; j < jmaxLocal + 1; j++) {
+            for (int i = 1; i < imaxLocal + 1; i++) {
+                wg[idx++] = (G(w, i, j, k) + G(w, i, j, k - 1)) / 2.0;
+            }
+        }
+    }
+#endif
+}
+
+void commPrintConfig(Comm* c)
+{
+#if defined(_MPI)
+    fflush(stdout);
+    MPI_Barrier(MPI_COMM_WORLD);
+    if (commIsMaster(c)) {
+        printf("Communication setup:\n");
+    }
+
+    for (int i = 0; i < c->size; i++) {
+        if (i == c->rank) {
+            printf("\tRank %d of %d\n", c->rank, c->size);
+            printf("\tNeighbours (front, back, bottom, top, left, right): %d, %d, %d, "
+                   "%d, %d, %d\n",
+                c->neighbours[FRONT],
+                c->neighbours[BACK],
+                c->neighbours[BOTTOM],
+                c->neighbours[TOP],
+                c->neighbours[LEFT],
+                c->neighbours[RIGHT]);
+            printf("\tCoordinates (k,j,i) %d %d %d\n",
+                c->coords[KCORD],
+                c->coords[JCORD],
+                c->coords[ICORD]);
+            printf("\tLocal domain size (k,j,i) %dx%dx%d\n",
+                c->kmaxLocal,
+                c->jmaxLocal,
+                c->imaxLocal);
+            fflush(stdout);
+        }
+    }
+    MPI_Barrier(MPI_COMM_WORLD);
+#endif
+}
+
+void commInit(Comm* c, int argc, char** argv)
+{
+#if defined(_MPI)
+    MPI_Init(&argc, &argv);
+    MPI_Comm_rank(MPI_COMM_WORLD, &(c->rank));
+    MPI_Comm_size(MPI_COMM_WORLD, &(c->size));
+#else
+    c->rank = 0;
+    c->size = 1;
+#endif
+}
+
+void commPartition(Comm* c, int kmax, int jmax, int imax)
+{
+#if defined(_MPI)
+    int dims[NDIMS]    = { 0, 0, 0 };
+    int periods[NDIMS] = { 0, 0, 0 };
+    MPI_Dims_create(c->size, NDIMS, dims);
+    MPI_Cart_create(MPI_COMM_WORLD, NCORDS, dims, periods, 0, &c->comm);
+    MPI_Cart_shift(c->comm, ICORD, 1, &c->neighbours[LEFT], &c->neighbours[RIGHT]);
+    MPI_Cart_shift(c->comm, JCORD, 1, &c->neighbours[BOTTOM], &c->neighbours[TOP]);
+    MPI_Cart_shift(c->comm, KCORD, 1, &c->neighbours[FRONT], &c->neighbours[BACK]);
+    MPI_Cart_get(c->comm, NCORDS, c->dims, periods, c->coords);
+
+    c->imaxLocal = sizeOfRank(c->coords[KDIM], dims[ICORD], imax);
+    c->jmaxLocal = sizeOfRank(c->coords[JDIM], dims[JCORD], jmax);
+    c->kmaxLocal = sizeOfRank(c->coords[IDIM], dims[KCORD], kmax);
+
+    // setup buffer types for communication
+    setupCommunication(c, LEFT, BULK);
+    setupCommunication(c, LEFT, HALO);
+    setupCommunication(c, RIGHT, BULK);
+    setupCommunication(c, RIGHT, HALO);
+    setupCommunication(c, BOTTOM, BULK);
+    setupCommunication(c, BOTTOM, HALO);
+    setupCommunication(c, TOP, BULK);
+    setupCommunication(c, TOP, HALO);
+    setupCommunication(c, FRONT, BULK);
+    setupCommunication(c, FRONT, HALO);
+    setupCommunication(c, BACK, BULK);
+    setupCommunication(c, BACK, HALO);
+#else
+    c->imaxLocal = imax;
+    c->jmaxLocal = jmax;
+    c->kmaxLocal = kmax;
+#endif
+}
+
+void commFinalize(Comm* c)
+{
+#if defined(_MPI)
+    for (int i = 0; i < NDIRS; i++) {
+        MPI_Type_free(&c->sbufferTypes[i]);
+        MPI_Type_free(&c->rbufferTypes[i]);
+    }
+
+    MPI_Finalize();
+#endif
+}
+
+void commUpdateDatatypes(
+    Comm* oldcomm, Comm* newcomm, int imaxLocal, int jmaxLocal, int kmaxLocal)
+{
+#if defined _MPI
+
+    int result = MPI_Comm_dup(oldcomm->comm, &newcomm->comm);
+
+    if (result == MPI_ERR_COMM) {
+        printf("\nNull communicator. Duplication failed !!\n");
+    }
+
+    newcomm->rank = oldcomm->rank;
+    newcomm->size = oldcomm->size;
+
+    newcomm->imaxLocal = imaxLocal / 2;
+    newcomm->jmaxLocal = jmaxLocal / 2;
+    newcomm->kmaxLocal = kmaxLocal / 2;
+
+    setupCommunication(newcomm, LEFT, BULK);
+    setupCommunication(newcomm, LEFT, HALO);
+    setupCommunication(newcomm, RIGHT, BULK);
+    setupCommunication(newcomm, RIGHT, HALO);
+    setupCommunication(newcomm, BOTTOM, BULK);
+    setupCommunication(newcomm, BOTTOM, HALO);
+    setupCommunication(newcomm, TOP, BULK);
+    setupCommunication(newcomm, TOP, HALO);
+    setupCommunication(newcomm, FRONT, BULK);
+    setupCommunication(newcomm, FRONT, HALO);
+    setupCommunication(newcomm, BACK, BULK);
+    setupCommunication(newcomm, BACK, HALO);
+#else
+    newcomm->imaxLocal = imaxLocal;
+    newcomm->jmaxLocal = jmaxLocal;
+    newcomm->kmaxLocal = kmaxLocal;
+#endif
+}
+
+void commFreeCommunicator(Comm* comm)
+{
+#ifdef _MPI
+    MPI_Comm_free(&comm->comm);
+#endif
+}

EnhancedSolver/3D-mpi/src/comm.h

@@ -0,0 +1,66 @@
+/*
+ * 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 __COMM_H_
+#define __COMM_H_
+#if defined(_MPI)
+#include <mpi.h>
+#endif
+/*
+ * Spatial directions:
+ * ICORD (0) from 0 (LEFT) to imax (RIGHT)
+ * JCORD (1) from 0 (BOTTOM) to jmax (TOP)
+ * KCORD (2) from 0 (FRONT) to kmax (BACK)
+ * All derived Subarray types are in C ordering
+ * with indices KDIM (0), JDIM(1) and IDIM(2)
+ * */
+typedef enum direction { LEFT = 0, RIGHT, BOTTOM, TOP, FRONT, BACK, NDIRS } Direction;
+typedef enum coordinates { ICORD = 0, JCORD, KCORD, NCORDS } Coordinates;
+typedef enum dimension { KDIM = 0, JDIM, IDIM, NDIMS } Dimension;
+enum layer { HALO = 0, BULK };
+enum op { MAX = 0, SUM };
+
+
+typedef struct {
+    int rank;
+    int size;
+#if defined(_MPI)
+    MPI_Comm comm;
+    MPI_Datatype sbufferTypes[NDIRS];
+    MPI_Datatype rbufferTypes[NDIRS];
+#endif
+    int neighbours[NDIRS];
+    int coords[NDIMS], dims[NDIMS];
+    int imaxLocal, jmaxLocal, kmaxLocal;
+} Comm;
+
+extern void commInit(Comm* c, int argc, char** argv);
+extern void commPartition(Comm* c, int kmax, int jmax, int imax);
+extern void commFinalize(Comm* comm);
+extern void commPrintConfig(Comm*);
+extern void commExchange(Comm*, double*);
+extern void commShift(Comm* c, double* f, double* g, double* h);
+extern void commReduction(double* v, int op);
+extern int commIsBoundary(Comm* c, Direction direction);
+extern void commGetOffsets(Comm* c, int offsets[], int kmax, int jmax, int imax);
+extern void commFreeCommunicator(Comm* comm);
+extern void commUpdateDatatypes(
+    Comm* oldcomm, Comm* newcomm, int imaxLocal, int jmaxLocal, int kmaxLocal);
+extern void commCollectResult(Comm* c,
+    double* ug,
+    double* vg,
+    double* wg,
+    double* pg,
+    double* u,
+    double* v,
+    double* w,
+    double* p,
+    int kmax,
+    int jmax,
+    int imax);
+
+static inline int commIsMaster(Comm* c) { return c->rank == 0; }
+#endif // __COMM_H_

EnhancedSolver/3D-mpi/src/discretization.h

@@ -0,0 +1,92 @@
+/*
+ * 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 "comm.h"
+#include "grid.h"
+#include "parameter.h"
+#include "util.h"
+
+enum OBJECTBOUNDARY {
+    FLUID = 0,
+    /* Front Corners */
+    FRONTTOPLEFTCORNER,
+    FRONTTOPRIGHTCORNER,
+    FRONTBOTTOMLEFTCORNER,
+    FRONTBOTTOMRIGHTCORNER,
+    /* Back Corners */
+    BACKTOPLEFTCORNER,
+    BACKTOPRIGHTCORNER,
+    BACKBOTTOMLEFTCORNER,
+    BACKBOTTOMRIGHTCORNER,
+    /* Faces */
+    FRONTFACE,
+    BACKFACE,
+    LEFTFACE,
+    RIGHTFACE,
+    TOPFACE,
+    BOTTOMFACE,
+    /* Front Lines remaining after Corners and Faces */
+    FRONTLEFTLINE,
+    FRONTRIGHTLINE,
+    FRONTTOPLINE,
+    FRONTBOTTOMLINE,
+    /* Bottom Lines remaining after Corners and Faces */
+    BACKLEFTLINE,
+    BACKRIGHTLINE,
+    BACKTOPLINE,
+    BACKBOTTOMLINE,
+    /* Mid Lines remaining after Corners and Faces */
+    MIDTOPLEFTLINE,
+    MIDTOPRIGHTLINE,
+    MIDBOTTOMLEFTLINE,
+    MIDBOTTOMRIGHTLINE,
+    /* Local where its an object but not a boundary */
+    OBSTACLE,
+    /*Ghost cells boundary */
+    OUTSIDEBOUNDARY
+};
+
+enum BC { NOSLIP = 1, SLIP, OUTFLOW, PERIODIC };
+
+enum SHAPE { NOSHAPE = 0, RECT, CIRCLE };
+
+typedef struct {
+    /* geometry and grid information */
+    Grid grid;
+    /* arrays */
+    double *p, *rhs;
+    double *f, *g, *h;
+    double *u, *v, *w;
+    double* s;
+    /* parameters */
+    double eps, omega;
+    double re, tau, gamma;
+    double gx, gy, gz;
+    /* time stepping */
+    int itermax;
+    double dt, te;
+    double dtBound;
+    char* problem;
+    double xLocal, yLocal, zLocal, xOffset, yOffset, zOffset, xOffsetEnd, yOffsetEnd,
+        zOffsetEnd;
+    int bcLeft, bcRight, bcBottom, bcTop, bcFront, bcBack;
+    Comm comm;
+} Discretization;
+
+extern void initDiscretization(Discretization*, Parameter*);
+extern void computeRHS(Discretization*);
+extern void normalizePressure(Discretization*);
+extern void computeTimestep(Discretization*);
+extern void setBoundaryConditions(Discretization*);
+extern void setObjectBoundaryCondition(Discretization*);
+extern void setSpecialBoundaryCondition(Discretization*);
+extern void computeFG(Discretization*);
+extern void adaptUV(Discretization*);
+extern void printGrid(Discretization*, double*);
+#endif

EnhancedSolver/3D-mpi/src/grid.h

@@ -0,0 +1,16 @@
+/*
+ * 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 __GRID_H_
+#define __GRID_H_
+
+typedef struct {
+    double dx, dy, dz;
+    int imax, jmax, kmax;
+    double xlength, ylength, zlength;
+} Grid;
+
+#endif // __GRID_H_

EnhancedSolver/3D-mpi/src/likwid-marker.h

@@ -0,0 +1,54 @@
+/*
+ * =======================================================================================
+ *
+ *      Author:   Jan Eitzinger (je), jan.eitzinger@fau.de
+ *      Copyright (c) 2020 RRZE, University Erlangen-Nuremberg
+ *
+ *      Permission is hereby granted, free of charge, to any person obtaining a copy
+ *      of this software and associated documentation files (the "Software"), to deal
+ *      in the Software without restriction, including without limitation the rights
+ *      to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ *      copies of the Software, and to permit persons to whom the Software is
+ *      furnished to do so, subject to the following conditions:
+ *
+ *      The above copyright notice and this permission notice shall be included in all
+ *      copies or substantial portions of the Software.
+ *
+ *      THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ *      IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ *      FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ *      AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ *      LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ *      OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+ *      SOFTWARE.
+ *
+ * =======================================================================================
+ */
+#ifndef LIKWID_MARKERS_H
+#define LIKWID_MARKERS_H
+
+#ifdef LIKWID_PERFMON
+#include <likwid.h>
+#define LIKWID_MARKER_INIT                likwid_markerInit()
+#define LIKWID_MARKER_THREADINIT          likwid_markerThreadInit()
+#define LIKWID_MARKER_SWITCH              likwid_markerNextGroup()
+#define LIKWID_MARKER_REGISTER(regionTag) likwid_markerRegisterRegion(regionTag)
+#define LIKWID_MARKER_START(regionTag)    likwid_markerStartRegion(regionTag)
+#define LIKWID_MARKER_STOP(regionTag)     likwid_markerStopRegion(regionTag)
+#define LIKWID_MARKER_CLOSE               likwid_markerClose()
+#define LIKWID_MARKER_RESET(regionTag)    likwid_markerResetRegion(regionTag)
+#define LIKWID_MARKER_GET(regionTag, nevents, events, time, count)                       \
+    likwid_markerGetRegion(regionTag, nevents, events, time, count)
+#else /* LIKWID_PERFMON */
+#define LIKWID_MARKER_INIT
+#define LIKWID_MARKER_THREADINIT
+#define LIKWID_MARKER_SWITCH
+#define LIKWID_MARKER_REGISTER(regionTag)
+#define LIKWID_MARKER_START(regionTag)
+#define LIKWID_MARKER_STOP(regionTag)
+#define LIKWID_MARKER_CLOSE
+#define LIKWID_MARKER_GET(regionTag, nevents, events, time, count)
+#define LIKWID_MARKER_RESET(regionTag)
+#endif /* LIKWID_PERFMON */
+
+#endif /*LIKWID_MARKERS_H*/

EnhancedSolver/3D-mpi/src/main.c

@@ -0,0 +1,151 @@
+/*
+ * 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.
+ */
+#include <stdio.h>
+#include <stdlib.h>
+#include <unistd.h>
+
+#include "allocate.h"
+#include "discretization.h"
+#include "parameter.h"
+#include "particletracing.h"
+#include "progress.h"
+#include "solver.h"
+#include "timing.h"
+#include "vtkWriter.h"
+
+int main(int argc, char** argv)
+{
+    double timeStart, timeStop;
+    Parameter p;
+    Solver s;
+    Discretization d;
+    ParticleTracer particletracer;
+
+    commInit(&d.comm, argc, argv);
+    initParameter(&p);
+    FILE* fp;
+    if (commIsMaster(&d.comm)) fp = initResidualWriter();
+
+    if (argc != 2) {
+        printf("Usage: %s <configFile>\n", argv[0]);
+        exit(EXIT_SUCCESS);
+    }
+
+    readParameter(&p, argv[1]);
+    commPartition(&d.comm, p.kmax, p.jmax, p.imax);
+
+    if (commIsMaster(&d.comm)) {
+        printParameter(&p);
+    }
+
+    initDiscretization(&d, &p);
+    initSolver(&s, &d, &p);
+    initParticleTracer(&particletracer, &p, &d);
+
+
+#ifndef VERBOSE
+    initProgress(d.te);
+#endif
+
+    double tau = d.tau;
+    double te  = d.te;
+    double t   = 0.0;
+    int nt     = 0;
+    double res = 0.0;
+
+    timeStart = getTimeStamp();
+    while (t <= te) {
+        if (tau > 0.0) computeTimestep(&d);
+        setBoundaryConditions(&d);
+        setSpecialBoundaryCondition(&d);
+        setObjectBoundaryCondition(&d);
+        computeFG(&d);
+        computeRHS(&d);
+
+        if (nt % 100 == 0) normalizePressure(&d);
+        res = solve(&s, d.p, d.rhs);
+        adaptUV(&d);
+        trace(&particletracer, d.u, d.v, d.w, d.s, t);
+
+        if (commIsMaster(&d.comm)) writeResidual(fp, t, res);
+
+        t += d.dt;
+        nt++;
+
+#ifdef VERBOSE
+        if (commIsMaster(s.comm)) {
+            printf("TIME %f , TIMESTEP %f\n", t, d.dt);
+        }
+#else
+        printProgress(t);
+#endif
+    }
+    timeStop = getTimeStamp();
+
+#ifndef VERBOSE
+    stopProgress();
+#endif
+
+    if (commIsMaster(s.comm)) {
+        printf("Solution took %.2fs\n", timeStop - timeStart);
+    }
+
+    timeStart = getTimeStamp();
+#ifdef _VTK_WRITER_MPI
+    VtkOptions opts = { .grid = s.grid, .comm = s.comm };
+    vtkOpen(&opts, s.problem);
+    vtkScalar(&opts, "pressure", d.p);
+    vtkVector(&opts, "velocity", (VtkVector) { d.u, d.v, d.w });
+    vtkClose(&opts);
+#else
+
+    if (commIsMaster(&d.comm)) fclose(fp);
+    freeParticles(&particletracer);
+
+    double *pg, *ug, *vg, *wg;
+
+    if (commIsMaster(s.comm)) {
+        size_t bytesize = s.grid->imax * s.grid->jmax * s.grid->kmax * sizeof(double);
+
+        pg = allocate(64, bytesize);
+        ug = allocate(64, bytesize);
+        vg = allocate(64, bytesize);
+        wg = allocate(64, bytesize);
+    }
+
+    commCollectResult(s.comm,
+        ug,
+        vg,
+        wg,
+        pg,
+        d.u,
+        d.v,
+        d.w,
+        d.p,
+        s.grid->kmax,
+        s.grid->jmax,
+        s.grid->imax);
+
+    if (commIsMaster(s.comm)) {
+        VtkOptions opts = { .grid = s.grid };
+        vtkOpen(&opts, s.problem);
+        vtkScalar(&opts, "pressure", pg);
+        vtkVector(&opts, "velocity", (VtkVector) { ug, vg, wg });
+        vtkClose(&opts);
+    }
+
+#endif
+
+    timeStop = getTimeStamp();
+
+    if (commIsMaster(s.comm)) {
+        printf("Result output took %.2fs\n", timeStop - timeStart);
+    }
+
+    commFinalize(s.comm);
+    return EXIT_SUCCESS;
+}

EnhancedSolver/3D-mpi/src/parameter.c

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

EnhancedSolver/3D-mpi/src/parameter.h

@@ -0,0 +1,35 @@
+/*
+ * 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 __PARAMETER_H_
+#define __PARAMETER_H_
+
+typedef struct {
+    int imax, jmax, kmax;
+    double xlength, ylength, zlength;
+    int itermax;
+    double eps, omg;
+    double re, tau, gamma;
+    double te, dt;
+    double gx, gy, gz;
+    char* name;
+    int bcLeft, bcRight, bcBottom, bcTop, bcFront, bcBack;
+    double u_init, v_init, w_init, p_init;
+    int levels, presmooth, postsmooth;
+
+    int numberOfParticles;
+    double startTime, injectTimePeriod, writeTimePeriod;
+
+    double x1, y1, z1, x2, y2, z2;
+
+    int shape;
+    double xCenter, yCenter, zCenter, xRectLength, yRectLength, zRectLength, circleRadius;
+} Parameter;
+
+void initParameter(Parameter*);
+void readParameter(Parameter*, const char*);
+void printParameter(Parameter*);
+#endif

EnhancedSolver/3D-mpi/src/particletracing.c

@@ -0,0 +1,749 @@
+/*
+ * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
+ * All rights reserved. This file is part of nusif-solver.
+ * Use of this source code is governed by a MIT style
+ * license that can be found in the LICENSE file.
+ */
+#include <float.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "particletracing.h"
+#include "solver.h"
+
+static int ts     = 0;
+unsigned int seed = 32767;
+
+#define XOFFSET    0
+#define YOFFSET    1
+#define ZOFFSET    2
+#define XOFFSETEND 3
+#define YOFFSETEND 4
+#define ZOFFSETEND 5
+
+static double sum(int* sizes, int size)
+{
+    double sum = 0;
+
+    for (int i = 0; i < size; ++i) {
+        sum += sizes[i];
+    }
+
+    return sum;
+}
+
+static int maxElement(int* sizes, int size)
+{
+
+    int maxValue = sizes[0];
+
+    for (size_t i = 1; i < size; ++i) {
+        if (sizes[i] > maxValue) {
+            maxValue = sizes[i];
+        }
+    }
+    return maxValue;
+}
+
+static double sumOffset(double* sizes, int init, int offset, int coord)
+{
+    double sum = 0;
+
+    for (int i = init - offset; coord > 0; i -= offset, --coord) {
+        sum += sizes[i];
+    }
+
+    return sum;
+}
+
+void printParticles(ParticleTracer* particletracer)
+{
+    for (int i = 0; i < particletracer->totalParticles; ++i) {
+        printf("Rank : %d Particle position X : %.2f, Y : %.2f, flag : %d, total pt : "
+               "%d, pointer : %d, xOffset : %.2f, yOffset : %.2f, xOffsetEnd : %.2f, "
+               "yOffsetEnd : %.2f\n",
+            particletracer->rank,
+            particletracer->particlePool[i].x,
+            particletracer->particlePool[i].y,
+            particletracer->particlePool[i].flag,
+            particletracer->totalParticles,
+            particletracer->pointer,
+            particletracer->xOffset,
+            particletracer->yOffset,
+            particletracer->xOffsetEnd,
+            particletracer->yOffsetEnd);
+    }
+}
+void injectParticles(ParticleTracer* particletracer, double* restrict s)
+{
+    compress(particletracer);
+    particleRandomizer(particletracer);
+
+    double x, y, z;
+    int imaxLocal = particletracer->imaxLocal;
+    int jmaxLocal = particletracer->jmaxLocal;
+    int kmaxLocal = particletracer->kmaxLocal;
+    for (int i = 0; i < particletracer->numberOfParticles; ++i) {
+        x = particletracer->linSpaceLine[i].x;
+        y = particletracer->linSpaceLine[i].y;
+        z = particletracer->linSpaceLine[i].z;
+        if (x >= particletracer->xOffset && y >= particletracer->yOffset &&
+            z >= particletracer->zOffset && x <= particletracer->xOffsetEnd &&
+            y <= particletracer->yOffsetEnd && y <= particletracer->zOffsetEnd) {
+
+            particletracer->particlePool[particletracer->pointer].x = x;
+            particletracer->particlePool[particletracer->pointer].y = y;
+            particletracer->particlePool[particletracer->pointer].z = z;
+
+            int i = particletracer->particlePool[particletracer->pointer].x /
+                    particletracer->dx;
+            int j = particletracer->particlePool[particletracer->pointer].y /
+                    particletracer->dy;
+            int k = particletracer->particlePool[particletracer->pointer].z /
+                    particletracer->dz;
+
+            int iOffset = particletracer->xOffset / particletracer->dx,
+                jOffset = particletracer->yOffset / particletracer->dy,
+                kOffset = particletracer->zOffset / particletracer->dz;
+
+            if (S(i - iOffset, j - jOffset, k - kOffset) == FLUID) {
+                particletracer->particlePool[particletracer->pointer].flag = true;
+                ++(particletracer->pointer);
+                ++(particletracer->totalParticles);
+            }
+        }
+    }
+}
+
+void advanceParticles(ParticleTracer* particletracer,
+    double* restrict u,
+    double* restrict v,
+    double* restrict w,
+    double* restrict s,
+    double time)
+{
+    int imax      = particletracer->imax;
+    int jmax      = particletracer->jmax;
+    int kmax      = particletracer->kmax;
+    int imaxLocal = particletracer->imaxLocal;
+    int jmaxLocal = particletracer->jmaxLocal;
+    int kmaxLocal = particletracer->kmaxLocal;
+
+    double dx = particletracer->dx;
+    double dy = particletracer->dy;
+    double dz = particletracer->dz;
+
+    double xlength = particletracer->xlength;
+    double ylength = particletracer->ylength;
+    double zlength = particletracer->zlength;
+
+    Particle buff[particletracer->size][(particletracer->estimatedNumParticles / particletracer->size)];
+    memset(buff, 0, sizeof(buff));
+    Particle recvbuff[particletracer->size][(particletracer->estimatedNumParticles / particletracer->size)];
+    memset(buff, 0, sizeof(recvbuff));
+
+    int particleBufIndex[particletracer->size],
+        recvparticleBufIndex[particletracer->size];
+
+    memset(particleBufIndex, 0, sizeof(particleBufIndex));
+    memset(recvparticleBufIndex, 0, sizeof(recvparticleBufIndex));
+
+    for (int i = 0; i < particletracer->totalParticles; ++i) {
+        if (particletracer->particlePool[i].flag == true) {
+            double xTemp = particletracer->particlePool[i].x;
+            double yTemp = particletracer->particlePool[i].y;
+            double zTemp = particletracer->particlePool[i].z;
+
+            double x = xTemp - particletracer->xOffset;
+            double y = yTemp - particletracer->yOffset;
+            double z = zTemp - particletracer->zOffset;
+
+            int iCoord = (int)(x / dx) + 1;
+            int jCoord = (int)((y + 0.5 * dy) / dy) + 1;
+            int kCoord = (int)((z + 0.5 * dz) / dz) + 1;
+
+            double x1 = (double)(iCoord - 1) * dx;
+            double y1 = ((double)(jCoord - 1) - 0.5) * dy;
+            double z1 = ((double)(kCoord - 1) - 0.5) * dz;
+
+            double x2 = (double)iCoord * dx;
+            double y2 = ((double)jCoord - 0.5) * dy;
+            double z2 = ((double)kCoord - 0.5) * dz;
+
+            double u_n =
+                (1.0 / (dx * dy * dz)) *
+                ((x2 - x) * (y2 - y) * (z2 - z) * U(iCoord - 1, jCoord - 1, kCoord - 1) +
+                    (x - x1) * (y2 - y) * (z2 - z) * U(iCoord, jCoord - 1, kCoord - 1) +
+                    (x2 - x) * (y - y1) * (z2 - z) * U(iCoord - 1, jCoord, kCoord - 1) +
+                    (x - x1) * (y - y1) * (z2 - z) * U(iCoord, jCoord, kCoord - 1) +
+                    (x2 - x) * (y2 - y) * (z - z1) * U(iCoord - 1, jCoord - 1, kCoord) +
+                    (x - x1) * (y2 - y) * (z - z1) * U(iCoord, jCoord - 1, kCoord) +
+                    (x2 - x) * (y - y1) * (z - z1) * U(iCoord - 1, jCoord, kCoord) +
+                    (x - x1) * (y - y1) * (z - z1) * U(iCoord, jCoord, kCoord));
+
+            double new_x = (x + particletracer->xOffset) + particletracer->dt * u_n;
+            particletracer->particlePool[i].x = new_x;
+
+            iCoord = (int)((x + 0.5 * dx) / dx) + 1;
+            jCoord = (int)(y / dy) + 1;
+            kCoord = (int)((z + 0.5 * dz) / dz) + 1;
+
+            x1 = ((double)(iCoord - 1) - 0.5) * dx;
+            y1 = (double)(jCoord - 1) * dy;
+            z1 = ((double)(kCoord - 1) - 0.5) * dz;
+
+            x2 = ((double)iCoord - 0.5) * dx;
+            y2 = (double)jCoord * dy;
+            z2 = ((double)kCoord - 0.5) * dz;
+
+            double v_n =
+                (1.0 / (dx * dy * dz)) *
+                ((x2 - x) * (y2 - y) * (z2 - z) * V(iCoord - 1, jCoord - 1, kCoord - 1) +
+                    (x - x1) * (y2 - y) * (z2 - z) * V(iCoord, jCoord - 1, kCoord - 1) +
+                    (x2 - x) * (y - y1) * (z2 - z) * V(iCoord - 1, jCoord, kCoord - 1) +
+                    (x - x1) * (y - y1) * (z2 - z) * V(iCoord, jCoord, kCoord - 1) +
+                    (x2 - x) * (y2 - y) * (z - z1) * V(iCoord - 1, jCoord - 1, kCoord) +
+                    (x - x1) * (y2 - y) * (z - z1) * V(iCoord, jCoord - 1, kCoord) +
+                    (x2 - x) * (y - y1) * (z - z1) * V(iCoord - 1, jCoord, kCoord) +
+                    (x - x1) * (y - y1) * (z - z1) * V(iCoord, jCoord, kCoord));
+
+            double new_y = (y + particletracer->yOffset) + particletracer->dt * v_n;
+            particletracer->particlePool[i].y = new_y;
+
+            iCoord = (int)((x + 0.5 * dx) / dx) + 1;
+            jCoord = (int)((y + 0.5 * dy) / dy) + 1;
+            kCoord = (int)(z / dz) + 1;
+
+            x1 = ((double)(iCoord - 1) - 0.5) * dx;
+            y1 = ((double)(jCoord - 1) - 0.5) * dy;
+            z1 = (double)(kCoord - 1) * dz;
+
+            x2 = ((double)iCoord - 0.5) * dx;
+            y2 = ((double)jCoord - 0.5) * dy;
+            z2 = (double)kCoord * dz;
+
+            double w_n =
+                (1.0 / (dx * dy * dz)) *
+                ((x2 - x) * (y2 - y) * (z2 - z) * W(iCoord - 1, jCoord - 1, kCoord - 1) +
+                    (x - x1) * (y2 - y) * (z2 - z) * W(iCoord, jCoord - 1, kCoord - 1) +
+                    (x2 - x) * (y - y1) * (z2 - z) * W(iCoord - 1, jCoord, kCoord - 1) +
+                    (x - x1) * (y - y1) * (z2 - z) * W(iCoord, jCoord, kCoord - 1) +
+                    (x2 - x) * (y2 - y) * (z - z1) * W(iCoord - 1, jCoord - 1, kCoord) +
+                    (x - x1) * (y2 - y) * (z - z1) * W(iCoord, jCoord - 1, kCoord) +
+                    (x2 - x) * (y - y1) * (z - z1) * W(iCoord - 1, jCoord, kCoord) +
+                    (x - x1) * (y - y1) * (z - z1) * W(iCoord, jCoord, kCoord));
+
+            double new_z = (z + particletracer->zOffset) + particletracer->dt * w_n;
+            particletracer->particlePool[i].z = new_z;
+
+            if (((new_x < particletracer->xOffset) ||
+                    (new_x >= particletracer->xOffsetEnd) ||
+                    (new_y < particletracer->yOffset) ||
+                    (new_y >= particletracer->yOffsetEnd) ||
+                    (new_z < particletracer->zOffset) ||
+                    (new_z >= particletracer->zOffsetEnd))) {
+                // New logic to transfer particles to neighbouring ranks or discard the
+                // particle.
+
+#ifdef _MPI
+                for (int i = 0; i < particletracer->size; ++i) {
+                    if ((new_x >=
+                            particletracer->offset[i + particletracer->size * XOFFSET]) &&
+                        (new_x <= particletracer
+                                      ->offset[i + particletracer->size * XOFFSETEND]) &&
+                        (new_y >=
+                            particletracer->offset[i + particletracer->size * YOFFSET]) &&
+                        (new_y <= particletracer
+                                      ->offset[i + particletracer->size * YOFFSETEND]) &&
+                        (new_z >=
+                            particletracer->offset[i + particletracer->size * ZOFFSET]) &&
+                        (new_z <= particletracer
+                                      ->offset[i + particletracer->size * ZOFFSETEND]) &&
+                        i != particletracer->rank) {
+                        buff[i][particleBufIndex[i]].x    = new_x;
+                        buff[i][particleBufIndex[i]].y    = new_y;
+                        buff[i][particleBufIndex[i]].z    = new_z;
+                        buff[i][particleBufIndex[i]].flag = true;
+                        ++particleBufIndex[i];
+                    }
+                }
+#endif
+
+                particletracer->particlePool[i].flag = false;
+                particletracer->removedParticles++;
+            }
+
+            int i_new = new_x / dx, j_new = new_y / dy, k_new = new_z / dz;
+
+            int iOffset = particletracer->xOffset / dx,
+                jOffset = particletracer->yOffset / dy,
+                kOffset = particletracer->zOffset / dz;
+
+            if (S(i_new - iOffset, j_new - jOffset, k_new - kOffset) != FLUID) {
+                particletracer->particlePool[i].flag = false;
+                particletracer->removedParticles++;
+            }
+        }
+    }
+
+#ifdef _MPI
+    for (int i = 0; i < particletracer->size; ++i) {
+        if (i != particletracer->rank) {
+            MPI_Send(&particleBufIndex[i], 1, MPI_INT, i, 0, particletracer->comm);
+        }
+    }
+    for (int i = 0; i < particletracer->size; ++i) {
+        if (i != particletracer->rank) {
+            MPI_Recv(&recvparticleBufIndex[i],
+                1,
+                MPI_INT,
+                i,
+                0,
+                particletracer->comm,
+                MPI_STATUS_IGNORE);
+
+            // if (0 !=recvparticleBufIndex[i]) {
+            //     printf("Rank %d will receive %d particles from rank %d\n",
+            //         particletracer->rank,
+            //         recvparticleBufIndex[i],
+            //         i);
+            // }
+        }
+    }
+
+    for (int i = 0; i < particletracer->size; ++i) {
+        if (i != particletracer->rank) {
+            MPI_Send(buff[i],
+                particleBufIndex[i],
+                particletracer->mpi_particle,
+                i,
+                0,
+                particletracer->comm);
+        }
+    }
+    for (int i = 0; i < particletracer->size; ++i) {
+        if (i != particletracer->rank) {
+            MPI_Recv(recvbuff[i],
+                recvparticleBufIndex[i],
+                particletracer->mpi_particle,
+                i,
+                0,
+                particletracer->comm,
+                MPI_STATUS_IGNORE);
+        }
+    }
+    for (int i = 0; i < particletracer->size; ++i) {
+        if (i != particletracer->rank) {
+            for (int j = 0; j < recvparticleBufIndex[i]; ++j) {
+                particletracer->particlePool[particletracer->pointer].x = recvbuff[i][j]
+                                                                              .x;
+                particletracer->particlePool[particletracer->pointer].y = recvbuff[i][j]
+                                                                              .y;
+                particletracer->particlePool[particletracer->pointer].z = recvbuff[i][j]
+                                                                              .z;
+                particletracer->particlePool[particletracer->pointer].flag = true;
+                ++(particletracer->pointer);
+                ++(particletracer->totalParticles);
+            }
+        }
+    }
+#endif
+}
+
+void freeParticles(ParticleTracer* particletracer)
+{
+    free(particletracer->particlePool);
+    free(particletracer->linSpaceLine);
+    free(particletracer->offset);
+}
+
+void writeParticles(ParticleTracer* particletracer)
+{
+    compress(particletracer);
+
+#ifdef _MPI
+    int collectedBuffIndex[particletracer->size];
+
+    MPI_Gather(&particletracer->totalParticles,
+        1,
+        MPI_INT,
+        collectedBuffIndex,
+        1,
+        MPI_INT,
+        0,
+        particletracer->comm);
+
+    if (particletracer->rank != 0) {
+        Particle buff[particletracer->totalParticles];
+        for (int i = 0; i < particletracer->totalParticles; ++i) {
+            buff[i].x    = particletracer->particlePool[i].x;
+            buff[i].y    = particletracer->particlePool[i].y;
+            buff[i].z    = particletracer->particlePool[i].z;
+            buff[i].flag = particletracer->particlePool[i].flag;
+            // printf("Rank : %d sending to rank 0 X : %.2f, Y : %.2f with totalpt :
+            // %d\n", particletracer->rank, buff[i].x, buff[i].y,
+            // particletracer->totalParticles);
+        }
+        MPI_Send(buff,
+            particletracer->totalParticles,
+            particletracer->mpi_particle,
+            0,
+            1,
+            particletracer->comm);
+    }
+#endif
+
+    if (particletracer->rank == 0) {
+        char filename[50];
+        FILE* fp;
+
+        snprintf(filename, 50, "vtk_files/particles_%d.vtk", ts);
+        fp = fopen(filename, "w");
+
+        if (fp == NULL) {
+            printf("Error!\n");
+            exit(EXIT_FAILURE);
+        }
+        fprintf(fp, "# vtk DataFile Version 3.0\n");
+        fprintf(fp, "PAMPI cfd solver particle tracing file\n");
+        fprintf(fp, "ASCII\n");
+
+        fprintf(fp, "DATASET UNSTRUCTURED_GRID\n");
+        fprintf(fp, "FIELD FieldData 2\n");
+        fprintf(fp, "TIME 1 1 double\n");
+        fprintf(fp, "%d\n", ts);
+        fprintf(fp, "CYCLE 1 1 int\n");
+        fprintf(fp, "1\n");
+
+#ifdef _MPI
+        int overallTotalParticles = sum(collectedBuffIndex, particletracer->size);
+
+        fprintf(fp, "POINTS %d float\n", overallTotalParticles);
+
+        // printf("Total particles : %d\n", overallTotalParticles);
+
+        for (int i = 1; i < particletracer->size; ++i) {
+            Particle recvBuff[collectedBuffIndex[i]];
+            MPI_Recv(&recvBuff,
+                collectedBuffIndex[i],
+                particletracer->mpi_particle,
+                i,
+                1,
+                particletracer->comm,
+                MPI_STATUS_IGNORE);
+
+            for (int j = 0; j < collectedBuffIndex[i]; ++j) {
+                double x = recvBuff[j].x;
+                double y = recvBuff[j].y;
+                double z = recvBuff[j].z;
+                fprintf(fp, "%f %f %f\n", x, y, z);
+            }
+        }
+#else
+        int overallTotalParticles = particletracer->totalParticles;
+        fprintf(fp, "POINTS %d float\n", overallTotalParticles);
+
+        // printf("Total particles : %d\n", overallTotalParticles);
+#endif
+        for (int i = 0; i < particletracer->totalParticles; ++i) {
+            double x = particletracer->particlePool[i].x;
+            double y = particletracer->particlePool[i].y;
+            double z = particletracer->particlePool[i].z;
+            fprintf(fp, "%f %f %f\n", x, y, z);
+        }
+
+        fprintf(fp, "CELLS %d %d\n", overallTotalParticles, 2 * overallTotalParticles);
+
+        for (int i = 0; i < overallTotalParticles; ++i) {
+            fprintf(fp, "1 %d\n", i);
+        }
+
+        fprintf(fp, "CELL_TYPES %d\n", overallTotalParticles);
+
+        for (int i = 0; i < overallTotalParticles; ++i) {
+            fprintf(fp, "1\n");
+        }
+
+        fclose(fp);
+    }
+
+    ++ts;
+}
+
+void initParticleTracer(
+    ParticleTracer* particletracer, Parameter* params, Discretization* d)
+{
+    /* initializing local properties from params */
+
+    particletracer->numberOfParticles = params->numberOfParticles;
+    particletracer->startTime         = params->startTime;
+    particletracer->injectTimePeriod  = params->injectTimePeriod;
+    particletracer->writeTimePeriod   = params->writeTimePeriod;
+
+    particletracer->rank = d->comm.rank;
+    particletracer->size = d->comm.size;
+
+    particletracer->dt = params->dt;
+    particletracer->dx = params->xlength / params->imax;
+    particletracer->dy = params->ylength / params->jmax;
+    particletracer->dz = params->zlength / params->kmax;
+
+    particletracer->xlength = params->xlength;
+    particletracer->ylength = params->ylength;
+    particletracer->zlength = params->zlength;
+
+    particletracer->x1 = params->x1;
+    particletracer->y1 = params->y1;
+    particletracer->z1 = params->z1;
+    particletracer->x2 = params->x2;
+    particletracer->y2 = params->y2;
+    particletracer->z2 = params->z2;
+
+    particletracer->lastInjectTime = params->startTime;
+    particletracer->lastUpdateTime = params->startTime;
+    particletracer->lastWriteTime  = params->startTime;
+
+    particletracer->pointer          = 0;
+    particletracer->totalParticles   = 0;
+    particletracer->removedParticles = 0;
+
+    particletracer->imax = params->imax;
+    particletracer->jmax = params->jmax;
+    particletracer->kmax = params->kmax;
+
+    particletracer->imaxLocal = d->comm.imaxLocal;
+    particletracer->jmaxLocal = d->comm.jmaxLocal;
+    particletracer->kmaxLocal = d->comm.kmaxLocal;
+
+    particletracer->estimatedNumParticles = (particletracer->imaxLocal *
+                                             particletracer->jmaxLocal *
+                                             particletracer->kmaxLocal);
+
+    particletracer->particlePool = malloc(
+        sizeof(Particle) * particletracer->estimatedNumParticles);
+
+    memset(particletracer->particlePool,
+        0,
+        sizeof(Particle) * particletracer->estimatedNumParticles);
+
+    particletracer->linSpaceLine = malloc(
+        sizeof(Particle) * particletracer->numberOfParticles);
+
+    memset(particletracer->linSpaceLine,
+        0,
+        sizeof(Particle) * particletracer->numberOfParticles);
+
+    particletracer->offset = (double*)malloc(sizeof(double) * 6 * particletracer->size);
+    double offset[6][particletracer->size];
+
+#ifdef _MPI
+    /* duplicating communication from solver */
+    MPI_Comm_dup(d->comm.comm, &particletracer->comm);
+
+    memcpy(particletracer->dims, d->comm.dims, sizeof(d->comm.dims));
+
+    memcpy(particletracer->coords, d->comm.coords, sizeof(d->comm.coords));
+
+#endif
+
+    particletracer->xOffset    = d->xOffset;
+    particletracer->yOffset    = d->yOffset;
+    particletracer->zOffset    = d->zOffset;
+    particletracer->xOffsetEnd = d->xOffsetEnd;
+    particletracer->yOffsetEnd = d->yOffsetEnd;
+    particletracer->zOffsetEnd = d->zOffsetEnd;
+
+#ifdef _MPI
+    MPI_Allgather(&particletracer->xOffset,
+        1,
+        MPI_DOUBLE,
+        offset[0],
+        1,
+        MPI_DOUBLE,
+        particletracer->comm);
+    MPI_Allgather(&particletracer->yOffset,
+        1,
+        MPI_DOUBLE,
+        offset[1],
+        1,
+        MPI_DOUBLE,
+        particletracer->comm);
+    MPI_Allgather(&particletracer->zOffset,
+        1,
+        MPI_DOUBLE,
+        offset[2],
+        1,
+        MPI_DOUBLE,
+        particletracer->comm);
+    MPI_Allgather(&particletracer->xOffsetEnd,
+        1,
+        MPI_DOUBLE,
+        offset[3],
+        1,
+        MPI_DOUBLE,
+        particletracer->comm);
+    MPI_Allgather(&particletracer->yOffsetEnd,
+        1,
+        MPI_DOUBLE,
+        offset[4],
+        1,
+        MPI_DOUBLE,
+        particletracer->comm);
+    MPI_Allgather(&particletracer->zOffsetEnd,
+        1,
+        MPI_DOUBLE,
+        offset[5],
+        1,
+        MPI_DOUBLE,
+        particletracer->comm);
+
+#endif
+    memcpy(particletracer->offset, offset, sizeof(offset));
+
+    particleRandomizer(particletracer);
+
+#ifdef _MPI
+    // Create the mpi_particle datatype
+    MPI_Datatype mpi_particle;
+    int lengths[4] = { 1, 1, 1, 1 };
+
+    MPI_Aint displacements[4];
+    Particle dummy_particle;
+    MPI_Aint base_address;
+    MPI_Get_address(&dummy_particle, &base_address);
+    MPI_Get_address(&dummy_particle.x, &displacements[0]);
+    MPI_Get_address(&dummy_particle.y, &displacements[1]);
+    MPI_Get_address(&dummy_particle.z, &displacements[2]);
+    MPI_Get_address(&dummy_particle.flag, &displacements[3]);
+    displacements[0]      = MPI_Aint_diff(displacements[0], base_address);
+    displacements[1]      = MPI_Aint_diff(displacements[1], base_address);
+    displacements[2]      = MPI_Aint_diff(displacements[2], base_address);
+    displacements[3]      = MPI_Aint_diff(displacements[3], base_address);
+    MPI_Datatype types[4] = { MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_C_BOOL };
+    MPI_Type_create_struct(4,
+        lengths,
+        displacements,
+        types,
+        &particletracer->mpi_particle);
+    MPI_Type_commit(&particletracer->mpi_particle);
+#endif
+}
+
+void printParticleTracerParameters(ParticleTracer* particletracer)
+{
+    printf("Particle Tracing data:\n");
+    printf("Rank : %d\n", particletracer->rank);
+    printf("\tNumber of particles : %d being injected for every period of %.2f\n",
+        particletracer->numberOfParticles,
+        particletracer->injectTimePeriod);
+    printf("\tstartTime : %.2f\n", particletracer->startTime);
+    printf("\t(Line along which the particles are to be injected) \n\tx1 : %.2f, y1 : "
+           "%.2f, z1 : %.2f, x2 : %.2f, y2 : %.2f,  z2 : %.2f\n",
+        particletracer->x1,
+        particletracer->y1,
+        particletracer->z1,
+        particletracer->x2,
+        particletracer->y2,
+        particletracer->z2);
+    printf("\tPointer : %d, TotalParticles : %d\n",
+        particletracer->pointer,
+        particletracer->totalParticles);
+    printf("\tdt : %.2f, dx : %.2f, dy : %.2f, dz : %.2f\n",
+        particletracer->dt,
+        particletracer->dx,
+        particletracer->dy,
+        particletracer->dz);
+    printf("\tcoord[0] : %d, coord[1] : %d, coord[2] : %d\n",
+        particletracer->coords[IDIM],
+        particletracer->coords[JDIM],
+        particletracer->coords[KDIM]);
+    printf("\txOffset : %.2f, yOffset : %.2f, zOffset : %.2f\n",
+        particletracer->xOffset,
+        particletracer->yOffset,
+        particletracer->zOffset);
+    printf("\txOffsetEnd : %.2f, yOffsetEnd : %.2f, zOffsetEnd : %.2f\n",
+        particletracer->xOffsetEnd,
+        particletracer->yOffsetEnd,
+        particletracer->zOffsetEnd);
+    printf("\txLocal : %.2f, yLocal : %.2f, zLocal : %.2f\n",
+        particletracer->xLocal,
+        particletracer->yLocal,
+        particletracer->zLocal);
+}
+
+void trace(ParticleTracer* particletracer,
+    double* restrict u,
+    double* restrict v,
+    double* restrict w,
+    double* restrict s,
+    double time)
+{
+    if (time >= particletracer->startTime) {
+        if ((time - particletracer->lastInjectTime) >= particletracer->injectTimePeriod) {
+            injectParticles(particletracer, s);
+            particletracer->lastInjectTime = time;
+        }
+        if ((time - particletracer->lastWriteTime) >= particletracer->writeTimePeriod) {
+            writeParticles(particletracer);
+            particletracer->lastWriteTime = time;
+        }
+        advanceParticles(particletracer, u, v, w, s, time);
+
+        if (particletracer->removedParticles > (particletracer->totalParticles * 0.2)) {
+            compress(particletracer);
+        }
+        particletracer->lastUpdateTime = time;
+    }
+}
+
+void compress(ParticleTracer* particletracer)
+{
+    Particle* memPool = particletracer->particlePool;
+    Particle tempPool[particletracer->totalParticles];
+
+    memset(tempPool, 0, sizeof(Particle) * particletracer->totalParticles);
+
+    int totalParticles = 0;
+
+    for (int i = 0; i < particletracer->totalParticles; ++i) {
+        if (memPool[i].flag == true) {
+            tempPool[totalParticles].x    = memPool[i].x;
+            tempPool[totalParticles].y    = memPool[i].y;
+            tempPool[totalParticles].z    = memPool[i].z;
+            tempPool[totalParticles].flag = true;
+            ++totalParticles;
+        }
+    }
+
+    particletracer->totalParticles   = totalParticles;
+    particletracer->removedParticles = 0;
+    particletracer->pointer          = totalParticles;
+
+    memcpy(particletracer->particlePool, tempPool, totalParticles * sizeof(Particle));
+}
+
+void particleRandomizer(ParticleTracer* particletracer)
+{
+    memset(particletracer->linSpaceLine,
+        0,
+        sizeof(Particle) * particletracer->numberOfParticles);
+
+    for (int i = 0; i < particletracer->numberOfParticles; ++i) {
+
+        particletracer->linSpaceLine[i].x = (((double)rand() / RAND_MAX) *
+                                                (particletracer->x2 -
+                                                    particletracer->x1)) +
+                                            particletracer->x1;
+        particletracer->linSpaceLine[i].y = (((double)rand() / RAND_MAX) *
+                                                (particletracer->y2 -
+                                                    particletracer->y1)) +
+                                            particletracer->y1;
+        particletracer->linSpaceLine[i].z = (((double)rand() / RAND_MAX) *
+                                                (particletracer->z2 -
+                                                    particletracer->z1)) +
+                                            particletracer->z1;
+        particletracer->linSpaceLine[i].flag = true;
+    }
+}

EnhancedSolver/3D-mpi/src/particletracing.h

@@ -0,0 +1,65 @@
+/*
+ * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
+ * All rights reserved. This file is part of nusif-solver.
+ * Use of this source code is governed by a MIT style
+ * license that can be found in the LICENSE file.
+ */
+#ifndef __PARTICLETRACING_H_
+#define __PARTICLETRACING_H_
+#include "allocate.h"
+#include "parameter.h"
+#include "discretization.h"
+#include <mpi.h>
+#include <stdbool.h>
+
+#define NDIMS 3
+
+typedef enum COORD { X = 0, Y, NCOORD } COORD;
+typedef struct {
+    double x, y, z;
+    bool flag;
+} Particle;
+
+typedef struct {
+    int numberOfParticles, totalParticles;
+    double startTime, injectTimePeriod, writeTimePeriod, lastInjectTime, lastUpdateTime,
+        lastWriteTime;
+
+    int estimatedNumParticles, removedParticles;
+
+    double dx, dy, dz, dt;
+    Particle* linSpaceLine;
+    Particle* particlePool;
+
+    int pointer;
+
+    double imax, jmax, kmax, xlength, ylength, zlength, imaxLocal, jmaxLocal, kmaxLocal;
+
+    double x1, y1, z1, x2, y2, z2;
+
+#ifdef _MPI
+    MPI_Comm comm;
+    MPI_Datatype mpi_particle;
+#endif
+
+    int rank, size;
+    int coords[NDIMS], dims[NDIMS];
+
+    double xLocal, yLocal, zLocal, xOffset, yOffset, zOffset, xOffsetEnd, yOffsetEnd,
+        zOffsetEnd;
+
+    double* offset;
+
+} ParticleTracer;
+
+void initParticleTracer(ParticleTracer*, Parameter*, Discretization*);
+void injectParticles(ParticleTracer*, double*);
+void advanceParticles(ParticleTracer*, double*, double*, double*, double*, double);
+void freeParticles(ParticleTracer*);
+void writeParticles(ParticleTracer*);
+void printParticleTracerParameters(ParticleTracer*);
+void printParticles(ParticleTracer*);
+void trace(ParticleTracer*, double*, double*, double*, double*, double);
+void compress(ParticleTracer*);
+extern void particleRandomizer(ParticleTracer*);
+#endif

EnhancedSolver/3D-mpi/src/progress.c

@@ -0,0 +1,65 @@
+/*
+ * 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 "progress.h"
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+static double _end;
+static int _current;
+
+void initProgress(double end)
+{
+    _end     = end;
+    _current = 0;
+
+    printf("[          ]");
+    fflush(stdout);
+}
+
+void printProgress(double current)
+{
+    int new = (int)rint((current / _end) * 10.0);
+
+    if (new > _current) {
+        char progress[11];
+        _current    = new;
+        progress[0] = 0;
+
+        for (int i = 0; i < 10; i++) {
+            if (i < _current) {
+                sprintf(progress + strlen(progress), "#");
+            } else {
+                sprintf(progress + strlen(progress), " ");
+            }
+        }
+        printf("\r[%s]", progress);
+    }
+    fflush(stdout);
+}
+
+void stopProgress()
+{
+    printf("\n");
+    fflush(stdout);
+}
+
+FILE* initResidualWriter()
+{
+    FILE* fp;
+    fp = fopen("residual.dat", "w");
+
+    if (fp == NULL) {
+        printf("Error!\n");
+        exit(EXIT_FAILURE);
+    }
+
+    return fp;
+}
+
+void writeResidual(FILE* fp, double ts, double res) { fprintf(fp, "%f, %f\n", ts, res); }

EnhancedSolver/3D-mpi/src/progress.h

@@ -0,0 +1,17 @@
+/*
+ * 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.
+ */
+
+#include <stdio.h>
+
+#ifndef __PROGRESS_H_
+#define __PROGRESS_H_
+extern void initProgress(double);
+extern void printProgress(double);
+extern void stopProgress(void);
+extern FILE* initResidualWriter(void);
+extern void writeResidual(FILE*, double, double);
+#endif

EnhancedSolver/3D-mpi/src/solver-mg.c

@@ -0,0 +1,389 @@
+/*
+ * 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)
+static void restrictMG(Solver* s, int level, Comm* comm)
+{
+    double* r   = s->r[level + 1];
+    double* old = s->r[level];
+
+    int imaxLocal = comm->imaxLocal;
+    int jmaxLocal = comm->jmaxLocal;
+    int kmaxLocal = comm->kmaxLocal;
+
+    commExchange(comm, old);
+
+    for (int k = 1; k < comm->kmaxLocal + 1; k++) {
+        for (int j = 1; j < comm->jmaxLocal + 1; j++) {
+            for (int i = 1; i < comm->imaxLocal + 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, Comm* comm)
+{
+    double* old = s->r[level + 1];
+    double* e   = s->r[level];
+
+    int imaxLocal = comm->imaxLocal;
+    int jmaxLocal = comm->jmaxLocal;
+    int kmaxLocal = comm->kmaxLocal;
+
+    for (int k = 2; k < kmaxLocal + 1; k += 2) {
+        for (int j = 2; j < jmaxLocal + 1; j += 2) {
+            for (int i = 2; i < imaxLocal + 1; i += 2) {
+                E(i, j, k) = OLD(i / 2, j / 2, k / 2);
+            }
+        }
+    }
+}
+
+static void correct(Solver* s, double* p, int level, Comm* comm)
+{
+    double* e = s->e[level];
+
+    int imaxLocal = comm->imaxLocal;
+    int jmaxLocal = comm->jmaxLocal;
+    int kmaxLocal = comm->kmaxLocal;
+
+    for (int k = 1; k < kmaxLocal + 1; ++k) {
+        for (int j = 1; j < jmaxLocal + 1; ++j) {
+            for (int i = 1; i < imaxLocal + 1; ++i) {
+                P(i, j, k) += E(i, j, k);
+            }
+        }
+    }
+}
+
+static void setBoundaryCondition(
+    Solver* s, double* p, int imaxLocal, int jmaxLocal, int kmaxLocal)
+{
+#ifdef _MPI
+    if (commIsBoundary(s->comm, FRONT)) {
+        for (int j = 1; j < jmaxLocal + 1; j++) {
+            for (int i = 1; i < imaxLocal + 1; i++) {
+                P(i, j, 0) = P(i, j, 1);
+            }
+        }
+    }
+
+    if (commIsBoundary(s->comm, BACK)) {
+        for (int j = 1; j < jmaxLocal + 1; j++) {
+            for (int i = 1; i < imaxLocal + 1; i++) {
+                P(i, j, kmaxLocal + 1) = P(i, j, kmaxLocal);
+            }
+        }
+    }
+
+    if (commIsBoundary(s->comm, BOTTOM)) {
+        for (int k = 1; k < kmaxLocal + 1; k++) {
+            for (int i = 1; i < imaxLocal + 1; i++) {
+                P(i, 0, k) = P(i, 1, k);
+            }
+        }
+    }
+
+    if (commIsBoundary(s->comm, TOP)) {
+        for (int k = 1; k < kmaxLocal + 1; k++) {
+            for (int i = 1; i < imaxLocal + 1; i++) {
+                P(i, jmaxLocal + 1, k) = P(i, jmaxLocal, k);
+            }
+        }
+    }
+
+    if (commIsBoundary(s->comm, LEFT)) {
+        for (int k = 1; k < kmaxLocal + 1; k++) {
+            for (int j = 1; j < jmaxLocal + 1; j++) {
+                P(0, j, k) = P(1, j, k);
+            }
+        }
+    }
+
+    if (commIsBoundary(s->comm, RIGHT)) {
+        for (int k = 1; k < kmaxLocal + 1; k++) {
+            for (int j = 1; j < jmaxLocal + 1; j++) {
+                P(imaxLocal + 1, j, k) = P(imaxLocal, j, k);
+            }
+        }
+    }
+#else
+    for (int j = 1; j < jmaxLocal + 1; j++) {
+        for (int i = 1; i < imaxLocal + 1; i++) {
+            P(i, j, 0)             = P(i, j, 1);
+            P(i, j, kmaxLocal + 1) = P(i, j, kmaxLocal);
+        }
+    }
+
+    for (int k = 1; k < kmaxLocal + 1; k++) {
+        for (int i = 1; i < imaxLocal + 1; i++) {
+            P(i, 0, k)             = P(i, 1, k);
+            P(i, jmaxLocal + 1, k) = P(i, jmaxLocal, k);
+        }
+    }
+
+    for (int k = 1; k < kmaxLocal + 1; k++) {
+        for (int j = 1; j < jmaxLocal + 1; j++) {
+            P(0, j, k)             = P(1, j, k);
+            P(imaxLocal + 1, j, k) = P(imaxLocal, j, k);
+        }
+    }
+#endif
+}
+
+static double smooth(Solver* s, double* p, double* rhs, int level, Comm* comm)
+{
+    int imaxLocal = comm->imaxLocal;
+    int jmaxLocal = comm->jmaxLocal;
+    int kmaxLocal = comm->kmaxLocal;
+
+    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;
+
+        commExchange(comm, p);
+
+        for (int k = 1; k < kmaxLocal + 1; k++) {
+            isw = jsw;
+            for (int j = 1; j < jmaxLocal + 1; j++) {
+                for (int i = isw; i < imaxLocal + 1; i += 2) {
+
+                    P(i, j, k) -=
+                        factor *
+                        (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));
+                }
+                isw = 3 - isw;
+            }
+            jsw = 3 - jsw;
+        }
+        ksw = 3 - ksw;
+    }
+}
+
+static double calculateResidual(Solver* s, double* p, double* rhs, int level, Comm* comm)
+{
+
+    int imaxLocal = comm->imaxLocal;
+    int jmaxLocal = comm->jmaxLocal;
+    int kmaxLocal = comm->kmaxLocal;
+
+    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;
+
+        commExchange(comm, p);
+
+        for (int k = 1; k < kmaxLocal + 1; k++) {
+            isw = jsw;
+            for (int j = 1; j < jmaxLocal + 1; j++) {
+                for (int i = isw; i < imaxLocal + 1; i += 2) {
+
+                    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));
+
+                    res += (R(i, j, k) * R(i, j, k));
+                }
+                isw = 3 - isw;
+            }
+            jsw = 3 - jsw;
+        }
+        ksw = 3 - ksw;
+    }
+
+    commReduction(&res, SUM);
+
+    res = res / (double)(imaxLocal * jmaxLocal * kmaxLocal);
+#ifdef DEBUG
+    if (commIsMaster(s->comm)) {
+        printf("%d Residuum: %e\n", it, res);
+    }
+#endif
+    return res;
+}
+
+static double multiGrid(Solver* s, double* p, double* rhs, int level, Comm* comm)
+{
+    int imaxLocal = comm->imaxLocal;
+    int jmaxLocal = comm->jmaxLocal;
+    int kmaxLocal = comm->kmaxLocal;
+
+    double res = 0.0;
+
+    // coarsest level
+    if (level == COARSEST_LEVEL) {
+        for (int i = 0; i < 5; i++) {
+            smooth(s, p, rhs, level, comm);
+        }
+        return res;
+    }
+
+    // pre-smoothing
+    for (int i = 0; i < s->presmooth; i++) {
+        smooth(s, p, rhs, level, comm);
+        if (level == FINEST_LEVEL)
+            setBoundaryCondition(s, p, imaxLocal, jmaxLocal, kmaxLocal);
+    }
+
+    res = calculateResidual(s, p, rhs, level, comm);
+
+    // restrict
+    restrictMG(s, level, comm);
+
+    Comm newcomm;
+    commUpdateDatatypes(s->comm, &newcomm, imaxLocal, jmaxLocal, kmaxLocal);
+
+    // MGSolver on residual and error.
+    multiGrid(s, s->e[level + 1], s->r[level + 1], level + 1, &newcomm);
+
+    commFreeCommunicator(&newcomm);
+
+    // prolongate
+    prolongate(s, level, comm);
+
+    // correct p on finer level using residual
+    correct(s, p, level, comm);
+    if (level == FINEST_LEVEL)
+        setBoundaryCondition(s, p, imaxLocal, jmaxLocal, kmaxLocal);
+
+    // post-smoothing
+    for (int i = 0; i < s->postsmooth; i++) {
+        smooth(s, p, rhs, level, comm);
+        if (level == FINEST_LEVEL)
+            setBoundaryCondition(s, p, imaxLocal, jmaxLocal, kmaxLocal);
+    }
+
+    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;
+    s->presmooth  = p->presmooth;
+    s->postsmooth = p->postsmooth;
+    s->comm       = &d->comm;
+    s->problem    = p->name;
+
+    int imax   = s->grid->imax;
+    int jmax   = s->grid->jmax;
+    int kmax   = s->grid->kmax;
+    int levels = s->levels;
+    if (commIsMaster(s->comm)) 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;
+        }
+    }
+}
+
+double solve(Solver* s, double* p, double* rhs)
+{
+    double res = multiGrid(s, p, rhs, 0, s->comm);
+#ifdef VERBOSE
+    if (commIsMaster(s->comm)) {
+        printf("Residuum: %.6f\n", res);
+    }
+#endif
+
+    return res;
+}

EnhancedSolver/3D-mpi/src/solver-rb.c

@@ -0,0 +1,175 @@
+/*
+ * 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 "comm.h"
+#include "parameter.h"
+#include "solver.h"
+#include "util.h"
+
+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;
+    s->presmooth  = p->presmooth;
+    s->postsmooth = p->postsmooth;
+    s->comm       = &d->comm;
+    s->problem    = p->name;
+}
+
+double solve(Solver* s, double* p, double* rhs)
+{
+    int imaxLocal = s->comm->imaxLocal;
+    int jmaxLocal = s->comm->jmaxLocal;
+    int kmaxLocal = s->comm->kmaxLocal;
+
+    int imax = s->grid->imax;
+    int jmax = s->grid->jmax;
+    int kmax = s->grid->kmax;
+
+    double eps  = s->eps;
+    int itermax = s->itermax;
+    double dx2  = s->grid->dx * s->grid->dx;
+    double dy2  = s->grid->dy * s->grid->dy;
+    double dz2  = s->grid->dz * s->grid->dz;
+    double idx2 = 1.0 / dx2;
+    double idy2 = 1.0 / dy2;
+    double idz2 = 1.0 / dz2;
+
+    double factor = s->omega * 0.5 * (dx2 * dy2 * dz2) /
+                    (dy2 * dz2 + dx2 * dz2 + dx2 * dy2);
+    double epssq = eps * eps;
+    int it       = 0;
+    double res   = 1.0;
+    int pass, ksw, jsw, isw;
+
+    while ((res >= epssq) && (it < itermax)) {
+        ksw = 1;
+
+        for (pass = 0; pass < 2; pass++) {
+            jsw = ksw;
+            commExchange(s->comm, p);
+
+            for (int k = 1; k < kmaxLocal + 1; k++) {
+                isw = jsw;
+                for (int j = 1; j < jmaxLocal + 1; j++) {
+                    for (int i = isw; i < imaxLocal + 1; i += 2) {
+
+                        double r =
+                            RHS(i, j, k) -
+                            ((P(i + 1, j, k) - 2.0 * P(i, j, k) + P(i - 1, j, k)) * idx2 +
+                                (P(i, j + 1, k) - 2.0 * P(i, j, k) + P(i, j - 1, k)) *
+                                    idy2 +
+                                (P(i, j, k + 1) - 2.0 * P(i, j, k) + P(i, j, k - 1)) *
+                                    idz2);
+
+                        P(i, j, k) -= (factor * r);
+                        res += (r * r);
+                    }
+                    isw = 3 - isw;
+                }
+                jsw = 3 - jsw;
+            }
+            ksw = 3 - ksw;
+        }
+#ifdef _MPI
+        if (commIsBoundary(s->comm, FRONT)) {
+            for (int j = 1; j < jmaxLocal + 1; j++) {
+                for (int i = 1; i < imaxLocal + 1; i++) {
+                    P(i, j, 0) = P(i, j, 1);
+                }
+            }
+        }
+
+        if (commIsBoundary(s->comm, BACK)) {
+            for (int j = 1; j < jmaxLocal + 1; j++) {
+                for (int i = 1; i < imaxLocal + 1; i++) {
+                    P(i, j, kmaxLocal + 1) = P(i, j, kmaxLocal);
+                }
+            }
+        }
+
+        if (commIsBoundary(s->comm, BOTTOM)) {
+            for (int k = 1; k < kmaxLocal + 1; k++) {
+                for (int i = 1; i < imaxLocal + 1; i++) {
+                    P(i, 0, k) = P(i, 1, k);
+                }
+            }
+        }
+
+        if (commIsBoundary(s->comm, TOP)) {
+            for (int k = 1; k < kmaxLocal + 1; k++) {
+                for (int i = 1; i < imaxLocal + 1; i++) {
+                    P(i, jmaxLocal + 1, k) = P(i, jmaxLocal, k);
+                }
+            }
+        }
+
+        if (commIsBoundary(s->comm, LEFT)) {
+            for (int k = 1; k < kmaxLocal + 1; k++) {
+                for (int j = 1; j < jmaxLocal + 1; j++) {
+                    P(0, j, k) = P(1, j, k);
+                }
+            }
+        }
+
+        if (commIsBoundary(s->comm, RIGHT)) {
+            for (int k = 1; k < kmaxLocal + 1; k++) {
+                for (int j = 1; j < jmaxLocal + 1; j++) {
+                    P(imaxLocal + 1, j, k) = P(imaxLocal, j, k);
+                }
+            }
+        }
+#else
+        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);
+            }
+        }
+#endif
+        commReduction(&res, SUM);
+        res = res / (double)(imax * jmax * kmax);
+#ifdef DEBUG
+        if (commIsMaster(&s->comm)) {
+            printf("%d Residuum: %e\n", it, res);
+        }
+#endif
+        commExchange(s->comm, p);
+        it++;
+    }
+
+#ifdef VERBOSE
+    if (commIsMaster(s->comm)) {
+        printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
+    }
+#endif
+
+return res;
+}

EnhancedSolver/3D-mpi/src/solver.h

@@ -0,0 +1,39 @@
+/*
+ * 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 __SOLVER_H_
+#define __SOLVER_H_
+#include "comm.h"
+#include "grid.h"
+#include "parameter.h"
+#include "discretization.h"
+
+typedef struct {
+    /* geometry and grid information */
+    Grid* grid;
+    /* arrays */
+    double *p, *rhs;
+    double *f, *g, *h;
+    double *u, *v, *w;
+    /* parameters */
+    double eps, omega;
+    double re, tau, gamma;
+    double gx, gy, gz;
+    /* time stepping */
+    int itermax;
+    double dt, te;
+    double dtBound;
+    char* problem;
+    int bcLeft, bcRight, bcBottom, bcTop, bcFront, bcBack;
+    /* communication */
+    double **r, **e;
+    int levels, presmooth, postsmooth;
+    Comm* comm;
+} Solver;
+
+extern double solve(Solver* , double* , double* );
+extern void initSolver(Solver*, Discretization*, Parameter*);
+#endif

EnhancedSolver/3D-mpi/src/test.c

@@ -0,0 +1,129 @@
+/*
+ * 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.
+ */
+#include <stdio.h>
+#include <string.h>
+
+#include "test.h"
+
+#define X(v, i, j, k)                                                                    \
+    v[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)]
+
+void testInit(Solver* s)
+{
+    int imaxLocal = s->comm->imaxLocal;
+    int jmaxLocal = s->comm->jmaxLocal;
+    int kmaxLocal = s->comm->kmaxLocal;
+    int myrank    = s->comm->rank;
+    double* p     = s->p;
+    double* f     = s->f;
+    double* g     = s->g;
+    double* h     = s->h;
+
+    for (int k = 0; k < kmaxLocal + 2; k++) {
+        for (int j = 0; j < jmaxLocal + 2; j++) {
+            for (int i = 0; i < imaxLocal + 2; i++) {
+                X(p, i, j, k) = 10.0;
+                X(f, i, j, k) = myrank + 1.0;
+                X(g, i, j, k) = myrank + 1.0;
+                X(h, i, j, k) = myrank + 1.0;
+            }
+        }
+    }
+
+    for (int k = 1; k < kmaxLocal + 1; k++) {
+        for (int j = 1; j < jmaxLocal + 1; j++) {
+            for (int i = 1; i < imaxLocal + 1; i++) {
+                X(p, i, j, k) = myrank + 1.0;
+                X(f, i, j, k) = myrank + 1.0;
+                X(g, i, j, k) = myrank + 1.0;
+                X(h, i, j, k) = myrank + 1.0;
+            }
+        }
+    }
+}
+
+static char* direction2String(Direction dir)
+{
+    switch (dir) {
+    case LEFT:
+        return "left";
+        break;
+    case RIGHT:
+        return "right";
+        break;
+    case BOTTOM:
+        return "bottom";
+        break;
+    case TOP:
+        return "top";
+        break;
+    case FRONT:
+        return "front";
+        break;
+    case BACK:
+        return "back";
+        break;
+    case NDIRS:
+        return "ERROR";
+    default:
+        return "ERROR";
+    }
+}
+
+static void printPlane(Solver* s, double* a, int ymax, int xmax, Direction dir)
+{
+    int imaxLocal = s->comm->imaxLocal;
+    int jmaxLocal = s->comm->jmaxLocal;
+    int kmaxLocal = s->comm->kmaxLocal;
+    char filename[50];
+    snprintf(filename, 50, "halo-%s-r%d.txt", direction2String(dir), s->comm->rank);
+    FILE* fh = fopen(filename, "w");
+
+    for (int y = 0; y < ymax; y++) {
+        for (int x = 0; x < xmax; x++) {
+            switch (dir) {
+            case LEFT:
+                fprintf(fh, "%12.8f ", X(a, 0, x, y));
+                break;
+            case RIGHT:
+                fprintf(fh, "%12.8f ", X(a, imaxLocal + 1, x, y));
+                break;
+            case BOTTOM:
+                fprintf(fh, "%12.8f ", X(a, x, 0, y));
+                break;
+            case TOP:
+                fprintf(fh, "%12.8f ", X(a, x, jmaxLocal + 1, y));
+                break;
+            case FRONT:
+                fprintf(fh, "%12.8f ", X(a, x, y, 0));
+                break;
+            case BACK:
+                fprintf(fh, "%12.8f ", X(a, x, y, kmaxLocal + 1));
+                break;
+            case NDIRS:
+                printf("ERROR\n");
+                break;
+            }
+        }
+        fprintf(fh, "\n");
+    }
+    fclose(fh);
+}
+
+void testPrintHalo(Solver* s, double* a)
+{
+    int imaxLocal = s->comm->imaxLocal;
+    int jmaxLocal = s->comm->jmaxLocal;
+    int kmaxLocal = s->comm->kmaxLocal;
+
+    printPlane(s, a, kmaxLocal + 2, imaxLocal + 2, BOTTOM);
+    printPlane(s, a, kmaxLocal + 2, imaxLocal + 2, TOP);
+    printPlane(s, a, kmaxLocal + 2, jmaxLocal + 2, LEFT);
+    printPlane(s, a, kmaxLocal + 2, jmaxLocal + 2, RIGHT);
+    printPlane(s, a, jmaxLocal + 2, imaxLocal + 2, FRONT);
+    printPlane(s, a, jmaxLocal + 2, imaxLocal + 2, BACK);
+}

EnhancedSolver/3D-mpi/src/test.h

@@ -0,0 +1,13 @@
+/*
+ * 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 __TEST_H_
+#define __TEST_H_
+#include "solver.h"
+
+extern void testInit(Solver* s);
+extern void testPrintHalo(Solver* s, double* a);
+#endif // __TEST_H_

EnhancedSolver/3D-mpi/src/timing.c

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

EnhancedSolver/3D-mpi/src/timing.h

@@ -0,0 +1,13 @@
+/*
+ * 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 __TIMING_H_
+#define __TIMING_H_
+
+extern double getTimeStamp(void);
+extern double getTimeResolution(void);
+
+#endif // __TIMING_H_

EnhancedSolver/3D-mpi/src/util.h

@@ -0,0 +1,35 @@
+/*
+ * 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 __UTIL_H_
+#define __UTIL_H_
+#define HLINE                                                                            \
+    "----------------------------------------------------------------------------\n"
+
+#ifndef MIN
+#define MIN(x, y) ((x) < (y) ? (x) : (y))
+#endif
+#ifndef MAX
+#define MAX(x, y) ((x) > (y) ? (x) : (y))
+#endif
+#ifndef ABS
+#define ABS(a) ((a) >= 0 ? (a) : -(a))
+#endif
+
+#define P(i, j, k)   p[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)]
+#define F(i, j, k)   f[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)]
+#define G(i, j, k)   g[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)]
+#define H(i, j, k)   h[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)]
+#define U(i, j, k)   u[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)]
+#define V(i, j, k)   v[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)]
+#define W(i, j, k)   w[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)]
+#define S(i, j, k)   s[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)]
+#define E(i, j, k)   e[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)]
+#define R(i, j, k)   r[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)]
+#define OLD(i, j, k) old[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)]
+#define RHS(i, j, k) rhs[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)]
+
+#endif // __UTIL_H_

EnhancedSolver/3D-mpi/src/vtkWriter-mpi.c

@@ -0,0 +1,227 @@
+/*
+ * 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 <mpi.h>
+#include <stdbool.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "allocate.h"
+#include "comm.h"
+#include "vtkWriter.h"
+
+// reset fileview for output of string headers
+static void resetFileview(VtkOptions* o)
+{
+    MPI_Offset disp;
+    MPI_File_sync(o->fh);
+    MPI_Barrier(o->comm.comm);
+    MPI_File_get_size(o->fh, &disp);
+    MPI_File_set_view(o->fh, disp, MPI_CHAR, MPI_CHAR, "native", MPI_INFO_NULL);
+}
+
+static void writeVersion(VtkOptions* o)
+{
+    char header[50] = "# vtk DataFile Version 3.0\n";
+    // always overwrite exiting files
+    MPI_File_set_view(o->fh, 0, MPI_CHAR, MPI_CHAR, "native", MPI_INFO_NULL);
+
+    if (commIsMaster(&o->comm)) {
+        MPI_File_write(o->fh, header, (int)strlen(header), MPI_CHAR, MPI_STATUS_IGNORE);
+    }
+}
+
+static void writeHeader(VtkOptions* o)
+{
+    char header[400];
+    char* cursor = header;
+
+    cursor += sprintf(cursor, "PAMPI cfd solver output\n");
+    cursor += sprintf(cursor, "BINARY\n");
+
+    cursor += sprintf(cursor, "DATASET STRUCTURED_POINTS\n");
+    cursor += sprintf(cursor,
+        "DIMENSIONS %d %d %d\n",
+        o->grid->imax,
+        o->grid->jmax,
+        o->grid->kmax);
+    cursor += sprintf(cursor,
+        "ORIGIN %f %f %f\n",
+        o->grid->dx * 0.5,
+        o->grid->dy * 0.5,
+        o->grid->dz * 0.5);
+    cursor += sprintf(cursor, "SPACING %f %f %f\n", o->grid->dx, o->grid->dy, o->grid->dz);
+    cursor += sprintf(cursor,
+        "POINT_DATA %d\n",
+        o->grid->imax * o->grid->jmax * o->grid->kmax);
+
+    if (commIsMaster(&o->comm)) {
+        MPI_File_write(o->fh, header, (int)strlen(header), MPI_CHAR, MPI_STATUS_IGNORE);
+    }
+}
+
+void vtkOpen(VtkOptions* o, char* problem)
+{
+    char filename[50];
+    snprintf(filename, 50, "%s-p%d.vtk", problem, o->comm.size);
+    MPI_File_open(o->comm.comm,
+        filename,
+        MPI_MODE_WRONLY | MPI_MODE_CREATE,
+        MPI_INFO_NULL,
+        &o->fh);
+
+    if (commIsMaster(&o->comm)) {
+        printf("Writing VTK output for %s\n", problem);
+    }
+
+    writeVersion(o);
+    writeHeader(o);
+}
+
+void vtkScalar(VtkOptions* o, char* name, double* s)
+{
+    resetFileview(o);
+    if (commIsMaster(&o->comm)) printf("Register scalar %s\n", name);
+
+    char header[100];
+    char* cursor = header;
+
+    cursor += sprintf(cursor, "SCALARS %s double\n", name);
+
+    if (commIsMaster(&o->comm)) {
+        MPI_File_write(o->fh, header, (int)strlen(header), MPI_CHAR, MPI_STATUS_IGNORE);
+    }
+
+    int offsets[NDIMS];
+    commGetOffsets(&o->comm, offsets, o->grid->kmax, o->grid->jmax, o->grid->imax);
+
+    // set global view in file
+    MPI_Offset disp;
+    MPI_Datatype fileViewType;
+    MPI_File_sync(o->fh);
+    MPI_Barrier(o->comm.comm);
+    MPI_File_get_size(o->fh, &disp);
+
+    MPI_Type_create_subarray(NDIMS,
+        (int[NDIMS]) { o->grid->kmax, o->grid->jmax, o->grid->imax },
+        (int[NDIMS]) { o->comm.kmaxLocal, o->comm.jmaxLocal, o->comm.imaxLocal },
+        offsets,
+        MPI_ORDER_C,
+        MPI_DOUBLE,
+        &fileViewType);
+    MPI_Type_commit(&fileViewType);
+    MPI_File_set_view(o->fh, disp, MPI_DOUBLE, fileViewType, "external32", MPI_INFO_NULL);
+
+#ifdef VERBOSE
+    printf("Rank: %d, Disp: %lld,  Size(k,j,i): %d %d %d, Offset(k,j,i): %d %d %d\n",
+        o->comm.rank,
+        disp,
+        o->comm.kmaxLocal,
+        o->comm.jmaxLocal,
+        o->comm.imaxLocal,
+        offsets[KDIM],
+        offsets[JDIM],
+        offsets[IDIM]);
+#endif
+
+    // create local bulk type
+    MPI_Datatype bulkType;
+
+    MPI_Type_create_subarray(NDIMS,
+        (int[NDIMS]) { o->comm.kmaxLocal + 2,
+            o->comm.jmaxLocal + 2,
+            o->comm.imaxLocal + 2 }, // oldsizes
+        (int[NDIMS]) { o->comm.kmaxLocal,
+            o->comm.jmaxLocal,
+            o->comm.imaxLocal },  // newsizes
+        (int[NDIMS]) { 1, 1, 1 }, // offsets
+        MPI_ORDER_C,
+        MPI_DOUBLE,
+        &bulkType);
+    MPI_Type_commit(&bulkType);
+
+    MPI_File_write(o->fh, s, 1, bulkType, MPI_STATUS_IGNORE);
+    MPI_Type_free(&bulkType);
+    MPI_Type_free(&fileViewType);
+
+    // Binary segment must be terminated with newline character
+    resetFileview(o);
+    if (commIsMaster(&o->comm)) {
+        MPI_File_write(o->fh, "\n", 1, MPI_CHAR, MPI_STATUS_IGNORE);
+    }
+}
+
+#define G(v, i, j, k)                                                                    \
+    v[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)]
+
+void vtkVector(VtkOptions* o, char* name, VtkVector vec)
+{
+    int imaxLocal = o->comm.imaxLocal;
+    int jmaxLocal = o->comm.jmaxLocal;
+    int kmaxLocal = o->comm.kmaxLocal;
+
+    if (commIsMaster(&o->comm)) printf("Register vector %s\n", name);
+    const size_t MAX_HEADER = 100;
+
+    char* header = (char*)malloc(MAX_HEADER);
+    sprintf(header, "VECTORS %s double\n", name);
+
+    resetFileview(o);
+    if (commIsMaster(&o->comm)) {
+        MPI_File_write(o->fh, header, (int)strlen(header), MPI_CHAR, MPI_STATUS_IGNORE);
+    }
+
+    int offsets[NDIMS];
+    commGetOffsets(&o->comm, offsets, o->grid->kmax, o->grid->jmax, o->grid->imax);
+
+    // set global view in file
+    MPI_Offset disp;
+    MPI_Datatype fileViewType, vectorType;
+    MPI_File_sync(o->fh);
+    MPI_Barrier(o->comm.comm);
+    MPI_File_get_size(o->fh, &disp);
+
+    MPI_Type_contiguous(NDIMS, MPI_DOUBLE, &vectorType);
+    MPI_Type_commit(&vectorType);
+
+    MPI_Type_create_subarray(NDIMS,
+        (int[NDIMS]) { o->grid->kmax, o->grid->jmax, o->grid->imax },
+        (int[NDIMS]) { kmaxLocal, jmaxLocal, imaxLocal },
+        offsets,
+        MPI_ORDER_C,
+        vectorType,
+        &fileViewType);
+    MPI_Type_commit(&fileViewType);
+    MPI_File_set_view(o->fh, disp, MPI_DOUBLE, fileViewType, "external32", MPI_INFO_NULL);
+
+    size_t cnt  = imaxLocal * jmaxLocal * kmaxLocal;
+    double* tmp = allocate(64, cnt * NDIMS * sizeof(double));
+    int idx     = 0;
+
+    for (int k = 1; k < kmaxLocal + 1; k++) {
+        for (int j = 1; j < jmaxLocal + 1; j++) {
+            for (int i = 1; i < imaxLocal + 1; i++) {
+                tmp[idx++] = (G(vec.u, i, j, k) + G(vec.u, i - 1, j, k)) / 2.0;
+                tmp[idx++] = (G(vec.v, i, j, k) + G(vec.v, i, j - 1, k)) / 2.0;
+                tmp[idx++] = (G(vec.w, i, j, k) + G(vec.w, i, j, k - 1)) / 2.0;
+            }
+        }
+    }
+
+    if (commIsMaster(&o->comm)) printf("Write %d vectors\n", (int)cnt);
+    MPI_File_write(o->fh, tmp, cnt, vectorType, MPI_STATUS_IGNORE);
+    MPI_Type_free(&fileViewType);
+    MPI_Type_free(&vectorType);
+
+    // Binary segment must be terminated with newline character
+    resetFileview(o);
+    if (commIsMaster(&o->comm)) {
+        MPI_File_write(o->fh, "\n", 1, MPI_CHAR, MPI_STATUS_IGNORE);
+    }
+}
+
+void vtkClose(VtkOptions* o) { MPI_File_close(&o->fh); }

EnhancedSolver/3D-mpi/src/vtkWriter-seq.c

@@ -0,0 +1,147 @@
+/*
+ * 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 <stdbool.h>
+#include <stdio.h>
+#include <string.h>
+
+#include "vtkWriter.h"
+#define Z(v, i, j, k) v[(k)*imax * jmax + (j)*imax + (i)]
+
+static double floatSwap(double f)
+{
+    union {
+        double f;
+        char b[8];
+    } dat1, dat2;
+
+    dat1.f    = f;
+    dat2.b[0] = dat1.b[7];
+    dat2.b[1] = dat1.b[6];
+    dat2.b[2] = dat1.b[5];
+    dat2.b[3] = dat1.b[4];
+    dat2.b[4] = dat1.b[3];
+    dat2.b[5] = dat1.b[2];
+    dat2.b[6] = dat1.b[1];
+    dat2.b[7] = dat1.b[0];
+    return dat2.f;
+}
+
+static void writeHeader(VtkOptions* o)
+{
+    fprintf(o->fh, "# vtk DataFile Version 3.0\n");
+    fprintf(o->fh, "PAMPI cfd solver output\n");
+    if (o->fmt == ASCII) {
+        fprintf(o->fh, "ASCII\n");
+    } else if (o->fmt == BINARY) {
+        fprintf(o->fh, "BINARY\n");
+    }
+
+    fprintf(o->fh, "DATASET STRUCTURED_POINTS\n");
+    fprintf(o->fh, "DIMENSIONS %d %d %d\n", o->grid->imax, o->grid->jmax, o->grid->kmax);
+    fprintf(o->fh,
+        "ORIGIN %f %f %f\n",
+        o->grid->dx * 0.5,
+        o->grid->dy * 0.5,
+        o->grid->dz * 0.5);
+    fprintf(o->fh, "SPACING %f %f %f\n", o->grid->dx, o->grid->dy, o->grid->dz);
+    fprintf(o->fh, "POINT_DATA %d\n", o->grid->imax * o->grid->jmax * o->grid->kmax);
+}
+
+void vtkOpen(VtkOptions* o, char* problem)
+{
+    char filename[50];
+
+    snprintf(filename, 50, "%s.vtk", problem);
+    o->fh = fopen(filename, "w");
+    writeHeader(o);
+
+    printf("Writing VTK output for %s\n", problem);
+}
+
+static void writeScalar(VtkOptions* o, double* s)
+{
+    int imax = o->grid->imax;
+    int jmax = o->grid->jmax;
+    int kmax = o->grid->kmax;
+
+    for (int k = 0; k < kmax; k++) {
+        for (int j = 0; j < jmax; j++) {
+            for (int i = 0; i < imax; i++) {
+                if (o->fmt == ASCII) {
+                    fprintf(o->fh, "%f\n", Z(s, i, j, k));
+                } else if (o->fmt == BINARY) {
+                    fwrite((double[1]) { floatSwap(Z(s, i, j, k)) },
+                        sizeof(double),
+                        1,
+                        o->fh);
+                }
+            }
+        }
+    }
+    if (o->fmt == BINARY) fprintf(o->fh, "\n");
+}
+
+static bool isInitialized(FILE* ptr)
+{
+    if (ptr == NULL) {
+        printf("vtkWriter not initialize! Call vtkOpen first!\n");
+        return false;
+    }
+    return true;
+}
+
+void vtkScalar(VtkOptions* o, char* name, double* s)
+{
+    printf("Register scalar %s\n", name);
+    if (!isInitialized(o->fh)) return;
+    fprintf(o->fh, "SCALARS %s double 1\n", name);
+    fprintf(o->fh, "LOOKUP_TABLE default\n");
+    writeScalar(o, s);
+}
+
+static void writeVector(VtkOptions* o, VtkVector vec)
+{
+    int imax = o->grid->imax;
+    int jmax = o->grid->jmax;
+    int kmax = o->grid->kmax;
+
+    for (int k = 0; k < kmax; k++) {
+        for (int j = 0; j < jmax; j++) {
+            for (int i = 0; i < imax; i++) {
+                if (o->fmt == ASCII) {
+                    fprintf(o->fh,
+                        "%f %f %f\n",
+                        Z(vec.u, i, j, k),
+                        Z(vec.v, i, j, k),
+                        Z(vec.w, i, j, k));
+                } else if (o->fmt == BINARY) {
+                    fwrite((double[3]) { floatSwap(Z(vec.u, i, j, k)),
+                               floatSwap(Z(vec.v, i, j, k)),
+                               floatSwap(Z(vec.w, i, j, k)) },
+                        sizeof(double),
+                        3,
+                        o->fh);
+                }
+            }
+        }
+    }
+    if (o->fmt == BINARY) fprintf(o->fh, "\n");
+}
+
+void vtkVector(VtkOptions* o, char* name, VtkVector vec)
+{
+    printf("Register vector %s\n", name);
+    if (!isInitialized(o->fh)) return;
+    fprintf(o->fh, "VECTORS %s double\n", name);
+    writeVector(o, vec);
+}
+
+void vtkClose(VtkOptions* o)
+{
+    fclose(o->fh);
+    o->fh = NULL;
+}

EnhancedSolver/3D-mpi/src/vtkWriter.h

@@ -0,0 +1,35 @@
+/*
+ * 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 __VTKWRITER_H_
+#define __VTKWRITER_H_
+#include <stdio.h>
+
+#include "comm.h"
+#include "grid.h"
+
+typedef enum VtkFormat { ASCII = 0, BINARY } VtkFormat;
+
+typedef struct VtkOptions {
+    Grid* grid;
+#ifdef _VTK_WRITER_MPI
+    MPI_File fh;
+#else
+    FILE* fh;
+    VtkFormat fmt;
+#endif // _VTK_WRITER_MPI
+    Comm comm;
+} VtkOptions;
+
+typedef struct VtkVector {
+    double *u, *v, *w;
+} VtkVector;
+
+extern void vtkOpen(VtkOptions* opts, char* filename);
+extern void vtkVector(VtkOptions* opts, char* name, VtkVector vec);
+extern void vtkScalar(VtkOptions* opts, char* name, double* p);
+extern void vtkClose(VtkOptions* opts);
+#endif // __VTKWRITER_H_

EnhancedSolver/3D-seq/Makefile

@@ -0,0 +1,88 @@
+#=======================================================================================
+# 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.
+#=======================================================================================
+
+#CONFIGURE BUILD SYSTEM
+TARGET	   = exe-$(TAG)
+BUILD_DIR  = ./$(TAG)
+SRC_DIR    = ./src
+MAKE_DIR   = ./
+Q         ?= @
+
+#DO NOT EDIT BELOW
+include $(MAKE_DIR)/config.mk
+include $(MAKE_DIR)/include_$(TAG).mk
+INCLUDES  += -I$(SRC_DIR) -I$(BUILD_DIR)
+
+VPATH     = $(SRC_DIR)
+SRC       = $(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)
+
+${TARGET}: $(BUILD_DIR) $(OBJ)
+	$(info ===>  LINKING  $(TARGET))
+	$(Q)${LINKER} ${LFLAGS} -o $(TARGET) $(OBJ) $(LIBS)
+
+$(BUILD_DIR)/%.o:  %.c $(MAKE_DIR)/include_$(TAG).mk $(MAKE_DIR)/config.mk
+	$(info ===>  COMPILE  $@)
+	$(CC) -c $(CPPFLAGS) $(CFLAGS) $< -o $@
+	$(Q)$(GCC) $(CPPFLAGS) -MT $(@:.d=.o) -MM  $< > $(BUILD_DIR)/$*.d
+
+$(BUILD_DIR)/%.s:  %.c
+	$(info ===>  GENERATE ASM  $@)
+	$(CC) -S $(CPPFLAGS) $(CFLAGS) $< -o $@
+
+.PHONY: clean distclean tags info asm format
+
+vis:
+	$(info ===>  GENERATE VISUALIZATION)
+	@gnuplot -e "filename='residual.dat'" ./residual.plot
+
+vis_clean:
+	$(info ===>  CLEAN VISUALIZATION)
+	@rm -f *.dat
+	@rm -f vtk_files/*.vtk
+	@rm -f *.png
+	@rm -f *.vtk
+
+
+clean: vis_clean
+
+clean:
+	$(info ===>  CLEAN)
+	@rm -rf $(BUILD_DIR)
+	@rm -f tags
+
+distclean: clean
+	$(info ===>  DIST CLEAN)
+	@rm -f $(TARGET)
+	@rm -f *.dat
+	@rm -f *.png
+
+info:
+	$(info $(CFLAGS))
+	$(Q)$(CC) $(VERSION)
+
+asm:  $(BUILD_DIR) $(ASM)
+
+tags:
+	$(info ===>  GENERATE TAGS)
+	$(Q)ctags -R
+
+format:
+	@for src in $(SOURCES) ; do \
+		echo "Formatting $$src" ; \
+		clang-format -i $$src ; \
+	done
+	@echo "Done"
+
+$(BUILD_DIR):
+	@mkdir $(BUILD_DIR)
+
+-include $(OBJ:.o=.d)

EnhancedSolver/3D-seq/README.md

@@ -0,0 +1,78 @@
+# C source skeleton
+
+## Build
+
+1. Configure the toolchain and additional options in `config.mk`:
+```
+# Supported: GCC, CLANG, ICC
+TAG ?= GCC
+ENABLE_OPENMP ?= false
+
+OPTIONS +=  -DARRAY_ALIGNMENT=64
+#OPTIONS +=  -DVERBOSE
+#OPTIONS +=  -DVERBOSE_AFFINITY
+#OPTIONS +=  -DVERBOSE_DATASIZE
+#OPTIONS +=  -DVERBOSE_TIMER
+```
+
+The verbosity options enable detailed output about solver, affinity settings, allocation sizes and timer resolution.
+For debugging you may want to set the VERBOSE option:
+```
+# Supported: GCC, CLANG, ICC
+TAG ?= GCC
+ENABLE_OPENMP ?= false
+
+OPTIONS +=  -DARRAY_ALIGNMENT=64
+OPTIONS +=  -DVERBOSE
+#OPTIONS +=  -DVERBOSE_AFFINITY
+#OPTIONS +=  -DVERBOSE_DATASIZE
+#OPTIONS +=  -DVERBOSE_TIMER
+`
+
+2. Build with:
+```
+make
+```
+
+You can build multiple toolchains in the same directory, but notice that the Makefile is only acting on the one currently set.
+Intermediate build results are located in the `<TOOLCHAIN>` directory.
+
+To output the executed commands use:
+```
+make Q=
+```
+
+3. Clean up with:
+```
+make clean
+```
+to clean intermediate build results.
+
+```
+make distclean
+```
+to clean intermediate build results and binary.
+
+4. (Optional) Generate assembler:
+```
+make asm
+```
+The assembler files will also be located in the `<TOOLCHAIN>` directory.
+
+## Usage
+
+You have to provide a parameter file describing the problem you want to solve:
+```
+./exe-CLANG  dcavity.par
+```
+
+Examples are given in in dcavity (a lid driven cavity test case) and canal (simulating a empty canal).
+
+You can plot the resulting velocity and pressure fields using gnuplot:
+```
+gnuplot vector.plot
+```
+and for the pressure:
+```
+gnuplot surface.plot
+```

EnhancedSolver/3D-seq/backstep.par

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

EnhancedSolver/3D-seq/canal.par

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

EnhancedSolver/3D-seq/config.mk

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

EnhancedSolver/3D-seq/dcavity.par

@@ -0,0 +1,88 @@
+#==============================================================================
+#                              Driven Cavity
+#==============================================================================
+
+# Problem specific Data:
+# ---------------------
+
+name dcavity        # name of flow setup
+
+bcLeft    1			#  flags for boundary conditions
+bcRight   1			#  1 = no-slip      3 = outflow
+bcBottom  1			#  2 = free-slip    4 = periodic
+bcTop     1			#
+bcFront   1			#
+bcBack    1			#
+
+gx    0.0			# Body forces (e.g. gravity)
+gy    0.0			#
+gz    0.0			#
+
+re    1000.0		    # Reynolds number
+
+u_init    0.0		# initial value for velocity in x-direction
+v_init    0.0		# initial value for velocity in y-direction
+w_init    0.0		# initial value for velocity in z-direction
+p_init    0.0		# initial value for pressure
+
+# Geometry Data:
+# -------------
+
+xlength    1.0		# domain size in x-direction
+ylength    1.0		# domain size in y-direction
+zlength    1.0		# domain size in z-direction
+imax       80		# number of interior cells in x-direction
+jmax       80		# number of interior cells in y-direction
+kmax       80		# number of interior cells in z-direction
+
+# Time Data:
+# ---------
+
+te       10.0		# final time
+dt       0.02	    # time stepsize
+tau      0.5		# safety factor for time stepsize control (<0 constant delt)
+
+# Pressure Iteration Data:
+# -----------------------
+
+itermax  1000		# maximal number of pressure iteration in one time step
+eps      0.001		# stopping tolerance for pressure iteration
+rho      0.5
+omg      1.7		# relaxation parameter for SOR iteration
+gamma    0.9		# upwind differencing factor gamma
+
+# Multigrid data:
+# ---------
+
+levels        2         # Multigrid levels
+presmooth     20        # Pre-smoothning iterations
+postsmooth    5         # Post-smoothning iterations
+
+# Particle Tracing Data:
+# -----------------------
+
+numberOfParticles   200
+startTime           0
+injectTimePeriod    1.0
+writeTimePeriod     0.1
+
+x1                  0.0
+y1                  0.1
+z1                  0.0
+x2                  1.0
+y2                  0.1
+z2                  1.0
+
+# Obstacle Geometry Data:
+# -----------------------
+# Shape 0 disable, 1 Rectangle/Square, 2 Circle
+
+shape               0 
+xCenter             0.5
+yCenter             0.5
+zCenter             0.5
+xRectLength         0.25
+yRectLength         0.25
+zRectLength         0.25
+circleRadius        0.25
+#===============================================================================

EnhancedSolver/3D-seq/include_CLANG.mk

@@ -0,0 +1,19 @@
+CC   = clang
+GCC  = cc
+LINKER = $(CC)
+
+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
+LFLAGS   = $(OPENMP) -lm
+DEFINES  = -D_GNU_SOURCE
+INCLUDES =

EnhancedSolver/3D-seq/include_GCC.mk

@@ -0,0 +1,14 @@
+CC   = gcc
+GCC  = gcc
+LINKER = $(CC)
+
+ifeq ($(ENABLE_OPENMP),true)
+OPENMP   = -fopenmp
+endif
+
+VERSION  = --version
+CFLAGS   = -Ofast -ffreestanding -std=c99 $(OPENMP)
+LFLAGS   = $(OPENMP)
+DEFINES  = -D_GNU_SOURCE
+INCLUDES =
+LIBS     =

EnhancedSolver/3D-seq/include_ICC.mk

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

EnhancedSolver/3D-seq/karman.par

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

EnhancedSolver/3D-seq/residual.plot

@@ -0,0 +1,9 @@
+set terminal png size 1800,768 enhanced font ,12
+set output 'residual.png'
+set datafile separator whitespace
+set xlabel "Timestep"
+set ylabel "Residual"
+
+set logscale y 2
+
+plot 'residual.dat' using 1:2 title "Residual"

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

@@ -0,0 +1,38 @@
+/*
+ * 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.
+ */
+#include <errno.h>
+#include <stddef.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "allocate.h"
+
+void* allocate(size_t alignment, size_t bytesize)
+{
+    int errorCode;
+    void* ptr;
+
+    errorCode = posix_memalign(&ptr, alignment, bytesize);
+
+    if (errorCode) {
+        if (errorCode == EINVAL) {
+            fprintf(stderr, "Error: Alignment parameter is not a power of two\n");
+            exit(EXIT_FAILURE);
+        }
+        if (errorCode == ENOMEM) {
+            fprintf(stderr, "Error: Insufficient memory to fulfill the request\n");
+            exit(EXIT_FAILURE);
+        }
+    }
+
+    if (ptr == NULL) {
+        fprintf(stderr, "Error: posix_memalign failed!\n");
+        exit(EXIT_FAILURE);
+    }
+
+    return ptr;
+}

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

@@ -0,0 +1,13 @@
+/*
+ * 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 __ALLOCATE_H_
+#define __ALLOCATE_H_
+#include <stdlib.h>
+
+extern void* allocate(size_t alignment, size_t bytesize);
+
+#endif

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

@@ -0,0 +1,86 @@
+/*
+ * 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 OBJECTBOUNDARY {
+    FLUID = 0,
+    /* Front Corners */
+    FRONTTOPLEFTCORNER,
+    FRONTTOPRIGHTCORNER,
+    FRONTBOTTOMLEFTCORNER,
+    FRONTBOTTOMRIGHTCORNER,
+    /* Back Corners */
+    BACKTOPLEFTCORNER,
+    BACKTOPRIGHTCORNER,
+    BACKBOTTOMLEFTCORNER,
+    BACKBOTTOMRIGHTCORNER,
+    /* Faces */
+    FRONTFACE,
+    BACKFACE,
+    LEFTFACE,
+    RIGHTFACE,
+    TOPFACE,
+    BOTTOMFACE,
+    /* Front Lines remaining after Corners and Faces */
+    FRONTLEFTLINE,
+    FRONTRIGHTLINE,
+    FRONTTOPLINE,
+    FRONTBOTTOMLINE,
+    /* Bottom Lines remaining after Corners and Faces */
+    BACKLEFTLINE,
+    BACKRIGHTLINE,
+    BACKTOPLINE,
+    BACKBOTTOMLINE,
+    /* Mid Lines remaining after Corners and Faces */
+    MIDTOPLEFTLINE,
+    MIDTOPRIGHTLINE,
+    MIDBOTTOMLEFTLINE,
+    MIDBOTTOMRIGHTLINE,
+    /* Local where its an object but not a boundary */
+    OBSTACLE,
+    /*Ghost cells boundary */
+    OUTSIDEBOUNDARY
+};
+
+enum BC { NOSLIP = 1, SLIP, OUTFLOW, PERIODIC };
+
+
+enum SHAPE { NOSHAPE = 0, RECT, CIRCLE };
+
+typedef struct {
+    /* geometry and grid information */
+    Grid grid;
+    /* arrays */
+    double *p, *rhs;
+    double *f, *g, *h;
+    double *u, *v, *w;
+    double* s;
+    /* 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 setObjectBoundaryCondition(Discretization*);
+extern void setSpecialBoundaryCondition(Discretization*);
+extern void computeFG(Discretization*);
+extern void adaptUV(Discretization*);
+#endif

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

@@ -0,0 +1,16 @@
+/*
+ * 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 __GRID_H_
+#define __GRID_H_
+
+typedef struct {
+    double dx, dy, dz;
+    int imax, jmax, kmax;
+    double xlength, ylength, zlength;
+} Grid;
+
+#endif // __GRID_H_

EnhancedSolver/3D-seq/src/likwid-marker.h

@@ -0,0 +1,54 @@
+/*
+ * =======================================================================================
+ *
+ *      Author:   Jan Eitzinger (je), jan.eitzinger@fau.de
+ *      Copyright (c) 2020 RRZE, University Erlangen-Nuremberg
+ *
+ *      Permission is hereby granted, free of charge, to any person obtaining a copy
+ *      of this software and associated documentation files (the "Software"), to deal
+ *      in the Software without restriction, including without limitation the rights
+ *      to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ *      copies of the Software, and to permit persons to whom the Software is
+ *      furnished to do so, subject to the following conditions:
+ *
+ *      The above copyright notice and this permission notice shall be included in all
+ *      copies or substantial portions of the Software.
+ *
+ *      THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ *      IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ *      FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ *      AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ *      LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ *      OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+ *      SOFTWARE.
+ *
+ * =======================================================================================
+ */
+#ifndef LIKWID_MARKERS_H
+#define LIKWID_MARKERS_H
+
+#ifdef LIKWID_PERFMON
+#include <likwid.h>
+#define LIKWID_MARKER_INIT                likwid_markerInit()
+#define LIKWID_MARKER_THREADINIT          likwid_markerThreadInit()
+#define LIKWID_MARKER_SWITCH              likwid_markerNextGroup()
+#define LIKWID_MARKER_REGISTER(regionTag) likwid_markerRegisterRegion(regionTag)
+#define LIKWID_MARKER_START(regionTag)    likwid_markerStartRegion(regionTag)
+#define LIKWID_MARKER_STOP(regionTag)     likwid_markerStopRegion(regionTag)
+#define LIKWID_MARKER_CLOSE               likwid_markerClose()
+#define LIKWID_MARKER_RESET(regionTag)    likwid_markerResetRegion(regionTag)
+#define LIKWID_MARKER_GET(regionTag, nevents, events, time, count)                       \
+    likwid_markerGetRegion(regionTag, nevents, events, time, count)
+#else /* LIKWID_PERFMON */
+#define LIKWID_MARKER_INIT
+#define LIKWID_MARKER_THREADINIT
+#define LIKWID_MARKER_SWITCH
+#define LIKWID_MARKER_REGISTER(regionTag)
+#define LIKWID_MARKER_START(regionTag)
+#define LIKWID_MARKER_STOP(regionTag)
+#define LIKWID_MARKER_CLOSE
+#define LIKWID_MARKER_GET(regionTag, nevents, events, time, count)
+#define LIKWID_MARKER_RESET(regionTag)
+#endif /* LIKWID_PERFMON */
+
+#endif /*LIKWID_MARKERS_H*/

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

@@ -0,0 +1,157 @@
+/*
+ * 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.
+ */
+#include <stdio.h>
+#include <stdlib.h>
+#include <unistd.h>
+
+#include "allocate.h"
+#include "discretization.h"
+#include "parameter.h"
+#include "particletracing.h"
+#include "progress.h"
+#include "solver.h"
+#include "timing.h"
+#include "vtkWriter.h"
+
+#define G(v, i, j, k) v[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)]
+
+static void createBulkArrays(
+    Discretization* s, double* pg, double* ug, double* vg, double* wg)
+{
+    int imax = s->grid.imax;
+    int jmax = s->grid.jmax;
+    int kmax = s->grid.kmax;
+    int idx  = 0;
+
+    for (int k = 1; k < kmax + 1; k++) {
+        for (int j = 1; j < jmax + 1; j++) {
+            for (int i = 1; i < imax + 1; i++) {
+                pg[idx++] = G(s->p, i, j, k);
+            }
+        }
+    }
+
+    idx = 0;
+
+    for (int k = 1; k < kmax + 1; k++) {
+        for (int j = 1; j < jmax + 1; j++) {
+            for (int i = 1; i < imax + 1; i++) {
+                ug[idx++] = (G(s->u, i, j, k) + G(s->u, i - 1, j, k)) / 2.0;
+            }
+        }
+    }
+
+    idx = 0;
+
+    for (int k = 1; k < kmax + 1; k++) {
+        for (int j = 1; j < jmax + 1; j++) {
+            for (int i = 1; i < imax + 1; i++) {
+                vg[idx++] = (G(s->v, i, j, k) + G(s->v, i, j - 1, k)) / 2.0;
+            }
+        }
+    }
+
+    idx = 0;
+
+    for (int k = 1; k < kmax + 1; k++) {
+        for (int j = 1; j < jmax + 1; j++) {
+            for (int i = 1; i < imax + 1; i++) {
+                wg[idx++] = (G(s->w, i, j, k) + G(s->w, i, j, k - 1)) / 2.0;
+            }
+        }
+    }
+}
+
+int main(int argc, char** argv)
+{
+    double timeStart, timeStop;
+    Parameter p;
+    Discretization d;
+    Solver s;
+    ParticleTracer particletracer;
+    initParameter(&p);
+
+    FILE* fp;
+    fp = initResidualWriter();
+
+    if (argc != 2) {
+        printf("Usage: %s <configFile>\n", argv[0]);
+        exit(EXIT_SUCCESS);
+    }
+
+    readParameter(&p, argv[1]);
+    printParameter(&p);
+    initDiscretization(&d, &p);
+    initSolver(&s, &d, &p);
+    initParticleTracer(&particletracer, &p);
+    printParticleTracerParameters(&particletracer);
+
+#ifndef VERBOSE
+    initProgress(d.te);
+#endif
+
+    double tau = d.tau;
+    double te  = d.te;
+    double t   = 0.0;
+    int nt     = 0;
+    double res = 0.0;
+
+    timeStart = getTimeStamp();
+    while (t <= te) {
+        if (tau > 0.0) computeTimestep(&d);
+        setBoundaryConditions(&d);
+        setSpecialBoundaryCondition(&d);
+        setObjectBoundaryCondition(&d);
+
+        computeFG(&d);
+        computeRHS(&d);
+        if (nt % 100 == 0) normalizePressure(&d);
+        res = solve(&s, d.p, d.rhs);
+        adaptUV(&d);
+
+        trace(&particletracer, d.u, d.v, d.w, d.s, t);
+
+        writeResidual(fp, t, res);
+
+        t += d.dt;
+        nt++;
+
+#ifdef VERBOSE
+        printf("TIME %f , TIMESTEP %f\n", t, d.dt);
+#else
+        printProgress(t);
+#endif
+    }
+    timeStop = getTimeStamp();
+    stopProgress();
+    printf("Solution took %.2fs\n", timeStop - timeStart);
+
+    timeStart = getTimeStamp();
+    double *pg, *ug, *vg, *wg;
+
+    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);
+
+    fclose(fp);
+    freeParticles(&particletracer);
+    
+    createBulkArrays(&d, pg, ug, vg, wg);
+    VtkOptions opts = { .grid = d.grid };
+    vtkOpen(&opts, d.problem);
+    vtkScalar(&opts, "pressure", pg);
+    vtkVector(&opts, "velocity", (VtkVector) { ug, vg, wg });
+    vtkClose(&opts);
+
+    timeStop = getTimeStamp();
+    printf("Result output took %.2fs\n", timeStop - timeStart);
+
+    return EXIT_SUCCESS;
+}

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

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

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

@@ -0,0 +1,35 @@
+/*
+ * 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 __PARAMETER_H_
+#define __PARAMETER_H_
+
+typedef struct {
+    int imax, jmax, kmax;
+    double xlength, ylength, zlength;
+    int itermax, levels;
+    double eps, omg, rho;
+    double re, tau, gamma;
+    double te, dt;
+    double gx, gy, gz;
+    char* name;
+    int bcLeft, bcRight, bcBottom, bcTop, bcFront, bcBack;
+    double u_init, v_init, w_init, p_init;
+    int presmooth, postsmooth;
+
+    int numberOfParticles;
+    double startTime, injectTimePeriod, writeTimePeriod;
+
+    double x1, y1, z1, x2, y2, z2;
+
+    int shape;
+    double xCenter, yCenter, zCenter, xRectLength, yRectLength, zRectLength, circleRadius;
+} Parameter;
+
+void initParameter(Parameter*);
+void readParameter(Parameter*, const char*);
+void printParameter(Parameter*);
+#endif

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

@@ -0,0 +1,313 @@
+/*
+ * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
+ * All rights reserved. This file is part of nusif-solver.
+ * Use of this source code is governed by a MIT style
+ * license that can be found in the LICENSE file.
+ */
+#include <float.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "discretization.h"
+#include "vtkWriter.h"
+
+#define U(i, j, k) u[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)]
+#define V(i, j, k) v[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)]
+#define W(i, j, k) w[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)]
+#define S(i, j, k) s[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)]
+
+static int ts     = 0;
+unsigned int seed = 32767;
+void printParticles(ParticleTracer* particletracer)
+{
+    for (int i = 0; i < particletracer->totalParticles; ++i) {
+        printf("Particle position X : %.2f, Y : %.2f, flag : %d\n",
+            particletracer->particlePool[i].x,
+            particletracer->particlePool[i].y,
+            particletracer->particlePool[i].flag);
+    }
+}
+void injectParticles(ParticleTracer* particletracer, double* s)
+{
+
+    int imax = particletracer->imax;
+    int jmax = particletracer->jmax;
+    int kmax = particletracer->kmax;
+
+    for (int i = 0; i < particletracer->numberOfParticles; ++i) {
+
+        particletracer->particlePool[particletracer->pointer].x =
+            (((double)rand() / RAND_MAX) * (particletracer->x2 - particletracer->x1)) +
+            particletracer->x1;
+        particletracer->particlePool[particletracer->pointer].y =
+            (((double)rand() / RAND_MAX) * (particletracer->y2 - particletracer->y1)) +
+            particletracer->y1;
+        particletracer->particlePool[particletracer->pointer].z =
+            (((double)rand() / RAND_MAX) * (particletracer->z2 - particletracer->z1)) +
+            particletracer->z1;
+
+        int i = particletracer->particlePool[particletracer->pointer].x /
+                particletracer->dx;
+        int j = particletracer->particlePool[particletracer->pointer].y /
+                particletracer->dy;
+        int k = particletracer->particlePool[particletracer->pointer].z /
+                particletracer->dz;
+
+        if (S(i, j, k) == FLUID) {
+            particletracer->particlePool[particletracer->pointer].flag = true;
+            ++(particletracer->pointer);
+            ++(particletracer->totalParticles);
+        } else
+            particletracer->particlePool[particletracer->pointer].flag = false;
+    }
+}
+
+void advanceParticles(ParticleTracer* particletracer,
+    double* restrict u,
+    double* restrict v,
+    double* restrict w,
+    double* restrict s,
+    double time)
+{
+    int imax = particletracer->imax;
+    int jmax = particletracer->jmax;
+    int kmax = particletracer->kmax;
+
+    double dx = particletracer->dx;
+    double dy = particletracer->dy;
+    double dz = particletracer->dz;
+
+    double xlength = particletracer->xlength;
+    double ylength = particletracer->ylength;
+    double zlength = particletracer->zlength;
+
+    for (int i = 0; i < particletracer->totalParticles; ++i) {
+        if (particletracer->particlePool[i].flag == true) {
+            double x = particletracer->particlePool[i].x;
+            double y = particletracer->particlePool[i].y;
+            double z = particletracer->particlePool[i].z;
+
+            int iCoord = (int)(x / dx) + 1;
+            int jCoord = (int)((y + 0.5 * dy) / dy) + 1;
+            int kCoord = (int)((z + 0.5 * dz) / dz) + 1;
+
+            double x1 = (double)(iCoord - 1) * dx;
+            double y1 = ((double)(jCoord - 1) - 0.5) * dy;
+            double z1 = ((double)(kCoord - 1) - 0.5) * dz;
+
+            double x2 = (double)iCoord * dx;
+            double y2 = ((double)jCoord - 0.5) * dy;
+            double z2 = ((double)kCoord - 0.5) * dz;
+
+            double u_n =
+                (1.0 / (dx * dy * dz)) *
+                ((x2 - x) * (y2 - y) * (z2 - z) * U(iCoord - 1, jCoord - 1, kCoord - 1) +
+                    (x - x1) * (y2 - y) * (z2 - z) * U(iCoord, jCoord - 1, kCoord - 1) +
+                    (x2 - x) * (y - y1) * (z2 - z) * U(iCoord - 1, jCoord, kCoord - 1) +
+                    (x - x1) * (y - y1) * (z2 - z) * U(iCoord, jCoord, kCoord - 1) +
+                    (x2 - x) * (y2 - y) * (z - z1) * U(iCoord - 1, jCoord - 1, kCoord) +
+                    (x - x1) * (y2 - y) * (z - z1) * U(iCoord, jCoord - 1, kCoord) +
+                    (x2 - x) * (y - y1) * (z - z1) * U(iCoord - 1, jCoord, kCoord) +
+                    (x - x1) * (y - y1) * (z - z1) * U(iCoord, jCoord, kCoord));
+
+            double new_x                      = x + particletracer->dt * u_n;
+            particletracer->particlePool[i].x = new_x;
+
+            iCoord = (int)((x + 0.5 * dx) / dx) + 1;
+            jCoord = (int)(y / dy) + 1;
+            kCoord = (int)((z + 0.5 * dz) / dz) + 1;
+
+            x1 = ((double)(iCoord - 1) - 0.5) * dx;
+            y1 = (double)(jCoord - 1) * dy;
+            z1 = ((double)(kCoord - 1) - 0.5) * dz;
+
+            x2 = ((double)iCoord - 0.5) * dx;
+            y2 = (double)jCoord * dy;
+            z2 = ((double)kCoord - 0.5) * dz;
+
+            double v_n =
+                (1.0 / (dx * dy * dz)) *
+                ((x2 - x) * (y2 - y) * (z2 - z) * V(iCoord - 1, jCoord - 1, kCoord - 1) +
+                    (x - x1) * (y2 - y) * (z2 - z) * V(iCoord, jCoord - 1, kCoord - 1) +
+                    (x2 - x) * (y - y1) * (z2 - z) * V(iCoord - 1, jCoord, kCoord - 1) +
+                    (x - x1) * (y - y1) * (z2 - z) * V(iCoord, jCoord, kCoord - 1) +
+                    (x2 - x) * (y2 - y) * (z - z1) * V(iCoord - 1, jCoord - 1, kCoord) +
+                    (x - x1) * (y2 - y) * (z - z1) * V(iCoord, jCoord - 1, kCoord) +
+                    (x2 - x) * (y - y1) * (z - z1) * V(iCoord - 1, jCoord, kCoord) +
+                    (x - x1) * (y - y1) * (z - z1) * V(iCoord, jCoord, kCoord));
+
+            double new_y                      = y + particletracer->dt * v_n;
+            particletracer->particlePool[i].y = new_y;
+
+            iCoord = (int)((x + 0.5 * dx) / dx) + 1;
+            jCoord = (int)((y + 0.5 * dy) / dy) + 1;
+            kCoord = (int)(z / dz) + 1;
+
+            x1 = ((double)(iCoord - 1) - 0.5) * dx;
+            y1 = ((double)(jCoord - 1) - 0.5) * dy;
+            z1 = (double)(kCoord - 1) * dz;
+
+            x2 = ((double)iCoord - 0.5) * dx;
+            y2 = ((double)jCoord - 0.5) * dy;
+            z2 = (double)kCoord * dz;
+
+            double w_n =
+                (1.0 / (dx * dy * dz)) *
+                ((x2 - x) * (y2 - y) * (z2 - z) * W(iCoord - 1, jCoord - 1, kCoord - 1) +
+                    (x - x1) * (y2 - y) * (z2 - z) * W(iCoord, jCoord - 1, kCoord - 1) +
+                    (x2 - x) * (y - y1) * (z2 - z) * W(iCoord - 1, jCoord, kCoord - 1) +
+                    (x - x1) * (y - y1) * (z2 - z) * W(iCoord, jCoord, kCoord - 1) +
+                    (x2 - x) * (y2 - y) * (z - z1) * W(iCoord - 1, jCoord - 1, kCoord) +
+                    (x - x1) * (y2 - y) * (z - z1) * W(iCoord, jCoord - 1, kCoord) +
+                    (x2 - x) * (y - y1) * (z - z1) * W(iCoord - 1, jCoord, kCoord) +
+                    (x - x1) * (y - y1) * (z - z1) * W(iCoord, jCoord, kCoord));
+
+            double new_z                      = z + particletracer->dt * w_n;
+            particletracer->particlePool[i].z = new_z;
+
+            if (((new_x < 0.0) || (new_x > xlength) || (new_y < 0.0) ||
+                    (new_y > ylength) || (new_z < 0.0) || (new_z > zlength))) {
+                particletracer->particlePool[i].flag = false;
+                particletracer->removedParticles++;
+            }
+            int i_new = new_x / dx, j_new = new_y / dy, k_new = new_z / dz;
+
+            if (S(i_new, j_new, k_new) != FLUID) {
+                particletracer->particlePool[i].flag = false;
+                particletracer->removedParticles++;
+            }
+        }
+    }
+}
+
+void freeParticles(ParticleTracer* particletracer) { free(particletracer->particlePool); }
+
+void writeParticles(ParticleTracer* particletracer)
+{
+    VtkOptions opts = { .particletracer = particletracer };
+    compress(particletracer);
+
+    char filename[50];
+    snprintf(filename, 50, "vtk_files/particles%d.vtk", ts);
+    vtkOpenPT(&opts, filename, ts);
+    vtkParticle(&opts, "particle");
+    vtkClose(&opts);
+    ++ts;
+}
+
+void initParticleTracer(ParticleTracer* particletracer, Parameter* params)
+{
+    particletracer->numberOfParticles = params->numberOfParticles;
+    particletracer->startTime         = params->startTime;
+    particletracer->injectTimePeriod  = params->injectTimePeriod;
+    particletracer->writeTimePeriod   = params->writeTimePeriod;
+
+    particletracer->dt = params->dt;
+    particletracer->dx = params->xlength / params->imax;
+    particletracer->dy = params->ylength / params->jmax;
+    particletracer->dz = params->zlength / params->kmax;
+
+    particletracer->xlength = params->xlength;
+    particletracer->ylength = params->ylength;
+    particletracer->zlength = params->zlength;
+
+    particletracer->x1 = params->x1;
+    particletracer->y1 = params->y1;
+    particletracer->z1 = params->z1;
+    particletracer->x2 = params->x2;
+    particletracer->y2 = params->y2;
+    particletracer->z2 = params->z2;
+
+    particletracer->lastInjectTime = params->startTime;
+    particletracer->lastUpdateTime = params->startTime;
+    particletracer->lastWriteTime  = params->startTime;
+
+    particletracer->pointer          = 0;
+    particletracer->totalParticles   = 0;
+    particletracer->removedParticles = 0;
+
+    particletracer->imax = params->imax;
+    particletracer->jmax = params->jmax;
+    particletracer->kmax = params->kmax;
+
+    particletracer->estimatedNumParticles = ((params->te - params->startTime) + 2) *
+                                            params->numberOfParticles;
+
+    particletracer->particlePool = malloc(
+        sizeof(Particle) * particletracer->estimatedNumParticles);
+}
+
+void printParticleTracerParameters(ParticleTracer* particletracer)
+{
+    printf("Particle Tracing data:\n");
+    printf("\tNumber of particles : %d being injected for every period of %.2f\n",
+        particletracer->numberOfParticles,
+        particletracer->injectTimePeriod);
+    printf("\tstartTime : %.2f\n", particletracer->startTime);
+    printf("\t(Line along which the particles are to be injected) \n\tx1 : %.2f, y1 : "
+           "%.2f, z1 : %.2f, x2 : %.2f, y2 : %.2f, z2 : %.2f\n",
+        particletracer->x1,
+        particletracer->y1,
+        particletracer->z1,
+        particletracer->x2,
+        particletracer->y2,
+        particletracer->z2);
+    printf("\tPointer : %d, TotalParticles : %d\n",
+        particletracer->pointer,
+        particletracer->totalParticles);
+    printf("\tdt : %.2f, dx : %.2f, dy : %.2f,  dz : %.2f\n",
+        particletracer->dt,
+        particletracer->dx,
+        particletracer->dy,
+        particletracer->dz);
+}
+
+void trace(ParticleTracer* particletracer,
+    double* u,
+    double* v,
+    double* w,
+    double* s,
+    double time)
+{
+    if (time >= particletracer->startTime) {
+        // printParticles(particletracer);
+        if ((time - particletracer->lastInjectTime) >= particletracer->injectTimePeriod) {
+            injectParticles(particletracer, s);
+            particletracer->lastInjectTime = time;
+        }
+        if ((time - particletracer->lastWriteTime) >= particletracer->writeTimePeriod) {
+            writeParticles(particletracer);
+            particletracer->lastWriteTime = time;
+        }
+        advanceParticles(particletracer, u, v, w, s, time);
+        if (particletracer->removedParticles > (particletracer->totalParticles * 0.2)) {
+            compress(particletracer);
+        }
+        particletracer->lastUpdateTime = time;
+    }
+}
+
+void compress(ParticleTracer* particletracer)
+{
+    Particle* memPool = particletracer->particlePool;
+    Particle tempPool[particletracer->totalParticles];
+    int totalParticles = 0;
+
+    for (int i = 0; i < particletracer->totalParticles; ++i) {
+        if (memPool[i].flag == true) {
+            tempPool[totalParticles].x    = memPool[i].x;
+            tempPool[totalParticles].y    = memPool[i].y;
+            tempPool[totalParticles].z    = memPool[i].z;
+            tempPool[totalParticles].flag = memPool[i].flag;
+            ++totalParticles;
+        }
+    }
+
+    particletracer->totalParticles   = totalParticles;
+    particletracer->removedParticles = 0;
+    particletracer->pointer          = totalParticles + 1;
+
+    memcpy(particletracer->particlePool, tempPool, totalParticles * sizeof(Particle));
+}

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

@@ -0,0 +1,46 @@
+/*
+ * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
+ * All rights reserved. This file is part of nusif-solver.
+ * Use of this source code is governed by a MIT style
+ * license that can be found in the LICENSE file.
+ */
+#ifndef __PARTICLETRACING_H_
+#define __PARTICLETRACING_H_
+#include "allocate.h"
+#include "parameter.h"
+
+#include <stdbool.h>
+
+typedef enum COORD { X = 0, Y, NCOORD } COORD;
+typedef struct {
+    double x, y, z;
+    bool flag;
+} Particle;
+
+typedef struct {
+    int numberOfParticles, totalParticles;
+    double startTime, injectTimePeriod, writeTimePeriod, lastInjectTime, lastUpdateTime,
+        lastWriteTime;
+
+    int estimatedNumParticles, activeParticles, removedParticles;
+
+    double dx, dy, dz, dt;
+    Particle* particlePool;
+
+    int pointer;
+
+    double imax, jmax, kmax, xlength, ylength, zlength;
+
+    double x1, y1, x2, y2, z1, z2;
+} ParticleTracer;
+
+extern void initParticleTracer(ParticleTracer*, Parameter*);
+extern void injectParticles(ParticleTracer*, double*);
+extern void advanceParticles(ParticleTracer*, double*, double*, double*, double*, double);
+extern void freeParticles(ParticleTracer*);
+extern void writeParticles(ParticleTracer*);
+extern void printParticleTracerParameters(ParticleTracer*);
+extern void printParticles(ParticleTracer*);
+extern void compress(ParticleTracer*);
+extern void trace(ParticleTracer*, double*, double*, double*, double*, double);
+#endif

EnhancedSolver/3D-seq/src/progress.c

@@ -0,0 +1,70 @@
+/*
+ * 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 <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "progress.h"
+
+static double _end;
+static int _current;
+
+void initProgress(double end)
+{
+    _end     = end;
+    _current = 0;
+
+    printf("[          ]");
+    fflush(stdout);
+}
+
+void printProgress(double current)
+{
+    int new = (int)rint((current / _end) * 10.0);
+
+    if (new > _current) {
+        char progress[11];
+        _current    = new;
+        progress[0] = 0;
+
+        for (int i = 0; i < 10; i++) {
+            if (i < _current) {
+                sprintf(progress + strlen(progress), "#");
+            } else {
+                sprintf(progress + strlen(progress), " ");
+            }
+        }
+        printf("\r[%s]", progress);
+    }
+    fflush(stdout);
+}
+
+void stopProgress()
+{
+    printf("\n");
+    fflush(stdout);
+}
+
+FILE* initResidualWriter()
+{
+    FILE* fp;
+    fp = fopen("residual.dat", "w");
+
+    if (fp == NULL) {
+        printf("Error!\n");
+        exit(EXIT_FAILURE);
+    }
+
+    return fp;
+
+}
+
+void writeResidual(FILE* fp, double ts, double res)
+{
+    fprintf(fp, "%f, %f\n", ts, res);
+}

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

@@ -0,0 +1,15 @@
+/*
+ * 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 __PROGRESS_H_
+#define __PROGRESS_H_
+
+extern void initProgress(double);
+extern void printProgress(double);
+extern void stopProgress(void);
+extern FILE* initResidualWriter(void);
+extern void writeResidual(FILE*, double, double);
+#endif

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

@@ -0,0 +1,304 @@
+/*
+ * 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) {
+
+                    P(i, j, k) -=
+                        factor *
+                        (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));
+                }
+                isw = 3 - isw;
+            }
+            jsw = 3 - jsw;
+        }
+        ksw = 3 - ksw;
+    }
+}
+
+static double calculateResidual(
+    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));
+
+                    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
+    if (level == COARSEST_LEVEL) {
+        for (int i = 0; i < 5; i++) {
+            smooth(s, p, rhs, level, imax, jmax, kmax);
+        }
+        return res;
+    }
+
+    // pre-smoothing
+    for (int i = 0; i < s->presmooth; i++) {
+        smooth(s, p, rhs, level, imax, jmax, kmax);
+        if (level == FINEST_LEVEL) setBoundaryCondition(p, imax, jmax, kmax);
+    }
+
+    res = calculateResidual(s, p, rhs, level, imax, jmax, kmax);
+
+    // restrict
+    restrictMG(s, level, imax, jmax, kmax);
+
+    // MGSolver on residual and error.
+    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 < s->postsmooth; i++) {
+        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;
+    s->presmooth  = p->presmooth;
+    s->postsmooth = p->postsmooth;
+
+    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;
+        }
+    }
+}
+
+double 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
+
+return res;
+}

EnhancedSolver/3D-seq/src/solver-rb.c

@@ -0,0 +1,101 @@
+/*
+ * 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;
+}
+
+double 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
+
+return res;
+}

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

@@ -0,0 +1,27 @@
+/*
+ * 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 __SOLVER_H_
+#define __SOLVER_H_
+#include "discretization.h"
+#include "grid.h"
+#include "parameter.h"
+
+typedef struct {
+    /* geometry and grid information */
+    Grid* grid;
+    /* parameters */
+    double eps, omega, rho;
+    int itermax;
+    int levels;
+    double **r, **e;
+    int presmooth, postsmooth;
+} Solver;
+
+extern void initSolver(Solver*, Discretization*, Parameter*);
+extern double solve(Solver*, double*, double*);
+
+#endif

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

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

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

@@ -0,0 +1,13 @@
+/*
+ * 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 __TIMING_H_
+#define __TIMING_H_
+
+extern double getTimeStamp(void);
+extern double getTimeResolution(void);
+
+#endif // __TIMING_H_

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

@@ -0,0 +1,32 @@
+/*
+ * 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 __UTIL_H_
+#define __UTIL_H_
+#define HLINE                                                                            \
+    "----------------------------------------------------------------------------\n"
+
+#ifndef MIN
+#define MIN(x, y) ((x) < (y) ? (x) : (y))
+#endif
+#ifndef MAX
+#define MAX(x, y) ((x) > (y) ? (x) : (y))
+#endif
+#ifndef ABS
+#define ABS(a) ((a) >= 0 ? (a) : -(a))
+#endif
+
+#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 S(i, j, k)   s[(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_