Skip to content

Commit ab545ae

Browse files
lucas56098lucas56098
committed
IO: added support for starting from snapshot (./ProteusGPU ./ics/param.txt 1 will start from last snapshot in output dir)
1 parent 0f9499c commit ab545ae

5 files changed

Lines changed: 227 additions & 23 deletions

File tree

‎src/begrun/begrun.cpp‎

Lines changed: 29 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -5,13 +5,14 @@
55
#include "../profiler/profiler.h"
66
#include <chrono>
77
#include <cmath>
8+
#include <cstdlib>
89
#include <iostream>
910
#include <thread>
1011

1112
namespace begrun {
1213

1314
// initalize Proteus
14-
void begrun(int argc, char* argv[]) {
15+
StartState begrun(int argc, char* argv[]) {
1516

1617
PROFILE_START("BEGRUN");
1718

@@ -37,12 +38,33 @@ namespace begrun {
3738
// load param.txt
3839
input = loadInputFiles(argc, argv);
3940

40-
// read IC file
41-
if (!input.readICFile(input.getParameter("ic_file"), icData)) { exit(EXIT_FAILURE); }
41+
// parse restart flag: ./ProteusGPU [param.txt] [0|1]
42+
// 0 = fresh start (default), 1 = restart from latest snapshot
43+
int restart_flag = (argc > 2) ? std::atoi(argv[2]) : 0;
44+
45+
StartState state = {0.0, 0};
46+
47+
if (restart_flag == 1) {
48+
// find and load latest snapshot
49+
std::string out_dir = input.getParameter("output_directory");
50+
int latest_n = InputHandler::findLatestSnapshot(out_dir);
51+
if (latest_n < 0) {
52+
std::cerr << "RESTART: Error! No snapshots found in " << out_dir << std::endl;
53+
exit(EXIT_FAILURE);
54+
}
55+
56+
std::string snap_path = out_dir + "snapshot_" + std::to_string(latest_n) + ".hdf5";
57+
std::cout << "RESTART: Loading snapshot " << snap_path << std::endl;
58+
59+
if (!input.readSnapshotFile(snap_path, icData, state.t_sim)) { exit(EXIT_FAILURE); }
60+
state.snap_num = latest_n + 1;
61+
} else {
62+
// fresh start: load IC file
63+
if (!input.readICFile(input.getParameter("ic_file"), icData)) { exit(EXIT_FAILURE); }
64+
}
4265

43-
// adapt buff for periodic bc to IC resolution (naive, will have to be improved later)
44-
buff = (1. / pow(icData.seedpos_dims[0], 1. / ((double)DIMENSION))) *
45-
4; // for approx unform grid this would be 4 layers of ghost points...
66+
// adapt buff for periodic bc to resolution
67+
buff = (1. / pow(icData.seedpos_dims[0], 1. / ((double)DIMENSION))) * 4;
4668
#ifdef MOVING_MESH
4769
// read mesh regularization parameters
4870
CellShapingSpeed = input.getParameterDouble("CellShapingSpeed");
@@ -56,6 +78,7 @@ namespace begrun {
5678
if (!output.initialize()) { exit(EXIT_FAILURE); }
5779

5880
PROFILE_END("BEGRUN");
81+
return state;
5982
}
6083

6184
void free_initial_conditions() {

‎src/begrun/begrun.h‎

Lines changed: 10 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -14,9 +14,16 @@
1414

1515
namespace begrun {
1616

17-
// called in main
18-
void begrun(int argc, char* argv[]);
19-
void free_initial_conditions();
17+
// simulation start state returned by begrun
18+
struct StartState {
19+
double t_sim;
20+
int snap_num;
21+
};
22+
23+
// called in main: ./ProteusGPU [param.txt] [restart_flag]
24+
// restart_flag: 0 = fresh start (default), 1 = restart from latest snapshot
25+
StartState begrun(int argc, char* argv[]);
26+
void free_initial_conditions();
2027

2128
// helpers
2229
void print_banner();

‎src/io/input.cpp‎

Lines changed: 159 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,7 @@
11
#include "input.h"
22
#include "../global/allvars.h"
33
#include <algorithm>
4+
#include <dirent.h>
45
#include <fstream>
56
#include <iostream>
67
#include <sstream>
@@ -109,7 +110,33 @@ bool InputHandler::getParameterBool(const std::string& key) const {
109110
"' to bool (expected: true/false/1/0/yes/no/on/off)");
110111
}
111112

113+
// check if a parameter exists
114+
bool InputHandler::hasParameter(const std::string& key) const {
115+
return parameters.find(key) != parameters.end();
116+
}
117+
112118
#ifdef USE_HDF5
119+
120+
// find the latest snapshot_N.hdf5 in a directory, return N (or -1 if none found)
121+
int InputHandler::findLatestSnapshot(const std::string& dir) {
122+
DIR* d = opendir(dir.c_str());
123+
if (!d) return -1;
124+
int max_num = -1;
125+
struct dirent* entry;
126+
while ((entry = readdir(d)) != NULL) {
127+
std::string name(entry->d_name);
128+
if (name.size() > 14 && name.substr(0, 9) == "snapshot_" && name.substr(name.size() - 5) == ".hdf5") {
129+
std::string num_str = name.substr(9, name.size() - 14);
130+
try {
131+
int num = std::stoi(num_str);
132+
if (num > max_num) max_num = num;
133+
} catch (...) {}
134+
}
135+
}
136+
closedir(d);
137+
return max_num;
138+
}
139+
113140
// opens IC.hdf5 file and reads initial conditions into ICData struct
114141
// will get extended if additional fields get added to IC file ...
115142
bool InputHandler::readICFile(const std::string& filename, ICData& icData) {
@@ -297,4 +324,136 @@ bool InputHandler::readICFile(const std::string& filename, ICData& icData) {
297324
std::cout << "INPUT: IC file loaded successfully!" << std::endl;
298325
return true;
299326
}
327+
328+
// read a snapshot file into ICData for restart (seeds from cells/seeds, hydro from hydro/)
329+
bool InputHandler::readSnapshotFile(const std::string& filename, ICData& icData, double& t_sim) {
330+
331+
std::ifstream f(filename);
332+
if (!f.good()) {
333+
std::cerr << "INPUT: Error! Snapshot file [" << filename << "] does not exist!" << std::endl;
334+
return false;
335+
}
336+
337+
hid_t file_id = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
338+
if (file_id < 0) {
339+
std::cerr << "INPUT: Error! Could not open snapshot file: " << filename << std::endl;
340+
return false;
341+
}
342+
343+
// read header/time
344+
hid_t header_group = H5Gopen(file_id, "header", H5P_DEFAULT);
345+
if (header_group < 0) {
346+
std::cerr << "INPUT: Error! Could not open header group in snapshot" << std::endl;
347+
H5Fclose(file_id);
348+
return false;
349+
}
350+
351+
hid_t attr_dim = H5Aopen(header_group, "dimension", H5P_DEFAULT);
352+
if (attr_dim >= 0) {
353+
H5Aread(attr_dim, H5T_NATIVE_INT, &icData.header.dimension);
354+
H5Aclose(attr_dim);
355+
}
356+
357+
#ifdef dim_2D
358+
if (icData.header.dimension != 2)
359+
#else
360+
if (icData.header.dimension != 3)
361+
#endif
362+
{
363+
std::cerr << "INPUT: Error! Snapshot dimension mismatch! Snapshot: " << icData.header.dimension
364+
<< "D, compiled: " << DIMENSION << "D" << std::endl;
365+
H5Gclose(header_group);
366+
H5Fclose(file_id);
367+
return false;
368+
}
369+
370+
hid_t attr_time = H5Aopen(header_group, "time", H5P_DEFAULT);
371+
if (attr_time >= 0) {
372+
H5Aread(attr_time, H5T_NATIVE_DOUBLE, &t_sim);
373+
H5Aclose(attr_time);
374+
} else {
375+
std::cerr << "INPUT: Error! Could not read time attribute from snapshot" << std::endl;
376+
H5Gclose(header_group);
377+
H5Fclose(file_id);
378+
return false;
379+
}
380+
H5Gclose(header_group);
381+
382+
// read cells/seeds -> icData.seedpos
383+
hid_t cells_group = H5Gopen(file_id, "cells", H5P_DEFAULT);
384+
if (cells_group < 0) {
385+
std::cerr << "INPUT: Error! Could not open cells group in snapshot" << std::endl;
386+
H5Fclose(file_id);
387+
return false;
388+
}
389+
390+
hid_t dataset_id = H5Dopen(cells_group, "seeds", H5P_DEFAULT);
391+
if (dataset_id < 0) {
392+
std::cerr << "INPUT: Error! Could not open seeds dataset in snapshot" << std::endl;
393+
H5Gclose(cells_group);
394+
H5Fclose(file_id);
395+
return false;
396+
}
397+
398+
hid_t dataspace_id = H5Dget_space(dataset_id);
399+
icData.seedpos_dims.resize(2);
400+
H5Sget_simple_extent_dims(dataspace_id, icData.seedpos_dims.data(), NULL);
401+
hsize_t totalElements = icData.seedpos_dims[0] * icData.seedpos_dims[1];
402+
icData.seedpos.resize(totalElements);
403+
H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, icData.seedpos.data());
404+
H5Sclose(dataspace_id);
405+
H5Dclose(dataset_id);
406+
H5Gclose(cells_group);
407+
408+
// read hydro/rho
409+
hid_t hydro_group = H5Gopen(file_id, "hydro", H5P_DEFAULT);
410+
if (hydro_group < 0) {
411+
std::cerr << "INPUT: Error! Could not open hydro group in snapshot" << std::endl;
412+
H5Fclose(file_id);
413+
return false;
414+
}
415+
416+
hsize_t n = icData.seedpos_dims[0];
417+
418+
dataset_id = H5Dopen(hydro_group, "rho", H5P_DEFAULT);
419+
if (dataset_id < 0) {
420+
std::cerr << "INPUT: Error! Could not open rho dataset in snapshot" << std::endl;
421+
H5Gclose(hydro_group);
422+
H5Fclose(file_id);
423+
return false;
424+
}
425+
icData.rho.resize(n);
426+
H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, icData.rho.data());
427+
H5Dclose(dataset_id);
428+
429+
// read hydro/vel
430+
dataset_id = H5Dopen(hydro_group, "vel", H5P_DEFAULT);
431+
if (dataset_id < 0) {
432+
std::cerr << "INPUT: Error! Could not open vel dataset in snapshot" << std::endl;
433+
H5Gclose(hydro_group);
434+
H5Fclose(file_id);
435+
return false;
436+
}
437+
icData.vel.resize(n * DIMENSION);
438+
H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, icData.vel.data());
439+
H5Dclose(dataset_id);
440+
441+
// read hydro/Energy
442+
dataset_id = H5Dopen(hydro_group, "Energy", H5P_DEFAULT);
443+
if (dataset_id < 0) {
444+
std::cerr << "INPUT: Error! Could not open Energy dataset in snapshot" << std::endl;
445+
H5Gclose(hydro_group);
446+
H5Fclose(file_id);
447+
return false;
448+
}
449+
icData.Energy.resize(n);
450+
H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, icData.Energy.data());
451+
H5Dclose(dataset_id);
452+
453+
H5Gclose(hydro_group);
454+
H5Fclose(file_id);
455+
456+
std::cout << "INPUT: Snapshot loaded successfully! (" << n << " cells, t = " << t_sim << ")" << std::endl;
457+
return true;
458+
}
300459
#endif

‎src/io/input.h‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -46,10 +46,17 @@ class InputHandler {
4646
int getParameterInt(const std::string& key) const;
4747
double getParameterDouble(const std::string& key) const;
4848
bool getParameterBool(const std::string& key) const;
49+
bool hasParameter(const std::string& key) const;
4950

5051
#ifdef USE_HDF5
5152
// read initial conditions from a HDF5 file
5253
bool readICFile(const std::string& filename, ICData& icData);
54+
55+
// read a snapshot file into ICData (for restart) and return the simulation time
56+
bool readSnapshotFile(const std::string& filename, ICData& icData, double& t_sim);
57+
58+
// find the latest snapshot_N.hdf5 in a directory, return N (or -1 if none found)
59+
static int findLatestSnapshot(const std::string& dir);
5360
#endif
5461
};
5562

‎src/main.cpp‎

Lines changed: 22 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -37,10 +37,12 @@ int main(int argc, char* argv[]) {
3737

3838
const auto wall_start = std::chrono::steady_clock::now();
3939

40-
// say hi and fill/prepare structs
41-
begrun::begrun(argc, argv);
40+
// say hi, load params, load IC or restart snapshot
41+
begrun::StartState state = begrun::begrun(argc, argv);
42+
double t_sim = state.t_sim;
43+
int snap_num = state.snap_num;
4244

43-
// init hydro values
45+
// init hydro values from icData (either fresh IC or restarted snapshot)
4446
primvars* primvar = hydro::init(icData.seedpos_dims[0]);
4547

4648
// compute voronoi mesh
@@ -50,21 +52,27 @@ int main(int argc, char* argv[]) {
5052
begrun::free_initial_conditions();
5153

5254
// start timestep loop
53-
std::cout << "HYDRO: started" << std::endl;
55+
if (t_sim > 0.0) {
56+
std::cout << "HYDRO: restarted from t = " << t_sim << " (snap_num = " << snap_num << ")" << std::endl;
57+
} else {
58+
std::cout << "HYDRO: started" << std::endl;
59+
}
5460

55-
double t_sim = 0.0;
56-
double t_end = input.getParameterDouble("time_end");
57-
double CFL = input.getParameterDouble("CFL_frac");
58-
int step = 0;
61+
const double t_start = t_sim;
62+
double t_end = input.getParameterDouble("time_end");
63+
double CFL = input.getParameterDouble("CFL_frac");
64+
int step = 0;
5965

6066
double output_dt = input.getParameterDouble("output_dt");
6167
double t_nextoutput = t_sim + output_dt;
62-
int snap_num = 0, next_log = 1;
68+
int next_log = 1;
6369

64-
// write first snapshot at t=0
70+
// write first snapshot at t=0 (only for fresh starts)
6571
#ifdef USE_HDF5
66-
output.snapshot(snap_num, mesh, primvar, icData.seedpos_dims[0], t_sim);
67-
snap_num += 1;
72+
if (snap_num == 0) {
73+
output.snapshot(snap_num, mesh, primvar, icData.seedpos_dims[0], t_sim);
74+
snap_num += 1;
75+
}
6876
#endif
6977

7078
PROFILE_START("HYDRO_MAIN");
@@ -84,8 +92,8 @@ int main(int argc, char* argv[]) {
8492
if (step >= next_log || t_sim >= t_end) {
8593
const double elapsed_s =
8694
std::chrono::duration<double>(std::chrono::steady_clock::now() - wall_start).count();
87-
std::cout << "HYDRO: Step " << step << " t = " << t_sim << " dt = " << dt
88-
<< " ETA = " << format_hms((t_sim > 0.0) ? elapsed_s * (t_end - t_sim) / t_sim : 0.0)
95+
std::cout << "HYDRO: Step " << step << " t = " << t_sim << " dt = " << dt << " ETA = "
96+
<< format_hms((t_sim > t_start) ? elapsed_s * (t_end - t_sim) / (t_sim - t_start) : 0.0)
8997
<< std::endl;
9098
const int guess = (elapsed_s > 1e-12) ? static_cast<int>(std::round(3.0 * step / elapsed_s)) : 1;
9199
const int bucket =

0 commit comments

Comments
 (0)