diff --git a/README.md b/README.md index d9858f2b..3b19ffd5 100644 --- a/README.md +++ b/README.md @@ -71,15 +71,20 @@ A single `.h5` containing particle and conjunction data from a full simulation r │ idA idB distanceSquared └── ParticleData └── - └── Particles - ├── (Dataset) IDs - ├── (Dataset) Positions - │ x y z - └── (Dataset) Velocities - x y z + │ └── Particles + │ ├── (Dataset) IDs + │ ├── (Dataset) Positions + │ │ x y z + │ └── (Dataset) Velocities + │ x y z + └─── (Dataset) ConstantProperties + id mass radius activityState ``` -Collision data is tracked every iteration, particle data only in intervals that are defined in the YAML file. To keep file size reasonable compression is supported. +Collision data is tracked every iteration, particle data only in intervals that are defined in the YAML file. +`ConstantProperties` contains properties of all particles that existed over the course of the simulation. +Due to burn ups or breakups the number of particles in any iteration might differ but `id`s are unique! +To keep file size reasonable compression is supported. ### CSV If HDF5 output is disabled entirely, collision data is written in a `.csv` file in ASCII layout. diff --git a/cfg/default_cfg.yaml b/cfg/default_cfg.yaml index a2467baa..75e3ed46 100644 --- a/cfg/default_cfg.yaml +++ b/cfg/default_cfg.yaml @@ -5,6 +5,8 @@ sim: minAltitude: 150 # Everything below this altitude above ground will be considered burning up [km] deltaT: 10.0 # [s] conjunctionThreshold: 0.5 # [km] + # TODO: should this go into the prop section? + # coefficientOfDrag: 2.2 # c_D for the drag component in the propagator TODO: more doc, where is it used to derive what? prop: # Which propagation model components should be applied useKEPComponent: true # Keplerian propagation diff --git a/cmake/modules/propagator.cmake b/cmake/modules/propagator.cmake index 20663521..30bdfe51 100644 --- a/cmake/modules/propagator.cmake +++ b/cmake/modules/propagator.cmake @@ -5,9 +5,9 @@ message(STATUS "Adding Propagator.") include(FetchContent) # Select https (default) or ssh path. -set(PropagatorRepoPath https://github.com/Wombatwarrior/BA_space_debris.git) +set(PropagatorRepoPath https://github.com/FG-TUM/OrbitPropagator.git) if (GIT_SUBMODULES_SSH) - set(PropagatorRepoPath git@github.com:Wombatwarrior/BA_space_debris.git) + set(PropagatorRepoPath git@github.com:FG-TUM/OrbitPropagator.git) endif () # disable stuff we do not need @@ -18,8 +18,8 @@ set(DebrisSim_Simulation OFF CACHE INTERNAL "") FetchContent_Declare( Propagatorfetch GIT_REPOSITORY ${PropagatorRepoPath} - # branch: IntegratorComparision; OpenMP parallelization - GIT_TAG 58ebb3a7d902a5adef9c9e4a9a62bb87367aa2b9 + # branch: main; latest feature: ActivityState + GIT_TAG 9f5be5c079b2ff9ff6e64e34751f1429b61c8b46 ) # Get Propagator source and binary directories from CMake project diff --git a/src/ladds/CollisionFunctor.cpp b/src/ladds/CollisionFunctor.cpp index bab04e74..42b3159d 100644 --- a/src/ladds/CollisionFunctor.cpp +++ b/src/ladds/CollisionFunctor.cpp @@ -11,8 +11,12 @@ #include -CollisionFunctor::CollisionFunctor(double cutoff, double dt, double minorCutoff) - : Functor(cutoff), _cutoffSquare(cutoff * cutoff), _dt(dt), _minorCutoffSquare(minorCutoff * minorCutoff) { +CollisionFunctor::CollisionFunctor(double cutoff, double dt, double minorCutoff, double minDetectionRadius) + : Functor(cutoff), + _cutoffSquare(cutoff * cutoff), + _dt(dt), + _minorCutoffSquare(minorCutoff * minorCutoff), + _minDetectionRadius(minDetectionRadius) { _threadData.resize(autopas::autopas_get_max_threads()); } @@ -26,10 +30,20 @@ void CollisionFunctor::AoSFunctor(Particle &i, Particle &j, bool newton3) { using autopas::utils::ArrayMath::mulScalar; using autopas::utils::ArrayMath::sub; - // skip interaction with deleted particles + // skip if interaction involves: + // - any deleted particles + // - two actively evasive satellites + // - one evasive and one of detectable size if (i.isDummy() or j.isDummy()) { return; } + const auto &iActivity = i.getActivityState(); + const auto &jActivity = j.getActivityState(); + if ((iActivity != Particle::ActivityState::passive and jActivity != Particle::ActivityState::passive) or + (iActivity != Particle::ActivityState::passive and j.getRadius() >= _minDetectionRadius) or + (jActivity != Particle::ActivityState::passive and i.getRadius() >= _minDetectionRadius)) { + return; + } // calculate distance between particles const auto dr = sub(i.getR(), j.getR()); diff --git a/src/ladds/CollisionFunctor.h b/src/ladds/CollisionFunctor.h index bee99211..5a9b4d7f 100644 --- a/src/ladds/CollisionFunctor.h +++ b/src/ladds/CollisionFunctor.h @@ -20,7 +20,15 @@ */ class CollisionFunctor final : public autopas::Functor { public: - explicit CollisionFunctor(double cutoff, double dt, double minorCutoff); + /** + * Constructor + * @param cutoff Distance for two particles to be considered for sub time step investigation (interpolation). + * @param dt time step width + * @param minorCutoff Distance where we consider particles to be colliding. + * @param minDetectionRadius All particles with a larger are assumed to be detectable by radar. + * Thus collisions with particles that are Particle::ActivityState::evasive will not be considered. + */ + CollisionFunctor(double cutoff, double dt, double minorCutoff, double minDetectionRadius); [[nodiscard]] bool isRelevantForTuning() final { return true; @@ -91,4 +99,9 @@ class CollisionFunctor final : public autopas::Functor>(autopas, csvFilename, OutputFile::CSV, selectedPropagatorComponents); } auto accumulator = std::make_unique>( - selectedPropagatorComponents, autopas, 0.0, *csvWriter); + selectedPropagatorComponents, autopas, 0.0, csvWriter.get()); auto deltaT = config.get("sim/deltaT"); auto integrator = std::make_unique>(autopas, *accumulator, deltaT); @@ -152,9 +152,10 @@ void Simulation::updateConstellation(AutoPas_t &autopas, std::tuple Simulation::collisionDetection(AutoPas_t &autopas, double deltaT, - double conjunctionThreshold) { + double conjunctionThreshold, + double minDetectionRadius) { // pairwise interaction - CollisionFunctor collisionFunctor(autopas.getCutoff(), deltaT, conjunctionThreshold); + CollisionFunctor collisionFunctor(autopas.getCutoff(), deltaT, conjunctionThreshold, minDetectionRadius); bool stillTuning = autopas.iteratePairwise(&collisionFunctor); return {collisionFunctor.getCollisions(), stillTuning}; } @@ -170,6 +171,7 @@ void Simulation::simulationLoop(AutoPas_t &autopas, const auto progressOutputFrequency = config.get("io/progressOutputFrequency", 50); const auto deltaT = config.get("sim/deltaT"); const auto conjunctionThreshold = config.get("sim/conjunctionThreshold"); + const auto minDetectionRadius = config.get("sim/minDetectionRadius", 0.1); const auto minAltitude = config.get("sim/minAltitude", 150.); std::vector delayedInsertion; @@ -228,7 +230,7 @@ void Simulation::simulationLoop(AutoPas_t &autopas, } timers.collisionDetection.start(); - auto [collisions, stillTuning] = collisionDetection(autopas, deltaT, conjunctionThreshold); + auto [collisions, stillTuning] = collisionDetection(autopas, deltaT, conjunctionThreshold, minDetectionRadius); timers.collisionDetection.stop(); if (tuningMode and not stillTuning) { diff --git a/src/ladds/Simulation.h b/src/ladds/Simulation.h index 1498ed7a..9854ee7d 100644 --- a/src/ladds/Simulation.h +++ b/src/ladds/Simulation.h @@ -82,7 +82,8 @@ class Simulation { */ std::tuple collisionDetection(AutoPas_t &autopas, double deltaT, - double conjunctionThreshold); + double conjunctionThreshold, + double minDetectionRadius); /** * Updates the configuration with the latest AutoPas configuration and writes it to a new YAML file. diff --git a/src/ladds/io/ConfigReader.cpp b/src/ladds/io/ConfigReader.cpp index d15aba58..52b15ec5 100644 --- a/src/ladds/io/ConfigReader.cpp +++ b/src/ladds/io/ConfigReader.cpp @@ -54,37 +54,6 @@ bool ConfigReader::defines(const std::string &valuePath, bool suppressWarning) { return true; } -/** - * Helper function for setValue recursively iterating through the yaml structure, - * inserting missing nodes and setting or updating the given value. - * @tparam T - * @tparam Iter - * @param node This is passed by copy as it works similar to a smart pointer. - * @param begin Start of the array containing the key path. - * @param end End of the array containing the key path. - * @param value - */ -template -void setValueAux(YAML::Node node, Iter begin, Iter end, T value) { - if (begin == end) { - return; - } - const auto &tag = *begin; - if (std::next(begin) == end) { - node[tag] = value; - return; - } - if (!node[tag]) { - node[tag] = YAML::Node(YAML::NodeType::Map); - } - setValueAux(node[tag], std::next(begin), end, value); -} - -void ConfigReader::setValue(const std::string &valuePath, const std::string &value) { - std::vector valuePathVec = autopas::utils::StringUtils::tokenize(valuePath, "/"); - setValueAux(config, valuePathVec.begin(), valuePathVec.end(), value); -} - void ConfigReader::dumpConfig(const std::string &filename) const { std::ofstream file(filename); file << config << std::endl; diff --git a/src/ladds/io/ConfigReader.h b/src/ladds/io/ConfigReader.h index 39d69af1..e9c9e200 100644 --- a/src/ladds/io/ConfigReader.h +++ b/src/ladds/io/ConfigReader.h @@ -38,7 +38,8 @@ class ConfigReader { * Primary function to retrieve values from the underlying YAML data. * @tparam T Type as which the retrieved data shall be read. * @param valuePath Path within the YAML file to the data field. - * @param fallbackValue Value to use if the field can not be found. If std::nullopt is used an exception is thrown. + * @param fallbackValue Value to use if the field can not be found. + * If std::nullopt is used an exception is thrown when no value was found. * @param suppressWarning If true no warning is logged when the default value is used. * @param fallbackToString If the fallback value is non trivial a function can be passed that provides a string * representation. @@ -68,6 +69,8 @@ class ConfigReader { fallbackValue.value()); } parsedValues.emplace(valuePath, fallbackToString(fallbackValue.value())); + // write to value to the datastructure so successive calls will return the same fallback. + setValue(valuePath, fallbackValue.value()); return fallbackValue.value(); } else { throw std::runtime_error("Config option not found: " + valuePath); @@ -102,12 +105,54 @@ class ConfigReader { /** * Update or insert a value for a given key / value path. * If any Node along the path does not exist it will be created. + * @tparam T Type of the value to write. Has to be std::string, bool, or convertible via std::to_string(). * @param valuePath * @param value */ - void setValue(const std::string &valuePath, const std::string &value); + template + void setValue(const std::string &valuePath, const T &value) { + std::vector valuePathVec = autopas::utils::StringUtils::tokenize(valuePath, "/"); + // try to convert value to a string if it isn't already std::string or string literal + if constexpr (std::is_same_v or std::is_same_v]>) { + // only set non-empty strings + if (not static_cast(value).empty()) { + setValueAux(config, valuePathVec.begin(), valuePathVec.end(), value); + } + } else if constexpr (std::is_same_v) { + const auto &v = static_cast(value); + setValueAux(config, valuePathVec.begin(), valuePathVec.end(), value ? "true" : "false"); + } else { + setValueAux(config, valuePathVec.begin(), valuePathVec.end(), std::to_string(value)); + } + } private: + /** + * Helper function for setValue recursively iterating through the yaml structure, + * inserting missing nodes and setting or updating the given value. + * @tparam T + * @tparam Iter + * @param node This is passed by copy as it works similar to a smart pointer. + * @param begin Start of the array containing the key path. + * @param end End of the array containing the key path. + * @param value + */ + template + void setValueAux(YAML::Node node, Iter begin, Iter end, T value) { + if (begin == end) { + return; + } + const auto &tag = *begin; + if (std::next(begin) == end) { + node[tag] = value; + return; + } + if (!node[tag]) { + node[tag] = YAML::Node(YAML::NodeType::Map); + } + setValueAux(node[tag], std::next(begin), end, value); + } + /** * Loads the config file, falls back to default cfg if not found. */ diff --git a/src/ladds/io/DatasetReader.cpp b/src/ladds/io/DatasetReader.cpp index e542528c..2caa1e08 100644 --- a/src/ladds/io/DatasetReader.cpp +++ b/src/ladds/io/DatasetReader.cpp @@ -6,18 +6,18 @@ #include "DatasetReader.h" -std::vector DatasetReader::readDataset(const std::string &csvFilepath) { - CSVReader DatasetReader::readDataset(const std::string &csvFilepath, double coefficientOfDrag) { + CSVReader csvReader{csvFilepath, true}; const auto parsedData = csvReader.getLines(); @@ -27,7 +27,13 @@ std::vector DatasetReader::readDataset(const std::string &csvFilepath) size_t particleId = 0; for (const auto &[id, name, mass, radius, activityState, rX, rY, rZ, vX, vY, vZ] : parsedData) { - particleCollection.emplace_back(std::array{rX, rY, rZ}, std::array{vX, vY, vZ}, id); + particleCollection.emplace_back(std::array{rX, rY, rZ}, + std::array{vX, vY, vZ}, + id, + activityState, + mass, + radius, + coefficientOfDrag); } return particleCollection; diff --git a/src/ladds/io/DatasetReader.h b/src/ladds/io/DatasetReader.h index 1ae5261e..66b63e4c 100644 --- a/src/ladds/io/DatasetReader.h +++ b/src/ladds/io/DatasetReader.h @@ -8,6 +8,9 @@ namespace DatasetReader { /** * Reads the passed position and velocity csv files. Returns a vector of particles. + * @param csvFilepath + * @param coefficientOfDrag c_D used to initialize all loaded particles. + * @return */ -std::vector readDataset(const std::string &csvFilepath); +std::vector readDataset(const std::string &csvFilepath, double coefficientOfDrag); } // namespace DatasetReader \ No newline at end of file diff --git a/src/ladds/io/SatelliteLoader.cpp b/src/ladds/io/SatelliteLoader.cpp index e5f650fa..37bfe812 100644 --- a/src/ladds/io/SatelliteLoader.cpp +++ b/src/ladds/io/SatelliteLoader.cpp @@ -14,6 +14,9 @@ void SatelliteLoader::loadSatellites(AutoPas_t &autopas, ConfigReader &config, const Logger &logger) { std::vector satellites; + // if the value is not set this is the only place that should define the default value as it is the first to access it + const auto coefficientOfDrag = config.get("sim/coefficientOfDrag", 2.2); + // load CSV ... const auto csvFileName = config.get("io/csv/fileName", ""); if (not csvFileName.empty()) { @@ -22,7 +25,7 @@ void SatelliteLoader::loadSatellites(AutoPas_t &autopas, ConfigReader &config, c SPDLOG_LOGGER_INFO(logger.get(), "Loading scenario from CSV: {}", csvFilePath); // Read in scenario - satellites = DatasetReader::readDataset(csvFilePath); + satellites = DatasetReader::readDataset(csvFilePath, coefficientOfDrag); SPDLOG_LOGGER_DEBUG(logger.get(), "Parsed {} satellites", satellites.size()); } else { // ... or load checkpoint ... @@ -34,7 +37,7 @@ void SatelliteLoader::loadSatellites(AutoPas_t &autopas, ConfigReader &config, c HDF5Reader hdfReader(checkpointPath); // either load the given iteration or fall back to the last iteration stored in the file auto iteration = config.get("io/checkpoint/iteration", hdfReader.readLastIterationNr()); - satellites = hdfReader.readParticles(iteration); + satellites = hdfReader.readParticles(iteration, coefficientOfDrag); } else { // ... or fail throw std::runtime_error("No valid input option found! Exiting..."); @@ -73,6 +76,7 @@ void SatelliteLoader::loadSatellites(AutoPas_t &autopas, ConfigReader &config, c std::vector SatelliteLoader::loadConstellations(ConfigReader &config, const Logger &logger) { std::vector constellations; + auto coefficientOfDrag = config.get("sim/coefficientOfDrag"); auto constellationList = config.get("io/constellationList", "", true); // altitudeSpread = 3 * sigma , altitudeDeviation = sigma (= standardDeviation) auto altitudeDeviation = config.get("io/altitudeSpread", 0.0) / 3.0; @@ -92,11 +96,12 @@ std::vector SatelliteLoader::loadConstellations(ConfigReader &con for (int i = 0; i < nConstellations; ++i) { unsigned long offset = constellationDataStr.find(';', 0); if (offset == 0) { - constellations.emplace_back(Constellation(constellationDataStr, insertionFrequency, altitudeDeviation)); + constellations.emplace_back( + Constellation(constellationDataStr, insertionFrequency, altitudeDeviation, coefficientOfDrag)); break; } else { - constellations.emplace_back( - Constellation(constellationDataStr.substr(0, offset), insertionFrequency, altitudeDeviation)); + constellations.emplace_back(Constellation( + constellationDataStr.substr(0, offset), insertionFrequency, altitudeDeviation, coefficientOfDrag)); constellationDataStr.erase(0, offset + 1); } } diff --git a/src/ladds/io/hdf5/HDF5Definitions.h b/src/ladds/io/hdf5/HDF5Definitions.h new file mode 100644 index 00000000..c89abb40 --- /dev/null +++ b/src/ladds/io/hdf5/HDF5Definitions.h @@ -0,0 +1,57 @@ +/** + * @file HDF5Definitions.h + * @author F. Gratl + * @date 19.01.22 + */ + +#pragma once + +namespace HDF5Definitions { +constexpr auto groupParticleData = "ParticleData/"; +constexpr auto datasetParticlePositions = "/Particles/Positions"; +constexpr auto datasetParticleVelocities = "/Particles/Velocities"; +constexpr auto datasetParticleIDs = "/Particles/IDs"; +constexpr auto tableParticleConstantProperties = "ConstantProperties"; +constexpr auto groupCollisionData = "CollisionData/"; +constexpr auto datasetCollisions = "/Collisions"; + +/** + * Type to which any floating point data will be cast before writing. + */ +using FloatType = float; +/** + * Type to which any integer data will be cast before writing. + */ +using IntType = unsigned int; + +/** + * This represents one line of 3D vector data in the HDF5 file. + */ +template +struct Vec3 { + T x, y, z; +}; + +/** + * Type for the information of a single collision. + */ +struct CollisionInfo { + unsigned int idA, idB; + float distanceSquared; + bool operator==(const CollisionInfo &rhs) const { + return idA == rhs.idA && idB == rhs.idB && distanceSquared == rhs.distanceSquared; + } + bool operator!=(const CollisionInfo &rhs) const { + return !(rhs == *this); + } +}; + +/** + * Type for information of a single particle that will stay constant throughout the simulation. + */ +struct ParticleConstantProperties { + IntType id; + FloatType mass, radius; + int activityState; +}; +} // namespace HDF5Definitions diff --git a/src/ladds/io/hdf5/HDF5Reader.cpp b/src/ladds/io/hdf5/HDF5Reader.cpp index 88f58f24..df24ae82 100644 --- a/src/ladds/io/hdf5/HDF5Reader.cpp +++ b/src/ladds/io/hdf5/HDF5Reader.cpp @@ -6,7 +6,8 @@ #include "HDF5Reader.h" -#include "HDF5Writer.h" +#include "HDF5Definitions.h" + HDF5Reader::HDF5Reader(const std::string &filename) #ifdef LADDS_HDF5 : file(filename, h5pp::FilePermission::READONLY) @@ -17,31 +18,51 @@ HDF5Reader::HDF5Reader(const std::string &filename) #endif } -std::vector HDF5Reader::readParticles(size_t iteration) const { +std::vector HDF5Reader::readParticles(size_t iteration, double coefficientOfDrag) const { std::vector particles{}; #ifdef LADDS_HDF5 - const auto pos = file.readDataset>>( - HDF5Writer::groupParticleData + std::to_string(iteration) + HDF5Writer::datasetParticlePositions); - const auto vel = file.readDataset>>( - HDF5Writer::groupParticleData + std::to_string(iteration) + HDF5Writer::datasetParticleVelocities); - const auto ids = file.readDataset>( - HDF5Writer::groupParticleData + std::to_string(iteration) + HDF5Writer::datasetParticleIDs); + const auto pos = file.readDataset>>( + HDF5Definitions::groupParticleData + std::to_string(iteration) + HDF5Definitions::datasetParticlePositions); + const auto vel = file.readDataset>>( + HDF5Definitions::groupParticleData + std::to_string(iteration) + HDF5Definitions::datasetParticleVelocities); + const auto ids = file.readDataset>( + HDF5Definitions::groupParticleData + std::to_string(iteration) + HDF5Definitions::datasetParticleIDs); + const auto pathParticleConstantProperties = + std::string(HDF5Definitions::groupParticleData) + HDF5Definitions::tableParticleConstantProperties; + const auto constantPropertiesVec = + file.readTableRecords>(pathParticleConstantProperties); + + // convert static data to a map for easier access + const auto constantPropertiesMap = [&]() { + // id is redundant in this datastructure but easier maintainable this way + std::unordered_map map; + map.reserve(constantPropertiesVec.size()); + for (const auto &data : constantPropertiesVec) { + map[data.id] = data; + } + return map; + }(); particles.reserve(pos.size()); for (size_t i = 0; i < pos.size(); ++i) { + const auto &constantProperties = constantPropertiesMap.at(ids[i]); particles.emplace_back(std::array{pos[i].x, pos[i].y, pos[i].z}, std::array{vel[i].x, vel[i].y, vel[i].z}, - ids[i]); + constantProperties.id, + static_cast(constantProperties.activityState), + constantProperties.mass, + constantProperties.radius, + coefficientOfDrag); } #endif return particles; } -std::vector HDF5Reader::readCollisions(size_t iteration) const { - std::vector collisions{}; +std::vector HDF5Reader::readCollisions(size_t iteration) const { + std::vector collisions{}; #ifdef LADDS_HDF5 - file.readDataset(collisions, - HDF5Writer::groupCollisionData + std::to_string(iteration) + HDF5Writer::datasetCollisions); + file.readDataset( + collisions, HDF5Definitions::groupCollisionData + std::to_string(iteration) + HDF5Definitions::datasetCollisions); #endif return collisions; } @@ -49,9 +70,9 @@ std::vector HDF5Reader::readCollisions(size_t iterati size_t HDF5Reader::readLastIterationNr() { std::vector allIterations; #ifdef LADDS_HDF5 - const auto groupStrLength = std::string(HDF5Writer::groupParticleData).size(); + const auto groupStrLength = std::string(HDF5Definitions::groupParticleData).size(); // gather all ID datasets. They contain the iteration number in their path - auto allIdDatasets = file.findDatasets(HDF5Writer::datasetParticleIDs); + auto allIdDatasets = file.findDatasets(HDF5Definitions::datasetParticleIDs); allIterations.reserve(allIdDatasets.size()); // extract all iteration numbers as size_t std::transform(allIdDatasets.begin(), allIdDatasets.end(), std::back_inserter(allIterations), [&](auto &datasetName) { diff --git a/src/ladds/io/hdf5/HDF5Reader.h b/src/ladds/io/hdf5/HDF5Reader.h index 6b947924..476ffd5b 100644 --- a/src/ladds/io/hdf5/HDF5Reader.h +++ b/src/ladds/io/hdf5/HDF5Reader.h @@ -12,15 +12,22 @@ #include -#include "HDF5Writer.h" +#include "HDF5Definitions.h" #include "ladds/particle/Particle.h" + class HDF5Reader { public: explicit HDF5Reader(const std::string &filename); - std::vector readParticles(size_t iteration) const; + /** + * Load a set of particles from the HDF5 file from a given iteration. + * @param iteration + * @param coefficientOfDrag c_D used to initialize all loaded particles. + * @return + */ + std::vector readParticles(size_t iteration, double coefficientOfDrag) const; - std::vector readCollisions(size_t iteration) const; + std::vector readCollisions(size_t iteration) const; /** * Finds the last iteration that is stored in the file. diff --git a/src/ladds/io/hdf5/HDF5Writer.cpp b/src/ladds/io/hdf5/HDF5Writer.cpp index 611d9ef1..9dc8ba2d 100644 --- a/src/ladds/io/hdf5/HDF5Writer.cpp +++ b/src/ladds/io/hdf5/HDF5Writer.cpp @@ -6,26 +6,48 @@ #include "HDF5Writer.h" +#include "HDF5Definitions.h" + HDF5Writer::HDF5Writer(const std::string &filename, unsigned int compressionLevel) #ifdef LADDS_HDF5 : _file(filename, h5pp::FilePermission::REPLACE), - collisionInfoH5Type(H5Tcreate(H5T_COMPOUND, sizeof(CollisionInfo))) + collisionInfoH5Type(H5Tcreate(H5T_COMPOUND, sizeof(HDF5Definitions::CollisionInfo))), + particleConstantPropertiesH5Type(H5Tcreate(H5T_COMPOUND, sizeof(HDF5Definitions::ParticleConstantProperties))) #endif { #ifdef LADDS_HDF5 _file.setCompressionLevel(compressionLevel); + // CollisionInfo H5Tinsert(collisionInfoH5Type, "idA", - HOFFSET(CollisionInfo, idA), - h5pp::type::getH5NativeType()); + HOFFSET(HDF5Definitions::CollisionInfo, idA), + h5pp::type::getH5NativeType()); H5Tinsert(collisionInfoH5Type, "idB", - HOFFSET(CollisionInfo, idB), - h5pp::type::getH5NativeType()); + HOFFSET(HDF5Definitions::CollisionInfo, idB), + h5pp::type::getH5NativeType()); H5Tinsert(collisionInfoH5Type, "distanceSquared", - HOFFSET(CollisionInfo, distanceSquared), - h5pp::type::getH5NativeType()); + HOFFSET(HDF5Definitions::CollisionInfo, distanceSquared), + h5pp::type::getH5NativeType()); + + // ParticleConstantProperties + H5Tinsert(particleConstantPropertiesH5Type, + "id", + HOFFSET(HDF5Definitions::ParticleConstantProperties, id), + h5pp::type::getH5NativeType()); + H5Tinsert(particleConstantPropertiesH5Type, + "mass", + HOFFSET(HDF5Definitions::ParticleConstantProperties, mass), + h5pp::type::getH5NativeType()); + H5Tinsert(particleConstantPropertiesH5Type, + "radius", + HOFFSET(HDF5Definitions::ParticleConstantProperties, radius), + h5pp::type::getH5NativeType()); + H5Tinsert(particleConstantPropertiesH5Type, + "activityState", + HOFFSET(HDF5Definitions::ParticleConstantProperties, activityState), + h5pp::type::getH5NativeType()); #else throw std::runtime_error("LADDS was compiled without HDF5 support, so the HDF5Writer can't do anything!"); #endif @@ -34,9 +56,11 @@ HDF5Writer::HDF5Writer(const std::string &filename, unsigned int compressionLeve void HDF5Writer::writeParticles(size_t iteration, const AutoPas_t &autopas) { #ifdef LADDS_HDF5 - std::vector> vecPos; - std::vector> vecVel; - std::vector vecId; + std::vector> vecPos; + std::vector> vecVel; + std::vector vecId; + HDF5Definitions::IntType maxParticleId{0}; + std::vector newConstantProperties; vecPos.reserve(autopas.getNumberOfParticles()); vecVel.reserve(autopas.getNumberOfParticles()); @@ -45,19 +69,47 @@ void HDF5Writer::writeParticles(size_t iteration, const AutoPas_t &autopas) { for (const auto &particle : autopas) { const auto &pos = particle.getR(); const auto &vel = particle.getV(); + const auto &id = static_cast(particle.getID()); // pack data and make sure it is of the correct type - vecPos.emplace_back>( - {static_cast(pos[0]), static_cast(pos[1]), static_cast(pos[2])}); - vecVel.emplace_back>( - {static_cast(vel[0]), static_cast(vel[1]), static_cast(vel[2])}); - vecId.emplace_back(static_cast(particle.getID())); + vecPos.emplace_back>( + {static_cast(pos[0]), + static_cast(pos[1]), + static_cast(pos[2])}); + vecVel.emplace_back>( + {static_cast(vel[0]), + static_cast(vel[1]), + static_cast(vel[2])}); + vecId.emplace_back(id); + // track the highest particle id that was written to the file + // All particles that have a higher id than the highest id from the last time something was written are new + // and their static properties need to be recorded. + maxParticleId = std::max(maxParticleId, id); + if (maxWrittenParticleID == 0 or id > maxWrittenParticleID) { + newConstantProperties.emplace_back( + HDF5Definitions::ParticleConstantProperties{id, + static_cast(particle.getMass()), + static_cast(particle.getRadius()), + static_cast(particle.getActivityState())}); + } } - const auto group = groupParticleData + std::to_string(iteration); + const auto group = HDF5Definitions::groupParticleData + std::to_string(iteration); const auto &compressionLvl = _file.getCompressionLevel(); - _file.writeDataset_compressed(vecPos, group + datasetParticlePositions, compressionLvl); - _file.writeDataset_compressed(vecVel, group + datasetParticleVelocities, compressionLvl); - _file.writeDataset_compressed(vecId, group + datasetParticleIDs, compressionLvl); + _file.writeDataset_compressed(vecPos, group + HDF5Definitions::datasetParticlePositions, compressionLvl); + _file.writeDataset_compressed(vecVel, group + HDF5Definitions::datasetParticleVelocities, compressionLvl); + _file.writeDataset_compressed(vecId, group + HDF5Definitions::datasetParticleIDs, compressionLvl); + + if (not newConstantProperties.empty()) { + const auto particleConstantPropertiesFullPath = + std::string(HDF5Definitions::groupParticleData) + HDF5Definitions::tableParticleConstantProperties; + // it table does not exist yet create it + if (not _file.linkExists(particleConstantPropertiesFullPath)) { + _file.createTable(particleConstantPropertiesH5Type, particleConstantPropertiesFullPath, "ConstantProperties"); + } + _file.appendTableRecords(newConstantProperties, particleConstantPropertiesFullPath); + } + + maxWrittenParticleID = maxParticleId; #endif } @@ -66,22 +118,16 @@ void HDF5Writer::writeConjunctions(size_t iteration, const CollisionFunctor::Col if (collisions.empty()) { return; } - std::vector data; + std::vector data; data.reserve(collisions.size()); for (const auto &[p1, p2, distanceSquare] : collisions) { - data.emplace_back({static_cast(p1->getID()), - static_cast(p2->getID()), - static_cast(distanceSquare)}); + data.emplace_back( + {static_cast(p1->getID()), + static_cast(p2->getID()), + static_cast(distanceSquare)}); } - const auto group = groupCollisionData + std::to_string(iteration); - _file.writeDataset(data, group + datasetCollisions, collisionInfoH5Type); + const auto group = HDF5Definitions::groupCollisionData + std::to_string(iteration); + _file.writeDataset(data, group + HDF5Definitions::datasetCollisions, collisionInfoH5Type); #endif } - -bool HDF5Writer::CollisionInfo::operator==(const HDF5Writer::CollisionInfo &rhs) const { - return idA == rhs.idA && idB == rhs.idB && distanceSquared == rhs.distanceSquared; -} -bool HDF5Writer::CollisionInfo::operator!=(const HDF5Writer::CollisionInfo &rhs) const { - return !(rhs == *this); -} diff --git a/src/ladds/io/hdf5/HDF5Writer.h b/src/ladds/io/hdf5/HDF5Writer.h index b0500c17..2ec4743b 100644 --- a/src/ladds/io/hdf5/HDF5Writer.h +++ b/src/ladds/io/hdf5/HDF5Writer.h @@ -12,6 +12,7 @@ #include +#include "HDF5Definitions.h" #include "ladds/TypeDefinitions.h" #include "ladds/io/ConjunctionWriterInterface.h" @@ -23,6 +24,7 @@ class HDF5Writer final : public ConjuctionWriterInterface { public: /** * Constructor setting up the file and creating the custom data type. + * @note Here the actual file is created on disk. * @param filename * @param compressionLevel */ @@ -44,53 +46,24 @@ class HDF5Writer final : public ConjuctionWriterInterface { */ void writeConjunctions(size_t iteration, const CollisionFunctor::CollisionCollectionT &collisions) override; + private: /** - * Type to which any floating point data will be cast before writing. - */ - using FloatType = float; - /** - * Type to which any integer data will be cast before writing. - */ - using IntType = unsigned int; - - /** - * This represents one line of 3D vector data in the HDF5 file. - */ - template - struct Vec3 { - T x, y, z; - }; - - static constexpr auto groupParticleData = "ParticleData/"; - static constexpr auto datasetParticlePositions = "/Particles/Positions"; - static constexpr auto datasetParticleVelocities = "/Particles/Velocities"; - static constexpr auto datasetParticleIDs = "/Particles/IDs"; - - static constexpr auto groupCollisionData = "CollisionData/"; - static constexpr auto datasetCollisions = "/Collisions"; - - /** - * Type for the information of a single collision. + * Highest partilce ID that was written in any previous iteration. + * For anything below this ID constant particle properties are already written. */ - struct CollisionInfo { - unsigned int idA, idB; - float distanceSquared; - bool operator==(const CollisionInfo &rhs) const; - bool operator!=(const CollisionInfo &rhs) const; - }; - - private: + HDF5Definitions::IntType maxWrittenParticleID{0}; #ifdef LADDS_HDF5 /** * Actual file that will be created. All of the data this writer gets ends up in this one file. */ h5pp::File _file; -#endif - -#ifdef LADDS_HDF5 /** * Object holding the info for the hdf5 compound type of the collision data. */ h5pp::hid::h5t collisionInfoH5Type; + /** + * Object holding the info for the hdf5 compound type of the constant particle properties data. + */ + h5pp::hid::h5t particleConstantPropertiesH5Type; #endif }; diff --git a/src/ladds/particle/BreakupWrapper.cpp b/src/ladds/particle/BreakupWrapper.cpp index 00cf9fa6..c81c030c 100644 --- a/src/ladds/particle/BreakupWrapper.cpp +++ b/src/ladds/particle/BreakupWrapper.cpp @@ -27,7 +27,7 @@ void BreakupWrapper::simulateBreakup(const CollisionFunctor::CollisionCollection // insert resulting debris into the simulation const auto &newSatellites = breakup->getResult(); for (const auto &newSat : newSatellites) { - autopas.addParticle(SatelliteToParticleConverter::convertSatelliteToParticle(newSat)); + autopas.addParticle(SatelliteToParticleConverter::convertSatelliteToParticle(newSat, coefficientOfDrag)); } } } diff --git a/src/ladds/particle/BreakupWrapper.h b/src/ladds/particle/BreakupWrapper.h index 3e7ec3f5..6dab73eb 100644 --- a/src/ladds/particle/BreakupWrapper.h +++ b/src/ladds/particle/BreakupWrapper.h @@ -23,6 +23,7 @@ class BreakupWrapper { BreakupWrapper(ConfigReader &config, AutoPas_t &autopas) : minLc{config.get("sim/breakup/minLc", 0.05)}, enforceMassConservation{config.get("sim/breakup/enforceMassConservation", true)}, + coefficientOfDrag{config.get("sim/coefficientOfDrag")}, autopas{autopas} { // find the highest existing particle id: // Particles are not sorted by id and might neither be starting by 0 nor be consecutive (e.g. due to burn-ups) @@ -37,6 +38,7 @@ class BreakupWrapper { private: const double minLc; const bool enforceMassConservation; + double coefficientOfDrag; size_t maxExistingParticleId{0ul}; AutoPas_t &autopas; }; diff --git a/src/ladds/particle/Constellation.cpp b/src/ladds/particle/Constellation.cpp index d1130751..a74cd8ac 100644 --- a/src/ladds/particle/Constellation.cpp +++ b/src/ladds/particle/Constellation.cpp @@ -14,7 +14,10 @@ std::mt19937 Constellation::generator{42}; -Constellation::Constellation(const std::string &constellation_data_str, size_t interval, double altitudeDeviation) +Constellation::Constellation(const std::string &constellation_data_str, + size_t interval, + double altitudeDeviation, + double coefficientOfDrag) : interval(interval), altitudeDeviation(altitudeDeviation), distribution(0., altitudeDeviation) { // split the 3 comma seperated arguments auto seperator1 = constellation_data_str.find(',', 0); @@ -31,7 +34,8 @@ Constellation::Constellation(const std::string &constellation_data_str, size_t i // set variables using 3 args std::vector sats = readDatasetConstellation(std::string(DATADIR) + argConstellationName + "/pos_" + argConstellationName + ".csv", - std::string(DATADIR) + argConstellationName + "/v_" + argConstellationName + ".csv"); + std::string(DATADIR) + argConstellationName + "/v_" + argConstellationName + ".csv", + coefficientOfDrag); // convert vector to deque constellationSize = sats.size(); @@ -122,7 +126,8 @@ size_t Constellation::getConstellationSize() const { } std::vector Constellation::readDatasetConstellation(const std::string &position_filepath, - const std::string &velocity_filepath) { + const std::string &velocity_filepath, + double coefficientOfDrag) { CSVReader pos_csvReader{position_filepath, false}; CSVReader vel_csvReader{velocity_filepath, false}; std::vector particleCollection; @@ -138,18 +143,20 @@ std::vector Constellation::readDatasetConstellation(const std::string particleCollection.reserve(positions.size()); size_t particleId = 0; - std::transform(positions.begin(), - positions.end(), - velocities.begin(), - std::back_insert_iterator>(particleCollection), - [&](const auto &pos, const auto &vel) { - const auto &[x, y, z] = pos; - const auto &[vx, vy, vz] = vel; - - const std::array posArray = {x, y, z}; - const std::array velArray = {vx, vy, vz}; - return Particle(posArray, velArray, particleId++); - }); + std::transform( + positions.begin(), + positions.end(), + velocities.begin(), + std::back_insert_iterator>(particleCollection), + [&](const auto &pos, const auto &vel) { + const auto &[x, y, z] = pos; + const auto &[vx, vy, vz] = vel; + + const std::array posArray = {x, y, z}; + const std::array velArray = {vx, vy, vz}; + return Particle( + posArray, velArray, particleId++, Particle::ActivityState::evasivePreserving, 1., 1., coefficientOfDrag); + }); return particleCollection; } diff --git a/src/ladds/particle/Constellation.h b/src/ladds/particle/Constellation.h index 4f57eab5..6969eece 100644 --- a/src/ladds/particle/Constellation.h +++ b/src/ladds/particle/Constellation.h @@ -29,8 +29,12 @@ class Constellation { * passed for internal logic * @param altitudeDeviation : used to create satellites with normally distributed * altitudes. Equals the standard deviation of a normal distribution + * @param coefficientOfDrag c_D used to initialize all satellites. */ - Constellation(const std::string &constellation_data_str, size_t interval, double altitudeDeviation); + Constellation(const std::string &constellation_data_str, + size_t interval, + double altitudeDeviation, + double coefficientOfDrag); /** * determines which satellites are being added to the simulation by adding each shell * within a time span proportional to the shells size. shells are added plane by plane @@ -53,9 +57,14 @@ class Constellation { /** * Reads the passed position and velocity csv files. Returns a vector of particles. + * @param position_filepath + * @param velocity_filepath + * @param coefficientOfDrag + * @return */ static std::vector readDatasetConstellation(const std::string &position_filepath, - const std::string &velocity_filepath); + const std::string &velocity_filepath, + double coefficientOfDrag); /** * changes the pos vector by adding a random, normal distributed offset to the altitude diff --git a/src/ladds/particle/Particle.cpp b/src/ladds/particle/Particle.cpp index 3a1caae4..d7e8886a 100644 --- a/src/ladds/particle/Particle.cpp +++ b/src/ladds/particle/Particle.cpp @@ -72,7 +72,66 @@ void Particle::setBcInv(const double bcInv) { bc_inv = bcInv; } +Particle::ActivityState Particle::getActivityState() const { + return activityState; +} + +void Particle::setActivityState(Particle::ActivityState activityState) { + Particle::activityState = activityState; +} + std::ostream &operator<<(std::ostream &os, const Particle &particle) { - os << particle.toString(); + // clang-format off + os << particle.toString() + << "\nMass : " << particle.getMass() + << "\nRadius : " << particle.getRadius() + << "\nActivityState : " << static_cast(particle.getActivityState()); + // clang-format on return os; } + +std::istream &operator>>(std::istream &s, Particle::ActivityState &state) { + std::string str; + s >> str; + if (str == "passive") { + state = Particle::ActivityState::passive; + } else if (str == "evasive") { + state = Particle::ActivityState::evasive; + } else if (str == "evasivePreserving") { + state = Particle::ActivityState::evasivePreserving; + } + return s; +} + +bool Particle::operator==(const Particle &rhs) const { + // clang-format off + return static_cast &>(*this) + == static_cast &>(rhs) and + acc_t0 == rhs.acc_t0 and + acc_t1 == rhs.acc_t1 and + aom == rhs.aom and + mass == rhs.mass and + radius == rhs.radius and + bc_inv == rhs.bc_inv and + activityState == rhs.activityState; + // clang-format on +} +bool Particle::operator!=(const Particle &rhs) const { + return !(rhs == *this); +} + +double Particle::getMass() const { + return mass; +} + +void Particle::setMass(double mass) { + Particle::mass = mass; +} + +double Particle::getRadius() const { + return radius; +} + +void Particle::setRadius(double radius) { + Particle::radius = radius; +} diff --git a/src/ladds/particle/Particle.h b/src/ladds/particle/Particle.h index e2fde120..19366a8b 100644 --- a/src/ladds/particle/Particle.h +++ b/src/ladds/particle/Particle.h @@ -16,9 +16,51 @@ */ class Particle final : public autopas::ParticleFP64 { public: - explicit Particle(std::array pos, std::array v, size_t debrisId) - : autopas::ParticleFP64(pos, v, debrisId) {} + /** + * Describes how active a particle behaves. This is relevant for the propagator to determine which forces are applied. + */ + enum class ActivityState : int { + /** + * Simply float around space and is subject to all external influences. + */ + passive, + + /** + * Can actively change its trajectory to avoid collisions. + */ + evasive, + + /** + * Same as evasive but additionally can actively maneuver to preserve their orbit. + * This is modelled only only applying Keplerian forces. + */ + evasivePreserving, + }; + + /** + * Constructor + * @param pos + * @param v + * @param debrisId + * @param activityState + */ + Particle(std::array pos, + std::array v, + size_t debrisId, + ActivityState activityState, + double mass, + double radius, + double coefficientOfDrag) + : autopas::ParticleFP64(pos, v, debrisId), + aom(M_PI * radius * radius * 1e-6 / mass), // convert m^2 -> km^2 + mass(mass), + radius(radius), + bc_inv(coefficientOfDrag * aom * 1e6), // convert km^2 -> m^2 + activityState(activityState) {} + /** + * Destructor. + */ ~Particle() final = default; /** @@ -98,117 +140,153 @@ class Particle final : public autopas::ParticleFP64 { } /** - * @brief Calculates distance from the origin of the coordinate frame + * Calculates distance from the origin of the coordinate frame [km] * * @return Euclidean norm of the #position vector */ [[nodiscard]] double getHeight() const; /** - * @brief Calculates speed of the debris + * Calculates speed of the debris * * @return Euclidean norm of the #velocity vector */ [[nodiscard]] double getSpeed() const; /** - * @brief Calculates the euclidean norm of the #acc_t0 + * Calculates the euclidean norm of the #acc_t0 * * @return Calculates the euclidean norm of the #acc_t0 */ [[nodiscard]] double getAccT0Norm() const; /** - * @brief Calculates the euclidean norm of the #acc_t1 + * Calculates the euclidean norm of the #acc_t1 * * @return Calculates the euclidean norm of the #acc_t1 */ [[nodiscard]] double getAccT1Norm() const; /** - * @brief Getter function for #position vector + * Getter function for #position vector [km] * * @return 3D vector representation of the debris #position */ [[nodiscard]] const std::array &getPosition() const; /** - * @brief Setter function for #position vector + * Setter function for #position vector [km] * * @param position 3D vector representation of the debris #position */ void setPosition(const std::array &position); /** - * @brief Getter function for #velocity vector + * Getter function for #velocity vector [km/s] * * @return 3D vector representation of the debris #velocity */ [[nodiscard]] const std::array &getVelocity() const; /** - * @brief Setter function for #velocity vector + * Setter function for #velocity vector * * @param velocity 3D vector representation of the debris #velocity */ void setVelocity(const std::array &velocity); /** - * @brief Getter function for #acc_t0 vector + * Getter function for #acc_t0 vector * * @return 3D vector representation of the debris #acc_t0 */ [[nodiscard]] const std::array &getAccT0() const; /** - * @brief Setter function for #acc_t0 vector + * Setter function for #acc_t0 vector * * @param accT0 3D vector representation of the debris #acc_t0 */ void setAccT0(const std::array &accT0); /** - * @brief Getter function for #acc_t1 vector + * Getter function for #acc_t1 vector * * @return 3D vector representation of the debris #acc_t1 */ [[nodiscard]] const std::array &getAccT1() const; /** - * @brief Setter function for #acc_t1 vector + * Setter function for #acc_t1 vector * * @param accT1 3D vector representation of the debris #acc_t1 */ void setAccT1(const std::array &accT1); /** - * @brief Getter function for #aom + * Getter function for #aom * * @return value of #aom */ [[nodiscard]] double getAom() const; /** - * @brief Setter function for #aom + * Setter function for #aom * * @param aom New value #aom */ void setAom(double aom); /** - * @brief Getter function for #bc_inv + * Getter function for #bc_inv * * @return Value of #bc_inv */ [[nodiscard]] double getBcInv() const; /** - * @brief Setter function for #bc_inv + * Setter function for #bc_inv * * @param bcInv New value of#bc_inv */ void setBcInv(double bcInv); + /** + * Getter for the activityState. + * @return + */ + ActivityState getActivityState() const; + + /** + * Setter for the activityState. + * @param activityState + */ + void setActivityState(ActivityState activityState); + + /** + * Getter for mass. + * @return + */ + double getMass() const; + + /** + * Setter for mass. + * @param mass + */ + void setMass(double mass); + + /** + * Getter for radius. + * @return + */ + double getRadius() const; + + /** + * Setter for radius. + * @param radius + */ + void setRadius(double radius); + /** * Stream operator * @param os @@ -217,22 +295,59 @@ class Particle final : public autopas::ParticleFP64 { */ friend std::ostream &operator<<(std::ostream &os, const Particle &particle); + /** + * Equality operator + * @param rhs + * @return + */ + bool operator==(const Particle &rhs) const; + + /** + * Inequality operator + * @param rhs + * @return + */ + bool operator!=(const Particle &rhs) const; + private: /** - * 3D vector representation of the debris acceleration at the last time step. + * 3D vector representation of the debris acceleration at the last time step. [km/s^2] */ std::array acc_t0{}; /** - * 3D vector representation of the debris acceleration at the current time step + * 3D vector representation of the debris acceleration at the current time step. [km/s^2] */ std::array acc_t1{}; /** - * Area to mass ration. + * Area to mass relation [km^2/kg]. + * @note Used in the Propagator in SRPComponent::apply(). + */ + double aom; + /** + * Mass of the object [kg]. + */ + double mass; + /** + * Radius of the object when approximating it as a ball [m]. */ - double aom{0.}; + double radius; /** - * C_cA)/m is the inverse of the ballistic coefficient. - * Used for Acceleration::DragComponent::apply(). + * c_D*A/m is the inverse of the ballistic coefficient. [kg m^2] + * @note Used in the Propagator in Acceleration::DragComponent::apply(). */ - double bc_inv{0.}; + double bc_inv; + + /** + * If a particle can actively influence its orbit. See ActivityState. + */ + ActivityState activityState; }; + +/** + * Input stream operator for the activity state. + * @note Needed for the CSVReader. + * @param s + * @param state + * @return + */ +std::istream &operator>>(std::istream &s, Particle::ActivityState &state); \ No newline at end of file diff --git a/src/ladds/particle/SatelliteToParticleConverter.h b/src/ladds/particle/SatelliteToParticleConverter.h index 4e59f8e4..ddb712af 100644 --- a/src/ladds/particle/SatelliteToParticleConverter.h +++ b/src/ladds/particle/SatelliteToParticleConverter.h @@ -10,19 +10,29 @@ namespace SatelliteToParticleConverter { /** * This function is used to convert the satellite data into particles. + * @note Resulting particles always have Particle::ActivityState::passive as we consider results of a collision to be + * broken. If a non-catastrophic collision should be simulated the Particle::ActivityState has to be changed on the + * returned object. * @param satellite: Satellite to convert. * @retval Converted particles. */ -[[nodiscard]] inline Particle convertSatelliteToParticle(const Satellite &satellite) { +[[nodiscard]] inline Particle convertSatelliteToParticle(const Satellite &satellite, double coefficientOfDrag) { // Convert all entries from meters to kilometers const auto &position = autopas::utils::ArrayMath::mulScalar(satellite.getPosition(), 1. / 1000.0); const auto &velocity = autopas::utils::ArrayMath::mulScalar(satellite.getVelocity(), 1. / 1000.0); - - return Particle(position, velocity, satellite.getId()); + const auto radius = std::sqrt(satellite.getArea() * M_1_PI); + return Particle{position, + velocity, + satellite.getId(), + Particle::ActivityState::passive, + satellite.getMass(), + radius, + coefficientOfDrag}; } /** * This function is used to convert the particle data into satellites. + * @note This conversion does not retain the information about the satellites radius and activity state! * @param particles: Particles to convert. * @retval Converted satellites. */ @@ -34,7 +44,7 @@ namespace SatelliteToParticleConverter { // Converting from km to meters .setPosition(autopas::utils::ArrayMath::mulScalar(particle.getPosition(), 1000.0)) .setVelocity(autopas::utils::ArrayMath::mulScalar(particle.getVelocity(), 1000.0)) - .setMass(1); + .setMass(particle.getMass()); return satelliteBuilder.getResult(); } diff --git a/tests/testladds/CollisionFunctorIntegrationTest.cpp b/tests/testladds/CollisionFunctorIntegrationTest.cpp index 30309fde..d60a8796 100644 --- a/tests/testladds/CollisionFunctorIntegrationTest.cpp +++ b/tests/testladds/CollisionFunctorIntegrationTest.cpp @@ -28,7 +28,11 @@ void CollisionFunctorIntegrationTest::SetUpTestSuite() { 0., 0., }, - 0}; + 0, + Particle::ActivityState::passive, + 1., + 1., + 2.2}; _debris.reserve(numDebris); // fix seed for randomPosition() @@ -76,7 +80,7 @@ TEST_P(CollisionFunctorIntegrationTest, testAutoPasAlgorithm) { GTEST_SKIP_("SoAFunctor currently not implemented!"); } - CollisionFunctor functor(_cutoff, 10.0, 0.1 * _cutoff); + CollisionFunctor functor(_cutoff, 10.0, 0.1 * _cutoff, 0.1); if (not functor.allowsNonNewton3() and newton3 == autopas::Newton3Option::disabled) { GTEST_SKIP_("Functor does not support Newton3==disabled!"); diff --git a/tests/testladds/CollisionFunctorTest.cpp b/tests/testladds/CollisionFunctorTest.cpp index 74357a34..97c15644 100644 --- a/tests/testladds/CollisionFunctorTest.cpp +++ b/tests/testladds/CollisionFunctorTest.cpp @@ -27,10 +27,14 @@ TEST(CollisionFunctorTest, ThreeParticles) { for (size_t i = 0; i < numDebris; ++i) { debris.emplace_back(std::array{static_cast(i), 0., 0.}, std::array{0., static_cast(i), 0.}, - i); + i, + Particle::ActivityState::passive, + 1., + 1., + 2.2); } - CollisionFunctor collisionFunctor(cutoff, 10.0, cutoff); + CollisionFunctor collisionFunctor(cutoff, 10.0, cutoff, 0.01); for (size_t i = 0; i < debris.size(); ++i) { for (size_t j = i + 1; j < debris.size(); ++j) { @@ -51,6 +55,140 @@ TEST(CollisionFunctorTest, ThreeParticles) { EXPECT_THAT(collisionFunctor.getCollisions(), ::testing::UnorderedElementsAreArray(expected)); } +/** + * Set up two particles per ActivityState, all very close to each other. + * We expect exactly the following collisions: + * passiveSmall x passiveSmall + * passiveSmall x passiveLarge + * passiveSmall x evasive + * passiveSmall x evasivePreserving + * passiveLarge x passiveLarge + * Thus the following collisions to be disregarded: + * evasive x passiveLarge + * evasive x evasive + * evasive x evadingPreserving + * evasivePreserving x passiveLarge + * evasivePreserving x evadingPreserving + */ +TEST(CollisionFunctorTest, MixActivityStates) { + constexpr double cutoff{1.}; + constexpr bool newton3{false}; + constexpr size_t numDebris{8}; + + std::vector debris; + debris.reserve(numDebris); + + // passive small 1 + debris.emplace_back(std::array{0., 0., 0.}, + std::array{0., 0., 0.}, + 0, + Particle::ActivityState::passive, + 1., + 0.001, + 2.2); + // passive small 2 + debris.emplace_back(std::array{0., 0., 0.1}, + std::array{0., 0., 0.}, + 0, + Particle::ActivityState::passive, + 1., + 0.001, + 2.2); + // passive large 1 + debris.emplace_back(std::array{0.1, 0., 0.}, + std::array{0., 0., 0.}, + 1, + Particle::ActivityState::passive, + 1., + 1., + 2.2); + // passive large 2 + debris.emplace_back(std::array{0.1, 0., 0.1}, + std::array{0., 0., 0.}, + 1, + Particle::ActivityState::passive, + 1., + 1., + 2.2); + // evasive 1 + debris.emplace_back(std::array{0., 0.1, 0.}, + std::array{0., 0., 0.}, + 2, + Particle::ActivityState::evasive, + 1., + 1., + 2.2); + // evasive 2 + debris.emplace_back(std::array{0., 0.1, 0.1}, + std::array{0., 0., 0.}, + 3, + Particle::ActivityState::evasive, + 1., + 1., + 2.2); + // evasivePreserving 1 + debris.emplace_back(std::array{0.1, 0.1, 0.}, + std::array{0., 0., 0.}, + 4, + Particle::ActivityState::evasivePreserving, + 1., + 1., + 2.2); + // evasivePreserving 2 + debris.emplace_back(std::array{0.1, 0.1, 0.1}, + std::array{0., 0., 0.}, + 5, + Particle::ActivityState::evasivePreserving, + 1., + 1., + 2.2); + + CollisionFunctor collisionFunctor(cutoff, 10.0, cutoff, 0.01); + + for (size_t i = 0; i < debris.size(); ++i) { + for (size_t j = i + 1; j < debris.size(); ++j) { + collisionFunctor.AoSFunctor(debris[i], debris[j], newton3); + } + } + + // needed to merge the functor's internal thread buffers + collisionFunctor.endTraversal(newton3); + + auto collisions = collisionFunctor.getCollisions(); + + // extract ID pairs from collisions + std::vector> collisionIdPairs; + collisionIdPairs.reserve(collisions.size()); + std::transform(collisions.begin(), collisions.end(), std::back_inserter(collisionIdPairs), [](const auto &collision) { + const auto &[d1, d2, dist] = collision; + return std::make_tuple(d1->getID(), d2->getID()); + }); + + // set up expectations + decltype(collisionIdPairs) expectations{ + // passive small 1 collides with everything + {debris[0].getID(), debris[1].getID()}, + {debris[0].getID(), debris[2].getID()}, + {debris[0].getID(), debris[3].getID()}, + {debris[0].getID(), debris[4].getID()}, + {debris[0].getID(), debris[5].getID()}, + {debris[0].getID(), debris[6].getID()}, + {debris[0].getID(), debris[7].getID()}, + // passive small 2 collides with everything (collsion with passive small 1 already covered above) + {debris[1].getID(), debris[2].getID()}, + {debris[1].getID(), debris[3].getID()}, + {debris[1].getID(), debris[4].getID()}, + {debris[1].getID(), debris[5].getID()}, + {debris[1].getID(), debris[6].getID()}, + {debris[1].getID(), debris[7].getID()}, + // passive large 1 and 2 collide + {debris[2].getID(), debris[3].getID()}, + }; + + EXPECT_EQ(collisions.size(), expectations.size()); + EXPECT_THAT(collisionIdPairs, ::testing::UnorderedElementsAreArray(expectations)); +} + TEST_P(CollisionFunctorTest, LinearInterpolationTest) { constexpr double cutoff{80.0}; constexpr bool newton3{false}; @@ -58,14 +196,14 @@ TEST_P(CollisionFunctorTest, LinearInterpolationTest) { const auto &[x1, x2, v1, v2, dt, expected_dist] = GetParam(); - CollisionFunctor collisionFunctor(cutoff, dt, 0.1 * cutoff); + CollisionFunctor collisionFunctor(cutoff, dt, 0.1 * cutoff, 0.1); std::vector debris; debris.reserve(numDebris); // Add two particles moving in the same direction on parallel lines - debris.emplace_back(x1, v1, 0.); - debris.emplace_back(x2, v2, 1.); + debris.emplace_back(x1, v1, 0., Particle::ActivityState::passive, 1., 1., 2.2); + debris.emplace_back(x2, v2, 1., Particle::ActivityState::passive, 1., 1., 2.2); for (size_t i = 0; i < debris.size(); ++i) { for (size_t j = i + 1; j < debris.size(); ++j) { diff --git a/tests/testladds/SimulationTest.cpp b/tests/testladds/SimulationTest.cpp index ed31c273..36422989 100644 --- a/tests/testladds/SimulationTest.cpp +++ b/tests/testladds/SimulationTest.cpp @@ -37,17 +37,28 @@ TEST_F(SimulationTest, testInsertionOverlap) { */ const std::vector expectedDelayedParticles{2, 1, 1, 0, 0, 0}; - const std::vector initialSatellites{Particle({6871., 0., 0.}, {0, 4.8958309146899, 5.83462408131549}, 1), - Particle({6870.99577848984, 4.89582991243564, 5.83462288687537}, - {-0.00844301944979238, 4.89582790671436, 5.83462049654983}, - 2)}; + const std::vector initialSatellites{ + Particle( + {6871., 0., 0.}, {0, 4.8958309146899, 5.83462408131549}, 1, Particle::ActivityState::passive, 1., 1., 2.2), + Particle({6870.99577848984, 4.89582991243564, 5.83462288687537}, + {-0.00844301944979238, 4.89582790671436, 5.83462049654983}, + 2, + Particle::ActivityState::passive, + 1., + 1., + 2.2)}; for (const auto &s : initialSatellites) { autopas->addParticle(s); } - const std::vector newSatellites{Particle({6871., 0., 0.}, {0, 4.8958309146899, 5.83462408131549}, 1), - Particle({6870.99577848984, 4.89582991243564, 5.83462288687537}, - {-0.00844301944979238, 4.89582790671436, 5.83462049654983}, - 4)}; + const std::vector newSatellites{ + Particle({6871., 0., 0.}, {0, 4.8958309146899, 5.83462408131549}, 1, Particle::ActivityState::passive, 0, 0, 2.2), + Particle({6870.99577848984, 4.89582991243564, 5.83462288687537}, + {-0.00844301944979238, 4.89582790671436, 5.83462049654983}, + 4, + Particle::ActivityState::passive, + 1., + 1., + 2.2)}; auto constellationCutoff = config["io"]["constellationCutoff"].as(); @@ -74,10 +85,12 @@ TEST_F(SimulationTest, testBurnUp) { auto [csvWriter, accumulator, integrator] = simulation.initIntegrator(*autopas, *configReader); const auto &minAltitude = Physics::R_EARTH + configReader->get("sim/minAltitude"); // initialize a particle 1km above burn up radius with a trajectory towards earth - autopas->addParticle(Particle({minAltitude + 1., 0., 0.}, {-10., 0., 0.}, 0)); + autopas->addParticle( + Particle({minAltitude + 1., 0., 0.}, {-10., 0., 0.}, 0, Particle::ActivityState::passive, 1., 1., 2.2)); // initialize a particle 1km above burn up radius with a trajectory away from earth // different position to avoid any interferences - autopas->addParticle(Particle({0., minAltitude + 1., 0.}, {0., 10., 0.}, 1)); + autopas->addParticle( + Particle({0., minAltitude + 1., 0.}, {0., 10., 0.}, 1, Particle::ActivityState::passive, 1., 1., 2.2)); ASSERT_EQ(autopas->getNumberOfParticles(), 2) << "Initial particles not found!"; // dummy because interface requires it @@ -98,10 +111,11 @@ TEST_F(SimulationTest, testBurnUp) { TEST_P(SimulationTest, testCheckedInsert) { const auto &[posTestParticle, positionIsSafe] = GetParam(); - autopas->addParticle(Particle(testCheckedInsertParticlePos, zeroVec, 0)); + autopas->addParticle( + Particle(testCheckedInsertParticlePos, zeroVec, 0, Particle::ActivityState::passive, 1., 1., 2.2)); // particle that will be inserted - Particle p1{posTestParticle, zeroVec, 1}; + Particle p1{posTestParticle, zeroVec, 1, Particle::ActivityState::passive, 1., 1., 2.2}; const auto escapedParticles = autopas->updateContainer(); ASSERT_TRUE(escapedParticles.empty()) << "Test setup faulty!"; diff --git a/tests/testladds/io/hdf5/HDF5WriterReaderTest.cpp b/tests/testladds/io/hdf5/HDF5WriterReaderTest.cpp index db1d8870..da3392d4 100644 --- a/tests/testladds/io/hdf5/HDF5WriterReaderTest.cpp +++ b/tests/testladds/io/hdf5/HDF5WriterReaderTest.cpp @@ -24,7 +24,13 @@ TEST_F(HDF5WriterReaderTest, WriteReadTestParticleData) { std::vector particles; particles.reserve(numParticles); for (size_t i = 0; i < numParticles; ++i) { - Particle p{{static_cast(i), static_cast(i), static_cast(i)}, {0., 0., 0.}, i}; + Particle p{{static_cast(i), static_cast(i), static_cast(i)}, + {1., 2., 3.}, + i, + Particle::ActivityState::evasive, + 1., + 1., + 2.2}; autopas.addParticle(p); particles.push_back(p); } @@ -39,12 +45,29 @@ TEST_F(HDF5WriterReaderTest, WriteReadTestParticleData) { HDF5Writer hdf5Writer(filename, 4); hdf5Writer.writeParticles(iterationNr, autopas); - // 3. read data + // 3. read data and check that read data is equal to generated data HDF5Reader hdf5Reader(filename); - auto particlesHDF5 = hdf5Reader.readParticles(iterationNr); + EXPECT_THAT(hdf5Reader.readParticles(iterationNr, 2.2), ::testing::UnorderedElementsAreArray(particles)) + << "Particle data of initial iteration is not correct!"; + + // 4. add more particles and check that they are added correctly to HDF5 + particles.emplace_back(std::array{0.5, 0.5, 0.5}, + std::array{1., 2., 3.}, + numParticles, + Particle::ActivityState::evasivePreserving, + 1., + 1., + 2.2); + autopas.addParticle(particles.back()); + + // 5. write data + hdf5Writer.writeParticles(iterationNr + 1, autopas); + + // 6. read data and check that read data is equal to generated data + EXPECT_THAT(hdf5Reader.readParticles(iterationNr + 1, 2.2), ::testing::UnorderedElementsAreArray(particles)) + << "Particle data of second iteration is not correct!"; + ; - // 4. check that read data is equal to generated data - EXPECT_THAT(particlesHDF5, ::testing::UnorderedElementsAreArray(particles)); // cleanup std::remove(filename); } @@ -55,7 +78,8 @@ TEST_F(HDF5WriterReaderTest, WriteReadTestCollisionData) { std::vector particles; particles.reserve(numParticles); for (size_t i = 0; i < numParticles; ++i) { - particles.push_back(Particle{{0., 0., static_cast(i)}, {0., 0., 0.}, i}); + particles.emplace_back( + {{0., 0., static_cast(i)}, {0., 0., 0.}, i, Particle::ActivityState::passive, 1., 1., 2.2}); } // These conjunctions are just randomly made up and have nothing to do with position data! @@ -82,10 +106,10 @@ TEST_F(HDF5WriterReaderTest, WriteReadTestCollisionData) { // 4. check that read data is equal to generated data EXPECT_EQ(conjunctions.size(), conjunctionsHDF5.size()); for (const auto &[ptrARef, ptrBRef, distRef] : conjunctions) { - const auto idARef = static_cast(ptrARef->getID()); - const auto idBRef = static_cast(ptrBRef->getID()); + const auto idARef = static_cast(ptrARef->getID()); + const auto idBRef = static_cast(ptrBRef->getID()); - HDF5Writer::CollisionInfo collisionInfo{idARef, idBRef, static_cast(distRef)}; + HDF5Definitions::CollisionInfo collisionInfo{idARef, idBRef, static_cast(distRef)}; EXPECT_THAT(conjunctionsHDF5, ::testing::Contains(collisionInfo)); } // cleanup diff --git a/tests/testladds/particle/BreakupWrapperTest.cpp b/tests/testladds/particle/BreakupWrapperTest.cpp index fc15a40b..4213c65b 100644 --- a/tests/testladds/particle/BreakupWrapperTest.cpp +++ b/tests/testladds/particle/BreakupWrapperTest.cpp @@ -12,12 +12,20 @@ * Crash two particles into each other and observe that new particles have higher IDs. */ TEST_F(BreakupWrapperTest, testSimulationLoop) { + configReader->setValue("sim/coefficientOfDrag", 2.2); auto [csvWriter, accumulator, integrator] = simulation.initIntegrator(*autopas, *configReader); // two particles 1000km above earth whose paths cross exactly at [R+1000, 0, 0] size_t highestIdBeforeCrash = 4; - autopas->addParticle(Particle({Physics::R_EARTH + 1000., -1., 0.}, {0., 2., 0.}, 1)); - autopas->addParticle(Particle({Physics::R_EARTH + 1000., 0., -1.}, {0., 0., 2.}, highestIdBeforeCrash)); + autopas->addParticle( + Particle({Physics::R_EARTH + 1000., -1., 0.}, {0., 2., 0.}, 1, Particle::ActivityState::passive, 1., 1., 2.2)); + autopas->addParticle(Particle({Physics::R_EARTH + 1000., 0., -1.}, + {0., 0., 2.}, + highestIdBeforeCrash, + Particle::ActivityState::passive, + 1., + 1., + 2.2)); // dummy std::vector constellations; diff --git a/tests/testladds/particle/SatelliteToParticleConverterTest.cpp b/tests/testladds/particle/SatelliteToParticleConverterTest.cpp index 4831210c..2f48dcfb 100644 --- a/tests/testladds/particle/SatelliteToParticleConverterTest.cpp +++ b/tests/testladds/particle/SatelliteToParticleConverterTest.cpp @@ -11,22 +11,24 @@ TEST(SatelliteToParticleConverterTest, SatelliteToParticle) { // Create a particle. const Satellite sat{"name", SatType::DEBRIS, {10., 10., 10.}}; - const auto particle = SatelliteToParticleConverter::convertSatelliteToParticle(sat); + const auto particle = SatelliteToParticleConverter::convertSatelliteToParticle(sat, 2.2); const auto &expectedPosition = autopas::utils::ArrayMath::mulScalar(sat.getPosition(), 1. / 1000.0); const auto &expectedVelocity = autopas::utils::ArrayMath::mulScalar(sat.getVelocity(), 1. / 1000.0); EXPECT_THAT(particle.getPosition(), ::testing::ElementsAreArray(expectedPosition)); EXPECT_THAT(particle.getVelocity(), ::testing::ElementsAreArray(expectedVelocity)); + EXPECT_EQ(particle.getMass(), sat.getMass()); const auto satConverted = SatelliteToParticleConverter::convertParticleToSatellite(particle); EXPECT_THAT(satConverted.getPosition(), ::testing::ElementsAreArray(sat.getPosition())); EXPECT_THAT(satConverted.getVelocity(), ::testing::ElementsAreArray(sat.getVelocity())); + EXPECT_EQ(satConverted.getMass(), sat.getMass()); } TEST(SatelliteToParticleConverterTest, ParticleToSatellite) { // Create a particle. - const Particle particle{{10., 11., 12}, {1., 2., 3.}, 42}; + const Particle particle{{10., 11., 12}, {1., 2., 3.}, 42, Particle::ActivityState::passive, 1., 1., 2.2}; const auto sat = SatelliteToParticleConverter::convertParticleToSatellite(particle); const auto &expectedPosition = autopas::utils::ArrayMath::mulScalar(particle.getPosition(), 1000.0); @@ -34,8 +36,10 @@ TEST(SatelliteToParticleConverterTest, ParticleToSatellite) { EXPECT_THAT(sat.getPosition(), ::testing::ElementsAreArray(expectedPosition)); EXPECT_THAT(sat.getVelocity(), ::testing::ElementsAreArray(expectedVelocity)); + EXPECT_EQ(sat.getMass(), particle.getMass()); - const auto particleConverted = SatelliteToParticleConverter::convertSatelliteToParticle(sat); + const auto particleConverted = SatelliteToParticleConverter::convertSatelliteToParticle(sat, 2.2); EXPECT_THAT(particleConverted.getPosition(), ::testing::ElementsAreArray(particle.getPosition())); EXPECT_THAT(particleConverted.getVelocity(), ::testing::ElementsAreArray(particle.getVelocity())); + EXPECT_EQ(particleConverted.getMass(), particle.getMass()); }