diff --git a/src/Domain/Creators/BinaryCompactObject.cpp b/src/Domain/Creators/BinaryCompactObject.cpp index ead7f92294ed..1c9cca7a782d 100644 --- a/src/Domain/Creators/BinaryCompactObject.cpp +++ b/src/Domain/Creators/BinaryCompactObject.cpp @@ -32,9 +32,9 @@ #include "Domain/CoordinateMaps/ProductMaps.hpp" #include "Domain/CoordinateMaps/ProductMaps.tpp" #include "Domain/CoordinateMaps/Wedge.hpp" -#include "Domain/Creators/BinaryCompactObjectHelpers.hpp" #include "Domain/Creators/DomainCreator.hpp" // IWYU pragma: keep #include "Domain/Creators/ExpandOverBlocks.hpp" +#include "Domain/Creators/TimeDependentOptions/BinaryCompactObject.hpp" #include "Domain/Domain.hpp" #include "Domain/DomainHelpers.hpp" #include "Domain/ExcisionSphere.hpp" @@ -50,6 +50,20 @@ struct BlockLogical; } // namespace Frame namespace domain::creators { +namespace bco { +std::unordered_map> +create_grid_anchors(const std::array& center_a, + const std::array& center_b) { + std::unordered_map> result{}; + result["Center" + get_output(ObjectLabel::A)] = + tnsr::I{center_a}; + result["Center" + get_output(ObjectLabel::B)] = + tnsr::I{center_b}; + result["Center"] = tnsr::I{std::array{0.0, 0.0, 0.0}}; + + return result; +} +} // namespace bco bool BinaryCompactObject::Object::is_excised() const { return inner_boundary_condition.has_value(); @@ -261,7 +275,7 @@ BinaryCompactObject::BinaryCompactObject( add_object_region("ObjectB", "Cube"); // 6 blocks first_outer_shell_block_ += 12; } - add_outer_region("Envelope"); // 10 blocks + add_outer_region("Envelope"); // 10 blocks first_outer_shell_block_ += 10; add_outer_region("OuterShell"); // 10 blocks diff --git a/src/Domain/Creators/BinaryCompactObject.hpp b/src/Domain/Creators/BinaryCompactObject.hpp index ecf97e107a7b..42214d2deece 100644 --- a/src/Domain/Creators/BinaryCompactObject.hpp +++ b/src/Domain/Creators/BinaryCompactObject.hpp @@ -20,8 +20,8 @@ #include "Domain/BoundaryConditions/GetBoundaryConditionsBase.hpp" #include "Domain/CoordinateMaps/CoordinateMap.hpp" #include "Domain/CoordinateMaps/Distribution.hpp" -#include "Domain/Creators/BinaryCompactObjectHelpers.hpp" #include "Domain/Creators/DomainCreator.hpp" +#include "Domain/Creators/TimeDependentOptions/BinaryCompactObject.hpp" #include "Domain/Domain.hpp" #include "Domain/Structure/DirectionMap.hpp" #include "Options/Auto.hpp" @@ -66,6 +66,22 @@ struct BlockLogical; namespace domain { namespace creators { +namespace bco { +/*! + * \brief Create a set of centers of objects for the binary domains. + * + * \details Will add the following centers to the set: + * + * - Center: The origin + * - CenterA: Center of object A + * - CenterB: Center of object B + * + * \return Object required by the DomainCreator%s + */ +std::unordered_map> +create_grid_anchors(const std::array& center_a, + const std::array& center_b); +} // namespace bco /*! * \ingroup ComputationalDomainGroup diff --git a/src/Domain/Creators/CMakeLists.txt b/src/Domain/Creators/CMakeLists.txt index e5a04729aab7..428692aa1261 100644 --- a/src/Domain/Creators/CMakeLists.txt +++ b/src/Domain/Creators/CMakeLists.txt @@ -10,7 +10,6 @@ spectre_target_sources( PRIVATE AlignedLattice.cpp BinaryCompactObject.cpp - BinaryCompactObjectHelpers.cpp BlockGroups.cpp Brick.cpp Cylinder.cpp @@ -24,9 +23,7 @@ spectre_target_sources( RotatedBricks.cpp RotatedIntervals.cpp RotatedRectangles.cpp - ShapeMapOptions.cpp Sphere.cpp - SphereTimeDependentMaps.cpp ) spectre_target_headers( @@ -35,7 +32,6 @@ spectre_target_headers( HEADERS AlignedLattice.hpp BinaryCompactObject.hpp - BinaryCompactObjectHelpers.hpp BlockGroups.hpp Brick.hpp Cylinder.hpp @@ -56,15 +52,14 @@ spectre_target_headers( RotatedBricks.hpp RotatedIntervals.hpp RotatedRectangles.hpp - ShapeMapOptions.hpp Sphere.hpp - SphereTimeDependentMaps.hpp ) target_link_libraries( ${LIBRARY} PRIVATE DomainStructure + H5 PUBLIC CoordinateMaps DataStructures @@ -82,3 +77,4 @@ target_link_libraries( add_subdirectory(Python) add_subdirectory(Tags) add_subdirectory(TimeDependence) +add_subdirectory(TimeDependentOptions) diff --git a/src/Domain/Creators/CylindricalBinaryCompactObject.cpp b/src/Domain/Creators/CylindricalBinaryCompactObject.cpp index 577c0527236d..d57a7798b553 100644 --- a/src/Domain/Creators/CylindricalBinaryCompactObject.cpp +++ b/src/Domain/Creators/CylindricalBinaryCompactObject.cpp @@ -20,8 +20,9 @@ #include "Domain/CoordinateMaps/UniformCylindricalFlatEndcap.hpp" #include "Domain/CoordinateMaps/UniformCylindricalSide.hpp" #include "Domain/CoordinateMaps/Wedge.hpp" -#include "Domain/Creators/BinaryCompactObjectHelpers.hpp" +#include "Domain/Creators/BinaryCompactObject.hpp" #include "Domain/Creators/ExpandOverBlocks.hpp" +#include "Domain/Creators/TimeDependentOptions/BinaryCompactObject.hpp" #include "Domain/DomainHelpers.hpp" #include "Domain/ExcisionSphere.hpp" #include "Domain/FunctionsOfTime/FixedSpeedCubic.hpp" diff --git a/src/Domain/Creators/CylindricalBinaryCompactObject.hpp b/src/Domain/Creators/CylindricalBinaryCompactObject.hpp index a70cb8c66127..28aef8061e45 100644 --- a/src/Domain/Creators/CylindricalBinaryCompactObject.hpp +++ b/src/Domain/Creators/CylindricalBinaryCompactObject.hpp @@ -16,8 +16,8 @@ #include "Domain/BoundaryConditions/BoundaryCondition.hpp" #include "Domain/BoundaryConditions/GetBoundaryConditionsBase.hpp" #include "Domain/CoordinateMaps/CoordinateMap.hpp" -#include "Domain/Creators/BinaryCompactObjectHelpers.hpp" #include "Domain/Creators/DomainCreator.hpp" +#include "Domain/Creators/TimeDependentOptions/BinaryCompactObject.hpp" #include "Domain/Domain.hpp" #include "Domain/Structure/DirectionMap.hpp" #include "Domain/Structure/ObjectLabel.hpp" diff --git a/src/Domain/Creators/Sphere.hpp b/src/Domain/Creators/Sphere.hpp index 91b986b65c25..9dcd1c27e845 100644 --- a/src/Domain/Creators/Sphere.hpp +++ b/src/Domain/Creators/Sphere.hpp @@ -16,8 +16,8 @@ #include "Domain/BoundaryConditions/GetBoundaryConditionsBase.hpp" #include "Domain/CoordinateMaps/Distribution.hpp" #include "Domain/Creators/DomainCreator.hpp" -#include "Domain/Creators/SphereTimeDependentMaps.hpp" #include "Domain/Creators/TimeDependence/TimeDependence.hpp" +#include "Domain/Creators/TimeDependentOptions/Sphere.hpp" #include "Domain/Domain.hpp" #include "Domain/Structure/DirectionMap.hpp" #include "Options/Auto.hpp" diff --git a/src/Domain/Creators/BinaryCompactObjectHelpers.cpp b/src/Domain/Creators/TimeDependentOptions/BinaryCompactObject.cpp similarity index 87% rename from src/Domain/Creators/BinaryCompactObjectHelpers.cpp rename to src/Domain/Creators/TimeDependentOptions/BinaryCompactObject.cpp index 1189a38397d9..a8782dc9683b 100644 --- a/src/Domain/Creators/BinaryCompactObjectHelpers.cpp +++ b/src/Domain/Creators/TimeDependentOptions/BinaryCompactObject.cpp @@ -1,7 +1,7 @@ // Distributed under the MIT License. // See LICENSE.txt for details. -#include "Domain/Creators/BinaryCompactObjectHelpers.hpp" +#include "Domain/Creators/TimeDependentOptions/BinaryCompactObject.hpp" #include #include @@ -16,7 +16,7 @@ #include "Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/ShapeMapTransitionFunction.hpp" #include "Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/SphereTransition.hpp" #include "Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/Wedge.hpp" -#include "Domain/Creators/ShapeMapOptions.hpp" +#include "Domain/Creators/TimeDependentOptions/ShapeMap.hpp" #include "Domain/FunctionsOfTime/FixedSpeedCubic.hpp" #include "Domain/FunctionsOfTime/FunctionOfTime.hpp" #include "Domain/FunctionsOfTime/PiecewisePolynomial.hpp" @@ -31,37 +31,25 @@ #include "Utilities/StdArrayHelpers.hpp" namespace domain::creators::bco { -std::unordered_map> -create_grid_anchors(const std::array& center_a, - const std::array& center_b) { - std::unordered_map> result{}; - result["Center" + get_output(ObjectLabel::A)] = - tnsr::I{center_a}; - result["Center" + get_output(ObjectLabel::B)] = - tnsr::I{center_b}; - result["Center"] = tnsr::I{std::array{0.0, 0.0, 0.0}}; - - return result; -} - template TimeDependentMapOptions::TimeDependentMapOptions( double initial_time, std::optional expansion_map_options, - std::optional rotation_options, + std::optional rotation_map_options, std::optional translation_map_options, std::optional> shape_options_A, std::optional> shape_options_B, const Options::Context& context) : initial_time_(initial_time), expansion_map_options_(expansion_map_options), - rotation_options_(rotation_options), - translation_options_(translation_map_options), + rotation_map_options_(rotation_map_options), + translation_map_options_(translation_map_options), shape_options_A_(shape_options_A), shape_options_B_(shape_options_B) { - if (not(expansion_map_options_.has_value() or rotation_options_.has_value() or - translation_options_.has_value() or shape_options_A_.has_value() or - shape_options_B_.has_value())) { + if (not(expansion_map_options_.has_value() or + rotation_map_options_.has_value() or + translation_map_options_.has_value() or + shape_options_A_.has_value() or shape_options_B_.has_value())) { PARSE_ERROR(context, "Time dependent map options were specified, but all options " "were 'None'. If you don't want time dependent maps, specify " @@ -136,37 +124,40 @@ TimeDependentMapOptions::create_functions_of_time( // (omega) to determine map parameters. In theory we could determine // each initial angle from the input axis-angle representation, but // we don't need to. - if (rotation_options_.has_value()) { + if (rotation_map_options_.has_value()) { result[rotation_name] = std::make_unique< FunctionsOfTime::QuaternionFunctionOfTime<3>>( initial_time_, std::array{DataVector{1.0, 0.0, 0.0, 0.0}}, std::array{ {{3, 0.0}, - {gsl::at(rotation_options_.value().initial_angular_velocity, 0), - gsl::at(rotation_options_.value().initial_angular_velocity, 1), - gsl::at(rotation_options_.value().initial_angular_velocity, 2)}, + {gsl::at(rotation_map_options_.value().initial_angular_velocity, + 0), + gsl::at(rotation_map_options_.value().initial_angular_velocity, + 1), + gsl::at(rotation_map_options_.value().initial_angular_velocity, + 2)}, {3, 0.0}, {3, 0.0}}}, expiration_times.at(rotation_name)); } // TranslationMap FunctionOfTime - if (translation_options_.has_value()) { - result[translation_name] = - std::make_unique>( - initial_time_, - std::array{ - {{gsl::at(translation_options_.value().initial_values, 0)[0], - gsl::at(translation_options_.value().initial_values, 0)[1], - gsl::at(translation_options_.value().initial_values, 0)[2]}, - {gsl::at(translation_options_.value().initial_values, 1)[0], - gsl::at(translation_options_.value().initial_values, 1)[1], - gsl::at(translation_options_.value().initial_values, 1)[2]}, - {gsl::at(translation_options_.value().initial_values, 2)[0], - gsl::at(translation_options_.value().initial_values, 2)[1], - gsl::at(translation_options_.value().initial_values, 2)[2]}}}, - expiration_times.at(translation_name)); + if (translation_map_options_.has_value()) { + result[translation_name] = std::make_unique< + FunctionsOfTime::PiecewisePolynomial<2>>( + initial_time_, + std::array{ + {{gsl::at(translation_map_options_.value().initial_values, 0)[0], + gsl::at(translation_map_options_.value().initial_values, 0)[1], + gsl::at(translation_map_options_.value().initial_values, 0)[2]}, + {gsl::at(translation_map_options_.value().initial_values, 1)[0], + gsl::at(translation_map_options_.value().initial_values, 1)[1], + gsl::at(translation_map_options_.value().initial_values, 1)[2]}, + {gsl::at(translation_map_options_.value().initial_values, 2)[0], + gsl::at(translation_map_options_.value().initial_values, 2)[1], + gsl::at(translation_map_options_.value().initial_values, 2)[2]}}}, + expiration_times.at(translation_name)); } // Size and Shape FunctionOfTime for objects A and B @@ -218,17 +209,17 @@ void TimeDependentMapOptions::build_maps( const std::optional>& object_B_radii, const double envelope_radius, const double domain_outer_radius) { - if (expansion_map_options_.has_value() or rotation_options_.has_value() or - translation_options_.has_value()) { + if (expansion_map_options_.has_value() or rotation_map_options_.has_value() or + translation_map_options_.has_value()) { rot_scale_trans_map_ = std::make_pair( RotScaleTrans{ expansion_map_options_.has_value() ? std::make_pair(expansion_name, expansion_outer_boundary_name) : std::optional>{}, - rotation_options_.has_value() ? rotation_name - : std::optional{}, - translation_options_.has_value() ? translation_name - : std::optional{}, + rotation_map_options_.has_value() ? rotation_name + : std::optional{}, + translation_map_options_.has_value() ? translation_name + : std::optional{}, envelope_radius, domain_outer_radius, domain::CoordinateMaps::TimeDependent::RotScaleTrans< 3>::BlockRegion::Inner}, @@ -236,10 +227,10 @@ void TimeDependentMapOptions::build_maps( expansion_map_options_.has_value() ? std::make_pair(expansion_name, expansion_outer_boundary_name) : std::optional>{}, - rotation_options_.has_value() ? rotation_name - : std::optional{}, - translation_options_.has_value() ? translation_name - : std::optional{}, + rotation_map_options_.has_value() ? rotation_name + : std::optional{}, + translation_map_options_.has_value() ? translation_name + : std::optional{}, envelope_radius, domain_outer_radius, domain::CoordinateMaps::TimeDependent::RotScaleTrans< 3>::BlockRegion::Transition}); diff --git a/src/Domain/Creators/BinaryCompactObjectHelpers.hpp b/src/Domain/Creators/TimeDependentOptions/BinaryCompactObject.hpp similarity index 95% rename from src/Domain/Creators/BinaryCompactObjectHelpers.hpp rename to src/Domain/Creators/TimeDependentOptions/BinaryCompactObject.hpp index 517969bbc06b..1fa2fe7968f9 100644 --- a/src/Domain/Creators/BinaryCompactObjectHelpers.hpp +++ b/src/Domain/Creators/TimeDependentOptions/BinaryCompactObject.hpp @@ -17,8 +17,8 @@ #include "Domain/CoordinateMaps/TimeDependent/Rotation.hpp" #include "Domain/CoordinateMaps/TimeDependent/Shape.hpp" #include "Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/ShapeMapTransitionFunction.hpp" -#include "Domain/Creators/ShapeMapOptions.hpp" -#include "Domain/Creators/SphereTimeDependentMaps.hpp" +#include "Domain/Creators/TimeDependentOptions/ShapeMap.hpp" +#include "Domain/Creators/TimeDependentOptions/Sphere.hpp" #include "Domain/FunctionsOfTime/FunctionOfTime.hpp" #include "Domain/Structure/ObjectLabel.hpp" #include "Options/Auto.hpp" @@ -39,22 +39,6 @@ struct Inertial; /// Namespace used to hold things used in both the BinaryCompactObject and /// CylindricalBinaryCompactObject domain creators. namespace domain::creators::bco { - -/*! - * \brief Create a set of centers of objects for the binary domains. - * - * \details Will add the following centers to the set: - * - * - Center: The origin - * - CenterA: Center of object A - * - CenterB: Center of object B - * - * \return Object required by the DomainCreator%s - */ -std::unordered_map> -create_grid_anchors(const std::array& center_a, - const std::array& center_b); - namespace detail { // Convenience type alias template @@ -259,8 +243,8 @@ struct TimeDependentMapOptions { TimeDependentMapOptions( double initial_time, std::optional expansion_map_options, - std::optional rotation_options, - std::optional translation_options, + std::optional rotation_map_options, + std::optional translation_map_options, std::optional> shape_options_A, std::optional> shape_options_B, const Options::Context& context = {}); @@ -384,8 +368,8 @@ struct TimeDependentMapOptions { double initial_time_{std::numeric_limits::signaling_NaN()}; std::optional expansion_map_options_{}; - std::optional rotation_options_{}; - std::optional translation_options_{}; + std::optional rotation_map_options_{}; + std::optional translation_map_options_{}; std::optional> shape_options_A_{}; std::optional> shape_options_B_{}; std::array, 2> inner_radii_{}; diff --git a/src/Domain/Creators/TimeDependentOptions/CMakeLists.txt b/src/Domain/Creators/TimeDependentOptions/CMakeLists.txt new file mode 100644 index 000000000000..1c9d0d4bbaf7 --- /dev/null +++ b/src/Domain/Creators/TimeDependentOptions/CMakeLists.txt @@ -0,0 +1,21 @@ +# Distributed under the MIT License. +# See LICENSE.txt for details. + +spectre_target_sources( + ${LIBRARY} + PRIVATE + BinaryCompactObject.cpp + FromVolumeFile.cpp + ShapeMap.cpp + Sphere.cpp +) + +spectre_target_headers( + ${LIBRARY} + 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 000000000000..f3513764c8af --- /dev/null +++ b/src/Domain/Creators/TimeDependentOptions/FromVolumeFile.cpp @@ -0,0 +1,211 @@ +// 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 "Domain/Structure/ObjectLabel.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]) { + using ::operator<<; + 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 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); + + bool is_quat_fot = false; + const auto set_values = [this, &function_of_time, &time, &is_quat_fot, + &context]() { + 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); + auto all_angle_values = + quat_function_of_time->angle_func_and_all_derivs(time); + if (all_angle_values.size() > angle_values.size()) { + PARSE_ERROR(context, + "FromVolumeFile bad size of angle values (needed at most " + << angle_values.size() << ", got " + << all_angle_values.size() << ")"); + } + angle_values = make_array<4>(DataVector{3, 0.0}); + for (size_t i = 0; i < all_angle_values.size(); i++) { + gsl::at(angle_values, i) = std::move(all_angle_values[i]); + } + is_quat_fot = true; + } + }; + + set_values.operator()<2>(); + set_values.operator()<3>(); + + // The FoT isn't a QuaternionFunctionOfTime so we must compute the angle + // derivatives ourselves + if (not is_quat_fot) { + 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)); + } + + angle_values = make_array<4>(DataVector{4, 0.0}); + // We can't compute the angle so set it to zero + // We'd need more quaternion derivatives to compute the higher angle + // derivatives. + // dtq = 0.5 * q * w => w = 2.0 * conj(q) * dtq + // d2tq = 0.5 * (dtq * w + q * dtw) + // => dtw = conj(q) * (2.0 * d2tq - dtq * w) + 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]))); + // The above DataVectors have size 4, so we need to make them have size 3. + for (size_t i = 0; i < angle_values.size(); 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]}; + } + } +} + +template +FromVolumeFile>::FromVolumeFile( + const std::string& h5_filename, const std::string& subfile_name, + const double time, const Options::Context& context) { + const std::string object = domain::name(Object); + const std::string shape_function_of_time_name{"Shape" + object}; + const std::string size_function_of_time_name{"Size" + object}; + + 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); + auto all_size_values = size_function_of_time->func_and_all_derivs(time); + if (all_size_values.size() > size_values.size()) { + PARSE_ERROR(context, + "FromVolumeFile bad size of Size values (needed at most " + << size_values.size() << ", got " << all_size_values.size() + << ")"); + } + size_values = make_array<4>(DataVector{0.0}); + for (size_t i = 0; i < all_size_values.size(); i++) { + gsl::at(size_values, i) = std::move(all_size_values[i]); + } +} + +template struct FromVolumeFile; +template struct FromVolumeFile>; +template struct FromVolumeFile>; +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 000000000000..6b3bb6cd1211 --- /dev/null +++ b/src/Domain/Creators/TimeDependentOptions/FromVolumeFile.hpp @@ -0,0 +1,122 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include +#include + +#include "DataStructures/DataVector.hpp" +#include "Domain/Structure/ObjectLabel.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 {}; +template +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 { + FromVolumeFile() = default; + FromVolumeFile(const std::string& h5_filename, + const std::string& subfile_name, double time, + 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/src/Domain/Creators/ShapeMapOptions.cpp b/src/Domain/Creators/TimeDependentOptions/ShapeMap.cpp similarity index 98% rename from src/Domain/Creators/ShapeMapOptions.cpp rename to src/Domain/Creators/TimeDependentOptions/ShapeMap.cpp index bd2d7eec0cff..8f11d859ed19 100644 --- a/src/Domain/Creators/ShapeMapOptions.cpp +++ b/src/Domain/Creators/TimeDependentOptions/ShapeMap.cpp @@ -1,7 +1,7 @@ // Distributed under the MIT License. // See LICENSE.txt for details. -#include "Domain/Creators/ShapeMapOptions.hpp" +#include "Domain/Creators/TimeDependentOptions/ShapeMap.hpp" #include #include diff --git a/src/Domain/Creators/ShapeMapOptions.hpp b/src/Domain/Creators/TimeDependentOptions/ShapeMap.hpp similarity index 100% rename from src/Domain/Creators/ShapeMapOptions.hpp rename to src/Domain/Creators/TimeDependentOptions/ShapeMap.hpp diff --git a/src/Domain/Creators/SphereTimeDependentMaps.cpp b/src/Domain/Creators/TimeDependentOptions/Sphere.cpp similarity index 98% rename from src/Domain/Creators/SphereTimeDependentMaps.cpp rename to src/Domain/Creators/TimeDependentOptions/Sphere.cpp index fba8cf1eacbe..84ff57e457d3 100644 --- a/src/Domain/Creators/SphereTimeDependentMaps.cpp +++ b/src/Domain/Creators/TimeDependentOptions/Sphere.cpp @@ -1,7 +1,7 @@ // Distributed under the MIT License. // See LICENSE.txt for details. -#include "Domain/Creators/SphereTimeDependentMaps.hpp" +#include "Domain/Creators/TimeDependentOptions/Sphere.hpp" #include #include @@ -17,7 +17,7 @@ #include "Domain/CoordinateMaps/CoordinateMap.tpp" #include "Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/ShapeMapTransitionFunction.hpp" #include "Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/SphereTransition.hpp" -#include "Domain/Creators/ShapeMapOptions.hpp" +#include "Domain/Creators/TimeDependentOptions/ShapeMap.hpp" #include "Domain/FunctionsOfTime/FixedSpeedCubic.hpp" #include "Domain/FunctionsOfTime/FunctionOfTime.hpp" #include "Domain/FunctionsOfTime/PiecewisePolynomial.hpp" diff --git a/src/Domain/Creators/SphereTimeDependentMaps.hpp b/src/Domain/Creators/TimeDependentOptions/Sphere.hpp similarity index 99% rename from src/Domain/Creators/SphereTimeDependentMaps.hpp rename to src/Domain/Creators/TimeDependentOptions/Sphere.hpp index 430688ab0b21..8b9af3abb9e0 100644 --- a/src/Domain/Creators/SphereTimeDependentMaps.hpp +++ b/src/Domain/Creators/TimeDependentOptions/Sphere.hpp @@ -16,7 +16,7 @@ #include "Domain/CoordinateMaps/TimeDependent/RotScaleTrans.hpp" #include "Domain/CoordinateMaps/TimeDependent/Shape.hpp" #include "Domain/CoordinateMaps/TimeDependent/Translation.hpp" -#include "Domain/Creators/ShapeMapOptions.hpp" +#include "Domain/Creators/TimeDependentOptions/ShapeMap.hpp" #include "Domain/FunctionsOfTime/FunctionOfTime.hpp" #include "Domain/Structure/ObjectLabel.hpp" #include "Options/Auto.hpp" diff --git a/src/Domain/FunctionsOfTime/FunctionOfTime.hpp b/src/Domain/FunctionsOfTime/FunctionOfTime.hpp index cc03790d1bc9..017fdef6ee51 100644 --- a/src/Domain/FunctionsOfTime/FunctionOfTime.hpp +++ b/src/Domain/FunctionsOfTime/FunctionOfTime.hpp @@ -81,6 +81,19 @@ class FunctionOfTime : public PUP::able { /// The DataVector can be of any size virtual std::array func_and_2_derivs(double t) const = 0; + /// \brief All derivatives a function of time has to offer (because it can be + /// more than 2) + /// + /// \details Defaults to the return values from `func_and_2_derivs` since some + /// functions of time only go up to 2 derivatives by design + virtual std::vector func_and_all_derivs(double t) const { + std::array tmp_func_and_2_derivs = func_and_2_derivs(t); + + return std::vector{std::move(tmp_func_and_2_derivs[0]), + std::move(tmp_func_and_2_derivs[1]), + std::move(tmp_func_and_2_derivs[2])}; + } + WRAPPED_PUPable_abstract(FunctionOfTime); // NOLINT }; } // namespace FunctionsOfTime diff --git a/src/Domain/FunctionsOfTime/PiecewisePolynomial.cpp b/src/Domain/FunctionsOfTime/PiecewisePolynomial.cpp index 102598df8ace..5b893f1efec6 100644 --- a/src/Domain/FunctionsOfTime/PiecewisePolynomial.cpp +++ b/src/Domain/FunctionsOfTime/PiecewisePolynomial.cpp @@ -100,6 +100,19 @@ PiecewisePolynomial::func_and_derivs(const double t) const { return result; } +template +std::vector PiecewisePolynomial::func_and_all_derivs( + double t) const { + auto tmp_func_and_max_deriv = func_and_derivs(t); + std::vector result{}; + result.resize(tmp_func_and_max_deriv.size()); + for (size_t i = 0; i < tmp_func_and_max_deriv.size(); i++) { + result[i] = std::move(gsl::at(tmp_func_and_max_deriv, i)); + } + + return result; +} + template void PiecewisePolynomial::store_entry( const double time_of_update, diff --git a/src/Domain/FunctionsOfTime/PiecewisePolynomial.hpp b/src/Domain/FunctionsOfTime/PiecewisePolynomial.hpp index 59998de98dc2..93871e71ca16 100644 --- a/src/Domain/FunctionsOfTime/PiecewisePolynomial.hpp +++ b/src/Domain/FunctionsOfTime/PiecewisePolynomial.hpp @@ -10,6 +10,7 @@ #include #include #include +#include #include "DataStructures/DataVector.hpp" #include "Domain/FunctionsOfTime/FunctionOfTime.hpp" @@ -66,6 +67,10 @@ class PiecewisePolynomial : public FunctionOfTime { return func_and_derivs<2>(t); } + /// Return the function and all derivs up to and including the `MaxDeriv` at + /// an arbitrary time `t`. + std::vector func_and_all_derivs(double t) const override; + /// Updates the `MaxDeriv`th derivative of the function at the given time. /// `updated_max_deriv` is a vector of the `MaxDeriv`ths for each component. /// `next_expiration_time` is the next expiration time. diff --git a/src/Domain/FunctionsOfTime/QuaternionFunctionOfTime.cpp b/src/Domain/FunctionsOfTime/QuaternionFunctionOfTime.cpp index 3a91c580810d..1f65badfd625 100644 --- a/src/Domain/FunctionsOfTime/QuaternionFunctionOfTime.cpp +++ b/src/Domain/FunctionsOfTime/QuaternionFunctionOfTime.cpp @@ -336,6 +336,42 @@ QuaternionFunctionOfTime::quat_func_and_2_derivs( quaternion_to_datavector(dt2quat)}; } +template +std::array +QuaternionFunctionOfTime::quat_func_and_3_derivs( + const double t) const { + boost::math::quaternion quat = setup_func(t); + + // Get angle and however many derivatives we need + std::vector angle_and_all_derivs = + angle_f_of_t_.func_and_all_derivs(t); + + if (angle_and_all_derivs.size() < 4) { + ERROR( + "Need more angle derivs to compute the third derivative of the " + "quaternion. Currently only have " + << MaxDeriv); + } + + boost::math::quaternion omega = + datavector_to_quaternion(angle_and_all_derivs[1]); + boost::math::quaternion dtomega = + datavector_to_quaternion(angle_and_all_derivs[2]); + + boost::math::quaternion dt2omega = + datavector_to_quaternion(angle_and_all_derivs[3]); + + boost::math::quaternion dtquat = 0.5 * quat * omega; + boost::math::quaternion dt2quat = + 0.5 * (dtquat * omega + quat * dtomega); + boost::math::quaternion dt3quat = + 0.5 * (dt2quat * omega + 2.0 * dtquat * dtomega + quat * dt2omega); + + return std::array{ + quaternion_to_datavector(quat), quaternion_to_datavector(dtquat), + quaternion_to_datavector(dt2quat), quaternion_to_datavector(dt3quat)}; +} + template bool operator==(const QuaternionFunctionOfTime& lhs, const QuaternionFunctionOfTime& rhs) { @@ -384,7 +420,7 @@ std::ostream& operator<<( return os; } -// do explicit instantiation of MaxDeriv = {2,3,4} +// do explicit instantiation of MaxDeriv = {2,3} #define DIM(data) BOOST_PP_TUPLE_ELEM(0, data) #define INSTANTIATE(_, data) \ @@ -398,7 +434,7 @@ std::ostream& operator<<( template std::ostream& operator<<( \ std::ostream& os, const QuaternionFunctionOfTime&); -GENERATE_INSTANTIATIONS(INSTANTIATE, (2, 3, 4)) +GENERATE_INSTANTIATIONS(INSTANTIATE, (2, 3)) #undef DIM #undef INSTANTIATE diff --git a/src/Domain/FunctionsOfTime/QuaternionFunctionOfTime.hpp b/src/Domain/FunctionsOfTime/QuaternionFunctionOfTime.hpp index 0c6199878c10..97bdf0f99175 100644 --- a/src/Domain/FunctionsOfTime/QuaternionFunctionOfTime.hpp +++ b/src/Domain/FunctionsOfTime/QuaternionFunctionOfTime.hpp @@ -10,6 +10,7 @@ #include #include #include +#include #include "DataStructures/DataVector.hpp" #include "Domain/FunctionsOfTime/FunctionOfTime.hpp" @@ -128,6 +129,10 @@ class QuaternionFunctionOfTime : public FunctionOfTime { /// time `t`. std::array quat_func_and_2_derivs(double t) const; + /// Returns the quaternion and the first three derivatives at an arbitrary + /// time `t`. + std::array quat_func_and_3_derivs(double t) const; + /// Returns stored angle at an arbitrary time `t`. std::array angle_func(const double t) const { return angle_f_of_t_.func(t); @@ -145,6 +150,10 @@ class QuaternionFunctionOfTime : public FunctionOfTime { return angle_f_of_t_.func_and_2_derivs(t); } + std::vector angle_func_and_all_derivs(const double t) const { + return angle_f_of_t_.func_and_all_derivs(t); + } + private: template friend bool operator==( // NOLINT(readability-redundant-declaration) diff --git a/src/Domain/FunctionsOfTime/RegisterDerivedWithCharm.cpp b/src/Domain/FunctionsOfTime/RegisterDerivedWithCharm.cpp index 3ac32e1358ea..549c3c5040de 100644 --- a/src/Domain/FunctionsOfTime/RegisterDerivedWithCharm.cpp +++ b/src/Domain/FunctionsOfTime/RegisterDerivedWithCharm.cpp @@ -29,7 +29,6 @@ void register_derived_with_charm() { FunctionsOfTime::PiecewisePolynomial<4>, FunctionsOfTime::QuaternionFunctionOfTime<2>, FunctionsOfTime::QuaternionFunctionOfTime<3>, - FunctionsOfTime::QuaternionFunctionOfTime<4>, FunctionsOfTime::SettleToConstant, FunctionsOfTime::SettleToConstantQuaternion>(); } diff --git a/src/Evolution/Ringdown/StrahlkorperCoefsInRingdownDistortedFrame.cpp b/src/Evolution/Ringdown/StrahlkorperCoefsInRingdownDistortedFrame.cpp index 3297c6c1e769..6733dcecfc4d 100644 --- a/src/Evolution/Ringdown/StrahlkorperCoefsInRingdownDistortedFrame.cpp +++ b/src/Evolution/Ringdown/StrahlkorperCoefsInRingdownDistortedFrame.cpp @@ -10,7 +10,7 @@ #include "DataStructures/Matrix.hpp" #include "Domain/CoordinateMaps/Distribution.hpp" #include "Domain/Creators/Sphere.hpp" -#include "Domain/Creators/SphereTimeDependentMaps.hpp" +#include "Domain/Creators/TimeDependentOptions/Sphere.hpp" #include "Domain/StrahlkorperTransformations.hpp" #include "IO/H5/Dat.hpp" #include "IO/H5/File.hpp" @@ -32,64 +32,63 @@ std::vector strahlkorper_coefs_in_ringdown_distorted_frame( ylm::read_surface_ylm( path_to_horizons_h5, surface_subfile_name, requested_number_of_times_from_end); - const size_t l_max{ahc_inertial_h5[0].l_max()}; + const size_t l_max{ahc_inertial_h5[0].l_max()}; - std::vector ahc_times{}; - { - // Read the AhC times from the H5 file - const h5::H5File ahc_h5_file{ - path_to_horizons_h5}; - const auto& dat = ahc_h5_file.get(surface_subfile_name); - const Matrix& coefs_for_times = dat.get_data_subset( - {0}, dat.get_dimensions()[0] - requested_number_of_times_from_end, - ahc_inertial_h5.size()); - for (size_t i = 0; i < coefs_for_times.rows(); ++i) { - ahc_times.push_back(coefs_for_times(i, 0)); - } + std::vector ahc_times{}; + { + // Read the AhC times from the H5 file + const h5::H5File ahc_h5_file{path_to_horizons_h5}; + const auto& dat = ahc_h5_file.get(surface_subfile_name); + const Matrix& coefs_for_times = dat.get_data_subset( + {0}, dat.get_dimensions()[0] - requested_number_of_times_from_end, + ahc_inertial_h5.size()); + for (size_t i = 0; i < coefs_for_times.rows(); ++i) { + ahc_times.push_back(coefs_for_times(i, 0)); } - // Create a time-dependent domain; only the the time-dependent map options - // matter; the domain is just a spherical shell with inner and outer - // radii chosen so any conceivable common horizon will fit between them. - const domain::creators::sphere::TimeDependentMapOptions::ShapeMapOptions - shape_map_options{l_max, std::nullopt}; - const domain::creators::sphere::TimeDependentMapOptions::ExpansionMapOptions - expansion_map_options{exp_func_and_2_derivs, settling_timescale, - exp_outer_bdry_func_and_2_derivs, - settling_timescale}; - const domain::creators::sphere::TimeDependentMapOptions::RotationMapOptions - rotation_map_options{rot_func_and_2_derivs, settling_timescale}; - const domain::creators::sphere::TimeDependentMapOptions - time_dependent_map_options{match_time, shape_map_options, - rotation_map_options, expansion_map_options, - std::nullopt}; - const domain::creators::Sphere domain_creator{ - 0.01, - 200.0, - // nullptr because no boundary condition - domain::creators::Sphere::Excision{nullptr}, - static_cast(0), - static_cast(5), - false, - std::nullopt, - {100.0}, - domain::CoordinateMaps::Distribution::Linear, - ShellWedges::All, - time_dependent_map_options}; + } + // Create a time-dependent domain; only the the time-dependent map options + // matter; the domain is just a spherical shell with inner and outer + // radii chosen so any conceivable common horizon will fit between them. + const domain::creators::sphere::TimeDependentMapOptions::ShapeMapOptions + shape_map_options{l_max, std::nullopt}; + const domain::creators::sphere::TimeDependentMapOptions::ExpansionMapOptions + expansion_map_options{exp_func_and_2_derivs, settling_timescale, + exp_outer_bdry_func_and_2_derivs, + settling_timescale}; + const domain::creators::sphere::TimeDependentMapOptions::RotationMapOptions + rotation_map_options{rot_func_and_2_derivs, settling_timescale}; + const domain::creators::sphere::TimeDependentMapOptions + time_dependent_map_options{match_time, shape_map_options, + rotation_map_options, expansion_map_options, + std::nullopt}; + const domain::creators::Sphere domain_creator{ + 0.01, + 200.0, + // nullptr because no boundary condition + domain::creators::Sphere::Excision{nullptr}, + static_cast(0), + static_cast(5), + false, + std::nullopt, + {100.0}, + domain::CoordinateMaps::Distribution::Linear, + ShellWedges::All, + time_dependent_map_options}; - const auto temporary_domain = domain_creator.create_domain(); - const auto functions_of_time = domain_creator.functions_of_time(); + const auto temporary_domain = domain_creator.create_domain(); + const auto functions_of_time = domain_creator.functions_of_time(); - // Loop over the selected horizons, transforming each to the - // ringdown distorted frame - std::vector ahc_ringdown_distorted_coefs{}; - ylm::Strahlkorper current_ahc; - for (size_t i = 0; i < requested_number_of_times_from_end; ++i) { - strahlkorper_in_different_frame( - make_not_null(¤t_ahc), gsl::at(ahc_inertial_h5, i), - temporary_domain, functions_of_time, gsl::at(ahc_times, i)); - ahc_ringdown_distorted_coefs.push_back(current_ahc.coefficients()); - } + // Loop over the selected horizons, transforming each to the + // ringdown distorted frame + std::vector ahc_ringdown_distorted_coefs{}; + ylm::Strahlkorper current_ahc; + for (size_t i = 0; i < requested_number_of_times_from_end; ++i) { + strahlkorper_in_different_frame( + make_not_null(¤t_ahc), gsl::at(ahc_inertial_h5, i), + temporary_domain, functions_of_time, gsl::at(ahc_times, i)); + ahc_ringdown_distorted_coefs.push_back(current_ahc.coefficients()); + } - return ahc_ringdown_distorted_coefs; + return ahc_ringdown_distorted_coefs; } } // namespace evolution::Ringdown diff --git a/tests/Unit/Domain/Creators/CMakeLists.txt b/tests/Unit/Domain/Creators/CMakeLists.txt index 527d78015a33..456426228b4f 100644 --- a/tests/Unit/Domain/Creators/CMakeLists.txt +++ b/tests/Unit/Domain/Creators/CMakeLists.txt @@ -1,15 +1,11 @@ # Distributed under the MIT License. # See LICENSE.txt for details. -add_subdirectory(Python) -add_subdirectory(TimeDependence) - set(LIBRARY "Test_DomainCreators") set(LIBRARY_SOURCES Test_AlignedLattice.cpp Test_BinaryCompactObject.cpp - Test_BinaryCompactObjectHelpers.cpp Test_BlockGroups.cpp Test_Brick.cpp Test_Cylinder.cpp @@ -23,11 +19,11 @@ set(LIBRARY_SOURCES Test_RotatedBricks.cpp Test_RotatedIntervals.cpp Test_RotatedRectangles.cpp - Test_ShapeMapOptions.cpp Test_Sphere.cpp - Test_SphereTimeDependentMaps.cpp ) +add_subdirectory(TimeDependentOptions) + add_test_library(${LIBRARY} "${LIBRARY_SOURCES}") target_link_libraries( @@ -42,3 +38,6 @@ target_link_libraries( SphericalHarmonics SphericalHarmonicsIO ) + +add_subdirectory(Python) +add_subdirectory(TimeDependence) diff --git a/tests/Unit/Domain/Creators/Test_BinaryCompactObject.cpp b/tests/Unit/Domain/Creators/Test_BinaryCompactObject.cpp index df24b83b8f3f..ba4448611878 100644 --- a/tests/Unit/Domain/Creators/Test_BinaryCompactObject.cpp +++ b/tests/Unit/Domain/Creators/Test_BinaryCompactObject.cpp @@ -27,7 +27,7 @@ #include "Domain/Creators/BinaryCompactObject.hpp" #include "Domain/Creators/DomainCreator.hpp" #include "Domain/Creators/OptionTags.hpp" -#include "Domain/Creators/ShapeMapOptions.hpp" +#include "Domain/Creators/TimeDependentOptions/ShapeMap.hpp" #include "Domain/Domain.hpp" #include "Domain/FunctionsOfTime/FixedSpeedCubic.hpp" #include "Domain/FunctionsOfTime/PiecewisePolynomial.hpp" diff --git a/tests/Unit/Domain/Creators/Test_CylindricalBinaryCompactObject.cpp b/tests/Unit/Domain/Creators/Test_CylindricalBinaryCompactObject.cpp index 044a78387d90..da2a7a690431 100644 --- a/tests/Unit/Domain/Creators/Test_CylindricalBinaryCompactObject.cpp +++ b/tests/Unit/Domain/Creators/Test_CylindricalBinaryCompactObject.cpp @@ -23,10 +23,10 @@ #include "DataStructures/Tensor/TypeAliases.hpp" // IWYU pragma: keep #include "Domain/Block.hpp" // IWYU pragma: keep #include "Domain/BoundaryConditions/BoundaryCondition.hpp" -#include "Domain/Creators/BinaryCompactObjectHelpers.hpp" #include "Domain/Creators/CylindricalBinaryCompactObject.hpp" #include "Domain/Creators/DomainCreator.hpp" #include "Domain/Creators/OptionTags.hpp" +#include "Domain/Creators/TimeDependentOptions/BinaryCompactObject.hpp" #include "Domain/Domain.hpp" #include "Domain/ExcisionSphere.hpp" #include "Domain/FunctionsOfTime/FixedSpeedCubic.hpp" diff --git a/tests/Unit/Domain/Creators/Test_Sphere.cpp b/tests/Unit/Domain/Creators/Test_Sphere.cpp index e025604e549e..d58f9edd6a65 100644 --- a/tests/Unit/Domain/Creators/Test_Sphere.cpp +++ b/tests/Unit/Domain/Creators/Test_Sphere.cpp @@ -33,11 +33,11 @@ #include "Domain/CoordinateMaps/Wedge.hpp" #include "Domain/Creators/DomainCreator.hpp" #include "Domain/Creators/OptionTags.hpp" -#include "Domain/Creators/ShapeMapOptions.hpp" #include "Domain/Creators/Sphere.hpp" #include "Domain/Creators/TimeDependence/None.hpp" #include "Domain/Creators/TimeDependence/RegisterDerivedWithCharm.hpp" #include "Domain/Creators/TimeDependence/UniformTranslation.hpp" +#include "Domain/Creators/TimeDependentOptions/ShapeMap.hpp" #include "Domain/Domain.hpp" #include "Domain/FunctionsOfTime/PiecewisePolynomial.hpp" #include "Domain/Structure/BlockNeighbor.hpp" diff --git a/tests/Unit/Domain/Creators/TimeDependentOptions/CMakeLists.txt b/tests/Unit/Domain/Creators/TimeDependentOptions/CMakeLists.txt new file mode 100644 index 000000000000..97123f3d6e07 --- /dev/null +++ b/tests/Unit/Domain/Creators/TimeDependentOptions/CMakeLists.txt @@ -0,0 +1,10 @@ +# Distributed under the MIT License. +# See LICENSE.txt for details. + +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/Test_BinaryCompactObjectHelpers.cpp b/tests/Unit/Domain/Creators/TimeDependentOptions/Test_BinaryCompactObject.cpp similarity index 99% rename from tests/Unit/Domain/Creators/Test_BinaryCompactObjectHelpers.cpp rename to tests/Unit/Domain/Creators/TimeDependentOptions/Test_BinaryCompactObject.cpp index 0075cc5d57f0..0a65626a5492 100644 --- a/tests/Unit/Domain/Creators/Test_BinaryCompactObjectHelpers.cpp +++ b/tests/Unit/Domain/Creators/TimeDependentOptions/Test_BinaryCompactObject.cpp @@ -10,7 +10,8 @@ #include #include "DataStructures/DataVector.hpp" -#include "Domain/Creators/BinaryCompactObjectHelpers.hpp" +#include "Domain/Creators/BinaryCompactObject.hpp" +#include "Domain/Creators/TimeDependentOptions/BinaryCompactObject.hpp" #include "Domain/FunctionsOfTime/FixedSpeedCubic.hpp" #include "Domain/FunctionsOfTime/FunctionOfTime.hpp" #include "Domain/FunctionsOfTime/PiecewisePolynomial.hpp" @@ -639,8 +640,9 @@ void test_errors() { } // namespace // [[TimeOut, 45]] -SPECTRE_TEST_CASE("Unit.Domain.Creators.BinaryCompactObjectHelpers", - "[Domain][Unit]") { +SPECTRE_TEST_CASE( + "Unit.Domain.Creators.TimeDependentOptions.BinaryCompactObject", + "[Domain][Unit]") { for (const auto& [include_expansion, include_rotation, include_translation, include_shape_a, include_shape_b] : cartesian_product(make_array(true, false), make_array(true, false), 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 000000000000..472dee16da9e --- /dev/null +++ b/tests/Unit/Domain/Creators/TimeDependentOptions/Test_FromVolumeFile.cpp @@ -0,0 +1,362 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "DataStructures/DataVector.hpp" +#include "Domain/Structure/ObjectLabel.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("Is a QuaternionFunctionOfTime"); + // 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\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}, + 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); + } + + // Write new function of time + { + auto quat_and_derivs = + functions_of_time.at(function_of_time_name)->func_and_2_derivs(0.0); + functions_of_time.clear(); + functions_of_time[function_of_time_name] = + std::make_unique>( + 0.0, std::move(quat_and_derivs), 100.0); + write_volume_data(filename, subfile_name, functions_of_time); + } + + { + INFO("Is not a QuaternionFunctionOfTime"); + // 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\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}, + 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); + } +} + +template <> +void test>() { + const std::string filename{"BlandHorseRadish.h5"}; + if (file_system::check_if_file_exists(filename)) { + file_system::rm(filename, true); + } + const std::string subfile_name{"VolumeData"}; + const std::string shape_name{"ShapeB"}; + const std::string size_name{"SizeB"}; + + std::unordered_map> + functions_of_time{}; + // For reading these in, we don't care how many components of the DataVector + // there are. Normally there'd be a lot more. + functions_of_time[shape_name] = + std::make_unique>( + 0.0, + std::array{DataVector{1, 0.0}, DataVector{1, 0.0}, + DataVector{1, 2.0}}, + 100.0); + functions_of_time[size_name] = + std::make_unique>( + 0.0, + std::array{DataVector{1, 0.1}, DataVector{1, 0.0}, DataVector{1, 0.3}, + DataVector{1, 0.0}}, + 100.0); + + write_volume_data(filename, subfile_name, functions_of_time); + + const auto from_volume_file = TestHelpers::test_creation< + domain::creators::time_dependent_options::FromVolumeFile< + domain::creators::time_dependent_options::names::ShapeSize< + domain::ObjectLabel::B>>>("H5Filename: " + filename + + "\nSubfileName: " + subfile_name + + "\nTime: 50.0\n"); + + std::array expected_shape_values{ + DataVector{2500.0}, DataVector{100.0}, DataVector{2.0}}; + std::array expected_size_values{ + DataVector{0.1 + 0.5 * 2500.0 * 0.3}, DataVector{50.0 * 0.3}, + DataVector{0.3}, DataVector{0.0}}; + + CHECK(from_volume_file.shape_values == expected_shape_values); + CHECK_ITERABLE_APPROX(from_volume_file.size_values, expected_size_values); + + 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>(); + test_errors(); +} +} // namespace diff --git a/tests/Unit/Domain/Creators/Test_ShapeMapOptions.cpp b/tests/Unit/Domain/Creators/TimeDependentOptions/Test_ShapeMap.cpp similarity index 98% rename from tests/Unit/Domain/Creators/Test_ShapeMapOptions.cpp rename to tests/Unit/Domain/Creators/TimeDependentOptions/Test_ShapeMap.cpp index 39e47566f55e..18b41797975f 100644 --- a/tests/Unit/Domain/Creators/Test_ShapeMapOptions.cpp +++ b/tests/Unit/Domain/Creators/TimeDependentOptions/Test_ShapeMap.cpp @@ -8,7 +8,7 @@ #include #include -#include "Domain/Creators/ShapeMapOptions.hpp" +#include "Domain/Creators/TimeDependentOptions/ShapeMap.hpp" #include "Domain/Structure/ObjectLabel.hpp" #include "Framework/TestCreation.hpp" #include "Framework/TestHelpers.hpp" @@ -302,7 +302,8 @@ void test_funcs(const gsl::not_null generator) { } } // namespace -SPECTRE_TEST_CASE("Unit.Domain.Creators.ShapeMapOptions", "[Domain][Unit]") { +SPECTRE_TEST_CASE("Unit.Domain.Creators.TimeDependentOptions.ShapeMap", + "[Domain][Unit]") { MAKE_GENERATOR(generator); test_kerr_schild_boyer_lindquist(); test_ylms_from_file(); diff --git a/tests/Unit/Domain/Creators/Test_SphereTimeDependentMaps.cpp b/tests/Unit/Domain/Creators/TimeDependentOptions/Test_Sphere.cpp similarity index 97% rename from tests/Unit/Domain/Creators/Test_SphereTimeDependentMaps.cpp rename to tests/Unit/Domain/Creators/TimeDependentOptions/Test_Sphere.cpp index ba20038c5eb9..e3b43a0cb3e2 100644 --- a/tests/Unit/Domain/Creators/Test_SphereTimeDependentMaps.cpp +++ b/tests/Unit/Domain/Creators/TimeDependentOptions/Test_Sphere.cpp @@ -11,8 +11,8 @@ #include #include "DataStructures/DataVector.hpp" -#include "Domain/Creators/ShapeMapOptions.hpp" -#include "Domain/Creators/SphereTimeDependentMaps.hpp" +#include "Domain/Creators/TimeDependentOptions/ShapeMap.hpp" +#include "Domain/Creators/TimeDependentOptions/Sphere.hpp" #include "Domain/FunctionsOfTime/FixedSpeedCubic.hpp" #include "Domain/FunctionsOfTime/FunctionOfTime.hpp" #include "Domain/FunctionsOfTime/PiecewisePolynomial.hpp" @@ -188,7 +188,7 @@ void test(const std::optional use_non_zero_shape) { } } // namespace -SPECTRE_TEST_CASE("Unit.Domain.Creators.SphereTimeDependentMaps", +SPECTRE_TEST_CASE("Unit.Domain.Creators.TimeDependentOptions.Sphere", "[Domain][Unit]") { test({true}); test({false}); diff --git a/tests/Unit/Domain/FunctionsOfTime/Test_PiecewisePolynomial.cpp b/tests/Unit/Domain/FunctionsOfTime/Test_PiecewisePolynomial.cpp index dea2c5d9a2a6..d32315b0e979 100644 --- a/tests/Unit/Domain/FunctionsOfTime/Test_PiecewisePolynomial.cpp +++ b/tests/Unit/Domain/FunctionsOfTime/Test_PiecewisePolynomial.cpp @@ -32,6 +32,16 @@ void test(const gsl::not_null f_of_t, *f_of_t_derived; CHECK(*f_of_t_derived == f_of_t_derived_copy); while (t < final_time) { + const auto lambdas_all = f_of_t->func_and_all_derivs(t); + CHECK(approx(lambdas_all[0][0]) == cube(t)); + CHECK(approx(lambdas_all[0][1]) == square(t)); + CHECK(approx(lambdas_all[1][0]) == 3.0 * square(t)); + CHECK(approx(lambdas_all[1][1]) == 2.0 * t); + CHECK(approx(lambdas_all[2][0]) == 6.0 * t); + CHECK(approx(lambdas_all[2][1]) == 2.0); + CHECK(approx(lambdas_all[3][0]) == 6.0); + CHECK(approx(lambdas_all[3][1]) == 0.0); + const auto lambdas0 = f_of_t->func_and_2_derivs(t); CHECK(approx(lambdas0[0][0]) == cube(t)); CHECK(approx(lambdas0[0][1]) == square(t)); @@ -79,6 +89,10 @@ void test_non_const_deriv( CHECK(*f_of_t_derived != f_of_t_derived_copy); } t *= 1.0 + std::numeric_limits::epsilon(); + const auto lambdas_all = f_of_t->func_and_all_derivs(t); + CHECK(approx(lambdas_all[0][0]) == 33.948); + CHECK(approx(lambdas_all[1][0]) == 19.56); + CHECK(approx(lambdas_all[2][0]) == 7.2); const auto lambdas0 = f_of_t->func_and_2_derivs(t); CHECK(approx(lambdas0[0][0]) == 33.948); CHECK(approx(lambdas0[1][0]) == 19.56); @@ -96,6 +110,15 @@ void test_non_const_deriv( template void test_within_roundoff(const FunctionsOfTime::FunctionOfTime& f_of_t) { + const auto lambdas_all = f_of_t.func_and_all_derivs(1.0 - 5.0e-16); + CHECK(approx(lambdas_all[0][0]) == 1.0); + CHECK(approx(lambdas_all[1][0]) == 3.0); + CHECK(approx(lambdas_all[2][0]) == 6.0); + CHECK(approx(lambdas_all[3][0]) == 6.0); + CHECK(approx(lambdas_all[0][1]) == 1.0); + CHECK(approx(lambdas_all[1][1]) == 2.0); + CHECK(approx(lambdas_all[2][1]) == 2.0); + CHECK(approx(lambdas_all[3][1]) == 0.0); const auto lambdas0 = f_of_t.func_and_2_derivs(1.0 - 5.0e-16); CHECK(approx(lambdas0[0][0]) == 1.0); CHECK(approx(lambdas0[1][0]) == 3.0); @@ -153,6 +176,14 @@ void check_func_and_derivs( CHECK(f_of_t.func(t) == q.func_and_derivs<0>(t)); CHECK(f_of_t.func_and_deriv(t) == q.func_and_derivs<1>(t)); CHECK(f_of_t.func_and_2_derivs(t) == q.func_and_derivs<2>(t)); + auto func_and_all_derivs = f_of_t.func_and_all_derivs(t); + auto func_and_2_derivs = q.func_and_derivs<2>(t); + CHECK(func_and_all_derivs.size() == MaxDeriv + 1); + for (size_t i = 0; + i < std::min(func_and_all_derivs.size(), func_and_2_derivs.size()); + i++) { + CHECK(func_and_all_derivs[i] == gsl::at(func_and_2_derivs, i)); + } test_copy_semantics(f_of_t); auto f_of_t_copy = f_of_t; diff --git a/tests/Unit/Domain/FunctionsOfTime/Test_QuaternionFunctionOfTime.cpp b/tests/Unit/Domain/FunctionsOfTime/Test_QuaternionFunctionOfTime.cpp index e009c7633030..9bbb9c5877fa 100644 --- a/tests/Unit/Domain/FunctionsOfTime/Test_QuaternionFunctionOfTime.cpp +++ b/tests/Unit/Domain/FunctionsOfTime/Test_QuaternionFunctionOfTime.cpp @@ -3,6 +3,8 @@ #include "Framework/TestingFramework.hpp" +#include +#include #include #include @@ -200,6 +202,12 @@ SPECTRE_TEST_CASE("Unit.Domain.FunctionsOfTime.QuaternionFunctionOfTime", CHECK_ITERABLE_APPROX(qfot_serialized_deserialized.func_and_deriv(2.0), qfot.func_and_deriv(2.0)); + CHECK_THROWS_WITH( + qfot.quat_func_and_3_derivs(2.0), + Catch::Matchers::ContainsSubstring( + "Need more angle derivs to compute the third derivative of the " + "quaternion. Currently only have 2")); + test_copy_semantics(qfot); test_move_semantics(std::move(qfot_serialized_deserialized), qfot); } @@ -392,7 +400,7 @@ SPECTRE_TEST_CASE("Unit.Domain.FunctionsOfTime.QuaternionFunctionOfTime", } CHECK(qfot.expiration_after(0.0) == time_step); - // Get the quaternion and 2 derivatives at a certain time. + // Get the quaternion and 3 derivatives at a certain time. double check_time = 5.398; const std::array quat_func_and_2_derivs = qfot.quat_func_and_2_derivs(check_time); @@ -413,6 +421,7 @@ SPECTRE_TEST_CASE("Unit.Domain.FunctionsOfTime.QuaternionFunctionOfTime", const double omega = 3 * fac1 * check_time * check_time + 2 * fac2 * check_time + fac3; const double dtomega = 6 * fac1 * check_time + 2 * fac2; + const double dt2omega = 6.0 * fac1; DataVector a_quat{{cos(0.5 * phi), 0.0, 0.0, sin(0.5 * phi)}}; DataVector a_dtquat{ @@ -420,6 +429,12 @@ SPECTRE_TEST_CASE("Unit.Domain.FunctionsOfTime.QuaternionFunctionOfTime", DataVector a_dt2quat{{-0.5 * (a_dtquat[3] * omega + a_quat[3] * dtomega), 0.0, 0.0, 0.5 * (a_dtquat[0] * omega + a_quat[0] * dtomega)}}; + DataVector a_dt3quat{ + {-0.5 * (a_dt2quat[3] * omega + 2.0 * a_dtquat[3] * dtomega + + a_quat[3] * dt2omega), + 0.0, 0.0, + 0.5 * (a_dt2quat[0] * omega + 2.0 * a_dtquat[0] * dtomega + + a_quat[0] * dt2omega)}}; // Compare analytic solution to numerical Approx custom_approx = Approx::custom().epsilon(5.0e-12).scale(1.0); @@ -438,6 +453,13 @@ SPECTRE_TEST_CASE("Unit.Domain.FunctionsOfTime.QuaternionFunctionOfTime", CHECK_ITERABLE_CUSTOM_APPROX(quat_func_and_2_derivs[2], a_dt2quat, custom_approx); } + { + INFO(" Compare third derivative of quaternion"); + const std::array quat_func_and_3_derivs = + qfot.quat_func_and_3_derivs(check_time); + CHECK_ITERABLE_CUSTOM_APPROX(quat_func_and_3_derivs[3], a_dt3quat, + custom_approx); + } } { INFO("QuaternionFunctionOfTime: No updates"); diff --git a/tests/Unit/Domain/Test_Domain.cpp b/tests/Unit/Domain/Test_Domain.cpp index abb1e048550c..a82241e13985 100644 --- a/tests/Unit/Domain/Test_Domain.cpp +++ b/tests/Unit/Domain/Test_Domain.cpp @@ -35,9 +35,9 @@ #include "Domain/CoordinateMaps/TimeDependent/Translation.hpp" #include "Domain/CoordinateMaps/Wedge.hpp" #include "Domain/Creators/BinaryCompactObject.hpp" -#include "Domain/Creators/BinaryCompactObjectHelpers.hpp" #include "Domain/Creators/ExpandOverBlocks.hpp" #include "Domain/Creators/RegisterDerivedWithCharm.hpp" +#include "Domain/Creators/TimeDependentOptions/BinaryCompactObject.hpp" #include "Domain/Domain.hpp" #include "Domain/DomainHelpers.hpp" #include "Domain/ExcisionSphere.hpp" diff --git a/tests/Unit/Evolution/Ringdown/Test_StrahlkorperCoefsInRingdownDistortedFrame.cpp b/tests/Unit/Evolution/Ringdown/Test_StrahlkorperCoefsInRingdownDistortedFrame.cpp index 7bd0f2dd13d9..cd3dcdd678ad 100644 --- a/tests/Unit/Evolution/Ringdown/Test_StrahlkorperCoefsInRingdownDistortedFrame.cpp +++ b/tests/Unit/Evolution/Ringdown/Test_StrahlkorperCoefsInRingdownDistortedFrame.cpp @@ -15,7 +15,7 @@ #include "DataStructures/DataVector.hpp" #include "Domain/CoordinateMaps/Distribution.hpp" #include "Domain/Creators/Sphere.hpp" -#include "Domain/Creators/SphereTimeDependentMaps.hpp" +#include "Domain/Creators/TimeDependentOptions/Sphere.hpp" #include "Domain/StrahlkorperTransformations.hpp" #include "Evolution/Ringdown/StrahlkorperCoefsInRingdownDistortedFrame.hpp" #include "Framework/TestHelpers.hpp" diff --git a/tests/Unit/ParallelAlgorithms/Interpolation/Test_ComputeDestVars.cpp b/tests/Unit/ParallelAlgorithms/Interpolation/Test_ComputeDestVars.cpp index 7efac1661819..eee09f7dc8ee 100644 --- a/tests/Unit/ParallelAlgorithms/Interpolation/Test_ComputeDestVars.cpp +++ b/tests/Unit/ParallelAlgorithms/Interpolation/Test_ComputeDestVars.cpp @@ -10,8 +10,8 @@ #include "DataStructures/Tensor/Tensor.hpp" #include "DataStructures/Variables.hpp" #include "Domain/Creators/Sphere.hpp" -#include "Domain/Creators/SphereTimeDependentMaps.hpp" #include "Domain/Creators/Tags/FunctionsOfTime.hpp" +#include "Domain/Creators/TimeDependentOptions/Sphere.hpp" #include "Domain/FunctionsOfTime/FunctionOfTime.hpp" #include "Domain/FunctionsOfTime/RegisterDerivedWithCharm.hpp" #include "Domain/Structure/ElementId.hpp"