From f7e04a19acaf75a0ceb764224c30df1f6eeff2a6 Mon Sep 17 00:00:00 2001 From: Jan Eitzinger Date: Mon, 6 Feb 2023 20:09:11 +0100 Subject: [PATCH] First draft. Compiles. Untested. --- BasicSolver/3D-mpi-io/src/comm.c | 202 ++--------------------- BasicSolver/3D-mpi-io/src/comm.h | 13 +- BasicSolver/3D-mpi-io/src/main.c | 33 +--- BasicSolver/3D-mpi-io/src/vtkWriter.c | 222 +++++++++++++++----------- 4 files changed, 147 insertions(+), 323 deletions(-) diff --git a/BasicSolver/3D-mpi-io/src/comm.c b/BasicSolver/3D-mpi-io/src/comm.c index 568b550..add6b20 100644 --- a/BasicSolver/3D-mpi-io/src/comm.c +++ b/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) \ - 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) +void commGetOffsets(Comm* c, int offsets[], int kmax, int jmax, int imax) { - int imaxLocal = c->imaxLocal; - int jmaxLocal = c->jmaxLocal; - int kmaxLocal = c->kmaxLocal; + int sum = 0; - 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, 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]); - } + for (int i = 0; i < c->coords[ICORD]; i++) { + sum += sizeOfRank(i, c->dims[ICORD], imax); } + offsets[IDIM] = sum; + sum = 0; - 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); - } - } + for (int i = 0; i < c->coords[JCORD]; i++) { + sum += sizeOfRank(i, c->dims[JCORD], jmax); } + offsets[JDIM] = sum; + sum = 0; - 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; - } - } + for (int i = 0; i < c->coords[KCORD]; i++) { + sum += sizeOfRank(i, c->dims[KCORD], kmax); } - - 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); + offsets[KDIM] = sum; } void commPrintConfig(Comm* c) diff --git a/BasicSolver/3D-mpi-io/src/comm.h b/BasicSolver/3D-mpi-io/src/comm.h index 243db14..fead94f 100644 --- a/BasicSolver/3D-mpi-io/src/comm.h +++ b/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, - double* ug, - double* vg, - double* wg, - double* pg, - double* u, - double* v, - double* w, - double* p, - int kmax, - int jmax, - int imax); +extern void commGetOffsets(Comm* c, int offsets[], int kmax, int jmax, int imax); static inline int commIsMaster(Comm* c) { return c->rank == 0; } #endif // __COMM_H_ diff --git a/BasicSolver/3D-mpi-io/src/main.c b/BasicSolver/3D-mpi-io/src/main.c index f53e2a1..231dec4 100644 --- a/BasicSolver/3D-mpi-io/src/main.c +++ b/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); - vtkVector(&opts, "velocity", (VtkVector) { ug, vg, wg }); + vtkScalar(&opts, "pressure", solver.p); + vtkVector(&opts, "velocity", (VtkVector) { solver.u, solver.v, solver.w }); vtkClose(&opts); commFree(&solver.comm); diff --git a/BasicSolver/3D-mpi-io/src/vtkWriter.c b/BasicSolver/3D-mpi-io/src/vtkWriter.c index d7aaff7..29b27f2 100644 --- a/BasicSolver/3D-mpi-io/src/vtkWriter.c +++ b/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 #include #include #include #include -#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; - 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", 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; + // reset fileview for output of string header + MPI_Offset disp; + MPI_File_get_size(o->fh, &disp); + MPI_File_set_view(o->fh, disp, MPI_CHAR, MPI_CHAR, "native", MPI_INFO_NULL); } 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; - if (o->mode == UNIX) { - if (commIsMaster(&o->comm)) { - if (!isInitialized(o->fh)) return; - fprintf(o->fh, "SCALARS %s float 1\n", name); - fprintf(o->fh, "LOOKUP_TABLE default\n"); - writeScalar(o, s); - } - } else if (o->mode == MPI) { + 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 (commIsMaster(&o->comm)) { + MPI_File_write(o->fh, header, strlen(header), MPI_CHAR, MPI_STATUS_IGNORE); + } + + 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) -{ - 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"); -} +#define G(v, i, j, k) \ + v[(k) * (imaxLocal + 2) * (jmaxLocal + 2) + (j) * (imaxLocal + 2) + (i)] void vtkVector(VtkOptions* o, char* name, VtkVector vec) { - if (commIsMaster(&o->comm)) printf("Register vector %s\n", name); + int imaxLocal = o->comm.imaxLocal; + int jmaxLocal = o->comm.jmaxLocal; + int kmaxLocal = o->comm.kmaxLocal; - if (o->mode == UNIX) { - if (commIsMaster(&o->comm)) { - if (!isInitialized(o->fh)) return; - fprintf(o->fh, "VECTORS %s float\n", name); - writeVector(o, vec); + 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 float\n", name); + + resetFileview(o); + if (commIsMaster(&o->comm)) { + MPI_File_write(o->fh, header, strlen(header), MPI_CHAR, MPI_STATUS_IGNORE); + } + + 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) { + } + + 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); } +void vtkClose(VtkOptions* o) { MPI_File_close(&o->fh); }