Merge 2D mpi versions. Cleanup.

BasicSolver/2D-mpi/Makefile

@@ -1,5 +1,5 @@
 #=======================================================================================
-# Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
+# 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.
@@ -18,9 +18,10 @@ include $(MAKE_DIR)/include_$(TAG).mk
 INCLUDES  += -I$(SRC_DIR) -I$(BUILD_DIR)
 
 VPATH     = $(SRC_DIR)
-SRC       = $(wildcard $(SRC_DIR)/*.c)
+SRC       = $(filter-out $(wildcard $(SRC_DIR)/*-*.c),$(wildcard $(SRC_DIR)/*.c))
 ASM       = $(patsubst $(SRC_DIR)/%.c, $(BUILD_DIR)/%.s, $(SRC))
 OBJ       = $(patsubst $(SRC_DIR)/%.c, $(BUILD_DIR)/%.o, $(SRC))
+OBJ      += $(BUILD_DIR)/comm-$(COMM_TYPE).o
 SOURCES   = $(SRC) $(wildcard $(SRC_DIR)/*.h)
 CPPFLAGS := $(CPPFLAGS) $(DEFINES) $(OPTIONS) $(INCLUDES)
 

BasicSolver/2D-mpi/config.mk

@@ -1,6 +1,8 @@
 # Supported: GCC, CLANG, ICC
 TAG ?= CLANG
+ENABLE_MPI ?= true
 ENABLE_OPENMP ?= false
+COMM_TYPE ?= v3
 
 #Feature options
 OPTIONS +=  -DARRAY_ALIGNMENT=64

BasicSolver/2D-mpi/include_CLANG.mk

@@ -1,4 +1,10 @@
+ifeq ($(ENABLE_MPI),true)
 CC   = mpicc
+DEFINES  = -D_MPI
+else
+CC   = cc
+endif
+
 GCC  = cc
 LINKER = $(CC)
 
@@ -9,8 +15,7 @@ LIBS     = # -lomp
 endif
 
 VERSION  = --version
-CFLAGS   = -Ofast -std=c99 $(OPENMP)
-#CFLAGS   = -Ofast -fnt-store=aggressive  -std=c99 $(OPENMP) #AMD CLANG
-LFLAGS   = $(OPENMP)
-DEFINES  = -D_GNU_SOURCE# -DDEBUG
-INCLUDES = -I/usr/local/include
+CFLAGS   = -Ofast -std=c17
+LFLAGS   = $(OPENMP) -lm
+DEFINES  += -D_GNU_SOURCE# -DDEBUG
+INCLUDES = -I/opt/homebrew/include

BasicSolver/2D-mpi/include_GCC.mk

@@ -1,4 +1,10 @@
+ifeq ($(ENABLE_MPI),true)
+CC   = mpicc
+DEFINES  = -D_MPI
+else
 CC   = gcc
+endif
+
 GCC  = gcc
 LINKER = $(CC)
 
@@ -9,6 +15,6 @@ endif
 VERSION  = --version
 CFLAGS   = -Ofast -ffreestanding -std=c99 $(OPENMP)
 LFLAGS   = $(OPENMP)
-DEFINES  = -D_GNU_SOURCE
+DEFINES  += -D_GNU_SOURCE
 INCLUDES =
 LIBS     =

BasicSolver/2D-mpi/include_ICC.mk

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

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

@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
+ * 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.

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

@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
+ * 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.

BasicSolver/2D-mpi/src/allocate.c

@@ -1,14 +1,17 @@
 /*
- * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
+ * 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>
 
-void* allocate(int alignment, size_t bytesize)
+#include "allocate.h"
+
+void* allocate(size_t alignment, size_t bytesize)
 {
     int errorCode;
     void* ptr;

BasicSolver/2D-mpi/src/allocate.h

@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
+ * 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.
@@ -8,6 +8,6 @@
 #define __ALLOCATE_H_
 #include <stdlib.h>
 
-extern void* allocate(int alignment, size_t bytesize);
+extern void* allocate(size_t alignment, size_t bytesize);
 
 #endif

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

@@ -0,0 +1,247 @@
+/*
+ * 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 "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 int sum(int* sizes, int position)
+{
+    int sum = 0;
+
+    for (int i = 0; i < position; i++) {
+        sum += sizes[i];
+    }
+
+    return sum;
+}
+
+static void gatherArray(
+    Comm* c, int cnt, int* rcvCounts, int* displs, double* src, double* dst)
+{
+    double* sendbuffer = src + (c->imaxLocal + 2);
+
+    if (c->rank == 0) {
+        sendbuffer = src;
+    }
+
+    MPI_Gatherv(sendbuffer,
+        cnt,
+        MPI_DOUBLE,
+        dst,
+        rcvCounts,
+        displs,
+        MPI_DOUBLE,
+        0,
+        MPI_COMM_WORLD);
+}
+#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, int direction)
+{
+#if defined(_MPI)
+    switch (direction) {
+    case LEFT:
+        return 1;
+        break;
+    case RIGHT:
+        return 1;
+        break;
+    case BOTTOM:
+        return c->rank == 0;
+        break;
+    case TOP:
+        return c->rank == (c->size - 1);
+        break;
+    }
+#endif
+
+    return 1;
+}
+
+void commExchange(Comm* c, double* grid)
+{
+#if defined(_MPI)
+    MPI_Request requests[4] = { MPI_REQUEST_NULL,
+        MPI_REQUEST_NULL,
+        MPI_REQUEST_NULL,
+        MPI_REQUEST_NULL };
+
+    /* exchange ghost cells with top neighbor */
+    if (c->rank + 1 < c->size) {
+        int top     = c->rank + 1;
+        double* src = grid + (c->jmaxLocal) * (c->imaxLocal + 2) + 1;
+        double* dst = grid + (c->jmaxLocal + 1) * (c->imaxLocal + 2) + 1;
+
+        MPI_Isend(src, c->imaxLocal, MPI_DOUBLE, top, 1, MPI_COMM_WORLD, &requests[0]);
+        MPI_Irecv(dst, c->imaxLocal, MPI_DOUBLE, top, 2, MPI_COMM_WORLD, &requests[1]);
+    }
+
+    /* exchange ghost cells with bottom neighbor */
+    if (c->rank > 0) {
+        int bottom  = c->rank - 1;
+        double* src = grid + (c->imaxLocal + 2) + 1;
+        double* dst = grid + 1;
+
+        MPI_Isend(src, c->imaxLocal, MPI_DOUBLE, bottom, 2, MPI_COMM_WORLD, &requests[2]);
+        MPI_Irecv(dst, c->imaxLocal, MPI_DOUBLE, bottom, 1, MPI_COMM_WORLD, &requests[3]);
+    }
+
+    MPI_Waitall(4, requests, MPI_STATUSES_IGNORE);
+#endif
+}
+
+void commShift(Comm* c, double* f, double* g)
+{
+#if defined(_MPI)
+    MPI_Request requests[2] = { MPI_REQUEST_NULL, MPI_REQUEST_NULL };
+
+    /* shift G */
+    /* receive ghost cells from bottom neighbor */
+    if (c->rank > 0) {
+        int bottom = c->rank - 1;
+        MPI_Irecv(g + 1,
+            c->imaxLocal,
+            MPI_DOUBLE,
+            bottom,
+            0,
+            MPI_COMM_WORLD,
+            &requests[0]);
+    }
+
+    if (c->rank + 1 < c->size) {
+        int top     = c->rank + 1;
+        double* buf = g + (c->jmaxLocal) * (c->imaxLocal + 2) + 1;
+        /* send ghost cells to top neighbor */
+        MPI_Isend(buf, c->imaxLocal, MPI_DOUBLE, top, 0, MPI_COMM_WORLD, &requests[1]);
+    }
+
+    MPI_Waitall(2, requests, MPI_STATUSES_IGNORE);
+#endif
+}
+
+// TODO Merge with seq
+void commCollectResult(Comm* c,
+    double* ug,
+    double* vg,
+    double* pg,
+    double* u,
+    double* v,
+    double* p,
+    int jmax,
+    int imax)
+{
+    int *rcvCounts, *displs;
+    int cnt = c->jmaxLocal * (imax + 2);
+
+    if (c->rank == 0) {
+        rcvCounts = (int*)malloc(c->size * sizeof(int));
+        displs    = (int*)malloc(c->size * sizeof(int));
+    }
+
+    if (c->rank == 0 && c->size == 1) {
+        cnt = (c->jmaxLocal + 2) * (imax + 2);
+    } else if (c->rank == 0 || c->rank == (c->size - 1)) {
+        cnt = (c->jmaxLocal + 1) * (imax + 2);
+    }
+
+    MPI_Gather(&cnt, 1, MPI_INTEGER, rcvCounts, 1, MPI_INTEGER, 0, MPI_COMM_WORLD);
+
+    if (c->rank == 0) {
+        displs[0]  = 0;
+        int cursor = rcvCounts[0];
+
+        for (int i = 1; i < c->size; i++) {
+            displs[i] = cursor;
+            cursor += rcvCounts[i];
+        }
+    }
+
+    gatherArray(c, cnt, rcvCounts, displs, p, pg);
+    gatherArray(c, cnt, rcvCounts, displs, u, ug);
+    gatherArray(c, cnt, rcvCounts, displs, v, vg);
+}
+
+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 (bottom, top, left, right): %d %d, %d, %d\n",
+                c->neighbours[BOTTOM],
+                c->neighbours[TOP],
+                c->neighbours[LEFT],
+                c->neighbours[RIGHT]);
+            printf("\tCoordinates (j,i) %d %d\n", c->coords[JDIM], c->coords[IDIM]);
+            printf("\tLocal domain size (j,i) %dx%d\n", 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 jmax, int imax)
+{
+#if defined(_MPI)
+    c->imaxLocal = imax;
+    c->jmaxLocal = sizeOfRank(c->rank, c->size, jmax);
+#else
+    c->imaxLocal = imax;
+    c->jmaxLocal = jmax;
+#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
+}

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

