First draft. Compiles. Untested.

BasicSolver/3D-mpi-io/src/comm.c

@@ -113,60 +113,6 @@ static void setupCommunication(Comm* c, Direction direction, int layer)
     }
 }
 
-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);
-}
-
 static int sum(int* sizes, int position)
 {
     int sum = 0;
@@ -284,148 +230,26 @@ void commShift(Comm* c, double* f, double* g, double* h)
     MPI_Waitall(6, requests, MPI_STATUSES_IGNORE);
 }
 
-#define G(v, i, j, k)                                                                    \
+void commGetOffsets(Comm* c, int offsets[], int kmax, int jmax, int imax)
-    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 sum = 0;
-    int jmaxLocal = c->jmaxLocal;
-    int kmaxLocal = c->kmaxLocal;
 
-    int offset[c->size * NDIMS];
+    for (int i = 0; i < c->coords[ICORD]; i++) {
-    int imaxLocalAll[c->size];
+        sum += sizeOfRank(i, c->dims[ICORD], imax);
-    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, coords[ICORD]);
-            offset[i * NDIMS + JDIM] = sum(jmaxLocalAll, coords[JCORD]);
-            offset[i * NDIMS + KDIM] = sum(kmaxLocalAll, 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]);
-        }
     }
+    offsets[IDIM] = sum;
+    sum           = 0;
 
-    size_t bytesize = imaxLocal * jmaxLocal * kmaxLocal * sizeof(double);
+    for (int i = 0; i < c->coords[JCORD]; i++) {
-    double* tmp     = allocate(64, bytesize);
+        sum += sizeOfRank(i, c->dims[JCORD], jmax);
-    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);
-            }
-        }
     }
+    offsets[JDIM] = sum;
+    sum           = 0;
 
-    assembleResult(c,
+    for (int i = 0; i < c->coords[KCORD]; i++) {
-        tmp,
+        sum += sizeOfRank(i, c->dims[KCORD], kmax);
-        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;
     }
-        }
+    offsets[KDIM] = sum;
-    }
-
-    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);
 }
 
 void commPrintConfig(Comm* c)

BasicSolver/3D-mpi-io/src/comm.h

@@ -40,18 +40,7 @@ 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 commCollectResult(Comm* c,
+extern void commGetOffsets(Comm* c, int offsets[], int kmax, int jmax, int imax);
-    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_

BasicSolver/3D-mpi-io/src/main.c

@@ -76,39 +76,10 @@ int main(int argc, char** argv)
         printf("Solution took %.2fs\n", timeStop - timeStart);
     }
 
-    // testInit(&solver);
-    // commExchange(&solver.comm, solver.p);
-    // testPrintHalo(&solver, solver.p);
-
-    double *pg, *ug, *vg, *wg;
-
-    if (commIsMaster(&solver.comm)) {
-        size_t bytesize = solver.grid.imax * solver.grid.jmax * solver.grid.kmax *
-                          sizeof(double);
-
-        pg = allocate(64, bytesize);
-        ug = allocate(64, bytesize);
-        vg = allocate(64, bytesize);
-        wg = allocate(64, bytesize);
-    }
-
-    commCollectResult(&solver.comm,
-        ug,
-        vg,
-        wg,
-        pg,
-        solver.u,
-        solver.v,
-        solver.w,
-        solver.p,
-        solver.grid.kmax,
-        solver.grid.jmax,
-        solver.grid.imax);
-
     VtkOptions opts = { .grid = solver.grid, .comm = solver.comm };
     vtkOpen(&opts, solver.problem);
-    vtkScalar(&opts, "pressure", pg);
+    vtkScalar(&opts, "pressure", solver.p);
-    vtkVector(&opts, "velocity", (VtkVector) { ug, vg, wg });
+    vtkVector(&opts, "velocity", (VtkVector) { solver.u, solver.v, solver.w });
     vtkClose(&opts);
 
     commFree(&solver.comm);

BasicSolver/3D-mpi-io/src/vtkWriter.c

@@ -4,29 +4,15 @@
  * 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 "mpi.h"
+#include "allocate.h"
+#include "comm.h"
 #include "vtkWriter.h"
-#define G(v, i, j, k) v[(k)*imax * jmax + (j)*imax + (i)]
-
-static float floatSwap(float f)
-{
-    union {
-        float f;
-        char b[4];
-    } dat1, dat2;
-
-    dat1.f    = f;
-    dat2.b[0] = dat1.b[3];
-    dat2.b[1] = dat1.b[2];
-    dat2.b[2] = dat1.b[1];
-    dat2.b[3] = dat1.b[0];
-    return dat2.f;
-}
 
 static void writeHeader(VtkOptions* o)
 {
@@ -54,8 +40,16 @@ static void writeHeader(VtkOptions* o)
     cursor += sprintf(cursor,
         "POINT_DATA %d\n",
         o->grid.imax * o->grid.jmax * o->grid.kmax);
+
+    if (commIsMaster(&o->comm)) {
+        MPI_File_write_at(o->fh, 0, header, strlen(header), MPI_CHAR, MPI_STATUS_IGNORE);
     }
 
+    o->disp = strlen(header);
+    free(header);
+}
+
+// TODO Check inputs for length
 void vtkOpen(VtkOptions* o, char* problem)
 {
     char filename[50];
@@ -67,98 +61,144 @@ void vtkOpen(VtkOptions* o, char* problem)
         MPI_INFO_NULL,
         &o->fh);
 
