xref: /honee/src/honee-file.c (revision 360b37c96ac5468b4453c10e4cd8c16bd6da1f1e)
193ca29b6SJames Wright // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors.
293ca29b6SJames Wright // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause
393ca29b6SJames Wright 
493ca29b6SJames Wright /// @file
593ca29b6SJames Wright /// Custom file I/O functions for HONEE
693ca29b6SJames Wright 
793ca29b6SJames Wright #include <honee-file.h>
893ca29b6SJames Wright 
9c9ff4f08SJames Wright const PetscInt32 HONEE_FILE_TOKEN    = 0xceedf00;  // for backwards compatibility
10c9ff4f08SJames Wright const PetscInt32 HONEE_FILE_TOKEN_32 = 0xceedf32;
11c9ff4f08SJames Wright const PetscInt32 HONEE_FILE_TOKEN_64 = 0xceedf64;
1293ca29b6SJames Wright 
13c9ff4f08SJames Wright // @brief Read in binary int based on it's data type
1493ca29b6SJames Wright static PetscErrorCode BinaryReadIntoInt(PetscViewer viewer, PetscInt *out, PetscDataType file_type) {
1593ca29b6SJames Wright   PetscFunctionBeginUser;
1693ca29b6SJames Wright   *out = -13;  // appease the overzealous GCC compiler warning Gods
1793ca29b6SJames Wright   if (file_type == PETSC_INT32) {
1893ca29b6SJames Wright     PetscInt32 val;
1993ca29b6SJames Wright     PetscCall(PetscViewerBinaryRead(viewer, &val, 1, NULL, PETSC_INT32));
2093ca29b6SJames Wright     *out = val;
2193ca29b6SJames Wright   } else if (file_type == PETSC_INT64) {
2293ca29b6SJames Wright     PetscInt64 val;
2393ca29b6SJames Wright     PetscCall(PetscViewerBinaryRead(viewer, &val, 1, NULL, PETSC_INT64));
2493ca29b6SJames Wright     *out = val;
2593ca29b6SJames Wright   } else {
2693ca29b6SJames Wright     PetscCall(PetscViewerBinaryRead(viewer, out, 1, NULL, PETSC_INT));
2793ca29b6SJames Wright   }
2893ca29b6SJames Wright   PetscFunctionReturn(PETSC_SUCCESS);
2993ca29b6SJames Wright }
3093ca29b6SJames Wright 
31c9ff4f08SJames Wright /**
32c9ff4f08SJames Wright   @brief Load vector from binary file, possibly with embedded solution time and step number
33c9ff4f08SJames Wright 
34c9ff4f08SJames Wright   Reads in Vec from binary file, possibly written by HONEE.
35c9ff4f08SJames Wright   If written by HONEE, will also load the solution time and timestep, otherwise not.
36c9ff4f08SJames Wright   Also handles case where file was written with different PetscInt size than is read with.
37c9ff4f08SJames Wright 
38*360b37c9SJames Wright   @param[in]  viewer      `PetscViewer` to read the vec from. Must be a binary viewer
39c9ff4f08SJames Wright   @param[out] Q           `Vec` to read the data into
40c9ff4f08SJames Wright   @param[out] time        Solution time the Vec was written at, or `NULL`. Set to 0 if legacy file format.
41c9ff4f08SJames Wright   @param[out] step_number Timestep number the Vec was written at, or `NULL`. Set to 0 if legacy file format.
42c9ff4f08SJames Wright **/
43*360b37c9SJames Wright PetscErrorCode HoneeLoadBinaryVec(PetscViewer viewer, Vec Q, PetscReal *time, PetscInt *step_number) {
4493ca29b6SJames Wright   PetscInt      file_step_number;
4593ca29b6SJames Wright   PetscInt32    token;
4693ca29b6SJames Wright   PetscReal     file_time;
4793ca29b6SJames Wright   PetscDataType file_type = PETSC_INT32;
48*360b37c9SJames Wright   MPI_Comm      comm      = PetscObjectComm((PetscObject)viewer);
4993ca29b6SJames Wright 
5093ca29b6SJames Wright   PetscFunctionBeginUser;
5193ca29b6SJames Wright   PetscCall(PetscViewerBinaryRead(viewer, &token, 1, NULL, PETSC_INT32));
52c9ff4f08SJames Wright   if (token == HONEE_FILE_TOKEN_32 || token == HONEE_FILE_TOKEN_64 ||
53c9ff4f08SJames Wright       token == HONEE_FILE_TOKEN) {  // New style format; we're reading a file with step number and time in the header
54c9ff4f08SJames Wright     if (token == HONEE_FILE_TOKEN_32) file_type = PETSC_INT32;
55c9ff4f08SJames Wright     else if (token == HONEE_FILE_TOKEN_64) file_type = PETSC_INT64;
5693ca29b6SJames Wright     PetscCall(BinaryReadIntoInt(viewer, &file_step_number, file_type));
5793ca29b6SJames Wright     PetscCall(PetscViewerBinaryRead(viewer, &file_time, 1, NULL, PETSC_REAL));
5893ca29b6SJames Wright   } else if (token == VEC_FILE_CLASSID) {  // Legacy format of just the vector, encoded as [VEC_FILE_CLASSID, length, ]
5993ca29b6SJames Wright     PetscInt length, N;
6093ca29b6SJames Wright     PetscCall(BinaryReadIntoInt(viewer, &length, file_type));
6193ca29b6SJames Wright     PetscCall(VecGetSize(Q, &N));
6293ca29b6SJames Wright     PetscCheck(length == N, comm, PETSC_ERR_ARG_INCOMP, "File Vec has length %" PetscInt_FMT " but DM has global Vec size %" PetscInt_FMT, length, N);
6393ca29b6SJames Wright     PetscCall(PetscViewerBinarySetSkipHeader(viewer, PETSC_TRUE));
64c9ff4f08SJames Wright     file_time        = 0;
65c9ff4f08SJames Wright     file_step_number = 0;
6693ca29b6SJames Wright   } else SETERRQ(comm, PETSC_ERR_FILE_UNEXPECTED, "Not a fluids header token or a PETSc Vec in file");
6793ca29b6SJames Wright 
68c9ff4f08SJames Wright   if (time) *time = file_time;
69c9ff4f08SJames Wright   if (step_number) *step_number = file_step_number;
7093ca29b6SJames Wright   PetscCall(VecLoad(Q, viewer));
7193ca29b6SJames Wright   PetscFunctionReturn(PETSC_SUCCESS);
7293ca29b6SJames Wright }
7393ca29b6SJames Wright 
74c9ff4f08SJames Wright /**
75c9ff4f08SJames Wright   @brief Open a PHASTA *.dat file, grabbing dimensions and file pointer
76c9ff4f08SJames Wright 
77c9ff4f08SJames Wright   This function opens the file specified by `path` using `PetscFOpen` and passes the file pointer in `fp`.
78c9ff4f08SJames Wright   It is not closed in this function, thus `fp` must be closed sometime after this function has been called (using `PetscFClose` for example).
79c9ff4f08SJames Wright 
80c9ff4f08SJames Wright   Assumes that the first line of the file has the number of rows and columns as the only two entries, separated by a single space.
81c9ff4f08SJames Wright 
82c9ff4f08SJames Wright   @param[in]  comm           MPI_Comm for the program
83c9ff4f08SJames Wright   @param[in]  path           Path to the file
84c9ff4f08SJames Wright   @param[in]  char_array_len Length of the character array that should contain each line
85c9ff4f08SJames Wright   @param[out] dims           Dimensions of the file, taken from the first line of the file
86c9ff4f08SJames Wright   @param[out] fp File        pointer to the opened file
87c9ff4f08SJames Wright **/
8893ca29b6SJames Wright PetscErrorCode PhastaDatFileOpen(const MPI_Comm comm, const char path[PETSC_MAX_PATH_LEN], const PetscInt char_array_len, PetscInt dims[2],
8993ca29b6SJames Wright                                  FILE **fp) {
9093ca29b6SJames Wright   int    ndims;
9193ca29b6SJames Wright   char   line[char_array_len];
9293ca29b6SJames Wright   char **array;
9393ca29b6SJames Wright 
9493ca29b6SJames Wright   PetscFunctionBeginUser;
9593ca29b6SJames Wright   PetscCall(PetscFOpen(comm, path, "r", fp));
9693ca29b6SJames Wright   PetscCall(PetscSynchronizedFGets(comm, *fp, char_array_len, line));
9793ca29b6SJames Wright   PetscCall(PetscStrToArray(line, ' ', &ndims, &array));
9893ca29b6SJames Wright   PetscCheck(ndims == 2, comm, PETSC_ERR_FILE_UNEXPECTED, "Found %d dimensions instead of 2 on the first line of %s", ndims, path);
9993ca29b6SJames Wright 
10093ca29b6SJames Wright   for (PetscInt i = 0; i < ndims; i++) dims[i] = atoi(array[i]);
10193ca29b6SJames Wright   PetscCall(PetscStrToArrayDestroy(ndims, array));
10293ca29b6SJames Wright   PetscFunctionReturn(PETSC_SUCCESS);
10393ca29b6SJames Wright }
10493ca29b6SJames Wright 
105c9ff4f08SJames Wright /**
106c9ff4f08SJames Wright   @brief Get the number of rows for the PHASTA file at path.
107c9ff4f08SJames Wright 
108c9ff4f08SJames Wright   Assumes that the first line of the file has the number of rows and columns as the only two entries, separated by a single space.
109c9ff4f08SJames Wright 
110c9ff4f08SJames Wright   @param[in]  comm  `MPI_Comm` for the program
111c9ff4f08SJames Wright   @param[in]  path  Path to the file
112c9ff4f08SJames Wright   @param[out] nrows Number of rows
113c9ff4f08SJames Wright **/
11493ca29b6SJames Wright PetscErrorCode PhastaDatFileGetNRows(const MPI_Comm comm, const char path[PETSC_MAX_PATH_LEN], PetscInt *nrows) {
11593ca29b6SJames Wright   const PetscInt char_array_len = 512;
11693ca29b6SJames Wright   PetscInt       dims[2];
11793ca29b6SJames Wright   FILE          *fp;
11893ca29b6SJames Wright 
11993ca29b6SJames Wright   PetscFunctionBeginUser;
12093ca29b6SJames Wright   PetscCall(PhastaDatFileOpen(comm, path, char_array_len, dims, &fp));
12193ca29b6SJames Wright   *nrows = dims[0];
12293ca29b6SJames Wright   PetscCall(PetscFClose(comm, fp));
12393ca29b6SJames Wright   PetscFunctionReturn(PETSC_SUCCESS);
12493ca29b6SJames Wright }
12593ca29b6SJames Wright 
126c9ff4f08SJames Wright /**
127c9ff4f08SJames Wright   @brief Read PetscReal values from PHASTA file
128c9ff4f08SJames Wright 
129c9ff4f08SJames Wright   @param[in]  comm  `MPI_Comm` for the program
130c9ff4f08SJames Wright   @param[in]  path  Path to the file
131c9ff4f08SJames Wright   @param[out] array Pointer to allocated array of correct size
132c9ff4f08SJames Wright **/
13393ca29b6SJames Wright PetscErrorCode PhastaDatFileReadToArrayReal(MPI_Comm comm, const char path[PETSC_MAX_PATH_LEN], PetscReal array[]) {
13493ca29b6SJames Wright   PetscInt       dims[2];
13593ca29b6SJames Wright   FILE          *fp;
13693ca29b6SJames Wright   const PetscInt char_array_len = 512;
13793ca29b6SJames Wright   char           line[char_array_len];
13893ca29b6SJames Wright 
13993ca29b6SJames Wright   PetscFunctionBeginUser;
14093ca29b6SJames Wright   PetscCall(PhastaDatFileOpen(comm, path, char_array_len, dims, &fp));
14193ca29b6SJames Wright 
14293ca29b6SJames Wright   for (PetscInt i = 0; i < dims[0]; i++) {
14393ca29b6SJames Wright     int    ndims;
14493ca29b6SJames Wright     char **row_array;
14593ca29b6SJames Wright 
14693ca29b6SJames Wright     PetscCall(PetscSynchronizedFGets(comm, fp, char_array_len, line));
14793ca29b6SJames Wright     PetscCall(PetscStrToArray(line, ' ', &ndims, &row_array));
14893ca29b6SJames Wright     PetscCheck(ndims == dims[1], comm, PETSC_ERR_FILE_UNEXPECTED,
14993ca29b6SJames Wright                "Line %" PetscInt_FMT " of %s does not contain enough columns (%d instead of %" PetscInt_FMT ")", i, path, ndims, dims[1]);
15093ca29b6SJames Wright 
15193ca29b6SJames Wright     for (PetscInt j = 0; j < dims[1]; j++) array[i * dims[1] + j] = (PetscReal)atof(row_array[j]);
15293ca29b6SJames Wright     PetscCall(PetscStrToArrayDestroy(ndims, row_array));
15393ca29b6SJames Wright   }
15493ca29b6SJames Wright 
15593ca29b6SJames Wright   PetscCall(PetscFClose(comm, fp));
15693ca29b6SJames Wright   PetscFunctionReturn(PETSC_SUCCESS);
15793ca29b6SJames Wright }
158