1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors.
2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause
3
4 /// @file
5 /// Custom file I/O functions for HONEE
6
7 #include <honee-file.h>
8
9 /**
10 @brief Check if a filename has a file extension
11
12 Developer's note: Could instead use PetscStrendswith
13
14 @param[in] comm `MPI_Comm` used for error handling
15 @param[in] filename Filename to check
16 @param[in] extension Extension to check for
17 @param[out] is_extension Whether the filename has the extension
18 **/
HoneeCheckFilenameExtension(MPI_Comm comm,const char filename[],const char extension[],PetscBool * is_extension)19 PetscErrorCode HoneeCheckFilenameExtension(MPI_Comm comm, const char filename[], const char extension[], PetscBool *is_extension) {
20 size_t len, ext_len;
21
22 PetscFunctionBeginUser;
23 PetscCall(PetscStrlen(filename, &len));
24 PetscCall(PetscStrlen(extension, &ext_len));
25 PetscCheck(ext_len, comm, PETSC_ERR_ARG_WRONG, "Zero-size extension: %s", extension);
26 if (len < ext_len) *is_extension = PETSC_FALSE;
27 else PetscCall(PetscStrncmp(filename + len - ext_len, extension, ext_len, is_extension));
28 PetscFunctionReturn(PETSC_SUCCESS);
29 }
30
31 const PetscInt32 HONEE_FILE_TOKEN = 0xceedf00; // for backwards compatibility
32 const PetscInt32 HONEE_FILE_TOKEN_32 = 0xceedf32;
33 const PetscInt32 HONEE_FILE_TOKEN_64 = 0xceedf64;
34
35 // @brief Read in binary int based on it's data type
BinaryReadIntoInt(PetscViewer viewer,PetscInt * out,PetscDataType file_type)36 static PetscErrorCode BinaryReadIntoInt(PetscViewer viewer, PetscInt *out, PetscDataType file_type) {
37 PetscFunctionBeginUser;
38 *out = -13; // appease the overzealous GCC compiler warning Gods
39 if (file_type == PETSC_INT32) {
40 PetscInt32 val;
41 PetscCall(PetscViewerBinaryRead(viewer, &val, 1, NULL, PETSC_INT32));
42 *out = val;
43 } else if (file_type == PETSC_INT64) {
44 PetscInt64 val;
45 PetscCall(PetscViewerBinaryRead(viewer, &val, 1, NULL, PETSC_INT64));
46 PetscCall(PetscIntCast(val, out));
47 } else {
48 PetscCall(PetscViewerBinaryRead(viewer, out, 1, NULL, PETSC_INT));
49 }
50 PetscFunctionReturn(PETSC_SUCCESS);
51 }
52
53 /**
54 @brief Load initial condition from file
55
56 @param[in] filename File to get data from (must be binary or CGNS file)
57 @param[out] solution_steps Number of timesteps that the initial condition was taken from
58 @param[out] solution_time Solution time that the initial condition was taken from
59 @param[out] Q `Vec` to hold the initial condition
60 **/
HoneeLoadInitialCondition(const char filename[],PetscInt * solution_steps,PetscReal * solution_time,Vec Q)61 PetscErrorCode HoneeLoadInitialCondition(const char filename[], PetscInt *solution_steps, PetscReal *solution_time, Vec Q) {
62 MPI_Comm comm;
63 PetscViewer viewer;
64 PetscBool isBin, isCGNS;
65
66 PetscFunctionBeginUser;
67 PetscCall(PetscObjectGetComm((PetscObject)Q, &comm));
68 PetscCall(HoneeCheckFilenameExtension(comm, filename, ".bin", &isBin));
69 PetscCall(HoneeCheckFilenameExtension(comm, filename, ".cgns", &isCGNS));
70
71 if (isBin) {
72 PetscCall(PetscViewerBinaryOpen(comm, filename, FILE_MODE_READ, &viewer));
73 PetscCall(HoneeLoadBinaryVec(viewer, Q, solution_time, solution_steps));
74 PetscCall(PetscViewerDestroy(&viewer));
75 } else if (isCGNS) {
76 #if defined(PETSC_HAVE_CGNS)
77 DM dm, dm_output;
78 Vec V_output, V_local;
79 PetscBool set;
80
81 PetscCall(VecGetDM(Q, &dm));
82 PetscCall(PetscViewerCGNSOpen(comm, filename, FILE_MODE_READ, &viewer));
83 PetscCall(DMGetOutputDM(dm, &dm_output));
84 PetscCall(DMCreateGlobalVector(dm_output, &V_output));
85 PetscCall(DMGetLocalVector(dm, &V_local));
86 PetscCall(VecLoad(V_output, viewer));
87 PetscCall(DMGlobalToLocal(dm_output, V_output, INSERT_VALUES, V_local));
88 PetscCall(DMLocalToGlobal(dm, V_local, INSERT_VALUES, Q));
89 PetscCall(VecDestroy(&V_output));
90 PetscCall(DMRestoreLocalVector(dm, &V_local));
91
92 PetscCall(PetscViewerCGNSGetSolutionTime(viewer, solution_time, &set));
93 if (!set) PetscCall(PetscPrintf(comm, "WARNING: Couldn't find solution time in file\n"));
94 PetscCall(PetscViewerCGNSGetSolutionIteration(viewer, solution_steps, &set));
95 if (!set) { // Based on assumption that solution name of has form "FlowSolutionXXX"
96 const char *name, flowsolution[] = "FlowSolution";
97 size_t flowsolutionlen;
98 PetscBool isFlowSolution;
99
100 PetscCall(PetscViewerCGNSGetSolutionName(viewer, &name));
101 PetscCall(PetscStrlen(flowsolution, &flowsolutionlen));
102 PetscCall(PetscStrncmp(flowsolution, name, flowsolutionlen, &isFlowSolution));
103 PetscCheck(isFlowSolution, comm, PETSC_ERR_FILE_UNEXPECTED,
104 "CGNS file %s does not have FlowIterations array and solution name have 'FlowSolution' (found '%s')", filename, name);
105 *solution_steps = atoi(&name[flowsolutionlen]);
106 }
107
108 PetscCall(PetscViewerDestroy(&viewer));
109 #else
110 SETERRQ(comm, PETSC_ERR_SUP, "Loading initial condition from CGNS requires PETSc to be built with support. Reconfigure using --with-cgns-dir");
111 #endif
112 } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "File does not have a valid extension, recieved '%s'", filename);
113 PetscFunctionReturn(PETSC_SUCCESS);
114 }
115
116 /**
117 @brief Load vector from binary file, possibly with embedded solution time and step number
118
119 Reads in Vec from binary file, possibly written by HONEE.
120 If written by HONEE, will also load the solution time and timestep, otherwise not.
121 Also handles case where file was written with different PetscInt size than is read with.
122
123 @param[in] viewer `PetscViewer` to read the vec from. Must be a binary viewer
124 @param[out] Q `Vec` to read the data into
125 @param[out] time Solution time the Vec was written at, or `NULL`. Set to 0 if legacy file format.
126 @param[out] step_number Timestep number the Vec was written at, or `NULL`. Set to 0 if legacy file format.
127 **/
HoneeLoadBinaryVec(PetscViewer viewer,Vec Q,PetscReal * time,PetscInt * step_number)128 PetscErrorCode HoneeLoadBinaryVec(PetscViewer viewer, Vec Q, PetscReal *time, PetscInt *step_number) {
129 PetscInt file_step_number;
130 PetscReal file_time;
131 PetscInt32 token;
132 MPI_Comm comm = PetscObjectComm((PetscObject)viewer);
133 const char *filename;
134
135 PetscFunctionBeginUser;
136 PetscCall(PetscViewerFileGetName(viewer, &filename));
137 PetscCall(PetscViewerBinaryRead(viewer, &token, 1, NULL, PETSC_INT32));
138 if (token == HONEE_FILE_TOKEN_32 || token == HONEE_FILE_TOKEN_64 || token == HONEE_FILE_TOKEN) {
139 // HONEE formatted file; we're reading a file with step number and time in the header
140 PetscDataType file_type = PETSC_INT32;
141 if (token == HONEE_FILE_TOKEN_32) file_type = PETSC_INT32;
142 else if (token == HONEE_FILE_TOKEN_64) file_type = PETSC_INT64;
143 PetscCall(BinaryReadIntoInt(viewer, &file_step_number, file_type));
144 PetscCall(PetscViewerBinaryRead(viewer, &file_time, 1, NULL, PETSC_REAL));
145 } else {
146 // Legacy format of just the vector, encoded as [VEC_FILE_CLASSID, length, ]
147 // NOTE: Only used for regression testing now
148 PetscInt length, N;
149
150 PetscCall(VecGetSize(Q, &N));
151 if (token == VEC_FILE_CLASSID) {
152 PetscCall(BinaryReadIntoInt(viewer, &length, PETSC_INT32));
153 } else { // See `VecLoad_Binary()` in PETSc for source code reference of how to handle Vecs written with 64-bit PETSc
154 PetscInt32 token2;
155
156 PetscCall(PetscViewerBinaryRead(viewer, &token2, 1, NULL, PETSC_INT32));
157 token = (PetscInt)((uint64_t)token << 32) + token2; // Reconstruct token
158 PetscCheck(token == VEC_FILE_CLASSID, comm, PETSC_ERR_FILE_UNEXPECTED, "Not a HONEE header token or a PETSc Vec in file '%s'", filename);
159 PetscCall(BinaryReadIntoInt(viewer, &length, PETSC_INT64));
160 }
161
162 PetscCheck(length == N, comm, PETSC_ERR_ARG_INCOMP, "File Vec '%s' has length %" PetscInt_FMT " but DM has global Vec size %" PetscInt_FMT,
163 filename, length, N);
164 PetscCall(PetscViewerBinarySetSkipHeader(viewer, PETSC_TRUE));
165 file_time = 0;
166 file_step_number = 0;
167 }
168
169 if (time) *time = file_time;
170 if (step_number) *step_number = file_step_number;
171 PetscCall(VecLoad(Q, viewer));
172 PetscFunctionReturn(PETSC_SUCCESS);
173 }
174
175 /**
176 @brief Write vector to binary file with solution time and step number
177
178 @param[in] viewer `PetscViewer` for binary file. Must be binary viewer and in write mode
179 @param[in] Q `Vec` of the solution
180 @param[in] time Solution time of the `Vec`
181 @param[in] step_number Timestep number of the Vec
182 **/
HoneeWriteBinaryVec(PetscViewer viewer,Vec Q,PetscReal time,PetscInt step_number)183 PetscErrorCode HoneeWriteBinaryVec(PetscViewer viewer, Vec Q, PetscReal time, PetscInt step_number) {
184 PetscInt32 token = PetscDefined(USE_64BIT_INDICES) ? HONEE_FILE_TOKEN_64 : HONEE_FILE_TOKEN_32;
185
186 PetscFunctionBeginUser;
187 { // Verify viewer is in correct state
188 PetscViewerType viewer_type;
189 PetscFileMode file_mode;
190 PetscBool is_binary_viewer;
191 MPI_Comm comm = PetscObjectComm((PetscObject)viewer);
192
193 PetscCall(PetscViewerGetType(viewer, &viewer_type));
194 PetscCall(PetscStrcmp(viewer_type, PETSCVIEWERBINARY, &is_binary_viewer));
195 PetscCheck(is_binary_viewer, comm, PETSC_ERR_ARG_WRONGSTATE, "Viewer must be binary type; instead got %s", viewer_type);
196 PetscCall(PetscViewerFileGetMode(viewer, &file_mode));
197 PetscCheck(file_mode == FILE_MODE_WRITE, comm, PETSC_ERR_ARG_WRONGSTATE, "Viewer must be binary type; instead got %s",
198 file_mode == -1 ? "UNDEFINED" : PetscFileModes[file_mode]);
199 }
200
201 PetscCall(PetscViewerBinaryWrite(viewer, &token, 1, PETSC_INT32));
202 PetscCall(PetscViewerBinaryWrite(viewer, &step_number, 1, PETSC_INT));
203 PetscCall(PetscViewerBinaryWrite(viewer, &time, 1, PETSC_REAL));
204 PetscCall(VecView(Q, viewer));
205 PetscFunctionReturn(PETSC_SUCCESS);
206 }
207
208 /**
209 @brief Open a PHASTA *.dat file, grabbing dimensions and file pointer
210
211 This function opens the file specified by `path` using `PetscFOpen` and passes the file pointer in `fp`.
212 It is not closed in this function, thus `fp` must be closed sometime after this function has been called (using `PetscFClose` for example).
213
214 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.
215
216 @param[in] comm MPI_Comm for the program
217 @param[in] path Path to the file
218 @param[in] char_array_len Length of the character array that should contain each line
219 @param[out] dims Dimensions of the file, taken from the first line of the file
220 @param[out] fp File pointer to the opened file
221 **/
PhastaDatFileOpen(const MPI_Comm comm,const char path[],const PetscInt char_array_len,PetscInt dims[2],FILE ** fp)222 PetscErrorCode PhastaDatFileOpen(const MPI_Comm comm, const char path[], const PetscInt char_array_len, PetscInt dims[2], FILE **fp) {
223 int ndims;
224 char line[char_array_len];
225 char **array;
226
227 PetscFunctionBeginUser;
228 PetscCall(PetscFOpen(comm, path, "r", fp));
229 PetscCall(PetscSynchronizedFGets(comm, *fp, char_array_len, line));
230 PetscCall(PetscStrToArray(line, ' ', &ndims, &array));
231 PetscCheck(ndims == 2, comm, PETSC_ERR_FILE_UNEXPECTED, "Found %d dimensions instead of 2 on the first line of %s", ndims, path);
232
233 for (PetscInt i = 0; i < ndims; i++) dims[i] = atoi(array[i]);
234 PetscCall(PetscStrToArrayDestroy(ndims, array));
235 PetscFunctionReturn(PETSC_SUCCESS);
236 }
237
238 /**
239 @brief Get the number of rows for the PHASTA file at path.
240
241 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.
242
243 @param[in] comm `MPI_Comm` for the program
244 @param[in] path Path to the file
245 @param[out] nrows Number of rows
246 **/
PhastaDatFileGetNRows(const MPI_Comm comm,const char path[],PetscInt * nrows)247 PetscErrorCode PhastaDatFileGetNRows(const MPI_Comm comm, const char path[], PetscInt *nrows) {
248 const PetscInt char_array_len = 512;
249 PetscInt dims[2];
250 FILE *fp;
251
252 PetscFunctionBeginUser;
253 PetscCall(PhastaDatFileOpen(comm, path, char_array_len, dims, &fp));
254 *nrows = dims[0];
255 PetscCall(PetscFClose(comm, fp));
256 PetscFunctionReturn(PETSC_SUCCESS);
257 }
258
259 /**
260 @brief Read PetscReal values from PHASTA file
261
262 @param[in] comm `MPI_Comm` for the program
263 @param[in] path Path to the file
264 @param[out] array Pointer to allocated array of correct size
265 **/
PhastaDatFileReadToArrayReal(MPI_Comm comm,const char path[],PetscReal array[])266 PetscErrorCode PhastaDatFileReadToArrayReal(MPI_Comm comm, const char path[], PetscReal array[]) {
267 PetscInt dims[2];
268 FILE *fp;
269 const PetscInt char_array_len = 512;
270 char line[char_array_len];
271
272 PetscFunctionBeginUser;
273 PetscCall(PhastaDatFileOpen(comm, path, char_array_len, dims, &fp));
274
275 for (PetscInt i = 0; i < dims[0]; i++) {
276 int ndims;
277 char **row_array;
278
279 PetscCall(PetscSynchronizedFGets(comm, fp, char_array_len, line));
280 PetscCall(PetscStrToArray(line, ' ', &ndims, &row_array));
281 PetscCheck(ndims == dims[1], comm, PETSC_ERR_FILE_UNEXPECTED,
282 "Line %" PetscInt_FMT " of %s does not contain enough columns (%d instead of %" PetscInt_FMT ")", i, path, ndims, dims[1]);
283
284 for (PetscInt j = 0; j < dims[1]; j++) array[i * dims[1] + j] = (PetscReal)atof(row_array[j]);
285 PetscCall(PetscStrToArrayDestroy(ndims, row_array));
286 }
287
288 PetscCall(PetscFClose(comm, fp));
289 PetscFunctionReturn(PETSC_SUCCESS);
290 }
291