-    if (commIsMaster(&o->comm)) printf("Writing VTK output for %s\n", problem);
+    if (commIsMaster(&o->comm)) {
+        printf("Writing VTK output for %s\n", problem);
+    }
     writeHeader(o);
 }
 
-static void writeScalar(VtkOptions* o, double* s)
+static void resetFileview(VtkOptions* o)
 {
-    int imax = o->grid.imax;
+    // reset fileview for output of string header
-    int jmax = o->grid.jmax;
+    MPI_Offset disp;
-    int kmax = o->grid.kmax;
+    MPI_File_get_size(o->fh, &disp);
-
+    MPI_File_set_view(o->fh, disp, MPI_CHAR, MPI_CHAR, "native", MPI_INFO_NULL);
-    for (int k = 0; k < kmax; k++) {
-        for (int j = 0; j < jmax; j++) {
-            for (int i = 0; i < imax; i++) {
-                if (o->fmt == ASCII) {
-                    fprintf(o->fh, "%f\n", G(s, i, j, k));
-                } else if (o->fmt == BINARY) {
-                    fwrite((float[1]) { floatSwap(G(s, i, j, k)) },
-                        sizeof(float),
-                        1,
-                        o->fh);
-                }
-            }
-        }
-    }
-    if (o->fmt == BINARY) fprintf(o->fh, "\n");
-}
-
-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)
 {
+    resetFileview(o);
+
     if (commIsMaster(&o->comm)) printf("Register scalar %s\n", name);
+    const size_t MAX_HEADER = 100;
+
+    char* header = (char*)malloc(MAX_HEADER);
+    char* cursor = header;
+
+    cursor += sprintf(cursor, "SCALARS %s float 1\n", name);
+    cursor += sprintf(cursor, "LOOKUP_TABLE default\n");
 
-    if (o->mode == UNIX) {
     if (commIsMaster(&o->comm)) {
-            if (!isInitialized(o->fh)) return;
+        MPI_File_write(o->fh, header, strlen(header), MPI_CHAR, MPI_STATUS_IGNORE);
-            fprintf(o->fh, "SCALARS %s float 1\n", name);
-            fprintf(o->fh, "LOOKUP_TABLE default\n");
-            writeScalar(o, s);
     }
-    } else if (o->mode == MPI) {
+
+    free(header);
+
+    int offsets[NDIMS];
+    commGetOffsets(&o->comm, offsets, o->grid.imax, o->grid.jmax, o->grid.kmax);
+
+    // set global view in file
+    MPI_Offset disp;
+    MPI_Datatype fileViewType;
+    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);
+
+    // 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);
     }
 }
 
-static void writeVector(VtkOptions* o, VtkVector vec)
+#define G(v, i, j, k)                                                                    \
-{
+    v[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)]
-    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",
-                        G(vec.u, i, j, k),
-                        G(vec.v, i, j, k),
-                        G(vec.w, i, j, k));
-                } else if (o->fmt == BINARY) {
-                    fwrite((float[3]) { floatSwap(G(vec.u, i, j, k)),
-                               floatSwap(G(vec.v, i, j, k)),
-                               floatSwap(G(vec.w, i, j, k)) },
-                        sizeof(float),
-                        3,
-                        o->fh);
-                }
-            }
-        }
-    }
-    if (o->fmt == BINARY) fprintf(o->fh, "\n");
-}
 
 void 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;
 
-    if (o->mode == UNIX) {
+    char* header = (char*)malloc(MAX_HEADER);
+    sprintf(header, "VECTORS %s float\n", name);
+
+    resetFileview(o);
     if (commIsMaster(&o->comm)) {
-            if (!isInitialized(o->fh)) return;
+        MPI_File_write(o->fh, header, strlen(header), MPI_CHAR, MPI_STATUS_IGNORE);
-            fprintf(o->fh, "VECTORS %s float\n", name);
+    }
-            writeVector(o, vec);
+
+    int offsets[NDIMS];
+    commGetOffsets(&o->comm, offsets, o->grid.imax, o->grid.jmax, o->grid.kmax);
+
+    // set global view in file
+    MPI_Offset disp;
+    MPI_Datatype fileViewType;
+    MPI_File_get_size(o->fh, &disp);
+
+    MPI_Type_create_subarray(NDIMS + 1,
+        (int[NDIMS + 1]) { o->grid.kmax, o->grid.jmax, o->grid.imax, NDIMS },
+        (int[NDIMS + 1]) { kmaxLocal, jmaxLocal, imaxLocal, NDIMS },
+        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);
+
+    size_t cnt  = imaxLocal * jmaxLocal * kmaxLocal * NDIMS;
+    double* tmp = allocate(64, cnt * 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;
             }
-    } else if (o->mode == MPI) {
         }
     }
 
-void vtkClose(VtkOptions* o) { MPI_File_close(o->fh); }
+    MPI_File_write(o->fh, tmp, cnt, MPI_DOUBLE, MPI_STATUS_IGNORE);
+    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);
+    }
+}
+
+void vtkClose(VtkOptions* o) { MPI_File_close(&o->fh); }