GCC Code Coverage Report


Directory: ./
File: src/output.cpp
Date: 2025-02-03 10:58:24
Exec Total Coverage
Lines: 48 346 13.9%
Functions: 8 26 30.8%
Branches: 28 102 27.5%

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