diff --git a/src/Domain/Creators/CMakeLists.txt b/src/Domain/Creators/CMakeLists.txt index ee60a1962b051..428692aa1261f 100644 --- a/src/Domain/Creators/CMakeLists.txt +++ b/src/Domain/Creators/CMakeLists.txt @@ -59,6 +59,7 @@ target_link_libraries( ${LIBRARY} PRIVATE DomainStructure + H5 PUBLIC CoordinateMaps DataStructures diff --git a/src/Domain/Creators/TimeDependentOptions/CMakeLists.txt b/src/Domain/Creators/TimeDependentOptions/CMakeLists.txt index 462672ff1414e..1c9d0d4bbaf71 100644 --- a/src/Domain/Creators/TimeDependentOptions/CMakeLists.txt +++ b/src/Domain/Creators/TimeDependentOptions/CMakeLists.txt @@ -5,6 +5,7 @@ spectre_target_sources( ${LIBRARY} PRIVATE BinaryCompactObject.cpp + FromVolumeFile.cpp ShapeMap.cpp Sphere.cpp ) @@ -14,6 +15,7 @@ spectre_target_headers( INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/src HEADERS BinaryCompactObject.hpp + FromVolumeFile.hpp ShapeMap.hpp Sphere.hpp ) diff --git a/src/Domain/Creators/TimeDependentOptions/FromVolumeFile.cpp b/src/Domain/Creators/TimeDependentOptions/FromVolumeFile.cpp new file mode 100644 index 0000000000000..c4e3c5afbcfa5 --- /dev/null +++ b/src/Domain/Creators/TimeDependentOptions/FromVolumeFile.cpp @@ -0,0 +1,189 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Domain/Creators/TimeDependentOptions/FromVolumeFile.hpp" + +#include +#include +#include +#include +#include +#include +#include + +#include "DataStructures/DataVector.hpp" +#include "Domain/FunctionsOfTime/FunctionOfTime.hpp" +#include "Domain/FunctionsOfTime/QuaternionFunctionOfTime.hpp" +#include "Domain/FunctionsOfTime/QuaternionHelpers.hpp" +#include "IO/H5/AccessType.hpp" +#include "IO/H5/File.hpp" +#include "IO/H5/VolumeData.hpp" +#include "Options/Context.hpp" +#include "Options/ParseError.hpp" +#include "Utilities/Algorithm.hpp" +#include "Utilities/PrettyType.hpp" +#include "Utilities/Serialization/Serialize.hpp" + +namespace domain::creators::time_dependent_options { + +namespace { +// Returns a clone of the FoT requested +std::unique_ptr get_function_of_time( + const std::string& function_of_time_name, const std::string& h5_filename, + const std::string& subfile_name, const double time, + const Options::Context& context) { + const h5::H5File h5_file{h5_filename}; + const auto& vol_file = h5_file.get(subfile_name); + + const std::vector obs_ids = vol_file.list_observation_ids(); + if (obs_ids.empty()) { + PARSE_ERROR(context, function_of_time_name + << ": There are no observation IDs in the subfile " + << subfile_name << " of H5 file " << h5_filename); + } + // Take last observation ID so we have all possible times available + std::optional> serialized_functions_of_time = + vol_file.get_functions_of_time(obs_ids[obs_ids.size() - 1]); + + if (not serialized_functions_of_time.has_value()) { + PARSE_ERROR(context, + function_of_time_name + << ": There are no functions of time in the subfile " + << subfile_name << " of the H5 file " << h5_filename + << ". Choose a different subfile or H5 file."); + } + + const auto functions_of_time = deserialize>>( + serialized_functions_of_time->data()); + + if (not functions_of_time.contains(function_of_time_name)) { + PARSE_ERROR(context, "No function of time named " + << function_of_time_name << " in the subfile " + << subfile_name << " of the H5 file " + << h5_filename); + } + + const auto& function_of_time = functions_of_time.at(function_of_time_name); + const std::array time_bounds = function_of_time->time_bounds(); + + if (time < time_bounds[0] or time > time_bounds[1]) { + PARSE_ERROR(context, function_of_time_name + << ": The requested time " << time + << " is out of the range of the function of time " + << time_bounds << " from the subfile " + << subfile_name << " of the H5 file " + << h5_filename); + } + + return function_of_time->get_clone(); +} +} // namespace + +template +FromVolumeFile::FromVolumeFile(const std::string& h5_filename, + const std::string& subfile_name, + const double time, + const Options::Context& context) { + const std::string function_of_time_name = pretty_type::name(); + + const auto function_of_time = get_function_of_time( + function_of_time_name, h5_filename, subfile_name, time, context); + + values = function_of_time->func_and_2_derivs(time); +} + +FromVolumeFile::FromVolumeFile( + const std::string& h5_filename, const std::string& subfile_name, + const double time, const Options::Context& context) { + const std::string expansion_function_of_time_name{"Expansion"}; + + const auto exp_function_of_time = + get_function_of_time(expansion_function_of_time_name, h5_filename, + subfile_name, time, context); + const auto exp_outer_boundary_function_of_time = + get_function_of_time(expansion_function_of_time_name + "OuterBoundary", + h5_filename, subfile_name, time, context); + + expansion_values = exp_function_of_time->func_and_2_derivs(time); + expansion_values_outer_boundary = + exp_outer_boundary_function_of_time->func_and_2_derivs(time); +} + +FromVolumeFile::FromVolumeFile( + const std::string& h5_filename, const std::string& subfile_name, + const double time, const bool is_quat_fot, + const Options::Context& context) { + const std::string function_of_time_name{"Rotation"}; + + const auto function_of_time = get_function_of_time( + function_of_time_name, h5_filename, subfile_name, time, context); + + if (is_quat_fot) { + bool filled_values = false; + const auto set_values = [this, &function_of_time, &time, + &filled_values]() { + const auto* const quat_function_of_time = + dynamic_cast*>( + function_of_time.get()); + + if (quat_function_of_time != nullptr) { + quaternions = quat_function_of_time->quat_func_and_2_derivs(time); + angle_values = quat_function_of_time->angle_func_and_2_derivs(time); + filled_values = true; + } + }; + + set_values.operator()<2>(); + set_values.operator()<3>(); + + if (not filled_values) { + PARSE_ERROR(context, + function_of_time_name + << ": The max deriv of the quaternion function of " + "time must be either 2 or 3."); + } + } else { + quaternions = function_of_time->func_and_2_derivs(time); + std::array, 3> quats; + for (size_t i = 0; i < quaternions.size(); i++) { + gsl::at(quats, i) = datavector_to_quaternion(gsl::at(quaternions, i)); + } + + // We can't compute the angle so set it to zero + angle_values[0] = DataVector{4, 0.0}; + angle_values[1] = quaternion_to_datavector(2.0 * conj(quats[0]) * quats[1]); + angle_values[2] = quaternion_to_datavector( + conj(quats[0]) * (2.0 * quats[2] - quats[1] * datavector_to_quaternion( + angle_values[0]))); + // We'd need more quaternion derivatives to compute the higher angle + // derivatives. Also the above DataVectors have size 4, so we need to make + // them have size 3. + angle_values[2] = DataVector{4, 0.0}; + for (size_t i = 0; i < 3; i++) { + gsl::at(angle_values, i) = + DataVector{gsl::at(angle_values, i)[1], gsl::at(angle_values, i)[2], + gsl::at(angle_values, i)[3]}; + } + + // angle_values = std::nullopt; + } +} + +FromVolumeFile::FromVolumeFile( + const std::string& h5_filename, const std::string& subfile_name, + const double time, const Options::Context& context) { + const std::string shape_function_of_time_name{"Shape"}; + const std::string size_function_of_time_name{"Size"}; + + const auto shape_function_of_time = get_function_of_time( + shape_function_of_time_name, h5_filename, subfile_name, time, context); + const auto size_function_of_time = get_function_of_time( + size_function_of_time_name, h5_filename, subfile_name, time, context); + + shape_values = shape_function_of_time->func_and_2_derivs(time); + size_values = size_function_of_time->func_and_2_derivs(time); +} + +template struct FromVolumeFile; +} // namespace domain::creators::time_dependent_options diff --git a/src/Domain/Creators/TimeDependentOptions/FromVolumeFile.hpp b/src/Domain/Creators/TimeDependentOptions/FromVolumeFile.hpp new file mode 100644 index 0000000000000..5abff350910fc --- /dev/null +++ b/src/Domain/Creators/TimeDependentOptions/FromVolumeFile.hpp @@ -0,0 +1,131 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include +#include +#include + +#include "DataStructures/DataVector.hpp" +#include "Options/Auto.hpp" +#include "Options/Context.hpp" +#include "Options/String.hpp" +#include "Utilities/TMPL.hpp" + +namespace domain::creators::time_dependent_options { +/*! + * \brief Structs meant to be used as template parameters for the + * `domain::creators::time_dependent_options::FromVolumeFile` classes. + */ +namespace names { +struct Translation {}; +struct Rotation {}; +struct Expansion {}; +struct ShapeSize {}; +} // namespace names + +namespace detail { +struct FromVolumeFileBase { + struct H5Filename { + using type = std::string; + static constexpr Options::String help{ + "Name of H5 file to read functions of time from."}; + }; + + struct SubfileName { + using type = std::string; + static constexpr Options::String help{ + "Subfile that holds the volume data. Must be an h5::VolumeData " + "subfile."}; + }; + + struct Time { + using type = double; + static constexpr Options::String help = + "Time in the H5File to get the coefficients at. Will likely be the " + "same as the initial time"; + }; + + using options = tmpl::list; + static constexpr Options::String help = + "Read function of time coefficients from a volume subfile of an H5 file."; + + FromVolumeFileBase() = default; +}; +} // namespace detail + +/// @{ +/*! + * \brief Read in FunctionOfTime coefficients from an H5 file and volume + * subfile. + * + * \details To use, template the class on one of the structs in + * `domain::creators::time_dependent_options::names`. The general struct will + * have one member, `values` that will hold the function of time and its first + * two derivatives. + * + * There are specializations for + * + * - `domain::creators::time_dependent_options::names::Rotation` because of + * quaternions + * - `domain::creators::time_dependent_options::name::Expansion` because it also + * has outer boundary values (a second function of time) + * - `domain::creators::time_dependent_options::names::ShapeSize` because it + * handles both the Shape and Size function of time. + */ +template +struct FromVolumeFile : public detail::FromVolumeFileBase { + FromVolumeFile() = default; + FromVolumeFile(const std::string& h5_filename, + const std::string& subfile_name, double time, + const Options::Context& context = {}); + + std::array values{}; +}; + +template <> +struct FromVolumeFile : public detail::FromVolumeFileBase { + FromVolumeFile() = default; + FromVolumeFile(const std::string& h5_filename, + const std::string& subfile_name, double time, + const Options::Context& context = {}); + + std::array expansion_values{}; + std::array expansion_values_outer_boundary{}; +}; + +template <> +struct FromVolumeFile : public detail::FromVolumeFileBase { + struct IsQuaternionFunctionOfTime { + using type = bool; + static constexpr Options::String help = + "Whether or not the function of time stored in the volume file is a " + "QuaternionFunctionOfTime or not. If the volume data came from an " + "inspiral, this should be set to 'True'."; + }; + + using options = tmpl::push_back; + + FromVolumeFile() = default; + FromVolumeFile(const std::string& h5_filename, + const std::string& subfile_name, double time, bool is_quat_fot, + const Options::Context& context = {}); + + std::array quaternions{}; + std::array angle_values{}; +}; + +template <> +struct FromVolumeFile : public detail::FromVolumeFileBase { + FromVolumeFile() = default; + FromVolumeFile(const std::string& h5_filename, + const std::string& subfile_name, double time, + const Options::Context& context = {}); + + std::array shape_values{}; + std::array size_values{}; +}; +/// @} +} // namespace domain::creators::time_dependent_options diff --git a/tests/Unit/Domain/Creators/TimeDependentOptions/CMakeLists.txt b/tests/Unit/Domain/Creators/TimeDependentOptions/CMakeLists.txt index b7814a48c2657..97123f3d6e079 100644 --- a/tests/Unit/Domain/Creators/TimeDependentOptions/CMakeLists.txt +++ b/tests/Unit/Domain/Creators/TimeDependentOptions/CMakeLists.txt @@ -4,6 +4,7 @@ set(LIBRARY_SOURCES ${LIBRARY_SOURCES} TimeDependentOptions/Test_BinaryCompactObject.cpp + TimeDependentOptions/Test_FromVolumeFile.cpp TimeDependentOptions/Test_ShapeMap.cpp TimeDependentOptions/Test_Sphere.cpp PARENT_SCOPE) diff --git a/tests/Unit/Domain/Creators/TimeDependentOptions/Test_FromVolumeFile.cpp b/tests/Unit/Domain/Creators/TimeDependentOptions/Test_FromVolumeFile.cpp new file mode 100644 index 0000000000000..d12e158508643 --- /dev/null +++ b/tests/Unit/Domain/Creators/TimeDependentOptions/Test_FromVolumeFile.cpp @@ -0,0 +1,290 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "DataStructures/DataVector.hpp" +#include "Framework/TestingFramework.hpp" + +#include +#include +#include +#include +#include + +#include "Domain/Creators/TimeDependentOptions/FromVolumeFile.hpp" +#include "Domain/FunctionsOfTime/FixedSpeedCubic.hpp" +#include "Domain/FunctionsOfTime/FunctionOfTime.hpp" +#include "Domain/FunctionsOfTime/PiecewisePolynomial.hpp" +#include "Domain/FunctionsOfTime/QuaternionFunctionOfTime.hpp" +#include "Domain/FunctionsOfTime/RegisterDerivedWithCharm.hpp" +#include "Framework/TestCreation.hpp" +#include "IO/H5/AccessType.hpp" +#include "IO/H5/File.hpp" +#include "IO/H5/TensorData.hpp" +#include "IO/H5/VolumeData.hpp" +#include "NumericalAlgorithms/Spectral/Basis.hpp" +#include "NumericalAlgorithms/Spectral/Quadrature.hpp" +#include "Utilities/ConstantExpressions.hpp" +#include "Utilities/FileSystem.hpp" +#include "Utilities/PrettyType.hpp" +#include "Utilities/Serialization/Serialize.hpp" + +namespace { +void write_volume_data( + const std::string& filename, const std::string& subfile_name, + const std::unordered_map< + std::string, std::unique_ptr>& + functions_of_time = {}) { + h5::H5File h5_file{filename, true}; + auto& vol_file = h5_file.insert(subfile_name); + + // We don't care about the volume data here, just the functions of time + vol_file.write_volume_data( + 0, 0.0, + {ElementVolumeData{"blah", + {TensorComponent{"RandomTensor", DataVector{3, 0.0}}}, + {3}, + {Spectral::Basis::Legendre}, + {Spectral::Quadrature::GaussLobatto}}}, + std::nullopt, + functions_of_time.empty() ? std::nullopt + : std::optional{serialize(functions_of_time)}); +} + +template +void test() { + const std::string filename{"HorseRadish.h5"}; + if (file_system::check_if_file_exists(filename)) { + file_system::rm(filename, true); + } + const std::string subfile_name{"VolumeData"}; + const std::string function_of_time_name = pretty_type::name(); + + std::unordered_map> + functions_of_time{}; + functions_of_time[function_of_time_name] = + std::make_unique>( + 0.0, + std::array{DataVector{3, 0.0}, DataVector{3, 1.0}, + DataVector{3, 0.0}}, + 100.0); + + write_volume_data(filename, subfile_name, functions_of_time); + + const double time = 50.0; + + const auto from_volume_file = TestHelpers::test_creation< + domain::creators::time_dependent_options::FromVolumeFile>( + "H5Filename: " + filename + "\nSubfileName: " + subfile_name + + "\nTime: 50.0\n"); + + std::array expected_values{ + DataVector{3, 1.0 * time}, DataVector{3, 1.0}, DataVector{3, 0.0}}; + + CHECK(from_volume_file.values == expected_values); + + if (file_system::check_if_file_exists(filename)) { + file_system::rm(filename, true); + } +} + +template <> +void test() { + const std::string filename{"SpicyHorseRadish.h5"}; + if (file_system::check_if_file_exists(filename)) { + file_system::rm(filename, true); + } + const std::string subfile_name{"VolumeData"}; + const std::string function_of_time_name{"Expansion"}; + + std::unordered_map> + functions_of_time{}; + functions_of_time[function_of_time_name] = + std::make_unique>( + 0.0, + std::array{DataVector{1, 0.0}, DataVector{1, 1.0}, + DataVector{1, 0.0}}, + 100.0); + const double velocity = -1e-5; + const double decay_timescale = 50.0; + functions_of_time[function_of_time_name + "OuterBoundary"] = + std::make_unique( + 1.0, 0.0, velocity, decay_timescale); + + write_volume_data(filename, subfile_name, functions_of_time); + + // Makes things easier + const double time = decay_timescale; + + const auto from_volume_file = TestHelpers::test_creation< + domain::creators::time_dependent_options::FromVolumeFile< + domain::creators::time_dependent_options::names::Expansion>>( + "H5Filename: " + filename + "\nSubfileName: " + subfile_name + + "\nTime: 50.0\n"); + + std::array expected_values{DataVector{1.0 * time}, + DataVector{1.0}, DataVector{0.0}}; + // Comes from the FixedSpeedCubic formula + std::array expected_values_outer_boundary{ + DataVector{1.0 + velocity * cube(time) / + (square(decay_timescale) + square(time))}, + DataVector{velocity}, DataVector{0.01 * velocity}}; + + CHECK(from_volume_file.expansion_values == expected_values); + CHECK_ITERABLE_APPROX(from_volume_file.expansion_values_outer_boundary, + expected_values_outer_boundary); + + if (file_system::check_if_file_exists(filename)) { + file_system::rm(filename, true); + } +} + +template <> +void test() { + const std::string filename{"ExtraSpicyHorseRadish.h5"}; + if (file_system::check_if_file_exists(filename)) { + file_system::rm(filename, true); + } + const std::string subfile_name{"VolumeData"}; + const std::string function_of_time_name{"Rotation"}; + + std::unordered_map> + functions_of_time{}; + functions_of_time[function_of_time_name] = + std::make_unique>( + 0.0, std::array{DataVector{1.0, 0.0, 0.0, 0.0}}, + std::array{DataVector{3, 2.0}, DataVector{3, 1.0}, DataVector{3, 0.0}, + DataVector{3, 0.0}}, + 100.0); + + write_volume_data(filename, subfile_name, functions_of_time); + + { + INFO("IsQuaternionFunctionOfTime = True"); + // Going at t=0 is easier for checking quaternions + const auto from_volume_file = TestHelpers::test_creation< + domain::creators::time_dependent_options::FromVolumeFile< + domain::creators::time_dependent_options::names::Rotation>>( + "H5Filename: " + filename + "\nSubfileName: " + subfile_name + + "\nTime: 0.0\nIsQuaternionFunctionOfTime: True\n"); + + // q + // dtq = 0.5 * q * omega + // d2tq = 0.5 * (dtq * omega + q * dtomega) + std::array expected_quaternion{ + DataVector{1.0, 0.0, 0.0, 0.0}, DataVector{0.0, 0.5, 0.5, 0.5}, + DataVector{-0.75, 0.0, 0.0, 0.0}}; + std::array expected_angle{ + DataVector{3, 2.0}, DataVector{3, 1.0}, DataVector{3, 0.0}}; + + CHECK(from_volume_file.quaternions == expected_quaternion); + CHECK(from_volume_file.angle_values == expected_angle); + } + + { + INFO("IsQuaternionFunctionOfTime = False"); + // Going at t=0 is easier for checking quaternions + const auto from_volume_file = TestHelpers::test_creation< + domain::creators::time_dependent_options::FromVolumeFile< + domain::creators::time_dependent_options::names::Rotation>>( + "H5Filename: " + filename + "\nSubfileName: " + subfile_name + + "\nTime: 0.0\nIsQuaternionFunctionOfTime: False\n"); + + // q + // dtq = 0.5 * q * omega + // d2tq = 0.5 * (dtq * omega + q * dtomega) + std::array expected_quaternion{ + DataVector{1.0, 0.0, 0.0, 0.0}, DataVector{0.0, 0.5, 0.5, 0.5}, + DataVector{-0.75, 0.0, 0.0, 0.0}}; + std::array expected_angle{ + DataVector{3, 0.0}, DataVector{3, 1.0}, DataVector{3, 0.0}}; + + CHECK(from_volume_file.quaternions == expected_quaternion); + CHECK(from_volume_file.angle_values == expected_angle); + } + + if (file_system::check_if_file_exists(filename)) { + file_system::rm(filename, true); + } +} + +void test_errors() { + const std::string filename{"HorseRadishErrors.h5"}; + if (file_system::check_if_file_exists(filename)) { + file_system::rm(filename, true); + } + std::string subfile_name{"VolumeData"}; + + { + h5::H5File h5_file{filename}; + h5_file.insert(subfile_name); + } + + using FromVolumeFile = + domain::creators::time_dependent_options::FromVolumeFile< + domain::creators::time_dependent_options::names::Expansion>; + + CHECK_THROWS_WITH( + (FromVolumeFile{filename, subfile_name, 0.0}), + Catch::Matchers::ContainsSubstring( + "Expansion: There are no observation IDs in the subfile ")); + + // Need new subfile to write to + subfile_name += "0"; + write_volume_data(filename, subfile_name); + + CHECK_THROWS_WITH( + (FromVolumeFile{filename, subfile_name, 0.0}), + Catch::Matchers::ContainsSubstring( + "Expansion: There are no functions of time in the subfile ")); + + subfile_name += "0"; + std::unordered_map> + functions_of_time{}; + functions_of_time["Translation"] = + std::make_unique>( + 0.0, + std::array{DataVector{3, 0.0}, DataVector{3, 1.0}, + DataVector{3, 0.0}}, + 100.0); + + write_volume_data(filename, subfile_name, functions_of_time); + + CHECK_THROWS_WITH((FromVolumeFile{filename, subfile_name, 0.1}), + Catch::Matchers::ContainsSubstring( + "No function of time named Expansion in the subfile ")); + + subfile_name += "0"; + functions_of_time.clear(); + functions_of_time["Expansion"] = + std::make_unique>( + 0.0, + std::array{DataVector{3, 0.0}, DataVector{3, 1.0}, + DataVector{3, 0.0}}, + 1.0); + + write_volume_data(filename, subfile_name, functions_of_time); + + CHECK_THROWS_WITH( + (FromVolumeFile{filename, subfile_name, 10.0}), + Catch::Matchers::ContainsSubstring("Expansion: The requested time") and + Catch::Matchers::ContainsSubstring( + "is out of the range of the function of time")); + + if (file_system::check_if_file_exists(filename)) { + file_system::rm(filename, true); + } +} + +SPECTRE_TEST_CASE("Unit.Domain.Creators.TimeDependentOptions.FromVolumeFile", + "[Unit][Domain]") { + domain::FunctionsOfTime::register_derived_with_charm(); + test(); + test(); + test(); + test_errors(); +} +} // namespace