Skip to content

Commit c0996da

Browse files
sanathkeshavsanathkeshav
authored
Use a persistent hdf5 file handler (#107)
1 parent 1893146 commit c0996da

8 files changed

Lines changed: 111 additions & 105 deletions

File tree

‎include/material_models/GBDiffusion.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -161,7 +161,7 @@ class GBDiffusion : public ThermalModel, public LinearModel<1, 3> {
161161
}
162162
}
163163

164-
void postprocess(Solver<1, 3> &solver, Reader &reader, const char *resultsFileName, int load_idx, int time_idx) override
164+
void postprocess(Solver<1, 3> &solver, Reader &reader, int load_idx, int time_idx) override
165165
{
166166
// Write GBnormals to HDF5 file if requested
167167
if (find(reader.resultsToWrite.begin(), reader.resultsToWrite.end(), "GBnormals") != reader.resultsToWrite.end()) {
@@ -174,7 +174,7 @@ class GBDiffusion : public ThermalModel, public LinearModel<1, 3> {
174174
GBnormals_field[element_idx * 3 + 2] = GBnormals[3 * mat_index + 2];
175175
}
176176
}
177-
reader.writeSlab("GBnormals", resultsFileName, solver.dataset_name, load_idx, time_idx, GBnormals_field, 3);
177+
reader.writeSlab("GBnormals", load_idx, time_idx, GBnormals_field, 3);
178178
FANS_free(GBnormals_field);
179179
}
180180
}