@@ -0,0 +1,355 @@
+/*
+ * 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 "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, int direction, int layer)
+{
+    MPI_Datatype type;
+    size_t dblsize = sizeof(double);
+    int imaxLocal  = c->imaxLocal;
+    int jmaxLocal  = c->jmaxLocal;
+    int sizes[NDIMS];
+    int subSizes[NDIMS];
+    int starts[NDIMS];
+    int offset = 0;
+}
+
+static void assembleResult(Comm* c,
+    double* src,
+    double* dst,
+    int imaxLocal[],
+    int jmaxLocal[],
+    int offset[],
+    int jmax,
+    int imax)
+{
+    MPI_Request* requests;
+    int numRequests = 1;
+
+    if (c->rank == 0) {
+        numRequests = c->size + 1;
+    } else {
+        numRequests = 1;
+    }
+
+    requests = (MPI_Request*)malloc(numRequests * sizeof(MPI_Request));
+
+    /* all ranks send their bulk array */
+    MPI_Datatype bulkType;
+    int oldSizes[NDIMS] = { c->jmaxLocal + 2, c->imaxLocal + 2 };
+    int newSizes[NDIMS] = { c->jmaxLocal, c->imaxLocal };
+    int starts[NDIMS]   = { 1, 1 };
+    MPI_Type_create_subarray(NDIMS,
+        oldSizes,
+        newSizes,
+        starts,
+        MPI_ORDER_C,
+        MPI_DOUBLE,
+        &bulkType);
+    MPI_Type_commit(&bulkType);
+
+    MPI_Isend(src, 1, bulkType, 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] = { jmax, imax };
+            int newSizes[NDIMS] = { jmaxLocal[i], imaxLocal[i] };
+            int starts[NDIMS]   = { 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_Waitall(numRequests, requests, MPI_STATUSES_IGNORE);
+}
+
+static int sum(int* sizes, int position)
+{
+    int sum = 0;
+
+    for (int i = 0; i < position; i++) {
+        sum += sizes[i];
+    }
+
+    return sum;
+}
+#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, int direction)
+{
+#if defined(_MPI)
+    switch (direction) {
+    case LEFT:
+        return c->coords[IDIM] == 0;
+        break;
+    case RIGHT:
+        return c->coords[IDIM] == (c->dims[IDIM] - 1);
+        break;
+    case BOTTOM:
+        return c->coords[JDIM] == 0;
+        break;
+    case TOP:
+        return c->coords[JDIM] == (c->dims[JDIM] - 1);
+        break;
+    }
+#endif
+
+    return 1;
+}
+
+void commExchange(Comm* c, double* grid)
+{
+#if defined(_MPI)
+    double* buf[8];
+    MPI_Request requests[8];
+    for (int i = 0; i < 8; i++)
+        requests[i] = MPI_REQUEST_NULL;
+
+    buf[0] = grid + 1;                                           // recv bottom
+    buf[1] = grid + (c->imaxLocal + 2) + 1;                      // send bottom
+    buf[2] = grid + (c->jmaxLocal + 1) * (c->imaxLocal + 2) + 1; // recv top
+    buf[3] = grid + (c->jmaxLocal) * (c->imaxLocal + 2) + 1;     // send top
+    buf[4] = grid + (c->imaxLocal + 2);                          // recv left
+    buf[5] = grid + (c->imaxLocal + 2) + 1;                      // send left
+    buf[6] = grid + (c->imaxLocal + 2) + (c->imaxLocal + 1);     // recv right
+    buf[7] = grid + (c->imaxLocal + 2) + (c->imaxLocal);         // send right
+
+    for (int i = 2; i < 4; i++) {
+        int tag = 0;
+        if (c->neighbours[i] != MPI_PROC_NULL) {
+            tag = c->neighbours[i];
+        }
+        /* exchange ghost cells with bottom/top neighbor */
+        MPI_Irecv(buf[i * 2],
+            1,
+            c->rbufferTypes[0],
+            c->neighbours[i],
+            tag,
+            c->comm,
+            &requests[i * 2]);
+        MPI_Isend(buf[(i * 2) + 1],
+            1,
+            c->rbufferTypes[0],
+            c->neighbours[i],
+            c->rank,
+            c->comm,
+            &requests[i * 2 + 1]);
+    }
+    for (int i = 0; i < 2; i++) {
+        int tag = 0;
+        if (c->neighbours[i] != MPI_PROC_NULL) {
+            tag = c->neighbours[i];
+        }
+        /* exchange ghost cells with left/right neighbor */
+        MPI_Irecv(buf[i * 2 + 4],
+            1,
+            c->sbufferTypes[0],
+            c->neighbours[i],
+            tag,
+            c->comm,
+            &requests[i * 2 + 4]);
+        MPI_Isend(buf[i * 2 + 5],
+            1,
+            c->sbufferTypes[0],
+            c->neighbours[i],
+            c->rank,
+            c->comm,
+            &requests[(i * 2) + 5]);
+    }
+
+    MPI_Waitall(8, requests, MPI_STATUSES_IGNORE);
+#endif
+}
+
+void commShift(Comm* c, double* f, double* g)
+{
+#if defined(_MPI)
+    MPI_Request requests[4] = { MPI_REQUEST_NULL,
+        MPI_REQUEST_NULL,
+        MPI_REQUEST_NULL,
+        MPI_REQUEST_NULL };
+
+    /* shift G */
+    double* buf = g + 1;
+    /* receive ghost cells from bottom neighbor */
+    MPI_Irecv(buf,
+        1,
+        c->rbufferTypes[0],
+        c->neighbours[BOTTOM],
+        0,
+        c->comm,
+        &requests[0]);
+
+    buf = g + (c->jmaxLocal) * (c->imaxLocal + 2) + 1;
+    /* send ghost cells to top neighbor */
+    MPI_Isend(buf, 1, c->rbufferTypes[0], c->neighbours[TOP], 0, c->comm, &requests[1]);
+
+    /* shift F */
+    buf = f + (c->imaxLocal + 2);
+    /* receive ghost cells from left neighbor */
+    MPI_Irecv(buf, 1, c->sbufferTypes[0], c->neighbours[LEFT], 1, c->comm, &requests[2]);
+
+    buf = f + (c->imaxLocal + 2) + (c->imaxLocal);
+    /* send ghost cells to right neighbor */
+    MPI_Isend(buf, 1, c->sbufferTypes[0], c->neighbours[RIGHT], 1, c->comm, &requests[3]);
+
+    MPI_Waitall(4, requests, MPI_STATUSES_IGNORE);
+#endif
+}
+
+void commCollectResult(Comm* c,
+    double* ug,
+    double* vg,
+    double* pg,
+    double* u,
+    double* v,
+    double* p,
+    int jmax,
+    int imax)
+{
+    int offset[c->size * NDIMS];
+    int imaxLocal[c->size];
+    int jmaxLocal[c->size];
+
+    MPI_Gather(&c->imaxLocal, 1, MPI_INT, imaxLocal, 1, MPI_INT, 0, MPI_COMM_WORLD);
+    MPI_Gather(&c->jmaxLocal, 1, MPI_INT, jmaxLocal, 1, MPI_INT, 0, MPI_COMM_WORLD);
+
+    if (c->rank == 0) {
+        for (int i = 0; i < c->size; i++) {
+            int coords[NDIMS];
+            MPI_Cart_coords(c->comm, i, NDIMS, coords);
+            offset[i * NDIMS + IDIM] = sum(imaxLocal, coords[IDIM]);
+            offset[i * NDIMS + JDIM] = sum(jmaxLocal, coords[JDIM]);
+            printf("Rank: %d, Coords(j,i): %d %d, Size(j,i): %d %d "
+                   "Offset(j,i): %d %d\n",
+                i,
+                coords[JDIM],
+                coords[IDIM],
+                jmaxLocal[i],
+                imaxLocal[i],
+                offset[i * NDIMS + JDIM],
+                offset[i * NDIMS + IDIM]);
+        }
+    }
+
+    /* collect P */
+    assembleResult(c, p, pg, imaxLocal, jmaxLocal, offset, jmax, imax);
+
+    /* collect U */
+    assembleResult(c, u, ug, imaxLocal, jmaxLocal, offset, jmax, imax);
+
+    /* collect V */
+    assembleResult(c, v, vg, imaxLocal, jmaxLocal, offset, jmax, imax);
+}
+
+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 (bottom, top, left, right): %d %d, %d, %d\n",
+                c->neighbours[BOTTOM],
+                c->neighbours[TOP],
+                c->neighbours[LEFT],
+                c->neighbours[RIGHT]);
+            printf("\tCoordinates (j,i) %d %d\n", c->coords[JDIM], c->coords[IDIM]);
+            printf("\tLocal domain size (j,i) %dx%d\n", 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 jmax, int imax)
+{
+#if defined(_MPI)
+    int dims[NDIMS]    = { 0, 0 };
+    int periods[NDIMS] = { 0, 0 };
+    MPI_Dims_create(c->size, NDIMS, dims);
+    MPI_Cart_create(MPI_COMM_WORLD, NDIMS, dims, periods, 0, &c->comm);
+    MPI_Cart_shift(c->comm, IDIM, 1, &c->neighbours[LEFT], &c->neighbours[RIGHT]);
+    MPI_Cart_shift(c->comm, JDIM, 1, &c->neighbours[BOTTOM], &c->neighbours[TOP]);
+    MPI_Cart_get(c->comm, NDIMS, c->dims, periods, c->coords);
+
+    c->imaxLocal = sizeOfRank(c->rank, dims[IDIM], imax);
+    c->jmaxLocal = sizeOfRank(c->rank, dims[JDIM], jmax);
+
+    MPI_Type_contiguous(c->imaxLocal, MPI_DOUBLE, &c->rbufferTypes[0]);
+    MPI_Type_commit(&c->rbufferTypes[0]);
+
+    MPI_Type_vector(c->jmaxLocal, 1, c->imaxLocal + 2, MPI_DOUBLE, &c->sbufferTypes[0]);
+    MPI_Type_commit(&c->sbufferTypes[0]);
+#else
+    c->imaxLocal = imax;
+    c->jmaxLocal = jmax;
+#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
+}

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

