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