‎include/material_models/J2Plasticity.h‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -118,7 +118,7 @@ class J2Plasticity : public SmallStrainMechModel {
118118
virtual double compute_q_trial(double psi_val, int mat_index) = 0;
119119
virtual double compute_gamma(double f_trial, int mat_index, int i, ptrdiff_t element_idx) = 0;
120120

121-
void postprocess(Solver<3, 6> &solver, Reader &reader, const char *resultsFileName, int load_idx, int time_idx) override;
121+
void postprocess(Solver<3, 6> &solver, Reader &reader, int load_idx, int time_idx) override;
122122

123123
protected:
124124
// Material properties
@@ -237,7 +237,7 @@ class J2ViscoPlastic_NonLinearIsotropicHardening : public J2Plasticity {
237237
vector<double> sigma_diff; // sqrt(2/3) * (sigma_inf - yield_stress)
238238
};
239239

240-
inline void J2Plasticity::postprocess(Solver<3, 6> &solver, Reader &reader, const char *resultsFileName, int load_idx, int time_idx)
240+
inline void J2Plasticity::postprocess(Solver<3, 6> &solver, Reader &reader, int load_idx, int time_idx)
241241
{
242242
int n_str = 6; // The plastic strain and stress vectors have 6 components each
243243
VectorXd mean_plastic_strain = VectorXd::Zero(solver.local_n0 * solver.n_y * solver.n_z * n_str);
@@ -251,9 +251,9 @@ inline void J2Plasticity::postprocess(Solver<3, 6> &solver, Reader &reader, cons
251251
mean_kinematic_hardening_variable.segment(n_str * elem_idx, n_str) = psi_bar_t[elem_idx].rowwise().mean();
252252
}
253253

254-
reader.writeSlab("plastic_strain", resultsFileName, solver.dataset_name, load_idx, time_idx, mean_plastic_strain.data(), n_str);
255-
reader.writeSlab("isotropic_hardening_variable", resultsFileName, solver.dataset_name, load_idx, time_idx, mean_isotropic_hardening_variable.data(), 1);
256-
reader.writeSlab("kinematic_hardening_variable", resultsFileName, solver.dataset_name, load_idx, time_idx, mean_kinematic_hardening_variable.data(), n_str);
254+
reader.writeSlab("plastic_strain", load_idx, time_idx, mean_plastic_strain.data(), n_str);
255+
reader.writeSlab("isotropic_hardening_variable", load_idx, time_idx, mean_isotropic_hardening_variable.data(), 1);
256+
reader.writeSlab("kinematic_hardening_variable", load_idx, time_idx, mean_kinematic_hardening_variable.data(), n_str);
257257
}
258258

259259
#endif // J2PLASTICITY_H

‎include/material_models/PseudoPlastic.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -47,14 +47,14 @@ class PseudoPlastic : public SmallStrainMechModel {
4747

4848
virtual void get_sigma(int i, int mat_index, ptrdiff_t element_idx) override = 0; // Pure virtual method
4949

50-
void postprocess(Solver<3, 6> &solver, Reader &reader, const char *resultsFileName, int load_idx, int time_idx) override
50+
void postprocess(Solver<3, 6> &solver, Reader &reader, int load_idx, int time_idx) override
5151
{
5252
VectorXf element_plastic_flag = VectorXf::Zero(solver.local_n0 * solver.n_y * solver.n_z);
5353
for (ptrdiff_t elem_idx = 0; elem_idx < solver.local_n0 * solver.n_y * solver.n_z; ++elem_idx) {
5454
element_plastic_flag(elem_idx) = plastic_flag[elem_idx].cast<float>().mean();
5555
}
5656

57-
reader.writeSlab("plastic_flag", resultsFileName, solver.dataset_name, load_idx, time_idx, element_plastic_flag.data(), 1);
57+
reader.writeSlab("plastic_flag", load_idx, time_idx, element_plastic_flag.data(), 1);
5858
}
5959

6060
protected:

‎include/matmodel.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -29,7 +29,7 @@ class Matmodel {
2929
void getStrainStress(double *strain, double *stress, Matrix<double, howmany * 8, 1> &ue, int mat_index, ptrdiff_t element_idx);
3030
void setGradient(vector<double> _g0);
3131

32-
virtual void postprocess(Solver<howmany, n_str> &solver, Reader &reader, const char resultsFileName[], int load_idx, int time_idx) {};
32+
virtual void postprocess(Solver<howmany, n_str> &solver, Reader &reader, int load_idx, int time_idx) {};
3333

3434
virtual void initializeInternalVariables(ptrdiff_t num_elements, int num_gauss_points) {}
3535
virtual void updateInternalVariables() {}

‎include/reader.h‎

Lines changed: 53 additions & 59 deletions
Original file line numberDiff line numberDiff line change
@@ -16,15 +16,15 @@ using namespace std;
1616
class Reader {
1717
public:
1818
// Default constructor
19-
Reader();
19+
Reader() = default;
2020

2121
// Destructor to free allocated memory
2222
~Reader();
2323

2424
// contents of input file:
25-
char ms_filename[4096]; // Name of Micro-structure hdf5 file
26-
char ms_datasetname[4096]; // Absolute path of Micro-structure in hdf5 file
27-
char results_prefix[4096];
25+
char ms_filename[4096]{}; // Name of Micro-structure hdf5 file
26+
char ms_datasetname[4096]{}; // Absolute path of Micro-structure in hdf5 file
27+
char results_prefix[4096]{};
2828
int n_mat;
2929
json materialProperties;
3030
double TOL;
@@ -35,15 +35,18 @@ class Reader {
3535
string problemType;
3636
string matmodel;
3737
string method;
38-
string strain_type; // "small" (default) or "large"
39-
string FE_type; // "HEX8" (default), "HEX8R", or "BBAR"
38+
string strain_type{"small"}; // "small" (default) or "large"
39+
string FE_type; // "HEX8" (default), "HEX8R", or "BBAR"
4040
vector<string> resultsToWrite;
41+
char results_filename[4096]{}; // Output HDF5 filename
42+
char dataset_name[8192]{}; // Base path for results in HDF5 file
43+
hid_t results_file_id = -1; // Open HDF5 file handle for results
4144

4245
// contents of microstructure file:
4346
vector<int> dims;
4447
vector<double> l_e;
4548
vector<double> L;
46-
unsigned short *ms; // Micro-structure
49+
unsigned short *ms{nullptr}; // Micro-structure
4750
bool is_zyx = true;
4851

4952
int world_rank;
@@ -56,30 +59,30 @@ class Reader {
5659
ptrdiff_t local_1_start;
5760

5861
// void Setup(ptrdiff_t howmany);
59-
void ReadInputFile(char fn[]);
62+
void ReadInputFile(char input_fn[]);
6063
void ReadMS(int hm);
6164
void ComputeVolumeFractions();
6265
// void ReadHDF5(char file_name[], char dset_name[]);
6366
void safe_create_group(hid_t file, const char *const name);
67+
void OpenResultsFile(const char *output_fn); // Open results file once
68+
void CloseResultsFile(); // Explicitly close results file
6469

6570
// Convenience methods to check if a result should be written and write it
6671
template <typename T>
67-
void writeData(const char *fieldName, const char *file_name, const char *dataset_name,
68-
int load_idx, int time_idx, T *data, hsize_t *dims, int ndims);
72+
void writeData(const char *fieldName, int load_idx, int time_idx, T *data, hsize_t *dims, int ndims);
6973

7074
template <typename T>
71-
void writeSlab(const char *fieldName, const char *file_name, const char *dataset_name,
72-
int load_idx, int time_idx, T *data, int size);
75+
void writeSlab(const char *fieldName, int load_idx, int time_idx, T *data, int size);
7376

7477
template <typename T>
75-
void WriteSlab(T *data, int _howmany, const char *file_name, const char *dset_name);
78+
void WriteSlab(T *data, int _howmany, const char *dset_name);
7679

7780
template <typename T>
78-
void WriteData(T *data, const char *file_name, const char *dset_name, hsize_t *dims, int rank);
81+
void WriteData(T *data, const char *dset_name, hsize_t *dims, int rank);
7982
};
8083

8184
template <typename T>
82-
void Reader::WriteData(T *data, const char *file_name, const char *dset_name, hsize_t *dims, int rank)
85+
void Reader::WriteData(T *data, const char *dset_name, hsize_t *dims, int rank)
8386
{
8487
hid_t data_type;
8588
if (std::is_same<T, double>::value) {
@@ -90,26 +93,10 @@ void Reader::WriteData(T *data, const char *file_name, const char *dset_name, hs
9093
throw std::invalid_argument("Conversion of this data type to H5 data type not yet implemented");
9194
}
9295

93-
hid_t plist_id;
94-
hid_t file_id;
95-
plist_id = H5Pcreate(H5P_FILE_ACCESS);
96-
97-
// TODO: refactor this into a general error handling method
98-
/* Save old error handler */
99-
herr_t (*old_func)(hid_t, void *);
100-
void *old_client_data;
101-
H5Eget_auto(H5E_DEFAULT, &old_func, &old_client_data);
102-
/* Turn off error handling */
103-
H5Eset_auto(H5E_DEFAULT, NULL, NULL);
104-
105-
file_id = H5Fopen(file_name, H5F_ACC_RDWR, plist_id);
106-
/* Restore previous error handler */
107-
H5Eset_auto(H5E_DEFAULT, old_func, old_client_data);
108-
96+
hid_t file_id = results_file_id;
10997
if (file_id < 0) {
110-
file_id = H5Fcreate(file_name, H5F_ACC_EXCL, H5P_DEFAULT, plist_id);
98+
throw std::runtime_error("WriteData: results file is not open");
11199
}
112-
H5Pclose(plist_id);
113100

114101
// Ensure all groups in the path are created
115102
safe_create_group(file_id, dset_name);
@@ -135,18 +122,36 @@ void Reader::WriteData(T *data, const char *file_name, const char *dset_name, hs
135122
throw std::runtime_error("Error creating dataset");
136123
}
137124

138-
// Write the data to the dataset
139-
if (H5Dwrite(dataset_id, data_type, H5S_ALL, H5S_ALL, H5P_DEFAULT, data) < 0) {
125+
// Perform the write: only rank 0 has meaningful data; other ranks use an empty memory selection
126+
hid_t memspace = H5Screate_simple(rank, dims, NULL);
127+
if (world_rank == 0) {
128+
/* root keeps full memspace */
129+
} else {
130+
/* non-root select none so they participate but write nothing */
131+
H5Sselect_none(memspace);
132+
/* also ensure file selection is none on non-root ranks */
133+
H5Sselect_none(dataspace_id);
134+
}
135+
136+
hid_t dxpl = H5Pcreate(H5P_DATASET_XFER);
137+
H5Pset_dxpl_mpio(dxpl, H5FD_MPIO_INDEPENDENT);
138+
139+
// Use dataspace_id as the FILE selection; on non-root ranks select none so
140+
// the number of elements selected in memspace and filespace match per rank.
141+
if (H5Dwrite(dataset_id, data_type, memspace, dataspace_id, dxpl, (world_rank == 0 ? data : nullptr)) < 0) {
142+
H5Pclose(dxpl);
143+
H5Sclose(memspace);
140144
H5Dclose(dataset_id);
141145
H5Sclose(dataspace_id);
142-
H5Fclose(file_id);
143146
throw std::runtime_error("Error writing data to dataset");
144147
}
145148

146-
// Close the dataset and the file
149+
H5Pclose(dxpl);
150+
H5Sclose(memspace);
151+
152+
// Close the dataset and dataspace
147153
H5Dclose(dataset_id);
148154
H5Sclose(dataspace_id);
149-
H5Fclose(file_id);
150155
}
151156

152157
// this function has to be here because of the template
@@ -165,7 +170,6 @@ template <typename T>
165170
void Reader::WriteSlab(
166171
T *data, // in: local slab, layout [X][Y][Z][k]
167172
int _howmany, // global size of the 4th axis (k)
168-
const char *file_name,
169173
const char *dset_name)
170174
{
171175
/*------------------------------------------------------------------*/
@@ -188,20 +192,13 @@ void Reader::WriteSlab(
188192
/*------------------------------------------------------------------*/
189193
/* 1. open or create the HDF5 file */
190194
/*------------------------------------------------------------------*/
191-
hid_t plist_id = H5Pcreate(H5P_FILE_ACCESS);
192-
H5Pset_fapl_mpio(plist_id, MPI_COMM_WORLD, MPI_INFO_NULL);
193-
194-
/* temporarily silence “file not found” during H5Fopen */
195+
hid_t file_id = results_file_id;
196+
hid_t plist_id;
197+
if (file_id < 0) {
198+
throw std::runtime_error("WriteSlab: results file is not open");
199+
}
195200
herr_t (*old_func)(hid_t, void *);
196201
void *old_client_data;
197-
H5Eget_auto(H5E_DEFAULT, &old_func, &old_client_data);
198-
H5Eset_auto(H5E_DEFAULT, nullptr, nullptr);
199-
200-
hid_t file_id = H5Fopen(file_name, H5F_ACC_RDWR, plist_id);
201-
H5Eset_auto(H5E_DEFAULT, old_func, old_client_data); /* restore */
202-
if (file_id < 0) /* create if absent */
203-
file_id = H5Fcreate(file_name, H5F_ACC_EXCL, H5P_DEFAULT, plist_id);
204-
H5Pclose(plist_id);
205202

206203
/*------------------------------------------------------------------*/
207204
/* 2. create the dataset (global dims = Z Y X k ) if necessary */
@@ -298,31 +295,28 @@ void Reader::WriteSlab(
298295
H5Sclose(filespace);
299296
H5Pclose(plist_id);
300297
H5Dclose(dset_id);
301-
H5Fclose(file_id);
302298
}
303299

304300
template <typename T>
305-
void Reader::writeData(const char *fieldName, const char *file_name, const char *dataset_name,
306-
int load_idx, int time_idx, T *data, hsize_t *dims, int ndims)
301+
void Reader::writeData(const char *fieldName, int load_idx, int time_idx, T *data, hsize_t *dims, int ndims)
307302
{
308-
if (world_rank != 0 || std::find(resultsToWrite.begin(), resultsToWrite.end(), fieldName) == resultsToWrite.end()) {
303+
if (std::find(resultsToWrite.begin(), resultsToWrite.end(), fieldName) == resultsToWrite.end()) {
309304
return;
310305
}
311306
char name[5096];
312307
snprintf(name, sizeof(name), "%s/load%i/time_step%i/%s", dataset_name, load_idx, time_idx, fieldName);
313-
WriteData(data, file_name, name, dims, ndims);
308+
WriteData(data, name, dims, ndims);
314309
}
315310

316311
template <typename T>
317-
void Reader::writeSlab(const char *fieldName, const char *file_name, const char *dataset_name,
318-
int load_idx, int time_idx, T *data, int size)
312+
void Reader::writeSlab(const char *fieldName, int load_idx, int time_idx, T *data, int size)
319313
{
320314
if (std::find(resultsToWrite.begin(), resultsToWrite.end(), fieldName) == resultsToWrite.end()) {
321315
return;
322316
}
323317
char name[5096];
324318
snprintf(name, sizeof(name), "%s/load%i/time_step%i/%s", dataset_name, load_idx, time_idx, fieldName);
325-
WriteSlab(data, size, file_name, name);
319+
WriteSlab(data, size, name);
326320
}
327321

328322
#endif

‎include/solver.h‎

Lines changed: 16 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -51,7 +51,7 @@ class Solver : private MixedBCController<howmany> {
5151
template <int padding>
5252
void compute_residual(RealArray &r_matrix, RealArray &u_matrix);
5353

54-
void postprocess(Reader &reader, const char resultsFileName[], int load_idx, int time_idx); //!< Computes Strain and stress
54+
void postprocess(Reader &reader, int load_idx, int time_idx); //!< Computes Strain and stress
5555

5656
void convolution();
5757
double compute_error(RealArray &r);
@@ -63,8 +63,6 @@ class Solver : private MixedBCController<howmany> {
6363
MatrixXd homogenized_tangent;
6464
MatrixXd get_homogenized_tangent(double pert_param);
6565

66-
char dataset_name[5096]; // Dataset name for postprocessing results
67-
6866
void enableMixedBC(const MixedBC &mbc, size_t step)
6967
{
7068
this->activate(*this, mbc, step);
@@ -134,9 +132,6 @@ Solver<howmany, n_str>::Solver(Reader &reader, Matmodel<howmany, n_str> *mat)
134132
matmodel->initializeInternalVariables(local_n0 * n_y * n_z, matmodel->n_gp);
135133

136134
computeFundamentalSolution();
137-
138-
// Set dataset name as member variable of the Solver class
139-
sprintf(dataset_name, "%s_results/%s", reader.ms_datasetname, reader.results_prefix);
140135
}
141136

142137
template <int howmany, int n_str>
@@ -438,7 +433,7 @@ double Solver<howmany, n_str>::compute_error(RealArray &r)
438433
}
439434

440435
template <int howmany, int n_str>
441-
void Solver<howmany, n_str>::postprocess(Reader &reader, const char resultsFileName[], int load_idx, int time_idx)
436+
void Solver<howmany, n_str>::postprocess(Reader &reader, int load_idx, int time_idx)
442437
{
443438
VectorXd strain = VectorXd::Zero(local_n0 * n_y * n_z * n_str);
444439
VectorXd stress = VectorXd::Zero(local_n0 * n_y * n_z * n_str);
@@ -568,29 +563,29 @@ void Solver<howmany, n_str>::postprocess(Reader &reader, const char resultsFileN
568563
}
569564

570565
hsize_t dims[1] = {static_cast<hsize_t>(n_str)};
571-
reader.writeData("stress_average", resultsFileName, dataset_name, load_idx, time_idx, stress_average.data(), dims, 1);
572-
reader.writeData("strain_average", resultsFileName, dataset_name, load_idx, time_idx, strain_average.data(), dims, 1);
566+
reader.writeData("stress_average", load_idx, time_idx, stress_average.data(), dims, 1);
567+
reader.writeData("strain_average", load_idx, time_idx, strain_average.data(), dims, 1);
573568
for (int mat_index = 0; mat_index < n_mat; ++mat_index) {
574569
char stress_name[512];
575570
char strain_name[512];
576571
sprintf(stress_name, "phase_stress_average_phase%d", mat_index);
577572
sprintf(strain_name, "phase_strain_average_phase%d", mat_index);
578-
reader.writeData(stress_name, resultsFileName, dataset_name, load_idx, time_idx, phase_stress_average[mat_index].data(), dims, 1);
579-
reader.writeData(strain_name, resultsFileName, dataset_name, load_idx, time_idx, phase_strain_average[mat_index].data(), dims, 1);
573+
reader.writeData(stress_name, load_idx, time_idx, phase_stress_average[mat_index].data(), dims, 1);
574+
reader.writeData(strain_name, load_idx, time_idx, phase_strain_average[mat_index].data(), dims, 1);
580575
}
581576
dims[0] = iter + 1;
582-
reader.writeData("absolute_error", resultsFileName, dataset_name, load_idx, time_idx, err_all.data(), dims, 1);
577+
reader.writeData("absolute_error", load_idx, time_idx, err_all.data(), dims, 1);
583578

584579
vector<int> rank_field(local_n0 * n_y * n_z, world_rank);
585-
reader.writeSlab("mpi_rank", resultsFileName, dataset_name, load_idx, time_idx, rank_field.data(), 1);
586-
reader.writeSlab("microstructure", resultsFileName, dataset_name, load_idx, time_idx, ms, 1);
587-
reader.writeSlab("displacement_fluctuation", resultsFileName, dataset_name, load_idx, time_idx, v_u, howmany);
588-
reader.writeSlab("displacement", resultsFileName, dataset_name, load_idx, time_idx, u_total.data(), howmany);
589-
reader.writeSlab("residual", resultsFileName, dataset_name, load_idx, time_idx, v_r, howmany);
590-
reader.writeSlab("strain", resultsFileName, dataset_name, load_idx, time_idx, strain.data(), n_str);
591-
reader.writeSlab("stress", resultsFileName, dataset_name, load_idx, time_idx, stress.data(), n_str);
580+
reader.writeSlab("mpi_rank", load_idx, time_idx, rank_field.data(), 1);
581+
reader.writeSlab("microstructure", load_idx, time_idx, ms, 1);
582+
reader.writeSlab("displacement_fluctuation", load_idx, time_idx, v_u, howmany);
583+
reader.writeSlab("displacement", load_idx, time_idx, u_total.data(), howmany);
584+
reader.writeSlab("residual", load_idx, time_idx, v_r, howmany);
585+
reader.writeSlab("strain", load_idx, time_idx, strain.data(), n_str);
586+
reader.writeSlab("stress", load_idx, time_idx, stress.data(), n_str);
592587

593-
matmodel->postprocess(*this, reader, resultsFileName, load_idx, time_idx);
588+
matmodel->postprocess(*this, reader, load_idx, time_idx);
594589

595590
// Compute homogenized tangent only if requested
596591
if (find(reader.resultsToWrite.begin(), reader.resultsToWrite.end(), "homogenized_tangent") != reader.resultsToWrite.end()) {
@@ -600,7 +595,7 @@ void Solver<howmany, n_str>::postprocess(Reader &reader, const char resultsFileN
600595
cout << "# Homogenized tangent: " << endl
601596
<< setprecision(12) << homogenized_tangent << endl
602597
<< endl;
603-
reader.writeData("homogenized_tangent", resultsFileName, dataset_name, load_idx, time_idx, homogenized_tangent.data(), dims, 2);
598+
reader.writeData("homogenized_tangent", load_idx, time_idx, homogenized_tangent.data(), dims, 2);
604599
}
605600
}
606601
}

0 commit comments

Comments
 (0)