@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
+ * 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.
@@ -10,6 +10,7 @@
 
 #include "comm.h"
 
+#if defined(_MPI)
 // subroutines local to this module
 static int sizeOfRank(int rank, int size, int N)
 {
@@ -146,19 +147,23 @@ static int sum(int* sizes, int position)
 
     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, int direction)
 {
+#if defined(_MPI)
     switch (direction) {
     case LEFT:
         return c->coords[IDIM] == 0;
@@ -173,12 +178,14 @@ int commIsBoundary(Comm* c, int direction)
         return c->coords[JDIM] == (c->dims[JDIM] - 1);
         break;
     }
+#endif
 
-    return 0;
+    return 1;
 }
 
 void commExchange(Comm* c, double* grid)
 {
+#if defined(_MPI)
     int counts[NDIRS]      = { 1, 1, 1, 1 };
     MPI_Aint displs[NDIRS] = { 0, 0, 0, 0 };
 
@@ -191,10 +198,12 @@ void commExchange(Comm* c, double* grid)
         displs,
         c->rbufferTypes,
         c->comm);
+#endif
 }
 
 void commShift(Comm* c, double* f, double* g)
 {
+#if defined(_MPI)
     MPI_Request requests[4] = { MPI_REQUEST_NULL,
         MPI_REQUEST_NULL,
         MPI_REQUEST_NULL,
@@ -227,8 +236,10 @@ void commShift(Comm* c, double* f, double* g)
         &requests[3]);
 
     MPI_Waitall(4, requests, MPI_STATUSES_IGNORE);
+#endif
 }
 
