From 1fcddc91b7fa5ad36173b651d46948f163a1e760 Mon Sep 17 00:00:00 2001
From: Jan Eitzinger <jan@moebiusband.org>
Date: Tue, 9 Jan 2024 14:36:42 +0000
Subject: [PATCH] Add VTK output for sequential case

---
 BasicSolver/3D-mpi/src/main.c      | 66 +++++++++++++++++++++++++++---
 BasicSolver/3D-mpi/src/vtkWriter.c | 52 +++++++----------------
 BasicSolver/3D-mpi/src/vtkWriter.h |  4 --
 3 files changed, 76 insertions(+), 46 deletions(-)

diff --git a/BasicSolver/3D-mpi/src/main.c b/BasicSolver/3D-mpi/src/main.c
index ae7b165..73d5feb 100644
--- a/BasicSolver/3D-mpi/src/main.c
+++ b/BasicSolver/3D-mpi/src/main.c
@@ -18,6 +18,56 @@
 #include "timing.h"
 #include "vtkWriter.h"
 
+#ifndef _MPI
+#define G(v, i, j, k) v[(k) * (imax + 2) * (jmax + 2) + (j) * (imax + 2) + (i)]
+
+static void createBulkArrays(Solver* s, double* pg, double* ug, double* vg, double* wg)
+{
+    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;
+            }
+        }
+    }
+}
+#endif /* ifndef _MPI */
+
 int main(int argc, char** argv)
 {
     double timeStart, timeStop;
@@ -84,6 +134,7 @@ int main(int argc, char** argv)
         wg = allocate(64, bytesize);
     }
 
+#ifdef _MPI
     commCollectResult(&s.comm,
         ug,
         vg,
@@ -96,12 +147,17 @@ int main(int argc, char** argv)
         s.grid.kmax,
         s.grid.jmax,
         s.grid.imax);
+#else
+    createBulkArrays(&s, pg, ug, vg, wg);
+#endif /* ifdef _MPI */
 
-    VtkOptions opts = { .grid = s.grid, .comm = s.comm };
-    vtkOpen(&opts, s.problem);
-    vtkScalar(&opts, "pressure", pg);
-    vtkVector(&opts, "velocity", (VtkVector) { ug, vg, wg });
-    vtkClose(&opts);
+    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);
+    }
 
     commFinalize(&s.comm);
     return EXIT_SUCCESS;
diff --git a/BasicSolver/3D-mpi/src/vtkWriter.c b/BasicSolver/3D-mpi/src/vtkWriter.c
index e96e773..3433b3a 100644
--- a/BasicSolver/3D-mpi/src/vtkWriter.c
+++ b/BasicSolver/3D-mpi/src/vtkWriter.c
@@ -55,16 +55,11 @@ void vtkOpen(VtkOptions* o, char* problem)
 {
     char filename[50];
 
-    if (o->mode == UNIX) {
-        if (commIsMaster(&o->comm)) {
-            snprintf(filename, 50, "%s-p%d.vtk", problem, o->comm.size);
-            o->fh = fopen(filename, "w");
-            writeHeader(o);
-        }
-    } else if (o->mode == MPI) {
-    }
+    snprintf(filename, 50, "%s.vtk", problem);
+    o->fh = fopen(filename, "w");
+    writeHeader(o);
 
-    if (commIsMaster(&o->comm)) printf("Writing VTK output for %s\n", problem);
+    printf("Writing VTK output for %s\n", problem);
 }
 
 static void writeScalar(VtkOptions* o, double* s)
@@ -101,17 +96,11 @@ static bool isInitialized(FILE* ptr)
 
 void vtkScalar(VtkOptions* o, char* name, double* s)
 {
-    if (commIsMaster(&o->comm)) printf("Register scalar %s\n", name);
-
-    if (o->mode == UNIX) {
-        if (commIsMaster(&o->comm)) {
-            if (!isInitialized(o->fh)) return;
-            fprintf(o->fh, "SCALARS %s double 1\n", name);
-            fprintf(o->fh, "LOOKUP_TABLE default\n");
-            writeScalar(o, s);
-        }
-    } else if (o->mode == MPI) {
-    }
+    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)
@@ -145,25 +134,14 @@ static void writeVector(VtkOptions* o, VtkVector vec)
 
 void vtkVector(VtkOptions* o, char* name, VtkVector vec)
 {
-    if (commIsMaster(&o->comm)) printf("Register vector %s\n", name);
-
-    if (o->mode == UNIX) {
-        if (commIsMaster(&o->comm)) {
-            if (!isInitialized(o->fh)) return;
-            fprintf(o->fh, "VECTORS %s double\n", name);
-            writeVector(o, vec);
-        }
-    } else if (o->mode == MPI) {
-    }
+    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)
 {
-    if (o->mode == UNIX) {
-        if (commIsMaster(&o->comm)) {
-            fclose(o->fh);
-            o->fh = NULL;
-        }
-    } else if (o->mode == MPI) {
-    }
+    fclose(o->fh);
+    o->fh = NULL;
 }
diff --git a/BasicSolver/3D-mpi/src/vtkWriter.h b/BasicSolver/3D-mpi/src/vtkWriter.h
index 7a91dbb..d166f62 100644
--- a/BasicSolver/3D-mpi/src/vtkWriter.h
+++ b/BasicSolver/3D-mpi/src/vtkWriter.h
@@ -8,18 +8,14 @@
 #define __VTKWRITER_H_
 #include <stdio.h>
 
-#include "comm.h"
 #include "grid.h"
 
 typedef enum VtkFormat { ASCII = 0, BINARY } VtkFormat;
-typedef enum VtkMode { UNIX = 0, MPI } VtkMode;
 
 typedef struct VtkOptions {
     VtkFormat fmt;
-    VtkMode mode;
     Grid grid;
     FILE* fh;
-    Comm comm;
 } VtkOptions;
 
 typedef struct VtkVector {