moebiusband/NuSiF-Solver
2
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Fix VTK MPI-IO Writer

Browse Source
This commit is contained in:
Jan Eitzinger 2023-02-07 09:48:24 +01:00
parent f7e04a19ac
commit 93a796af98
5 changed files with 77 additions and 39 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
BasicSolver/3D-mpi-io/src/main.c
View File

@ -76,6 +76,7 @@ int main(int argc, char** argv)
printf("Solution took %.2fs\n", timeStop - timeStart); printf("Solution took %.2fs\n", timeStop - timeStart);
} }
// testInit(&solver);
VtkOptions opts = { .grid = solver.grid, .comm = solver.comm }; VtkOptions opts = { .grid = solver.grid, .comm = solver.comm };
vtkOpen(&opts, solver.problem); vtkOpen(&opts, solver.problem);
vtkScalar(&opts, "pressure", solver.p); vtkScalar(&opts, "pressure", solver.p);

19
BasicSolver/3D-mpi-io/src/test.c
View File

@ -26,10 +26,21 @@ void testInit(Solver* s)
for (int k = 0; k < kmaxLocal + 2; k++) { for (int k = 0; k < kmaxLocal + 2; k++) {
for (int j = 0; j < jmaxLocal + 2; j++) { for (int j = 0; j < jmaxLocal + 2; j++) {
for (int i = 0; i < imaxLocal + 2; i++) { for (int i = 0; i < imaxLocal + 2; i++) {
G(p, i, j, k) = myrank; G(p, i, j, k) = 10.0;
G(f, i, j, k) = myrank; G(f, i, j, k) = myrank + 1.0;
G(g, i, j, k) = myrank; G(g, i, j, k) = myrank + 1.0;
G(h, i, j, k) = myrank; G(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++) {
G(p, i, j, k) = myrank + 1.0;
G(f, i, j, k) = myrank + 1.0;
G(g, i, j, k) = myrank + 1.0;
G(h, i, j, k) = myrank + 1.0;
} }
} }
} }

64
BasicSolver/3D-mpi-io/src/vtkWriter.c
View File

@ -14,6 +14,16 @@
#include "comm.h" #include "comm.h"
#include "vtkWriter.h" #include "vtkWriter.h"
static void resetFileview(VtkOptions* o)
{
// reset fileview for output of string header
MPI_Offset disp;
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);
// printf("Rank %d disp %lld\n", o->comm.rank, disp);
}
static void writeHeader(VtkOptions* o) static void writeHeader(VtkOptions* o)
{ {
const size_t MAX_HEADER = 200; const size_t MAX_HEADER = 200;
@ -42,10 +52,9 @@ static void writeHeader(VtkOptions* o)
o->grid.imax * o->grid.jmax * o->grid.kmax); o->grid.imax * o->grid.jmax * o->grid.kmax);
if (commIsMaster(&o->comm)) { if (commIsMaster(&o->comm)) {
MPI_File_write_at(o->fh, 0, header, strlen(header), MPI_CHAR, MPI_STATUS_IGNORE); MPI_File_write(o->fh, header, strlen(header), MPI_CHAR, MPI_STATUS_IGNORE);
} }
o->disp = strlen(header);
free(header); free(header);
} }
@ -57,22 +66,16 @@ void vtkOpen(VtkOptions* o, char* problem)
snprintf(filename, 50, "%s-p%d.vtk", problem, o->comm.size); snprintf(filename, 50, "%s-p%d.vtk", problem, o->comm.size);
MPI_File_open(o->comm.comm, MPI_File_open(o->comm.comm,
filename, filename,
MPI_MODE_WRONLY | MPI_MODE_CREATE, MPI_MODE_WRONLY | MPI_MODE_CREATE | MPI_MODE_EXCL,
MPI_INFO_NULL, MPI_INFO_NULL,
&o->fh); &o->fh);
if (commIsMaster(&o->comm)) { if (commIsMaster(&o->comm)) {
printf("Writing VTK output for %s\n", problem); printf("Writing VTK output for %s\n", problem);
} }
writeHeader(o);
}
static void resetFileview(VtkOptions* o) resetFileview(o);
{ writeHeader(o);
// 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) void vtkScalar(VtkOptions* o, char* name, double* s)
@ -85,7 +88,7 @@ void vtkScalar(VtkOptions* o, char* name, double* s)
char* header = (char*)malloc(MAX_HEADER); char* header = (char*)malloc(MAX_HEADER);
char* cursor = header; char* cursor = header;
cursor += sprintf(cursor, "SCALARS %s float 1\n", name); cursor += sprintf(cursor, "SCALARS %s double 1\n", name);
cursor += sprintf(cursor, "LOOKUP_TABLE default\n"); cursor += sprintf(cursor, "LOOKUP_TABLE default\n");
if (commIsMaster(&o->comm)) { if (commIsMaster(&o->comm)) {
@ -100,6 +103,7 @@ void vtkScalar(VtkOptions* o, char* name, double* s)
// set global view in file // set global view in file
MPI_Offset disp; MPI_Offset disp;
MPI_Datatype fileViewType; MPI_Datatype fileViewType;
MPI_Barrier(o->comm.comm);
MPI_File_get_size(o->fh, &disp); MPI_File_get_size(o->fh, &disp);
MPI_Type_create_subarray(NDIMS, MPI_Type_create_subarray(NDIMS,
@ -112,6 +116,18 @@ void vtkScalar(VtkOptions* o, char* name, double* s)
MPI_Type_commit(&fileViewType); MPI_Type_commit(&fileViewType);
MPI_File_set_view(o->fh, disp, MPI_DOUBLE, fileViewType, "external32", MPI_INFO_NULL); 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 // create local bulk type
MPI_Datatype bulkType; MPI_Datatype bulkType;
@ -152,7 +168,7 @@ void vtkVector(VtkOptions* o, char* name, VtkVector vec)
const size_t MAX_HEADER = 100; const size_t MAX_HEADER = 100;
char* header = (char*)malloc(MAX_HEADER); char* header = (char*)malloc(MAX_HEADER);
sprintf(header, "VECTORS %s float\n", name); sprintf(header, "VECTORS %s double\n", name);
resetFileview(o); resetFileview(o);
if (commIsMaster(&o->comm)) { if (commIsMaster(&o->comm)) {
@ -164,21 +180,25 @@ void vtkVector(VtkOptions* o, char* name, VtkVector vec)
// set global view in file // set global view in file
MPI_Offset disp; MPI_Offset disp;
MPI_Datatype fileViewType; MPI_Datatype fileViewType, vectorType;
MPI_Barrier(o->comm.comm);
MPI_File_get_size(o->fh, &disp); MPI_File_get_size(o->fh, &disp);
MPI_Type_create_subarray(NDIMS + 1, MPI_Type_contiguous(NDIMS, MPI_DOUBLE, &vectorType);
(int[NDIMS + 1]) { o->grid.kmax, o->grid.jmax, o->grid.imax, NDIMS }, MPI_Type_commit(&vectorType);
(int[NDIMS + 1]) { kmaxLocal, jmaxLocal, imaxLocal, NDIMS },
MPI_Type_create_subarray(NDIMS,
(int[NDIMS]) { o->grid.kmax, o->grid.jmax, o->grid.imax },
(int[NDIMS]) { kmaxLocal, jmaxLocal, imaxLocal },
offsets, offsets,
MPI_ORDER_C, MPI_ORDER_C,
MPI_DOUBLE, vectorType,
&fileViewType); &fileViewType);
MPI_Type_commit(&fileViewType); MPI_Type_commit(&fileViewType);
MPI_File_set_view(o->fh, disp, MPI_DOUBLE, fileViewType, "external32", MPI_INFO_NULL); MPI_File_set_view(o->fh, disp, MPI_DOUBLE, fileViewType, "external32", MPI_INFO_NULL);
size_t cnt = imaxLocal * jmaxLocal * kmaxLocal * NDIMS; size_t cnt = imaxLocal * jmaxLocal * kmaxLocal;
double* tmp = allocate(64, cnt * sizeof(double)); double* tmp = allocate(64, cnt * NDIMS * sizeof(double));
int idx = 0; int idx = 0;
for (int k = 1; k < kmaxLocal + 1; k++) { for (int k = 1; k < kmaxLocal + 1; k++) {
@ -191,8 +211,10 @@ void vtkVector(VtkOptions* o, char* name, VtkVector vec)
} }
} }
MPI_File_write(o->fh, tmp, cnt, MPI_DOUBLE, MPI_STATUS_IGNORE); 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(&fileViewType);
MPI_Type_free(&vectorType);
// Binary segment must be terminated with newline character // Binary segment must be terminated with newline character
resetFileview(o); resetFileview(o);

2
BasicSolver/3D-mpi/config.mk
View File

@ -1,5 +1,5 @@
# Supported: GCC, CLANG, ICC # Supported: GCC, CLANG, ICC
TAG ?= ICC TAG ?= CLANG
ENABLE_OPENMP ?= false ENABLE_OPENMP ?= false
#Feature options #Feature options

30
BasicSolver/3D-mpi/src/vtkWriter.c
View File

@ -11,18 +11,22 @@
#include "vtkWriter.h" #include "vtkWriter.h"
#define G(v, i, j, k) v[(k)*imax * jmax + (j)*imax + (i)] #define G(v, i, j, k) v[(k)*imax * jmax + (j)*imax + (i)]
static float floatSwap(float f) static double floatSwap(double f)
{ {
union { union {
float f; double f;
char b[4]; char b[8];
} dat1, dat2; } dat1, dat2;
dat1.f = f; dat1.f = f;
dat2.b[0] = dat1.b[3]; dat2.b[0] = dat1.b[7];
dat2.b[1] = dat1.b[2]; dat2.b[1] = dat1.b[6];
dat2.b[2] = dat1.b[1]; dat2.b[2] = dat1.b[5];
dat2.b[3] = dat1.b[0]; 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; return dat2.f;
} }
@ -75,8 +79,8 @@ static void writeScalar(VtkOptions* o, double* s)
if (o->fmt == ASCII) { if (o->fmt == ASCII) {
fprintf(o->fh, "%f\n", G(s, i, j, k)); fprintf(o->fh, "%f\n", G(s, i, j, k));
} else if (o->fmt == BINARY) { } else if (o->fmt == BINARY) {
fwrite((float[1]) { floatSwap(G(s, i, j, k)) }, fwrite((double[1]) { floatSwap(G(s, i, j, k)) },
sizeof(float), sizeof(double),
1, 1,
o->fh); o->fh);
} }
@ -102,7 +106,7 @@ void vtkScalar(VtkOptions* o, char* name, double* s)
if (o->mode == UNIX) { if (o->mode == UNIX) {
if (commIsMaster(&o->comm)) { if (commIsMaster(&o->comm)) {
if (!isInitialized(o->fh)) return; if (!isInitialized(o->fh)) return;
fprintf(o->fh, "SCALARS %s float 1\n", name); fprintf(o->fh, "SCALARS %s double 1\n", name);
fprintf(o->fh, "LOOKUP_TABLE default\n"); fprintf(o->fh, "LOOKUP_TABLE default\n");
writeScalar(o, s); writeScalar(o, s);
} }
@ -126,10 +130,10 @@ static void writeVector(VtkOptions* o, VtkVector vec)
G(vec.v, i, j, k), G(vec.v, i, j, k),
G(vec.w, i, j, k)); G(vec.w, i, j, k));
} else if (o->fmt == BINARY) { } else if (o->fmt == BINARY) {
fwrite((float[3]) { floatSwap(G(vec.u, i, j, k)), fwrite((double[3]) { floatSwap(G(vec.u, i, j, k)),
floatSwap(G(vec.v, i, j, k)), floatSwap(G(vec.v, i, j, k)),
floatSwap(G(vec.w, i, j, k)) }, floatSwap(G(vec.w, i, j, k)) },
sizeof(float), sizeof(double),
3, 3,
o->fh); o->fh);
} }
@ -146,7 +150,7 @@ void vtkVector(VtkOptions* o, char* name, VtkVector vec)
if (o->mode == UNIX) { if (o->mode == UNIX) {
if (commIsMaster(&o->comm)) { if (commIsMaster(&o->comm)) {
if (!isInitialized(o->fh)) return; if (!isInitialized(o->fh)) return;
fprintf(o->fh, "VECTORS %s float\n", name); fprintf(o->fh, "VECTORS %s double\n", name);
writeVector(o, vec); writeVector(o, vec);
} }
} else if (o->mode == MPI) { } else if (o->mode == MPI) {