+// TODO Merge with seq
 void commCollectResult(Comm* c,
     double* ug,
     double* vg,
@@ -276,6 +287,7 @@ void commCollectResult(Comm* c,
 
 void commPrintConfig(Comm* c)
 {
+#if defined(_MPI)
     fflush(stdout);
     MPI_Barrier(MPI_COMM_WORLD);
     if (commIsMaster(c)) {
@@ -296,13 +308,24 @@ void commPrintConfig(Comm* c)
         }
     }
     MPI_Barrier(MPI_COMM_WORLD);
+#endif
 }
 
-void commInit(Comm* c, int jmax, int imax)
+void commInit(Comm* c, int argc, char** argv)
 {
-    /* setup communication */
+#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 jmax, int imax)
+{
+#if defined(_MPI)
     int dims[NDIMS]    = { 0, 0 };
     int periods[NDIMS] = { 0, 0 };
     MPI_Dims_create(c->size, NDIMS, dims);
@@ -323,4 +346,20 @@ void commInit(Comm* c, int jmax, int imax)
     setupCommunication(c, BOTTOM, HALO);
     setupCommunication(c, TOP, BULK);
     setupCommunication(c, TOP, HALO);
+#else
+    c->imaxLocal = imax;
+    c->jmaxLocal = jmax;
+#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
 }

BasicSolver/2D-mpi/src/comm.h

@@ -1,12 +1,14 @@
 /*
- * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
+ * 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
 
 enum direction { LEFT = 0, RIGHT, BOTTOM, TOP, NDIRS };
 enum dimension { JDIM = 0, IDIM, NDIMS };
@@ -16,15 +18,19 @@ 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;
 } Comm;
 
-extern void commInit(Comm* c, int jmax, int imax);
+extern void commInit(Comm* c, int argc, char** argv);
+extern void commFinalize(Comm* c);
+extern void commPartition(Comm* c, int jmax, int imax);
 extern void commPrintConfig(Comm*);
 extern void commExchange(Comm*, double*);
 extern void commShift(Comm* c, double* f, double* g);

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

@@ -1,95 +1,84 @@
 /*
- * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
+ * 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 <float.h>
-#include <limits.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <unistd.h>
 
 #include "allocate.h"
+#include "comm.h"
 #include "parameter.h"
 #include "progress.h"
 #include "solver.h"
 #include "timing.h"
-#include <mpi.h>
 
 int main(int argc, char** argv)
 {
     int rank;
-    double S, E;
-    Parameter params;
-    Solver solver;
+    double timeStart, timeStop;
+    Parameter p;
+    Solver s;
 
-    MPI_Init(&argc, &argv);
-    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
-    initParameter(&params);
+    commInit(&s.comm, argc, argv);
+    initParameter(&p);
 
     if (argc != 2) {
         printf("Usage: %s <configFile>\n", argv[0]);
         exit(EXIT_SUCCESS);
     }
 
-    readParameter(&params, argv[1]);
-    if (rank == 0) {
-        printParameter(&params);
+    readParameter(&p, argv[1]);
+    commPartition(&s.comm, p.jmax, p.imax);
+    if (commIsMaster(&s.comm)) {
+        printParameter(&p);
     }
-    initSolver(&solver, &params);
-    /* debugExchange(&solver); */
-    /* exit(EXIT_SUCCESS); */
-    initProgress(solver.te);
+    initSolver(&s, &p);
+#ifndef VERBOSE
+    initProgress(s.te);
+#endif
 
-    double tau = solver.tau;
-    double te  = solver.te;
+    double tau = s.tau;
+    double te  = s.te;
     double t   = 0.0;
 
-    S = getTimeStamp();
+    timeStart = getTimeStamp();
     while (t <= te) {
-        if (tau > 0.0) {
-            computeTimestep(&solver);
-        }
-
-        setBoundaryConditions(&solver);
-        setSpecialBoundaryCondition(&solver);
-        computeFG(&solver);
-        computeRHS(&solver);
-        solve(&solver);
-        adaptUV(&solver);
-        t += solver.dt;
+        if (tau > 0.0) computeTimestep(&s);
+        setBoundaryConditions(&s);
+        setSpecialBoundaryCondition(&s);
+        computeFG(&s);
+        computeRHS(&s);
+        solve(&s);
+        adaptUV(&s);
+        t += s.dt;
 
 #ifdef VERBOSE
-        if (rank == 0) {
-            printf("TIME %f , TIMESTEP %f\n", t, solver.dt);
+        if (commIsMaster(&s.comm)) {
+            printf("TIME %f , TIMESTEP %f\n", t, s.dt);
         }
 #else
         printProgress(t);
 #endif
     }
-    E = getTimeStamp();
+    timeStop = getTimeStamp();
+#ifndef VERBOSE
     stopProgress();
-    if (rank == 0) {
-        printf("Solution took %.2fs\n", E - S);
+#endif
+    if (commIsMaster(&s.comm)) {
+        printf("Solution took %.2fs\n", timeStop - timeStart);
     }
-    size_t bytesize = solver.imax * solver.jmax * sizeof(double);
+    size_t bytesize = s.imax * s.jmax * sizeof(double);
 
     double* ug = allocate(64, bytesize);
     double* vg = allocate(64, bytesize);
     double* pg = allocate(64, bytesize);
 
-    commCollectResult(&solver.comm,
-        ug,
-        vg,
-        pg,
-        solver.u,
-        solver.v,
-        solver.p,
-        solver.jmax,
-        solver.imax);
-    writeResult(&solver, ug, vg, pg);
+    commCollectResult(&s.comm, ug, vg, pg, s.u, s.v, s.p, s.jmax, s.imax);
+    writeResult(&s, ug, vg, pg);
 
-    MPI_Finalize();
+    commFinalize(&s.comm);
     return EXIT_SUCCESS;
 }

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

@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
+ * 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.

BasicSolver/2D-mpi/src/parameter.h

@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
+ * 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.

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

@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
+ * 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.

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

@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
+ * 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.

BasicSolver/2D-mpi/src/solver.c

@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
+ * 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.
@@ -75,8 +75,6 @@ void initSolver(Solver* s, Parameter* params)
     s->tau      = params->tau;
     s->gamma    = params->gamma;
 
-    commInit(&s->comm, s->jmax, s->imax);
-
     /* allocate arrays */
     int imaxLocal = s->comm.imaxLocal;
     int jmaxLocal = s->comm.jmaxLocal;

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

@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
+ * 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.

BasicSolver/2D-mpi/src/timing.c

@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
+ * 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.
@@ -7,18 +7,16 @@
 #include <stdlib.h>
 #include <time.h>
 
-double getTimeStamp()
+double getTimeStamp(void)
 {
     struct timespec ts;
     clock_gettime(CLOCK_MONOTONIC, &ts);
     return (double)ts.tv_sec + (double)ts.tv_nsec * 1.e-9;
 }
 
-double getTimeResolution()
+double getTimeResolution(void)
 {
     struct timespec ts;
     clock_getres(CLOCK_MONOTONIC, &ts);
     return (double)ts.tv_sec + (double)ts.tv_nsec * 1.e-9;
 }
-
-double getTimeStamp_() { return getTimeStamp(); }

BasicSolver/2D-mpi/src/timing.h

@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
+ * 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.
@@ -7,8 +7,7 @@
 #ifndef __TIMING_H_
 #define __TIMING_H_
 
-extern double getTimeStamp();
-extern double getTimeResolution();
-extern double getTimeStamp_();
+extern double getTimeStamp(void);
+extern double getTimeResolution(void);
 
 #endif // __TIMING_H_

BasicSolver/2D-mpi/src/util.h

@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2022 NHR@FAU, University Erlangen-Nuremberg.
+ * 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.