Line | Branch | Exec | Source |
---|---|---|---|
1 | #include <H5Cpp.h> | ||
2 | #include <ctime> | ||
3 | #include <filesystem> | ||
4 | #include <fstream> | ||
5 | #include <iostream> | ||
6 | #include <limits> | ||
7 | #include <string> | ||
8 | |||
9 | #include "integrator.hpp" | ||
10 | #include "multiindex.hpp" | ||
11 | #include "observables.hpp" | ||
12 | #include "output.hpp" | ||
13 | #include "particle.hpp" | ||
14 | #include "simulation.hpp" | ||
15 | #include "system.hpp" | ||
16 | #include "timeseries.hpp" | ||
17 | #include "trajectories.hpp" | ||
18 | #include "version.hpp" | ||
19 | |||
20 | 1 | Teebuf::Teebuf(std::streambuf *sbuf1, std::streambuf *sbuf2) : sbuf1(sbuf1), sbuf2(sbuf2) {} | |
21 | |||
22 | 818 | int Teebuf::overflow(int c) { | |
23 |
1/2✓ Branch 0 taken 818 times.
✗ Branch 1 not taken.
|
818 | if (c != std::char_traits<char>::eof()) { |
24 | 818 | this->sbuf1->sputc(c); | |
25 | 818 | this->sbuf2->sputc(c); | |
26 | } | ||
27 | 818 | return std::char_traits<char>::not_eof(c); | |
28 | } | ||
29 | |||
30 | 3 | int Teebuf::sync() { | |
31 | 3 | const int r1 = sbuf1->pubsync(); | |
32 | 3 | const int r2 = sbuf2->pubsync(); | |
33 | 3 | return r1 == 0 && r2 == 0 ? 0 : 1; | |
34 | } | ||
35 | |||
36 | 1 | Output::Output(const std::filesystem::path &outputPath, bool dryrun) | |
37 |
1/2✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
|
1 | : outputPath(outputPath), dryrun(dryrun) { |
38 |
2/5✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 9 not taken.
|
1 | if (std::filesystem::exists(this->outputPath) && !std::filesystem::is_empty(this->outputPath)) { |
39 | ✗ | throw std::invalid_argument("Output path " + outputPath.string() + " already exists."); | |
40 | } | ||
41 | 1 | time_t walltimeBegin = time(0); | |
42 | 1 | char timestamp[20]; | |
43 | 1 | strftime(timestamp, sizeof(timestamp), "%Y-%m-%d_%H-%M-%S", localtime(&walltimeBegin)); | |
44 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
|
1 | if (this->outputPath.empty()) { |
45 | ✗ | this->outputPath = this->dryrun ? std::filesystem::temp_directory_path() : "results"; | |
46 | ✗ | this->outputPath /= std::string(timestamp); | |
47 | } | ||
48 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | std::filesystem::create_directories(this->outputPath); |
49 |
3/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
|
2 | this->logFile = std::ofstream(this->outputPath / "log.txt"); |
50 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
1 | this->tee = new Teebuf(this->logFile.rdbuf(), std::clog.rdbuf()); |
51 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | this->clogbuf = std::clog.rdbuf(tee); |
52 | 1 | std::cout << std::boolalpha; | |
53 | 1 | std::clog << std::boolalpha; | |
54 | // print version info after clog is tee'd to also have this info in log.txt | ||
55 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
2 | std::clog << programInfo(); |
56 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
1 | std::clog << "\nBegin of run: " << timestamp; |
57 |
3/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
|
1 | std::clog << "\nOutput path: " << this->outputPath.string(); |
58 | 1 | try { | |
59 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
1 | std::filesystem::remove("latest"); |
60 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
1 | std::filesystem::create_directory_symlink(this->outputPath, "latest"); |
61 | ✗ | } catch (const std::filesystem::filesystem_error &e) { | |
62 | ✗ | std::clog << "\nWARNING: Could not create 'latest' symlink"; | |
63 | } | ||
64 | 1 | } | |
65 | |||
66 | 2 | Output::~Output() { | |
67 |
1/2✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
|
2 | if (dryrun) { |
68 | 2 | std::clog | |
69 | 2 | << "\nWARNING: This was a dry-run. If no errors occured, you may run the actual simulation." | |
70 | 2 | << std::endl; | |
71 | 2 | std::filesystem::remove_all(outputPath); | |
72 | } | ||
73 | 2 | std::clog.flush(); | |
74 | 2 | std::clog.rdbuf(clogbuf); | |
75 |
1/2✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
|
2 | delete this->tee; |
76 | } | ||
77 | |||
78 | 1 | void Output::emitConfig(const YAML::Node &config) { | |
79 |
1/2✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
|
2 | std::ofstream out(outputPath / "config.yaml"); |
80 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | out << config; |
81 | 1 | } | |
82 | |||
83 | ✗ | void Output::writeParticles(const std::vector<Particle> &particles, const std::string &filename) { | |
84 | ✗ | std::ofstream stream(outputPath / filename); | |
85 | ✗ | stream << std::fixed << std::setprecision(10); | |
86 | ✗ | for (auto &p : particles) { | |
87 | ✗ | for (int i = 0; i < p.r.size(); ++i) { | |
88 | ✗ | stream << p.r[i] << " "; | |
89 | } | ||
90 | ✗ | for (int i = 0; i < p.v.size(); ++i) { | |
91 | ✗ | stream << p.v[i] << " "; | |
92 | } | ||
93 | ✗ | stream << p.sigma << " "; | |
94 | ✗ | stream << p.m << " "; | |
95 | ✗ | stream << p.gamma << " "; | |
96 | ✗ | stream << p.species << "\n"; | |
97 | } | ||
98 | } | ||
99 | |||
100 | ✗ | void Output::writeStats(Simulation *simulation, const std::string &filename) { | |
101 | ✗ | std::ofstream statsStream(outputPath / filename); | |
102 | ✗ | int col = 15; | |
103 | ✗ | statsStream << std::setw(5) << "stage" << std::setw(col) << "stepsMade" << std::setw(col) | |
104 | << "tBefore" << std::setw(col) << "tAfter" << std::setw(col) << "walltime" | ||
105 | ✗ | << "\n"; | |
106 | ✗ | for (size_t i = 0; i < simulation->stageStats.size(); ++i) { | |
107 | ✗ | auto &s = simulation->stageStats[i]; | |
108 | ✗ | int walltime = s.walltimeAfter ? difftime(s.walltimeAfter, s.walltimeBefore) : 0; | |
109 | ✗ | statsStream << std::setw(5) << i << std::setw(col) << s.stepsMade << std::setw(col) << s.tBefore | |
110 | ✗ | << std::setw(col) << s.tAfter << std::setw(col) << std::to_string(walltime) << "\n"; | |
111 | } | ||
112 | } | ||
113 | |||
114 | ✗ | void Output::write(Simulation *simulation) { | |
115 | ✗ | unsigned int stage = simulation->stageCounter; | |
116 | ✗ | for (auto &o : simulation->observables) { | |
117 | ✗ | writeObservables(o, stage); | |
118 | } | ||
119 | ✗ | if (simulation->timeseries) { | |
120 | ✗ | writeTimeseries(simulation->timeseries, stage); | |
121 | } | ||
122 | ✗ | writeParticles(simulation->system->particles, "particlesLast.dat"); | |
123 | ✗ | writeStats(simulation, "stats.dat"); | |
124 | } | ||
125 | |||
126 | ✗ | void Output::writeTrajectoryFrame(Trajectories *trajectories, System *system, int stage) { | |
127 | ✗ | throw std::invalid_argument("Not implemented. Use output format H5MD for particle trajectories."); | |
128 | } | ||
129 | |||
130 | 1 | H5MDOutput::H5MDOutput(const std::filesystem::path &outputPath, bool dryrun) | |
131 |
1/2✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
|
1 | : Output(outputPath, dryrun) {} |
132 | |||
133 | 4 | H5MDOutput::~H5MDOutput() { | |
134 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
|
2 | if (trajectoryWriter.joinable()) { |
135 | ✗ | trajectoryWriter.join(); | |
136 | } | ||
137 | 4 | } | |
138 | |||
139 | ✗ | static H5::Group createOrOpen(H5::Group &parent, const std::string &name) { | |
140 | ✗ | return parent.nameExists(name) ? parent.openGroup(name) : parent.createGroup(name); | |
141 | } | ||
142 | |||
143 | ✗ | static H5::DataSet createOrOpen(H5::Group &group, const std::string &name, const H5::DataType &type, | |
144 | const H5::DataSpace &dataspace) { | ||
145 | ✗ | return group.nameExists(name) ? group.openDataSet(name) | |
146 | ✗ | : group.createDataSet(name, type, dataspace); | |
147 | } | ||
148 | |||
149 | ✗ | void H5MDOutput::write(Simulation *simulation) { | |
150 | ✗ | initFile(outputPath, std::to_string(simulation->stageCounter) + ".h5md"); | |
151 | ✗ | Output::write(simulation); | |
152 | ✗ | file.flush(H5F_SCOPE_GLOBAL); | |
153 | } | ||
154 | |||
155 | ✗ | void H5MDOutput::writeMeta() { | |
156 | ✗ | H5::Group group = createOrOpen(file, "h5md"); | |
157 | // version | ||
158 | ✗ | H5::Attribute version = group.createAttribute("version", H5::PredType::NATIVE_INT, | |
159 | ✗ | H5::DataSpace(1, std::vector<hsize_t>{2}.data())); | |
160 | ✗ | version.write(version.getDataType(), std::vector<int>{1, 0}.data()); | |
161 | // author | ||
162 | ✗ | H5::Group author = createOrOpen(group, "author"); | |
163 | ✗ | std::string authorNameData = "Florian Sammueller"; | |
164 | ✗ | H5::Attribute authorName = | |
165 | ✗ | author.createAttribute("name", H5::StrType(H5::PredType::C_S1, H5T_VARIABLE), | |
166 | ✗ | H5::DataSpace(1, std::vector<hsize_t>{1}.data())); | |
167 | ✗ | authorName.write(authorName.getDataType(), authorNameData); | |
168 | ✗ | std::string emailData = "florian.sammueller@uni-bayreuth.de"; | |
169 | ✗ | H5::Attribute email = author.createAttribute( | |
170 | ✗ | "email", H5::StrType(H5::PredType::C_S1, H5T_VARIABLE), H5::DataSpace()); | |
171 | ✗ | email.write(email.getDataType(), emailData); | |
172 | // creator | ||
173 | ✗ | H5::Group creator = createOrOpen(group, "creator"); | |
174 | ✗ | std::string creatorNameData = "MBD"; | |
175 | ✗ | H5::Attribute creatorName = creator.createAttribute( | |
176 | ✗ | "name", H5::StrType(H5::PredType::C_S1, H5T_VARIABLE), H5::DataSpace()); | |
177 | ✗ | creatorName.write(creatorName.getDataType(), creatorNameData); | |
178 | ✗ | std::string creatorVersionData = programInfo(); | |
179 | ✗ | H5::Attribute creatorVersion = creator.createAttribute( | |
180 | ✗ | "version", H5::StrType(H5::PredType::C_S1, H5T_VARIABLE), H5::DataSpace()); | |
181 | ✗ | creatorVersion.write(creatorVersion.getDataType(), creatorVersionData); | |
182 | } | ||
183 | |||
184 | ✗ | void H5MDOutput::initFile(const std::filesystem::path &path, const std::string &filename) { | |
185 | ✗ | if (trajectoryWriter.joinable()) { | |
186 | ✗ | trajectoryWriter.join(); | |
187 | } | ||
188 | ✗ | if (!std::filesystem::exists(path / filename)) { | |
189 | ✗ | file.close(); // file could still point to an old one | |
190 | ✗ | file = H5::H5File(path / filename, H5F_ACC_TRUNC | H5F_ACC_SWMR_WRITE); | |
191 | ✗ | writeMeta(); // this is static info, so only write once | |
192 | } | ||
193 | } | ||
194 | |||
195 | ✗ | void H5MDOutput::writeObservables(Observables *observables, int stage) { | |
196 | ✗ | H5::Group observablesGroup = createOrOpen(file, "observables"); | |
197 | ✗ | H5::Group group = createOrOpen(observablesGroup, observables->type()); | |
198 | ✗ | std::vector<int> shapeTmp = observables->shape(); | |
199 | ✗ | std::vector<hsize_t> shape(shapeTmp.begin(), shapeTmp.end()); | |
200 | ✗ | std::vector<std::string> axisLabels = observables->axisLabels(); | |
201 | ✗ | H5::Group coordinatesGroup; | |
202 | ✗ | if (shape.size()) { | |
203 | ✗ | coordinatesGroup = createOrOpen(group, "coordinates"); | |
204 | } | ||
205 | ✗ | for (size_t axis = 0; axis < shape.size(); ++axis) { | |
206 | ✗ | std::vector<double> coordinates = observables->coordinates(axis); | |
207 | ✗ | H5::DataSet dataset = | |
208 | ✗ | createOrOpen(coordinatesGroup, std::to_string(axis), H5::PredType::NATIVE_DOUBLE, | |
209 | ✗ | H5::DataSpace(1, std::vector<hsize_t>{coordinates.size()}.data())); | |
210 | ✗ | dataset.write(coordinates.data(), dataset.getDataType()); | |
211 | ✗ | if (!dataset.attrExists("label")) { | |
212 | ✗ | H5::Attribute axisLabel = dataset.createAttribute( | |
213 | ✗ | "label", H5::StrType(H5::PredType::C_S1, H5T_VARIABLE), H5::DataSpace()); | |
214 | ✗ | axisLabel.write(axisLabel.getDataType(), axisLabels[axis]); | |
215 | } | ||
216 | } | ||
217 | ✗ | Multiindex<Eigen::Dynamic> speciesMultiindex(Eigen::Array<int, Eigen::Dynamic, 1>::Constant( | |
218 | ✗ | observables->numSpeciesIndices, observables->numSpecies)); | |
219 | ✗ | for (const auto &key : observables->keys) { | |
220 | ✗ | for (auto species : speciesMultiindex) { | |
221 | ✗ | H5::Group speciesGroup = group; | |
222 | ✗ | for (int s = 0; s < species.size(); ++s) { | |
223 | ✗ | speciesGroup = createOrOpen(speciesGroup, std::to_string(species[s])); | |
224 | } | ||
225 | ✗ | H5::DataSet dataset = createOrOpen(speciesGroup, key, H5::PredType::NATIVE_DOUBLE, | |
226 | ✗ | H5::DataSpace(shape.size(), shape.data())); | |
227 | ✗ | dataset.write(observables->data(key, species).data(), dataset.getDataType()); | |
228 | } | ||
229 | } | ||
230 | } | ||
231 | |||
232 | ✗ | void H5MDOutput::writeTimeseries(Timeseries *timeseries, int stage) { | |
233 | // According to H5MD spec, timeseries data shall also be in the group "observables" | ||
234 | ✗ | H5::Group observablesGroup = createOrOpen(file, "observables"); | |
235 | ✗ | std::vector<double> &time = timeseries->time; | |
236 | ✗ | if (!observablesGroup.nameExists("timeseries")) { // init extendable dataspace | |
237 | ✗ | H5::Group timeseriesGroup = observablesGroup.createGroup("timeseries"); | |
238 | ✗ | hsize_t dimsfChunk[1] = {1000}; // Might matter for performance | |
239 | ✗ | H5::DSetCreatPropList cparms; | |
240 | ✗ | cparms.setChunk(1, dimsfChunk); | |
241 | ✗ | double fillVal = std::numeric_limits<double>::quiet_NaN(); | |
242 | ✗ | cparms.setFillValue(H5::PredType::NATIVE_DOUBLE, &fillVal); | |
243 | ✗ | hsize_t dimsfInitial[1] = {0}; | |
244 | ✗ | hsize_t maxdimsf[1] = {H5S_UNLIMITED}; | |
245 | ✗ | timeseriesGroup.createDataSet("time", H5::PredType::NATIVE_DOUBLE, | |
246 | ✗ | H5::DataSpace(1, dimsfInitial, maxdimsf), cparms); | |
247 | ✗ | for (const auto &[key, _] : timeseries->data) { | |
248 | ✗ | H5::Group group = timeseriesGroup.createGroup(key); | |
249 | ✗ | group.link(H5L_TYPE_HARD, "/observables/timeseries/time", "time"); | |
250 | ✗ | group.createDataSet("value", H5::PredType::NATIVE_DOUBLE, | |
251 | ✗ | H5::DataSpace(1, dimsfInitial, maxdimsf), cparms); | |
252 | } | ||
253 | } | ||
254 | ✗ | H5::Group timeseriesGroup = observablesGroup.openGroup("timeseries"); | |
255 | ✗ | H5::DataSet datasetTime = timeseriesGroup.openDataSet("time"); | |
256 | ✗ | hsize_t dimsf[1] = {time.size()}; | |
257 | ✗ | hsize_t offset[1]; | |
258 | ✗ | datasetTime.getSpace().getSimpleExtentDims(offset); | |
259 | ✗ | hsize_t dimsfExtend[1] = {dimsf[0] - offset[0]}; | |
260 | ✗ | datasetTime.extend(dimsf); | |
261 | ✗ | H5::DataSpace mExtend(1, dimsfExtend); | |
262 | ✗ | H5::DataSpace fExtendTime = datasetTime.getSpace(); | |
263 | ✗ | fExtendTime.selectHyperslab(H5S_SELECT_SET, dimsfExtend, offset); | |
264 | ✗ | datasetTime.write(time.data() + offset[0], datasetTime.getDataType(), mExtend, fExtendTime); | |
265 | ✗ | for (const auto &[key, values] : timeseries->data) { | |
266 | ✗ | H5::Group group = timeseriesGroup.openGroup(key); | |
267 | ✗ | H5::DataSet datasetValue = group.openDataSet("value"); | |
268 | ✗ | datasetValue.extend(dimsf); | |
269 | ✗ | H5::DataSpace fExtendValue = datasetValue.getSpace(); | |
270 | ✗ | fExtendValue.selectHyperslab(H5S_SELECT_SET, dimsfExtend, offset); | |
271 | ✗ | datasetValue.write(values.data() + offset[0], datasetValue.getDataType(), mExtend, | |
272 | fExtendValue); | ||
273 | } | ||
274 | } | ||
275 | |||
276 | ✗ | void H5MDOutput::initParticlesGroup(Trajectories *trajectories, System *system) { | |
277 | ✗ | hsize_t dim = DIM; | |
278 | ✗ | hsize_t n = system->particles.size(); | |
279 | ✗ | H5::Group particlesGroup = file.createGroup("particles"); | |
280 | // box | ||
281 | ✗ | H5::Group boxGroup = particlesGroup.createGroup("box"); | |
282 | ✗ | H5::Attribute dimension = | |
283 | ✗ | boxGroup.createAttribute("dimension", H5::PredType::NATIVE_INT, H5::DataSpace()); | |
284 | ✗ | dimension.write(dimension.getDataType(), &dim); | |
285 | ✗ | std::vector<const char *> boundaryStrings = {"periodic", "periodic", "periodic"}; | |
286 | ✗ | H5::Attribute boundary = boxGroup.createAttribute( | |
287 | ✗ | "boundary", H5::StrType(H5::PredType::C_S1, H5T_VARIABLE), H5::DataSpace(1, &dim)); | |
288 | ✗ | boundary.write(boundary.getDataType(), boundaryStrings.data()); | |
289 | ✗ | H5::Attribute edges = | |
290 | ✗ | boxGroup.createAttribute("edges", H5::PredType::NATIVE_DOUBLE, H5::DataSpace(1, &dim)); | |
291 | ✗ | edges.write(edges.getDataType(), system->L.data()); | |
292 | // init trajectory props | ||
293 | // note that data pointer must be set before each use, see below in writeParticles() | ||
294 | ✗ | trajectoryQuantities = { | |
295 | // {name, {hdf5 datatype, dimension}} | ||
296 | {"position", {H5::PredType::NATIVE_DOUBLE, DIM}}, | ||
297 | {"image", {H5::PredType::NATIVE_INT, DIM}}, | ||
298 | {"velocity", {H5::PredType::NATIVE_DOUBLE, DIM}}, | ||
299 | {"force", {H5::PredType::NATIVE_DOUBLE, DIM}}, | ||
300 | {"Fext", {H5::PredType::NATIVE_DOUBLE, DIM}}, | ||
301 | {"Fint", {H5::PredType::NATIVE_DOUBLE, DIM}}, | ||
302 | {"rRand", {H5::PredType::NATIVE_DOUBLE, DIM}}, | ||
303 | {"mass", {H5::PredType::NATIVE_DOUBLE, 1}}, | ||
304 | {"species", {H5::PredType::NATIVE_INT, 1}}, | ||
305 | ✗ | }; | |
306 | // particle datasets | ||
307 | ✗ | hsize_t dimsfChunk[1] = {1}; | |
308 | ✗ | H5::DSetCreatPropList cparms; | |
309 | ✗ | cparms.setChunk(1, dimsfChunk); | |
310 | ✗ | hsize_t dimsfInitial[1] = {0}; | |
311 | ✗ | hsize_t maxdimsf[1] = {H5S_UNLIMITED}; | |
312 | ✗ | H5::DataSet datasetTime = particlesGroup.createDataSet( | |
313 | ✗ | "time", H5::PredType::NATIVE_DOUBLE, H5::DataSpace(1, dimsfInitial, maxdimsf), cparms); | |
314 | ✗ | for (const auto &qkey : trajectories->quantities) { | |
315 | ✗ | auto &props = trajectoryQuantities.at(qkey); | |
316 | ✗ | H5::Group group = particlesGroup.createGroup(qkey); | |
317 | ✗ | group.link(H5L_TYPE_HARD, "/particles/time", "time"); | |
318 | ✗ | hsize_t qdimsfChunk[3] = {dimsfChunk[0], n, props.dim}; | |
319 | ✗ | H5::DSetCreatPropList qcparms; | |
320 | ✗ | qcparms.setChunk(3, qdimsfChunk); | |
321 | ✗ | hsize_t qdimsfInitial[3] = {dimsfInitial[0], n, props.dim}; | |
322 | ✗ | hsize_t qmaxdimsf[3] = {H5S_UNLIMITED, H5S_UNLIMITED, props.dim}; | |
323 | ✗ | props.dataset = group.createDataSet("value", props.type, | |
324 | ✗ | H5::DataSpace(3, qdimsfInitial, qmaxdimsf), qcparms); | |
325 | ✗ | }; | |
326 | } | ||
327 | |||
328 | ✗ | void H5MDOutput::writeTrajectoryFrame(Trajectories *trajectories, System *system, int stage) { | |
329 | ✗ | hsize_t n = system->particles.size(); | |
330 | ✗ | if (trajectories->count == 0) { | |
331 | ✗ | initFile(outputPath, std::to_string(stage) + ".h5md"); | |
332 | ✗ | initParticlesGroup(trajectories, system); | |
333 | } | ||
334 | ✗ | if (trajectoryWriter.joinable()) { | |
335 | ✗ | trajectoryWriter.join(); | |
336 | } | ||
337 | // copy trajectory quantities into contiguous buffers | ||
338 | ✗ | trajectoryBuffers.r.reserve(n); | |
339 | ✗ | trajectoryBuffers.box.reserve(n); | |
340 | ✗ | trajectoryBuffers.v.reserve(n); | |
341 | ✗ | trajectoryBuffers.F.reserve(n); | |
342 | ✗ | trajectoryBuffers.Fext.reserve(n); | |
343 | ✗ | trajectoryBuffers.Fint.reserve(n); | |
344 | ✗ | trajectoryBuffers.rRand.reserve(n); | |
345 | ✗ | trajectoryBuffers.m.reserve(n); | |
346 | ✗ | trajectoryBuffers.species.reserve(n); | |
347 | ✗ | for (size_t i = 0; i < n; ++i) { | |
348 | ✗ | Particle &p = system->particles[i]; | |
349 | ✗ | trajectoryBuffers.r[i] = p.r; | |
350 | ✗ | trajectoryBuffers.box[i] = p.box; | |
351 | ✗ | trajectoryBuffers.v[i] = p.v; | |
352 | ✗ | trajectoryBuffers.F[i] = p.F(); | |
353 | ✗ | trajectoryBuffers.Fext[i] = p.Fext; | |
354 | ✗ | trajectoryBuffers.Fint[i] = p.Fint; | |
355 | ✗ | trajectoryBuffers.rRand[i] = p.rRand; | |
356 | ✗ | trajectoryBuffers.m[i] = p.m; | |
357 | ✗ | trajectoryBuffers.species[i] = p.species; | |
358 | } | ||
359 | ✗ | double t = system->t; | |
360 | // this may invalidate the pointers, so set them again | ||
361 | ✗ | trajectoryQuantities.at("position").data = trajectoryBuffers.r.data(); | |
362 | ✗ | trajectoryQuantities.at("image").data = trajectoryBuffers.box.data(); | |
363 | ✗ | trajectoryQuantities.at("velocity").data = trajectoryBuffers.v.data(); | |
364 | ✗ | trajectoryQuantities.at("force").data = trajectoryBuffers.F.data(); | |
365 | ✗ | trajectoryQuantities.at("Fext").data = trajectoryBuffers.Fext.data(); | |
366 | ✗ | trajectoryQuantities.at("Fint").data = trajectoryBuffers.Fint.data(); | |
367 | ✗ | trajectoryQuantities.at("rRand").data = trajectoryBuffers.rRand.data(); | |
368 | ✗ | trajectoryQuantities.at("mass").data = trajectoryBuffers.m.data(); | |
369 | ✗ | trajectoryQuantities.at("species").data = trajectoryBuffers.species.data(); | |
370 | // write to file | ||
371 | ✗ | trajectoryWriter = std::thread([this, trajectories, n, t]() { | |
372 | ✗ | H5::Group particlesGroup = file.openGroup("particles"); | |
373 | ✗ | hsize_t dimsf[1] = {(hsize_t)trajectories->count}; | |
374 | ✗ | hsize_t offset[1] = {dimsf[0] - 1}; | |
375 | ✗ | hsize_t dimsfExtend[1] = {1}; | |
376 | ✗ | H5::DataSpace mExtend(1, dimsfExtend); | |
377 | ✗ | H5::DataSet datasetTime = particlesGroup.openDataSet("time"); | |
378 | ✗ | datasetTime.extend(dimsf); | |
379 | ✗ | H5::DataSpace fExtendTime = datasetTime.getSpace(); | |
380 | ✗ | fExtendTime.selectHyperslab(H5S_SELECT_SET, dimsfExtend, offset); | |
381 | ✗ | datasetTime.write(&t, datasetTime.getDataType(), mExtend, fExtendTime); | |
382 | ✗ | for (const auto &qkey : trajectories->quantities) { | |
383 | ✗ | const auto &props = trajectoryQuantities.at(qkey); | |
384 | ✗ | hsize_t qdimsf[3] = {dimsf[0], n, props.dim}; | |
385 | ✗ | hsize_t qoffset[3] = {offset[0], 0, 0}; | |
386 | ✗ | hsize_t qdimsfExtend[3] = {1, n, props.dim}; | |
387 | ✗ | H5::DataSpace qmExtend(3, qdimsfExtend); | |
388 | ✗ | props.dataset.extend(qdimsf); | |
389 | ✗ | H5::DataSpace fExtendValue = props.dataset.getSpace(); | |
390 | ✗ | fExtendValue.selectHyperslab(H5S_SELECT_SET, qdimsfExtend, qoffset); | |
391 | ✗ | props.dataset.write(props.data, props.dataset.getDataType(), qmExtend, fExtendValue); | |
392 | } | ||
393 | ✗ | }); | |
394 | } | ||
395 | |||
396 | ✗ | ASCIIOutput::ASCIIOutput(const std::filesystem::path &outputPath, bool dryrun) | |
397 | ✗ | : Output(outputPath, dryrun) {} | |
398 | |||
399 | ✗ | void ASCIIOutput::writeObservables(Observables *observables, int stage) { | |
400 | ✗ | std::string filename = std::to_string(stage) + "_" + observables->type() + ".dat"; | |
401 | ✗ | std::vector<int> shape = observables->shape(); | |
402 | ✗ | std::vector<std::string> axisLabels = observables->axisLabels(); | |
403 | ✗ | std::ofstream coordFile; | |
404 | ✗ | if (shape.size()) { | |
405 | ✗ | coordFile = std::ofstream(outputPath / (filename + ".coord")); | |
406 | } | ||
407 | ✗ | for (size_t d = 0; d < shape.size(); ++d) { | |
408 | ✗ | std::vector<double> coords = observables->coordinates(d); | |
409 | ✗ | coordFile << axisLabels[d] << ": "; | |
410 | ✗ | for (const double r : coords) { | |
411 | ✗ | coordFile << r << " "; | |
412 | } | ||
413 | ✗ | coordFile << "\n"; | |
414 | } | ||
415 | ✗ | if (shape.size()) { | |
416 | ✗ | coordFile.close(); | |
417 | } | ||
418 | ✗ | std::ofstream file(outputPath / filename, std::ios::trunc); | |
419 | ✗ | Eigen::Array<int, Eigen::Dynamic, 1> speciesIndexShape(observables->numSpeciesIndices); | |
420 | ✗ | speciesIndexShape.setConstant(observables->numSpecies); | |
421 | ✗ | Multiindex<Eigen::Dynamic> speciesMultiindex(speciesIndexShape); | |
422 | ✗ | for (const auto &key : observables->keys) { | |
423 | ✗ | for (int linSpecies = 0; linSpecies < speciesMultiindex.linsize; ++linSpecies) { | |
424 | ✗ | Eigen::Array<int, Eigen::Dynamic, 1> species = speciesMultiindex.lin2multi(linSpecies); | |
425 | ✗ | file << key; | |
426 | ✗ | for (int s = 0; s < species.size(); ++s) { | |
427 | ✗ | file << " " << species[s]; | |
428 | } | ||
429 | ✗ | file << ":"; | |
430 | ✗ | std::vector<double> &data = observables->data(key, species); | |
431 | ✗ | int datasize = data.size(); | |
432 | ✗ | for (int i = 0; i < datasize; ++i) { | |
433 | ✗ | file << " " << data[i]; | |
434 | } | ||
435 | ✗ | file << "\n"; | |
436 | } | ||
437 | } | ||
438 | ✗ | file.close(); | |
439 | } | ||
440 | |||
441 | ✗ | static size_t countLines(const std::string &filename) { | |
442 | ✗ | size_t offset = 0; | |
443 | ✗ | std::ifstream infile(filename); | |
444 | ✗ | std::string dump; | |
445 | ✗ | while (std::getline(infile, dump)) { | |
446 | ✗ | ++offset; | |
447 | } | ||
448 | ✗ | return offset; | |
449 | } | ||
450 | |||
451 | ✗ | void ASCIIOutput::writeTimeseries(Timeseries *timeseries, int stage) { | |
452 | ✗ | std::string filename = std::to_string(stage) + "_timeseries.dat"; | |
453 | ✗ | size_t offset = 0; | |
454 | ✗ | std::ofstream file; | |
455 | ✗ | if (std::filesystem::exists(filename)) { | |
456 | ✗ | offset = countLines(outputPath / filename); | |
457 | ✗ | file.open(outputPath / filename, std::ios::app); | |
458 | } else { | ||
459 | ✗ | file.open(outputPath / filename, std::ios::out); | |
460 | ✗ | file << "time "; | |
461 | ✗ | for (const auto &[key, _] : timeseries->data) { | |
462 | ✗ | file << key << " "; | |
463 | } | ||
464 | ✗ | file << "\n"; | |
465 | } | ||
466 | ✗ | for (size_t i = offset; i < timeseries->time.size(); ++i) { | |
467 | ✗ | file << timeseries->time[i] << " "; | |
468 | ✗ | for (const auto &[_, data] : timeseries->data) { | |
469 | ✗ | file << data[i] << " "; | |
470 | } | ||
471 | ✗ | file << "\n"; | |
472 | } | ||
473 | ✗ | file.close(); | |
474 | } | ||
475 |