diff --git a/src/Domain/Creators/TimeDependentOptions/BinaryCompactObject.cpp b/src/Domain/Creators/TimeDependentOptions/BinaryCompactObject.cpp index 9796519ba1e2..b05d60301b3b 100644 --- a/src/Domain/Creators/TimeDependentOptions/BinaryCompactObject.cpp +++ b/src/Domain/Creators/TimeDependentOptions/BinaryCompactObject.cpp @@ -16,12 +16,16 @@ #include "Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/ShapeMapTransitionFunction.hpp" #include "Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/SphereTransition.hpp" #include "Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/Wedge.hpp" +#include "Domain/Creators/TimeDependentOptions/ExpansionMap.hpp" +#include "Domain/Creators/TimeDependentOptions/RotationMap.hpp" #include "Domain/Creators/TimeDependentOptions/ShapeMap.hpp" +#include "Domain/Creators/TimeDependentOptions/TranslationMap.hpp" #include "Domain/FunctionsOfTime/FixedSpeedCubic.hpp" #include "Domain/FunctionsOfTime/FunctionOfTime.hpp" #include "Domain/FunctionsOfTime/IntegratedFunctionOfTime.hpp" #include "Domain/FunctionsOfTime/PiecewisePolynomial.hpp" #include "Domain/FunctionsOfTime/QuaternionFunctionOfTime.hpp" +#include "Domain/Structure/ObjectLabel.hpp" #include "NumericalAlgorithms/SphericalHarmonics/Spherepack.hpp" #include "Options/ParseError.hpp" #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrHorizon.hpp" @@ -34,19 +38,18 @@ namespace domain::creators::bco { template TimeDependentMapOptions::TimeDependentMapOptions( - double initial_time, - std::optional expansion_map_options, - std::optional rotation_map_options, - std::optional translation_map_options, - std::optional> shape_options_A, - std::optional> shape_options_B, + double initial_time, ExpansionMapOptionType expansion_map_options, + RotationMapOptionType rotation_map_options, + TranslationMapOptionType translation_map_options, + ShapeMapOptionType shape_options_A, + ShapeMapOptionType shape_options_B, const Options::Context& context) : initial_time_(initial_time), - expansion_map_options_(expansion_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) { + expansion_map_options_(std::move(expansion_map_options)), + rotation_map_options_(std::move(rotation_map_options)), + translation_map_options_(std::move(translation_map_options)), + shape_options_A_(std::move(shape_options_A)), + shape_options_B_(std::move(shape_options_B)) { if (not(expansion_map_options_.has_value() or rotation_map_options_.has_value() or translation_map_options_.has_value() or @@ -60,10 +63,15 @@ TimeDependentMapOptions::TimeDependentMapOptions( const auto check_l_max = [&context](const auto& shape_option, const domain::ObjectLabel label) { - if (shape_option.has_value() and shape_option.value().l_max <= 1) { - PARSE_ERROR(context, "Initial LMax for object " - << label << " must be 2 or greater but is " - << shape_option.value().l_max << " instead."); + if (shape_option.has_value() and + time_dependent_options::l_max_from_shape_options( + shape_option.value()) <= 1) { + PARSE_ERROR(context, + "Initial LMax for object " + << label << " must be 2 or greater but is " + << time_dependent_options::l_max_from_shape_options( + shape_option.value()) + << " instead."); } }; @@ -82,37 +90,42 @@ TimeDependentMapOptions::create_worldtube_functions_of_time() std::unordered_map> result{}; + // The functions of time need to be valid only for the very first time step, // after that they need to be updated by the worldtube singleton. const double initial_expiration_time = initial_time_ + 1e-10; if (not expansion_map_options_.has_value()) { ERROR("Initial values for the expansion map need to be provided."); } + + const auto& expansion_map_opts = + std::get>( + expansion_map_options_.value()); result[expansion_name] = std::make_unique( initial_time_, - std::array{ - {gsl::at(expansion_map_options_.value().initial_values, 0), - gsl::at(expansion_map_options_.value().initial_values, 1)}}, + std::array{expansion_map_opts.initial_values[0][0], + expansion_map_opts.initial_values[1][0]}, initial_expiration_time, false); result[expansion_outer_boundary_name] = std::make_unique( 1.0, initial_time_, - expansion_map_options_.value().outer_boundary_velocity, - expansion_map_options_.value().outer_boundary_decay_time); + expansion_map_opts.asymptotic_velocity_outer_boundary.value(), + expansion_map_opts.decay_timescale_outer_boundary); + if (not rotation_map_options_.has_value()) { ERROR( "Initial values for the rotation map need to be provided when using " "the worldtube."); } + const auto& rotation_map_opts = + std::get>( + rotation_map_options_.value()); result[rotation_name] = std::make_unique( initial_time_, - std::array{ - 0., - gsl::at(rotation_map_options_.value().initial_angular_velocity, - 2)}, + std::array{0., rotation_map_opts.angles[1][2]}, initial_expiration_time, true); // Size and Shape FunctionOfTime for objects A and B. Only spherical excision @@ -122,15 +135,19 @@ TimeDependentMapOptions::create_worldtube_functions_of_time() "Initial size for both excision spheres need to be provided when using " "the worldtube."); } + const auto& shape_opts_A = std::get>(shape_options_A_.value()); + const auto& shape_opts_B = std::get>(shape_options_B_.value()); for (size_t i = 0; i < shape_names.size(); i++) { const auto make_initial_size_values = [](const auto& lambda_options) { return std::array{ - {gsl::at(lambda_options.value().initial_size_values.value(), 0), - gsl::at(lambda_options.value().initial_size_values.value(), 1)}}; + {gsl::at(lambda_options.initial_size_values.value(), 0), + gsl::at(lambda_options.initial_size_values.value(), 1)}}; }; const std::array initial_size_values = - i == 0 ? make_initial_size_values(shape_options_A_) - : make_initial_size_values(shape_options_B_); + i == 0 ? make_initial_size_values(shape_opts_A) + : make_initial_size_values(shape_opts_B); const size_t initial_l_max = 2; const DataVector shape_zeros{ ylm::Spherepack::spectral_size(initial_l_max, initial_l_max), 0.0}; @@ -186,96 +203,45 @@ TimeDependentMapOptions::create_functions_of_time( expiration_times[name] = expr_time; } - // ExpansionMap FunctionOfTime for the function \f$a(t)\f$ in the + // ExpansionMap FunctionOfTime for the function a(t) and b(t) in the // domain::CoordinateMaps::TimeDependent::RotScaleTrans map if (expansion_map_options_.has_value()) { - result[expansion_name] = - std::make_unique>( - initial_time_, - std::array{ - {{gsl::at(expansion_map_options_.value().initial_values, 0)}, - {gsl::at(expansion_map_options_.value().initial_values, 1)}, - {0.0}}}, - expiration_times.at(expansion_name)); - - // ExpansionMap FunctionOfTime for the function \f$b(t)\f$ in the - // domain::CoordinateMaps::TimeDependent::RotScaleTrans map - result[expansion_outer_boundary_name] = - std::make_unique( - 1.0, initial_time_, - expansion_map_options_.value().outer_boundary_velocity, - expansion_map_options_.value().outer_boundary_decay_time); + auto expansion_functions_of_time = time_dependent_options::get_expansion( + expansion_map_options_.value(), initial_time_, + expiration_times.at(expansion_name)); + + result.merge(expansion_functions_of_time); } // RotationMap FunctionOfTime for the rotation angles about each - // axis. The initial rotation angles don't matter as we never - // actually use the angles themselves. We only use their derivatives - // (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. + // axis. 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_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}}}, + result[rotation_name] = time_dependent_options::get_rotation( + rotation_map_options_.value(), initial_time_, expiration_times.at(rotation_name)); } // TranslationMap FunctionOfTime 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]}}}, + result[translation_name] = time_dependent_options::get_translation( + translation_map_options_.value(), initial_time_, expiration_times.at(translation_name)); } // Size and Shape FunctionOfTime for objects A and B - const auto build_shape_and_size_fot = - [&result, &expiration_times, this]( - const auto& shape_options, const double inner_radius, - const std::string& shape_name, const std::string& size_name) { - auto [shape_funcs, size_funcs] = - time_dependent_options::initial_shape_and_size_funcs(shape_options, - inner_radius); - - result[shape_name] = - std::make_unique>( - initial_time_, std::move(shape_funcs), - expiration_times.at(shape_name)); - result[size_name] = - std::make_unique>( - initial_time_, std::move(size_funcs), - expiration_times.at(size_name)); - }; - if (shape_options_A_.has_value()) { if (not deformed_radii_[0].has_value()) { ERROR( "A shape map was specified for object A, but no inner radius is " "available. The object must be enclosed by a sphere."); } - build_shape_and_size_fot(shape_options_A_.value(), *deformed_radii_[0], - shape_names[0], size_names[0]); + + auto shape_and_size = time_dependent_options::get_shape_and_size( + shape_options_A_.value(), initial_time_, + expiration_times.at(shape_names[0]), expiration_times.at(size_names[0]), + *deformed_radii_[0]); + + result.merge(shape_and_size); } if (shape_options_B_.has_value()) { if (not deformed_radii_[1].has_value()) { @@ -283,8 +249,13 @@ TimeDependentMapOptions::create_functions_of_time( "A shape map was specified for object B, but no inner radius is " "available. The object must be enclosed by a sphere."); } - build_shape_and_size_fot(shape_options_B_.value(), *deformed_radii_[1], - shape_names[1], size_names[1]); + + auto shape_and_size = time_dependent_options::get_shape_and_size( + shape_options_B_.value(), initial_time_, + expiration_times.at(shape_names[1]), expiration_times.at(size_names[1]), + *deformed_radii_[1]); + + result.merge(shape_and_size); } return result; @@ -354,8 +325,11 @@ void TimeDependentMapOptions::build_maps( // Store the inner radii for creating functions of time gsl::at(deformed_radii_, i) = filled ? radii[1] : radii[0]; - const size_t initial_l_max = i == 0 ? shape_options_A_.value().l_max - : shape_options_B_.value().l_max; + const size_t initial_l_max = + i == 0 ? time_dependent_options::l_max_from_shape_options( + shape_options_A_.value()) + : time_dependent_options::l_max_from_shape_options( + shape_options_B_.value()); std::unique_ptr @@ -385,8 +359,12 @@ void TimeDependentMapOptions::build_maps( const std::array axes{3, -3, 2, -2, 1, -1}; const bool transition_ends_at_cube = - i == 0 ? shape_options_A_->transition_ends_at_cube - : shape_options_B_->transition_ends_at_cube; + i == 0 ? time_dependent_options:: + transition_ends_at_cube_from_shape_options( + shape_options_A_.value()) + : time_dependent_options:: + transition_ends_at_cube_from_shape_options( + shape_options_B_.value()); // These centers must take in to account if we have an offset of the // center of the object and where the transition ends. The inner center @@ -469,8 +447,16 @@ TimeDependentMapOptions::distorted_to_inertial_map( } else { const bool transition_ends_at_cube = Object == domain::ObjectLabel::A - ? shape_options_A_->transition_ends_at_cube - : shape_options_B_->transition_ends_at_cube; + ? (shape_options_A_.has_value() + ? time_dependent_options:: + transition_ends_at_cube_from_shape_options( + shape_options_A_.value()) + : false) + : (shape_options_B_.has_value() + ? time_dependent_options:: + transition_ends_at_cube_from_shape_options( + shape_options_B_.value()) + : false); block_has_shape_map = include_distorted_map.has_value() and (transition_ends_at_cube or include_distorted_map.value() < 6); @@ -505,11 +491,15 @@ TimeDependentMapOptions::grid_to_distorted_map( if constexpr (IsCylindrical) { block_has_shape_map = block_has_shape_map and include_distorted_map; - } else { + } else if (block_has_shape_map) { const bool transition_ends_at_cube = Object == domain::ObjectLabel::A - ? shape_options_A_->transition_ends_at_cube - : shape_options_B_->transition_ends_at_cube; + ? time_dependent_options:: + transition_ends_at_cube_from_shape_options( + shape_options_A_.value()) + : time_dependent_options:: + transition_ends_at_cube_from_shape_options( + shape_options_B_.value()); block_has_shape_map = block_has_shape_map and include_distorted_map.has_value() and (transition_ends_at_cube or include_distorted_map.value() < 6); @@ -550,11 +540,15 @@ TimeDependentMapOptions::grid_to_inertial_map( if constexpr (IsCylindrical) { block_has_shape_map = block_has_shape_map and include_distorted_map; - } else { + } else if (block_has_shape_map) { const bool transition_ends_at_cube = Object == domain::ObjectLabel::A - ? shape_options_A_->transition_ends_at_cube - : shape_options_B_->transition_ends_at_cube; + ? time_dependent_options:: + transition_ends_at_cube_from_shape_options( + shape_options_A_.value()) + : time_dependent_options:: + transition_ends_at_cube_from_shape_options( + shape_options_B_.value()); block_has_shape_map = block_has_shape_map and include_distorted_map.has_value() and (transition_ends_at_cube or include_distorted_map.value() < 6); diff --git a/src/Domain/Creators/TimeDependentOptions/BinaryCompactObject.hpp b/src/Domain/Creators/TimeDependentOptions/BinaryCompactObject.hpp index 342261648ee2..72c15874e066 100644 --- a/src/Domain/Creators/TimeDependentOptions/BinaryCompactObject.hpp +++ b/src/Domain/Creators/TimeDependentOptions/BinaryCompactObject.hpp @@ -17,8 +17,10 @@ #include "Domain/CoordinateMaps/TimeDependent/Rotation.hpp" #include "Domain/CoordinateMaps/TimeDependent/Shape.hpp" #include "Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/ShapeMapTransitionFunction.hpp" +#include "Domain/Creators/TimeDependentOptions/ExpansionMap.hpp" +#include "Domain/Creators/TimeDependentOptions/RotationMap.hpp" #include "Domain/Creators/TimeDependentOptions/ShapeMap.hpp" -#include "Domain/Creators/TimeDependentOptions/Sphere.hpp" +#include "Domain/Creators/TimeDependentOptions/TranslationMap.hpp" #include "Domain/FunctionsOfTime/FunctionOfTime.hpp" #include "Domain/Structure/ObjectLabel.hpp" #include "Options/Auto.hpp" @@ -146,105 +148,30 @@ struct TimeDependentMapOptions { /// The outer boundary radius of the map is always set to /// the outer boundary of the Domain, so there is no option /// here to set the outer boundary radius. - struct ExpansionMapOptions { - using type = Options::Auto; - static std::string name() { return "ExpansionMap"; } - static constexpr Options::String help = { - "Options for the expansion map. Specify 'None' to not use this map."}; - struct InitialValues { - using type = std::array; - static constexpr Options::String help = { - "Initial value and deriv of expansion."}; - }; - struct AsymptoticVelocityOuterBoundary { - using type = double; - static constexpr Options::String help = { - "The asymptotic velocity of the outer boundary."}; - }; - struct DecayTimescaleOuterBoundaryVelocity { - using type = double; - static constexpr Options::String help = { - "The timescale for how fast the outer boundary velocity approaches " - "its asymptotic value."}; - }; - using options = tmpl::list; - ExpansionMapOptions() = default; - ExpansionMapOptions(std::array initial_values_in, - double outer_boundary_velocity_in, - double outer_boundary_decay_time_in) - : initial_values(initial_values_in), - outer_boundary_velocity(outer_boundary_velocity_in), - outer_boundary_decay_time(outer_boundary_decay_time_in) {} - - std::array initial_values{ - std::numeric_limits::signaling_NaN(), - std::numeric_limits::signaling_NaN()}; - double outer_boundary_velocity{ - std::numeric_limits::signaling_NaN()}; - double outer_boundary_decay_time{ - std::numeric_limits::signaling_NaN()}; - }; - - struct RotationMapOptions { - using type = Options::Auto; - static std::string name() { return "RotationMap"; } - static constexpr Options::String help = { - "Options for a time-dependent rotation map about an arbitrary axis. " - "Specify 'None' to not use this map."}; - - struct InitialAngularVelocity { - using type = std::array; - static constexpr Options::String help = {"The initial angular velocity."}; - }; + using ExpansionMapOptions = + domain::creators::time_dependent_options::ExpansionMapOptions; + using ExpansionMapOptionType = typename ExpansionMapOptions::type::value_type; - using options = tmpl::list; - - RotationMapOptions() = default; - explicit RotationMapOptions( - std::array initial_angular_velocity_in) - : initial_angular_velocity(initial_angular_velocity_in) {} - - std::array initial_angular_velocity{}; - }; + /// \brief Options for the rotation map + using RotationMapOptions = + domain::creators::time_dependent_options::RotationMapOptions; + using RotationMapOptionType = typename RotationMapOptions::type::value_type; /// \brief Options for the Translation Map, the outer radius is always set to /// the outer boundary of the Domain, so there's no option needed for outer /// boundary. - struct TranslationMapOptions { - using type = Options::Auto; - static std::string name() { return "TranslationMap"; } - static constexpr Options::String help = { - "Options for a time-dependent translation map. Specify 'None' to not " - "use this map."}; - - struct InitialValues { - using type = std::array, 3>; - static constexpr Options::String help = { - "Initial position, velocity and acceleration."}; - }; - - using options = tmpl::list; - TranslationMapOptions() = default; - explicit TranslationMapOptions( - std::array, 3> initial_values_in) - : initial_values(initial_values_in) {} + using TranslationMapOptions = + domain::creators::time_dependent_options::TranslationMapOptions<3>; + using TranslationMapOptionType = + typename TranslationMapOptions::type::value_type; - std::array, 3> initial_values{}; - }; - - // We use a type alias here instead of defining the ShapeMapOptions struct - // because there appears to be a bug in clang-10. If the definition of - // ShapeMapOptions is here inside TimeDependentMapOptions, on clang-10 there - // is a linking error that there is an undefined reference to - // Options::Option::parse_as> (and B). This doesn't - // show up for GCC. If we put the definition of ShapeMapOptions outside of - // TimeDependentMapOptions and just use a type alias here, the linking error - // goes away. + /// \brief Options for the shape map template using ShapeMapOptions = domain::creators::time_dependent_options::ShapeMapOptions< not IsCylindrical, Object>; + template + using ShapeMapOptionType = typename ShapeMapOptions::type::value_type; using options = tmpl::list expansion_map_options, - std::optional rotation_map_options, - std::optional translation_map_options, - std::optional> shape_options_A, - std::optional> shape_options_B, + double initial_time, ExpansionMapOptionType expansion_map_options, + RotationMapOptionType rotation_map_options, + TranslationMapOptionType translation_map_options, + ShapeMapOptionType shape_options_A, + ShapeMapOptionType shape_options_B, const Options::Context& context = {}); /*! @@ -407,11 +333,11 @@ struct TimeDependentMapOptions { static size_t get_index(domain::ObjectLabel object); double initial_time_{std::numeric_limits::signaling_NaN()}; - std::optional expansion_map_options_{}; - std::optional rotation_map_options_{}; - std::optional translation_map_options_{}; - std::optional> shape_options_A_{}; - std::optional> shape_options_B_{}; + ExpansionMapOptionType expansion_map_options_; + RotationMapOptionType rotation_map_options_; + TranslationMapOptionType translation_map_options_; + ShapeMapOptionType shape_options_A_{}; + ShapeMapOptionType shape_options_B_{}; std::array, 2> deformed_radii_{}; // Maps diff --git a/src/Domain/Creators/TimeDependentOptions/ExpansionMap.cpp b/src/Domain/Creators/TimeDependentOptions/ExpansionMap.cpp index 9b72e5eb321a..01da212fb929 100644 --- a/src/Domain/Creators/TimeDependentOptions/ExpansionMap.cpp +++ b/src/Domain/Creators/TimeDependentOptions/ExpansionMap.cpp @@ -4,87 +4,171 @@ #include "Domain/Creators/TimeDependentOptions/ExpansionMap.hpp" #include +#include #include #include #include "DataStructures/DataVector.hpp" #include "Domain/Creators/TimeDependentOptions/FromVolumeFile.hpp" +#include "Domain/FunctionsOfTime/FixedSpeedCubic.hpp" +#include "Domain/FunctionsOfTime/FunctionOfTime.hpp" +#include "Domain/FunctionsOfTime/PiecewisePolynomial.hpp" +#include "Domain/FunctionsOfTime/SettleToConstant.hpp" #include "Options/Context.hpp" #include "Options/ParseError.hpp" #include "Utilities/GenerateInstantiations.hpp" #include "Utilities/Gsl.hpp" +#include "Utilities/MakeArray.hpp" namespace domain::creators::time_dependent_options { -ExpansionMapOptions::ExpansionMapOptions( - const std::variant, FromVolumeFile>& - expansion_values, - const std::variant, FromVolumeFile>& - expansion_outer_boundary_values, - std::optional decay_timescale_outer_boundary_in, - std::optional decay_timescale_in, - std::optional asymptotic_velocity_outer_boundary_in, - const Options::Context& context) - : decay_timescale(std::move(decay_timescale_in)) { - const auto set_values = - [&, this](const gsl::not_null*> to_set, - const auto& input_values, const bool is_outer_boundary) { - if (std::holds_alternative>(input_values)) { - asymptotic_velocity_outer_boundary = - asymptotic_velocity_outer_boundary_in; - if (decay_timescale.has_value() == - asymptotic_velocity_outer_boundary.has_value()) { - PARSE_ERROR( - context, - "When specifying the ExpansionMap initial outer boundary " - "values directly, you must specify one of DecayTimescale or " - "AsymptoticVelocityOuterBoundary, but not both."); - } - auto& values = std::get>(input_values); - for (size_t i = 0; i < to_set->size(); i++) { - gsl::at(*to_set, i) = DataVector{1, gsl::at(values, i)}; - } - - if (is_outer_boundary) { - if (not decay_timescale_outer_boundary_in.has_value()) { - PARSE_ERROR(context, - "When specifying the ExpansionMap initial outer " - "boundary values directly, you must also specify a " - "'DecayTimescaleOuterBoundary'."); - } - - decay_timescale_outer_boundary = - decay_timescale_outer_boundary_in.value(); - } - - } else if (std::holds_alternative>( - input_values)) { - if (decay_timescale.has_value()) { - PARSE_ERROR( - context, - "When specifying the initial values from a volume file, " - "the decay timescale must be 'Auto'."); - } - auto& values_from_file = - std::get>(input_values); - *to_set = is_outer_boundary - ? values_from_file.expansion_values_outer_boundary - : values_from_file.expansion_values; - - if (is_outer_boundary) { - decay_timescale_outer_boundary = - decay_timescale_outer_boundary_in.value_or( - values_from_file.decay_timescale_outer_boundary); - asymptotic_velocity_outer_boundary = - asymptotic_velocity_outer_boundary_in.value_or( - values_from_file.velocity_outer_boundary); - } - } - }; - - // Expansion values - set_values(make_not_null(&initial_values), expansion_values, false); - // Outer boundary - set_values(make_not_null(&initial_values_outer_boundary), - expansion_outer_boundary_values, true); +#if defined(__GNUC__) && !defined(__clang__) +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wsuggest-attribute=noreturn" +#endif // defined(__GNUC__) && !defined(__clang__) +template +ExpansionMapOptions::ExpansionMapOptions( + const std::array& initial_values_in, + double decay_timescale_outer_boundary_in, + const std::array& initial_values_outer_boundary_in, + double decay_timescale_in, const Options::Context& context) + : decay_timescale_outer_boundary(decay_timescale_outer_boundary_in), + decay_timescale(decay_timescale_in) { + if constexpr (AllowSettleFoTs) { + initial_values = std::array{DataVector{initial_values_in[0]}, + DataVector{initial_values_in[1]}, + DataVector{initial_values_in[2]}}; + initial_values_outer_boundary = + std::array{DataVector{initial_values_outer_boundary_in[0]}, + DataVector{initial_values_outer_boundary_in[1]}, + DataVector{initial_values_outer_boundary_in[2]}}; + } else { + PARSE_ERROR(context, + "This class does not allow SettleToConst functions of time, " + "but the constructor for allowing SettleToConst functions of " + "time was used. Please use the other constructor."); + } } +#if defined(__GNUC__) && !defined(__clang__) +#pragma GCC diagnostic pop +#endif // defined(__GNUC__) && !defined(__clang__) + +template +ExpansionMapOptions::ExpansionMapOptions( + const std::array& initial_values_in, + double decay_timescale_outer_boundary_in, + double asymptotic_velocity_outer_boundary_in, + const Options::Context& /*context*/) + : decay_timescale_outer_boundary(decay_timescale_outer_boundary_in), + asymptotic_velocity_outer_boundary( + asymptotic_velocity_outer_boundary_in) { + initial_values = std::array{DataVector{initial_values_in[0]}, + DataVector{initial_values_in[1]}, + DataVector{initial_values_in[2]}}; + initial_values_outer_boundary = + std::array{DataVector{1.0}, DataVector{0.0}, DataVector{0.0}}; +} + +template +FunctionsOfTimeMap get_expansion( + const std::variant, FromVolumeFile>& + expansion_map_options, + const double initial_time, const double expiration_time) { + const std::string name{"Expansion"}; + const std::string name_outer_boundary{"ExpansionOuterBoundary"}; + FunctionsOfTimeMap result{}; + + if (std::holds_alternative(expansion_map_options)) { + const auto& from_vol_file = std::get(expansion_map_options); + const auto volume_fot = from_vol_file.retrieve_function_of_time( + {name, name_outer_boundary}, initial_time); + + // Expansion must be either a PiecewisePolynomial or a SettleToConstant + const auto* pp_volume_fot = + dynamic_cast*>( + volume_fot.at(name).get()); + const auto* settle_volume_fot = + dynamic_cast( + volume_fot.at(name).get()); + + if (UNLIKELY(pp_volume_fot == nullptr and settle_volume_fot == nullptr)) { + ERROR_NO_TRACE( + "Expansion function of time read from volume data is not a " + "PiecewisePolynomial<2> or a SettleToConstant. Cannot use it to " + "initialize the expansion map."); + } + + result[name] = + volume_fot.at(name)->create_at_time(initial_time, expiration_time); + + // Outer boundary must be either a FixedSpeedCubic or a SettleToConstant + const auto* outer_boundary_cubic_volume_fot = + dynamic_cast( + volume_fot.at(name_outer_boundary).get()); + const auto* outer_boundary_settle_volume_fot = + dynamic_cast( + volume_fot.at(name_outer_boundary).get()); + + if (UNLIKELY(outer_boundary_cubic_volume_fot == nullptr and + outer_boundary_settle_volume_fot == nullptr)) { + ERROR_NO_TRACE( + "ExpansionOuterBoundary function of time read from volume data is " + "not a FixedSpeedCubic or a SettleToConstant. Cannot use it to " + "initialize the expansion map."); + } + + ASSERT(outer_boundary_cubic_volume_fot != nullptr or AllowSettleFoTs, + "ExpansionOuterBoundary function of time in the volume file is a " + "SettleToConstant, but SettleToConstant functions of time aren't " + "allowed."); + + result[name_outer_boundary] = + volume_fot.at(name_outer_boundary)->get_clone(); + } else if (std::holds_alternative>( + expansion_map_options)) { + const auto& hard_coded_options = + std::get>(expansion_map_options); + + if (hard_coded_options.asymptotic_velocity_outer_boundary.has_value()) { + result[name] = + std::make_unique>( + initial_time, hard_coded_options.initial_values, expiration_time); + result[name_outer_boundary] = + std::make_unique( + hard_coded_options.initial_values_outer_boundary[0][0], + initial_time, + hard_coded_options.asymptotic_velocity_outer_boundary.value(), + hard_coded_options.decay_timescale_outer_boundary); + } else { + ASSERT(hard_coded_options.decay_timescale.has_value(), + "To construct an ExpansionMap SettleToConstant function of time, " + "a decay timescale must be supplied."); + result[name] = + std::make_unique( + hard_coded_options.initial_values, initial_time, + hard_coded_options.decay_timescale.value()); + result[name_outer_boundary] = + std::make_unique( + hard_coded_options.initial_values_outer_boundary, initial_time, + hard_coded_options.decay_timescale_outer_boundary); + } + } else { + ERROR("Unknown ExpansionMap."); + } + + return result; +} + +#define ALLOWSETTLE(data) BOOST_PP_TUPLE_ELEM(0, data) + +#define INSTANTIATE(_, data) \ + template struct ExpansionMapOptions; \ + template FunctionsOfTimeMap get_expansion( \ + const std::variant, \ + FromVolumeFile>& expansion_map_options, \ + double initial_time, double expiration_time); + +GENERATE_INSTANTIATIONS(INSTANTIATE, (true, false)) + +#undef INSTANTIATE +#undef ALLOWSETTLE } // namespace domain::creators::time_dependent_options diff --git a/src/Domain/Creators/TimeDependentOptions/ExpansionMap.hpp b/src/Domain/Creators/TimeDependentOptions/ExpansionMap.hpp index b02ceb964372..9587ed2d2235 100644 --- a/src/Domain/Creators/TimeDependentOptions/ExpansionMap.hpp +++ b/src/Domain/Creators/TimeDependentOptions/ExpansionMap.hpp @@ -5,104 +5,113 @@ #include #include +#include #include #include #include #include "DataStructures/DataVector.hpp" #include "Domain/Creators/TimeDependentOptions/FromVolumeFile.hpp" +#include "Domain/FunctionsOfTime/FunctionOfTime.hpp" #include "Options/Auto.hpp" #include "Options/Context.hpp" +#include "Options/Options.hpp" #include "Options/String.hpp" #include "Utilities/TMPL.hpp" namespace domain::creators::time_dependent_options { /*! - * \brief Class to be used as an option for initializing expansion map - * coefficients. + * \brief Class that holds hard coded expansion map options from the input file. + * + * \details This class can also be used as an option tag with the \p type type + * alias, `name()` function, and \p help string. */ +template struct ExpansionMapOptions { - using type = Options::Auto; + using type = Options::Auto< + std::variant, FromVolumeFile>, + Options::AutoLabel::None>; static std::string name() { return "ExpansionMap"; } static constexpr Options::String help = { "Options for a time-dependent expansion of the coordinates. Specify " "'None' to not use this map."}; struct InitialValues { - using type = - std::variant, FromVolumeFile>; - static constexpr Options::String help = { - "Initial values for the expansion map, its velocity and " - "acceleration."}; + using type = std::array; + static constexpr Options::String help = + "Initial values for the expansion map, its velocity and acceleration."; }; struct InitialValuesOuterBoundary { - using type = - std::variant, FromVolumeFile>; + using type = std::array; static constexpr Options::String help = { - "Initial values for the expansion map, its velocity and " - "acceleration at the outer boundary. Unless you are starting from a " - "checkpoint or continuing an evolution, this option should likely be " - "[1.0, 0.0, 0.0] at the start of an evolution"}; + "Initial values for the expansion map, its velocity and acceleration " + "at the outer boundary."}; }; struct DecayTimescaleOuterBoundary { - using type = Options::Auto; + using type = double; static constexpr Options::String help = { "A timescale for how fast the outer boundary expansion approaches its " - "asymptotic value. Can optionally specify 'Auto' when reading the " - "initial values 'FromVolumeFile' to use the decay timescale from the " - "function of time in the volume file. Cannot specify 'Auto' when " - "initial values are specified directly."}; + "asymptotic value."}; }; struct DecayTimescale { - using type = Options::Auto; + using type = double; static constexpr Options::String help = { - "If specified, a SettleToConstant function of time will be used for " - "the expansion map and this number will determine the timescale that " - "the expansion approaches its asymptotic value. If 'Auto' is " - "specified, a PiecewisePolynomial function of time will be used for " - "the expansion map. Note that if you are reading the initial values " - "from a volume file, you must specify 'Auto' for this option."}; + "A timescale for how fast the expansion approaches its asymptotic " + "value with a SettleToConstant function of time."}; }; struct AsymptoticVelocityOuterBoundary { - using type = Options::Auto; + using type = double; static constexpr Options::String help = { - "There are two choices for this option. If a value is specified, a " - "FixedSpeedCubic function of time will be used for the expansion map " - "at the outer boundary and this number will determine its velocity. If " - "'Auto' is specified, the behavior will depend on what is chosen for " - "'InitialValuesOuterBoundary'. If values are specified for " - "'InitialValuesOuterBoundary', then 'Auto' here means a " - "SettleToConstant function of time will be used for the expansion map " - "at the outer boundary. If 'FromVolumeFile' is specified for " - "'InitialValuesOuterBoundary', then a FixedSpeedCubic function of time " - "will be used and the velocity from the function of " - "time in the volume file will be used."}; + "The constant velocity of the outer boundary expansion."}; }; - using options = tmpl::list; + using common_options = tmpl::list; + using settle_options = + tmpl::push_back; + using non_settle_options = + tmpl::push_back; + + using options = tmpl::conditional_t< + AllowSettleFoTs, + tmpl::list>, + non_settle_options>; ExpansionMapOptions() = default; + // Constructor for SettleToConstant functions of time ExpansionMapOptions( - const std::variant, - FromVolumeFile>& expansion_values, - const std::variant, - FromVolumeFile>& - expansion_outer_boundary_values, - std::optional decay_timescale_outer_boundary_in, - std::optional decay_timescale_in, - std::optional asymptotic_velocity_outer_boundary_in, - const Options::Context& context = {}); + const std::array& initial_values_in, + double decay_timescale_outer_boundary_in, + const std::array& initial_values_outer_boundary_in, + double decay_timescale_in, const Options::Context& context = {}); + // Constructor for non SettleToConstant functions of time + ExpansionMapOptions(const std::array& initial_values_in, + double decay_timescale_outer_boundary_in, + double asymptotic_velocity_outer_boundary_in, + const Options::Context& context = {}); std::array initial_values{}; std::array initial_values_outer_boundary{}; double decay_timescale_outer_boundary{}; - std::optional decay_timescale{}; - std::optional asymptotic_velocity_outer_boundary{}; + std::optional decay_timescale; + std::optional asymptotic_velocity_outer_boundary; }; + +/*! + * \brief Helper functions that take the variant of the expansion map options, + * and return the fully constructed expansion functions of time. + * + * \details Even if the functions of time are read from a file, they will have a + * new \p initial_time and \p expiration_time (no expiration time for the outer + * boundary function of time though). + */ +template +FunctionsOfTimeMap get_expansion( + const std::variant, FromVolumeFile>& + expansion_map_options, + double initial_time, double expiration_time); } // namespace domain::creators::time_dependent_options diff --git a/src/Domain/Creators/TimeDependentOptions/FromVolumeFile.cpp b/src/Domain/Creators/TimeDependentOptions/FromVolumeFile.cpp index cee1142763ec..471bee911811 100644 --- a/src/Domain/Creators/TimeDependentOptions/FromVolumeFile.cpp +++ b/src/Domain/Creators/TimeDependentOptions/FromVolumeFile.cpp @@ -7,6 +7,7 @@ #include #include #include +#include #include #include #include @@ -23,36 +24,34 @@ #include "Options/Context.hpp" #include "Options/ParseError.hpp" #include "Utilities/Algorithm.hpp" +#include "Utilities/ErrorHandling/Error.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 std::string& subfile_name, const std::optional& time) { 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); + ERROR_NO_TRACE(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."); + ERROR_NO_TRACE(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 get_function_of_time( 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); + ERROR_NO_TRACE("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]) { + if (time.has_value() and + (time.value() < time_bounds[0] or time.value() > 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); + ERROR_NO_TRACE(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); - - const auto* fixed_speed_cubic = - dynamic_cast( - exp_outer_boundary_function_of_time.get()); - - if (fixed_speed_cubic == nullptr) { - PARSE_ERROR( - context, - "When reading the Expansion map parameters from a volume file, the " - "ExpansionOuterBoundary function of time must be a FixedSpeedCubic."); +FromVolumeFile::FromVolumeFile(std::string h5_filename, + std::string subfile_name) + : h5_filename_(std::move(h5_filename)), + subfile_name_(std::move(subfile_name)) {} + +FunctionsOfTimeMap FromVolumeFile::retrieve_function_of_time( + const std::unordered_set& function_of_time_names, + const std::optional& time) const { + FunctionsOfTimeMap result{}; + for (const std::string& name : function_of_time_names) { + result[name] = + get_function_of_time(name, h5_filename_, subfile_name_, time); } - velocity_outer_boundary = fixed_speed_cubic->velocity(); - decay_timescale_outer_boundary = fixed_speed_cubic->decay_timescale(); + return result; } - -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 index 7eb090387be0..f6d1d4eef8e4 100644 --- a/src/Domain/Creators/TimeDependentOptions/FromVolumeFile.hpp +++ b/src/Domain/Creators/TimeDependentOptions/FromVolumeFile.hpp @@ -4,9 +4,12 @@ #pragma once #include +#include +#include #include #include "DataStructures/DataVector.hpp" +#include "Domain/FunctionsOfTime/FunctionOfTime.hpp" #include "Domain/Structure/ObjectLabel.hpp" #include "Options/Auto.hpp" #include "Options/Context.hpp" @@ -14,20 +17,15 @@ #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. + * \brief Read in FunctionOfTime coefficients from an H5 file and volume + * subfile. + * + * \details The H5 file will only be accessed in the `retrieve_function_of_time` + * function. */ -namespace names { -struct Translation {}; -struct Rotation {}; -struct Expansion {}; -template -struct ShapeSize {}; -} // namespace names - -namespace detail { -struct FromVolumeFileBase { +struct FromVolumeFile { struct H5Filename { using type = std::string; static constexpr Options::String help{ @@ -41,84 +39,28 @@ struct FromVolumeFileBase { "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; + 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{}; - double velocity_outer_boundary{}; - double decay_timescale_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{}; + FromVolumeFile(std::string h5_filename, std::string subfile_name); + + /*! + * \brief Searches the last observation in the volume subfile and returns + * clones of all functions of time in \p function_of_time_names. + * + * \details If a value for \p time is specified, will ensure that \p time is + * within the `domain::FunctionsOfTime::FunctionOfTime::time_bounds()` of each + * function of time. + */ + FunctionsOfTimeMap retrieve_function_of_time( + const std::unordered_set& function_of_time_names, + const std::optional& time) const; + + private: + std::string h5_filename_; + std::string subfile_name_; }; /// @} } // namespace domain::creators::time_dependent_options diff --git a/src/Domain/Creators/TimeDependentOptions/RotationMap.cpp b/src/Domain/Creators/TimeDependentOptions/RotationMap.cpp index ed1d6e12ea53..d7051194696a 100644 --- a/src/Domain/Creators/TimeDependentOptions/RotationMap.cpp +++ b/src/Domain/Creators/TimeDependentOptions/RotationMap.cpp @@ -4,6 +4,7 @@ #include "Domain/Creators/TimeDependentOptions/RotationMap.hpp" #include +#include #include #include #include @@ -11,6 +12,8 @@ #include "DataStructures/DataVector.hpp" #include "Domain/Creators/TimeDependentOptions/FromVolumeFile.hpp" +#include "Domain/FunctionsOfTime/QuaternionFunctionOfTime.hpp" +#include "Domain/FunctionsOfTime/SettleToConstantQuaternion.hpp" #include "Options/Context.hpp" #include "Options/ParseError.hpp" #include "Utilities/ErrorHandling/Error.hpp" @@ -18,88 +21,104 @@ #include "Utilities/MakeArray.hpp" namespace domain::creators::time_dependent_options { -template -RotationMapOptions::RotationMapOptions( - std::variant>, - FromVolumeFile> - initial_quaternions, - std::optional>> initial_angles, - std::optional decay_timescale_in, const Options::Context& context) +template +void RotationMapOptions::initialize_angles_and_quats() { + quaternions = make_array<3, DataVector>(DataVector{4, 0.0}); + // Defautl to the identity quaternion + quaternions[0][0] = 1.0; + angles = make_array<4, DataVector>(DataVector{3, 0.0}); +} + +template +RotationMapOptions::RotationMapOptions( + const std::array& initial_angular_velocity, + const Options::Context& /*context*/) { + initialize_angles_and_quats(); + + angles[1] = DataVector{initial_angular_velocity}; +} + +template +RotationMapOptions::RotationMapOptions( + const std::vector>& initial_quaternions, + const double decay_timescale_in, const Options::Context& context) : decay_timescale(decay_timescale_in) { - quaternions = make_array(DataVector{4, 0.0}); - angles = make_array(DataVector{3, 0.0}); - - if (std::holds_alternative>>( - initial_quaternions)) { - auto& values = - std::get>>(initial_quaternions); - if (values.empty() or values.size() > quaternions.size()) { - PARSE_ERROR( - context, - "Must specify at least the value of the quaternion, and optionally " - "up to " - << NumDerivs << " time derivatives."); - } - for (size_t i = 0; i < values.size(); i++) { - gsl::at(quaternions, i) = - DataVector{values[i][0], values[i][1], values[i][2], values[i][3]}; - } + initialize_angles_and_quats(); - if (initial_angles.has_value()) { - auto& angle_values = initial_angles.value(); - if (angle_values.empty() or angle_values.size() > angles.size()) { - PARSE_ERROR( - context, - "When specifying the angle, you must specify at least the value, " - "and optionally up to " - << NumDerivs << " time derivatives."); - } - for (size_t i = 0; i < angle_values.size(); i++) { - gsl::at(angles, i) = DataVector{angle_values[i][0], angle_values[i][1], - angle_values[i][2]}; - } - } - } else if (std::holds_alternative>( - initial_quaternions)) { - if (decay_timescale.has_value()) { - PARSE_ERROR(context, - "When specifying the initial quaternions from a volume file, " - "the decay timescale must be 'Auto'."); - } - auto& values_from_file = - std::get>(initial_quaternions); + if (initial_quaternions.empty() or initial_quaternions.size() > 3) { + PARSE_ERROR( + context, + "Must specify at least the value of the quaternion, and optionally " + "up to 2 time derivatives."); + } + for (size_t i = 0; i < initial_quaternions.size(); i++) { + gsl::at(quaternions, i) = DataVector{initial_quaternions[i]}; + } +} + +template +std::unique_ptr get_rotation( + const std::variant, FromVolumeFile>& + rotation_map_options, + const double initial_time, const double expiration_time) { + const std::string name = "Rotation"; + std::unique_ptr result{}; - for (size_t i = 0; i < values_from_file.quaternions.size(); i++) { - gsl::at(quaternions, i) = gsl::at(values_from_file.quaternions, i); - gsl::at(angles, i) = gsl::at(values_from_file.angle_values, i); + if (std::holds_alternative(rotation_map_options)) { + const auto& from_vol_file = std::get(rotation_map_options); + const auto volume_fot = + from_vol_file.retrieve_function_of_time({name}, initial_time); + + // It must be either a QuaternionFoT or a SettleToConstant + const auto* quat_volume_fot = + dynamic_cast*>( + volume_fot.at(name).get()); + const auto* settle_volume_fot = + dynamic_cast( + volume_fot.at(name).get()); + + if (UNLIKELY(quat_volume_fot == nullptr and settle_volume_fot == nullptr)) { + ERROR_NO_TRACE( + "Rotation function of time read from volume data is not a " + "QuaternionFunctionOfTime<3> or a SettleToConstantQuaternion. Cannot " + "use it to initialize the rotation map."); } - if (initial_angles.has_value()) { - // Reset angle func so derivs that weren't specified are zero - angles = make_array(DataVector{3, 0.0}); - auto& angle_values = initial_angles.value(); - if (angle_values.empty() or angle_values.size() > angles.size()) { - PARSE_ERROR( - context, - "When specifying the angle, you must specify at least the value, " - "and optionally up to " - << NumDerivs << " time derivatives."); - } - for (size_t i = 0; i < angle_values.size(); i++) { - gsl::at(angles, i) = DataVector{angle_values[i][0], angle_values[i][1], - angle_values[i][2]}; - } + result = volume_fot.at(name)->create_at_time(initial_time, expiration_time); + } else if (std::holds_alternative>( + rotation_map_options)) { + const auto& hard_coded_options = + std::get>(rotation_map_options); + + if (hard_coded_options.decay_timescale.has_value()) { + result = + std::make_unique( + hard_coded_options.quaternions, initial_time, + hard_coded_options.decay_timescale.value()); + } else { + result = std::make_unique< + domain::FunctionsOfTime::QuaternionFunctionOfTime<3>>( + initial_time, std::array{hard_coded_options.quaternions[0]}, + hard_coded_options.angles, expiration_time); } + } else { + ERROR("Unknown RotationMap."); } + + return result; } -#define NUMDERIVS(data) BOOST_PP_TUPLE_ELEM(0, data) +#define ALLOWSETTLE(data) BOOST_PP_TUPLE_ELEM(0, data) -#define INSTANTIATE(_, data) \ - template struct RotationMapOptions; +#define INSTANTIATE(_, data) \ + template struct RotationMapOptions; \ + template std::unique_ptr \ + get_rotation(const std::variant, \ + FromVolumeFile>& rotation_map_options, \ + double initial_time, double expiration_time); -GENERATE_INSTANTIATIONS(INSTANTIATE, (2, 3)) +GENERATE_INSTANTIATIONS(INSTANTIATE, (true, false)) #undef INSTANTIATE -#undef NUMDERIVS +#undef ALLOWSETTLE } // namespace domain::creators::time_dependent_options diff --git a/src/Domain/Creators/TimeDependentOptions/RotationMap.hpp b/src/Domain/Creators/TimeDependentOptions/RotationMap.hpp index 4e198ef51110..fb2e87a460ae 100644 --- a/src/Domain/Creators/TimeDependentOptions/RotationMap.hpp +++ b/src/Domain/Creators/TimeDependentOptions/RotationMap.hpp @@ -15,68 +15,83 @@ #include "Domain/Creators/TimeDependentOptions/FromVolumeFile.hpp" #include "Options/Auto.hpp" #include "Options/Context.hpp" +#include "Options/Options.hpp" #include "Options/String.hpp" #include "Utilities/TMPL.hpp" namespace domain::creators::time_dependent_options { /*! - * \brief Class to be used as an option for initializing rotation map - * coefficients. + * \brief Class that holds hard coded rotation map options from the input file. + * + * \details This class can also be used as an option tag with the \p type type + * alias, `name()` function, and \p help string. */ -template +template struct RotationMapOptions { - using type = Options::Auto; + using type = Options::Auto< + std::variant, FromVolumeFile>, + Options::AutoLabel::None>; static std::string name() { return "RotationMap"; } static constexpr Options::String help = { "Options for a time-dependent rotation of the coordinates. Specify " "'None' to not use this map."}; struct InitialQuaternions { - using type = std::variant>, - FromVolumeFile>; + using type = std::vector>; static constexpr Options::String help = { "Initial values for the quaternion of the rotation map. You can " "optionally specify its first two time derivatives. If time " "derivatives aren't specified, zero will be used."}; }; - struct InitialAngles { - using type = Options::Auto>>; - static constexpr Options::String help = { - "Initial values for the angle of the rotation map. If 'Auto' is " - "specified, the behavior will depend on the 'InitialQuaternions' " - "option. If you are reading the quaternion from a volume file, 'Auto' " - "will use the angle values from the volume file. If you are simply " - "specifying the quaternion and (optionally) its time derivatives, then " - "'Auto' here will set the angle and its time derivatives to zero. If " - "values are specified for the angle and its time derivatives, then " - "those will override anything computed from the 'InitialQuaternions' " - "option."}; + struct InitialAngularVelocity { + using type = std::array; + static constexpr Options::String help = {"The initial angular velocity."}; }; struct DecayTimescale { - using type = Options::Auto; + using type = double; static constexpr Options::String help = { "The timescale for how fast the rotation approaches its asymptotic " - "value. If this is specified, a SettleToConstant function of time will " - "be used. If 'Auto' is specified, a PiecewisePolynomial function of " - "time will be used. Note that if you are reading the initial " - "quaternions from a volume file, then this option must be 'Auto'"}; + "value with a SettleToConstant function of time."}; }; - using options = tmpl::list; + using non_settle_options = tmpl::list; + using settle_options = tmpl::list; + + using options = tmpl::conditional_t< + AllowSettleFoTs, + tmpl::list>, + non_settle_options>; RotationMapOptions() = default; + // Constructor for non SettleToConstant functions of time + // NOLINTNEXTLINE(google-explicit-constructor) + RotationMapOptions(const std::array& initial_angular_velocity, + const Options::Context& context = {}); + // Constructor for SettleToConst functions of time RotationMapOptions( - std::variant>, - FromVolumeFile> - initial_quaternions, - std::optional>> initial_angles, - std::optional decay_timescale_in, - const Options::Context& context = {}); + const std::vector>& initial_quaternions, + double decay_timescale_in, const Options::Context& context = {}); + + std::array quaternions{}; + std::array angles{}; + std::optional decay_timescale; - std::array quaternions{}; - std::array angles{}; - std::optional decay_timescale{}; + private: + void initialize_angles_and_quats(); }; + +/*! + * \brief Helper function that takes the variant of the rotation map options, + * and returns the fully constructed rotation function of time. + * + * \details Even if the function of time is read from a file, it will have a + * new \p initial_time and \p expiration_time. + */ +template +std::unique_ptr get_rotation( + const std::variant, FromVolumeFile>& + rotation_map_options, + double initial_time, double expiration_time); } // namespace domain::creators::time_dependent_options diff --git a/src/Domain/Creators/TimeDependentOptions/ShapeMap.cpp b/src/Domain/Creators/TimeDependentOptions/ShapeMap.cpp index d0a7209cc8e6..382c0d91cfcd 100644 --- a/src/Domain/Creators/TimeDependentOptions/ShapeMap.cpp +++ b/src/Domain/Creators/TimeDependentOptions/ShapeMap.cpp @@ -8,6 +8,7 @@ #include #include #include +#include #include #include #include @@ -16,6 +17,8 @@ #include "DataStructures/DataVector.hpp" #include "DataStructures/ModalVector.hpp" #include "DataStructures/Tensor/Tensor.hpp" +#include "Domain/FunctionsOfTime/FunctionOfTime.hpp" +#include "Domain/FunctionsOfTime/PiecewisePolynomial.hpp" #include "Domain/Structure/ObjectLabel.hpp" #include "FromVolumeFile.hpp" #include "NumericalAlgorithms/SphericalHarmonics/Spherepack.hpp" @@ -55,26 +58,109 @@ YlmsFromSpEC::YlmsFromSpEC(std::string dat_filename_in, match_time_epsilon(match_time_epsilon_in), set_l1_coefs_to_zero(set_l1_coefs_to_zero_in) {} +template +FromVolumeFileShapeSize::FromVolumeFileShapeSize( + const bool transition_ends_at_cube_in, std::string h5_filename, + std::string subfile_name) + : FromVolumeFile(std::move(h5_filename), std::move(subfile_name)), + transition_ends_at_cube(transition_ends_at_cube_in) { + const auto shape_fot_map = + retrieve_function_of_time({"Shape" + name(Object)}, std::nullopt); + const auto& shape_fot = shape_fot_map.at("Shape" + name(Object)); + + const double initial_time = shape_fot->time_bounds()[0]; + const auto function = shape_fot->func(initial_time); + + // num_components = 2 * (l_max + 1)**2 if l_max == m_max which it is for the + // shape map. This is why we can divide by 2 and take the sqrt without + // worrying about odd numbers or non-perfect squares + l_max = -1 + sqrt(function[0].size() / 2); +} + +template +size_t l_max_from_shape_options( + const std::variant, + FromVolumeFileShapeSize>& shape_map_options) { + return std::visit([](auto variant) { return variant.l_max; }, + shape_map_options); +} + +template +bool transition_ends_at_cube_from_shape_options( + const std::variant, + FromVolumeFileShapeSize>& shape_map_options) { + return std::visit( + [](auto variant) { return variant.transition_ends_at_cube; }, + shape_map_options); +} + template -std::pair, std::array> -initial_shape_and_size_funcs( - const ShapeMapOptions& shape_options, - const double deformed_radius) { - const DataVector shape_zeros{ - ylm::Spherepack::spectral_size(shape_options.l_max, shape_options.l_max), - 0.0}; +FunctionsOfTimeMap get_shape_and_size( + const std::variant, + FromVolumeFileShapeSize>& shape_map_options, + const double initial_time, const double shape_expiration_time, + const double size_expiration_time, const double deformed_radius) { + const size_t l_max = l_max_from_shape_options(shape_map_options); + const DataVector shape_zeros{ylm::Spherepack::spectral_size(l_max, l_max), + 0.0}; + const std::string shape_name = "Shape" + name(Object); + const std::string size_name = "Size" + name(Object); + + FunctionsOfTimeMap result{}; + + if (std::holds_alternative>( + shape_map_options)) { + const auto& from_vol_file = + std::get>(shape_map_options); + const auto volume_fots = from_vol_file.retrieve_function_of_time( + {shape_name, size_name}, initial_time); + + const auto check_fot = [&](const std::string& name) { + // It must be a PiecewisePolynomial + if (UNLIKELY(dynamic_cast< + domain::FunctionsOfTime::PiecewisePolynomial*>( + volume_fots.at(name).get()) == nullptr)) { + ERROR_NO_TRACE(name << " function of time read from volume data is not " + "a PiecewisePolynomial<" + << MaxDeriv << ">. Cannot use it to initialize the " + << name << " map."); + } + }; + + check_fot.template operator()<2>(shape_name); + check_fot.template operator()<3>(size_name); + + result[shape_name] = + volume_fots.at(shape_name) + ->create_at_time(initial_time, shape_expiration_time); + result[size_name] = volume_fots.at(size_name)->create_at_time( + initial_time, size_expiration_time); + + return result; + } + + if (not std::holds_alternative< + ShapeMapOptions>( + shape_map_options)) { + ERROR("Unknown ShapeMap."); + } + + const auto& hard_coded_options = + std::get>( + shape_map_options); std::array shape_funcs = make_array<3, DataVector>(shape_zeros); std::array size_funcs = make_array<4, DataVector>(DataVector{1, 0.0}); - if (shape_options.initial_values.has_value()) { + if (hard_coded_options.initial_values.has_value()) { if (std::holds_alternative( - shape_options.initial_values.value())) { - const ylm::Spherepack ylm{shape_options.l_max, shape_options.l_max}; + hard_coded_options.initial_values.value())) { + const ylm::Spherepack ylm{hard_coded_options.l_max, + hard_coded_options.l_max}; const auto& mass_and_spin = std::get( - shape_options.initial_values.value()); + hard_coded_options.initial_values.value()); const DataVector radial_distortion = deformed_radius - get(gr::Solutions::kerr_schild_radius_from_boyer_lindquist( @@ -86,16 +172,15 @@ initial_shape_and_size_funcs( // Set l=0 for shape map to 0 because size control will adjust l=0 shape_funcs[0][0] = 0.0; } else if (std::holds_alternative( - shape_options.initial_values.value())) { + hard_coded_options.initial_values.value())) { const auto& files = - std::get(shape_options.initial_values.value()); + std::get(hard_coded_options.initial_values.value()); const std::string& h5_filename = files.h5_filename; const std::vector& subfile_names = files.subfile_names; const double match_time = files.match_time; const double match_time_epsilon = files.match_time_epsilon.value_or(1e-12); const bool set_l1_coefs_to_zero = files.set_l1_coefs_to_zero; - const size_t l_max = shape_options.l_max; ylm::SpherepackIterator iter{l_max, l_max}; for (size_t i = 0; i < subfile_names.size(); i++) { @@ -105,7 +190,7 @@ initial_shape_and_size_funcs( h5_filename, gsl::at(subfile_names, i), match_time, match_time_epsilon, files.check_frame); const ylm::Strahlkorper this_strahlkorper{ - shape_options.l_max, 1.0, std::array{0.0, 0.0, 0.0}}; + hard_coded_options.l_max, 1.0, std::array{0.0, 0.0, 0.0}}; // The coefficients in the shape map are stored as the negative // coefficients of the strahlkorper, so we need to multiply by -1 here. @@ -129,9 +214,9 @@ initial_shape_and_size_funcs( } } } else if (std::holds_alternative( - shape_options.initial_values.value())) { + hard_coded_options.initial_values.value())) { const auto& spec_option = - std::get(shape_options.initial_values.value()); + std::get(hard_coded_options.initial_values.value()); const std::string& dat_filename = spec_option.dat_filename; const double match_time = spec_option.match_time; const double match_time_epsilon = @@ -144,7 +229,7 @@ initial_shape_and_size_funcs( } std::string line{}; size_t total_col = 0; - std::optional l_max{}; + std::optional file_l_max{}; std::array center{}; ModalVector coefficients{}; // This will be actually set below @@ -173,7 +258,7 @@ initial_shape_and_size_funcs( continue; } - if (l_max.has_value()) { + if (file_l_max.has_value()) { ERROR("Found more than one time in the SpEC dat file " << dat_filename << " that is within a relative epsilon of " << match_time_epsilon << " of the time requested " << time); @@ -181,42 +266,43 @@ initial_shape_and_size_funcs( // Casting to an integer floors a double, so we add 0.5 before we take // the sqrt to avoid any rounding issues - const auto l_max_plus_one = + const auto file_l_max_plus_one = static_cast(sqrt(static_cast(total_col) + 0.5)); - if (l_max_plus_one == 0) { + if (file_l_max_plus_one == 0) { ERROR( "Invalid l_max from SpEC dat file. l_max + 1 was computed to be " "0"); } - l_max = l_max_plus_one - 1; + file_l_max = file_l_max_plus_one - 1; ss >> center[0]; ss >> center[1]; ss >> center[2]; - coefficients.destructive_resize( - ylm::Spherepack::spectral_size(l_max.value(), l_max.value())); + coefficients.destructive_resize(ylm::Spherepack::spectral_size( + file_l_max.value(), file_l_max.value())); - file_iter = ylm::SpherepackIterator{l_max.value(), l_max.value()}; + file_iter = + ylm::SpherepackIterator{file_l_max.value(), file_l_max.value()}; - for (int l = 0; l <= static_cast(l_max.value()); l++) { + for (int l = 0; l <= static_cast(file_l_max.value()); l++) { for (int m = -l; m <= l; m++) { ss >> coefficients[file_iter.set(static_cast(l), m)()]; } } } - if (not l_max.has_value()) { + if (not file_l_max.has_value()) { ERROR_NO_TRACE("Unable to find requested time " << time << " within an epsilon of " << match_time_epsilon << " in SpEC dat file " << dat_filename); } const ylm::Strahlkorper file_strahlkorper{ - l_max.value(), l_max.value(), coefficients, center}; + file_l_max.value(), file_l_max.value(), coefficients, center}; const ylm::Strahlkorper this_strahlkorper{ - shape_options.l_max, 1.0, std::array{0.0, 0.0, 0.0}}; - ylm::SpherepackIterator iter{shape_options.l_max, shape_options.l_max}; + l_max, 1.0, std::array{0.0, 0.0, 0.0}}; + ylm::SpherepackIterator iter{l_max, l_max}; shape_funcs[0] = -1.0 * file_strahlkorper.ylm_spherepack().prolong_or_restrict( @@ -231,46 +317,61 @@ initial_shape_and_size_funcs( shape_funcs[0][iter.set(1_st, m)()] = 0.0; } } - } else if (std::holds_alternative>>( - shape_options.initial_values.value())) { - const auto& volume_file_options = - std::get>>( - shape_options.initial_values.value()); - - shape_funcs = volume_file_options.shape_values; - size_funcs = volume_file_options.size_values; } } // If any size options were specified, those override the values from the // shape coefs - if (shape_options.initial_size_values.has_value()) { + if (hard_coded_options.initial_size_values.has_value()) { for (size_t i = 0; i < 3; i++) { gsl::at(size_funcs, i)[0] = - gsl::at(shape_options.initial_size_values.value(), i); + gsl::at(hard_coded_options.initial_size_values.value(), i); } } - return std::make_pair(std::move(shape_funcs), std::move(size_funcs)); + result[shape_name] = + std::make_unique>( + initial_time, std::move(shape_funcs), shape_expiration_time); + result[size_name] = std::make_unique>( + initial_time, std::move(size_funcs), size_expiration_time); + + return result; } -#define INCLUDETRANSITION(data) BOOST_PP_TUPLE_ELEM(0, data) -#define OBJECT(data) BOOST_PP_TUPLE_ELEM(1, data) +#define OBJECT(data) BOOST_PP_TUPLE_ELEM(0, data) +#define INCLUDETRANSITION(data) BOOST_PP_TUPLE_ELEM(1, data) + +#define INSTANTIATE(_, data) \ + template size_t l_max_from_shape_options( \ + const std::variant< \ + ShapeMapOptions, \ + FromVolumeFileShapeSize>& shape_map_options); \ + template bool transition_ends_at_cube_from_shape_options( \ + const std::variant< \ + ShapeMapOptions, \ + FromVolumeFileShapeSize>& shape_map_options); \ + template FunctionsOfTimeMap get_shape_and_size( \ + const std::variant< \ + ShapeMapOptions, \ + FromVolumeFileShapeSize>& shape_map_options, \ + double initial_time, double shape_expiration_time, \ + double size_expiration_time, double deformed_radius); + +GENERATE_INSTANTIATIONS(INSTANTIATE, + (domain::ObjectLabel::A, domain::ObjectLabel::B, + domain::ObjectLabel::None), + (true, false)) + +#undef INCLUDETRANSITION +#undef INSTANTIATE -#define INSTANTIATE(_, data) \ - template class ShapeMapOptions; \ - template std::pair, std::array> \ - initial_shape_and_size_funcs( \ - const ShapeMapOptions& \ - shape_options, \ - double deformed_radius); +#define INSTANTIATE(_, data) \ + template class FromVolumeFileShapeSize; -GENERATE_INSTANTIATIONS(INSTANTIATE, (true, false), +GENERATE_INSTANTIATIONS(INSTANTIATE, (domain::ObjectLabel::A, domain::ObjectLabel::B, domain::ObjectLabel::None)) -#undef INCLUDETRANSITION #undef OBJECT #undef INSTANTIATE - } // namespace domain::creators::time_dependent_options diff --git a/src/Domain/Creators/TimeDependentOptions/ShapeMap.hpp b/src/Domain/Creators/TimeDependentOptions/ShapeMap.hpp index 0d181fb2ae59..6b65787d5773 100644 --- a/src/Domain/Creators/TimeDependentOptions/ShapeMap.hpp +++ b/src/Domain/Creators/TimeDependentOptions/ShapeMap.hpp @@ -12,6 +12,7 @@ #include "DataStructures/DataVector.hpp" #include "Domain/Creators/TimeDependentOptions/FromVolumeFile.hpp" +#include "Domain/FunctionsOfTime/FunctionOfTime.hpp" #include "Domain/Structure/ObjectLabel.hpp" #include "NumericalAlgorithms/SphericalHarmonics/IO/ReadSurfaceYlm.hpp" #include "NumericalAlgorithms/SphericalHarmonics/Strahlkorper.hpp" @@ -122,10 +123,10 @@ struct YlmsFromFile { std::optional match_time_epsilon_in, bool set_l1_coefs_to_zero_in, bool check_frame_in = true); - std::string h5_filename{}; - std::vector subfile_names{}; + std::string h5_filename; + std::vector subfile_names; double match_time{}; - std::optional match_time_epsilon{}; + std::optional match_time_epsilon; bool set_l1_coefs_to_zero{}; bool check_frame{true}; }; @@ -172,15 +173,51 @@ struct YlmsFromSpEC { std::optional match_time_epsilon_in, bool set_l1_coefs_to_zero_in); - std::string dat_filename{}; + std::string dat_filename; double match_time{}; - std::optional match_time_epsilon{}; + std::optional match_time_epsilon; bool set_l1_coefs_to_zero{}; }; +namespace detail { +struct TransitionEndsAtCube { + using type = bool; + static constexpr Options::String help = { + "If 'true', the shape map transition function will be 0 at the cubical " + "boundary around the object. If 'false' the transition function will " + "be 0 at the outer radius of the inner sphere around the object"}; +}; +} // namespace detail + +/*! + * \brief Specialized version of `FromVolumeFile` for the shape map + * + * \details This is needed because the regular `FromVolumeFile` doesn't have + * options for domain settings like `TransitionEndsAtCube`. + */ +template +struct FromVolumeFileShapeSize : public FromVolumeFile { + using options = + tmpl::push_front; + + FromVolumeFileShapeSize() = default; + FromVolumeFileShapeSize(bool transition_ends_at_cube_in, + std::string h5_filename, std::string subfile_name); + + size_t l_max{}; + bool transition_ends_at_cube{}; + + private: + std::string h5_filename_; + std::string subfile_name_; +}; + /*! * \brief Class to be used as an option for initializing shape map coefficients. * + * \details This class can also be used as an option tag with the \p type type + * alias, `name()` function, and \p help string. + * * \tparam IncludeTransitionEndsAtCube This is mainly added for the * `domain::creators::BinaryCompactObject` domain. * \tparam Object Which object that this shape map represents. Use @@ -189,7 +226,10 @@ struct YlmsFromSpEC { */ template struct ShapeMapOptions { - using type = Options::Auto; + using type = Options::Auto< + std::variant, + FromVolumeFileShapeSize>, + Options::AutoLabel::None>; static std::string name() { return "ShapeMap" + get_output(Object); } static constexpr Options::String help = { "Options for a time-dependent distortion (shape) map about the " @@ -204,8 +244,7 @@ struct ShapeMapOptions { struct InitialValues { using type = Options::Auto< - std::variant>>, + std::variant, Spherical>; static constexpr Options::String help = { "Initial Ylm coefficients for the shape map. Specify 'Spherical' for " @@ -222,30 +261,20 @@ struct ShapeMapOptions { "set the radius of the sphere in the grid frame (before deformation)."}; }; - struct TransitionEndsAtCube { - using type = bool; - static constexpr Options::String help = { - "If 'true', the shape map transition function will be 0 at the cubical " - "boundary around the object. If 'false' the transition function will " - "be 0 at the outer radius of the inner sphere around the object"}; - }; - using common_options = tmpl::list; - using options = - tmpl::conditional_t, - common_options>; + using options = tmpl::conditional_t< + IncludeTransitionEndsAtCube, + tmpl::push_back, + common_options>; ShapeMapOptions() = default; - ShapeMapOptions( - size_t l_max_in, - std::optional< - std::variant>>> - initial_values_in, - std::optional> initial_size_values_in = - std::nullopt, - bool transition_ends_at_cube_in = false) + ShapeMapOptions(size_t l_max_in, + std::optional> + initial_values_in, + std::optional> initial_size_values_in = + std::nullopt, + bool transition_ends_at_cube_in = false) : l_max(l_max_in), initial_values(std::move(initial_values_in)), initial_size_values(initial_size_values_in), @@ -253,16 +282,42 @@ struct ShapeMapOptions { size_t l_max{}; std::optional< - std::variant>>> - initial_values{}; - std::optional> initial_size_values{}; + std::variant> + initial_values; + std::optional> initial_size_values; bool transition_ends_at_cube{false}; }; +/*! + * \brief Helper function to get LMax from the different variants that the shape + * map options could be. + */ +template +size_t l_max_from_shape_options( + const std::variant, + FromVolumeFileShapeSize>& shape_map_options); + +/*! + * \brief Helper function to get whether the shape map transition function ends + * at the cube from the different variants that the shape map options could be. + */ +template +bool transition_ends_at_cube_from_shape_options( + const std::variant, + FromVolumeFileShapeSize>& shape_map_options); + +/*! + * \brief Helper function that takes the variant of the shape map options, and + * returns the fully constructed shape and size functions of time. + * + * \details Even if the functions of time are read from a file, they will have a + * new \p initial_time, \p shape_expiration_time, and \p size_expiration_time. + * The \p deformed_radius is only used for the non-volume file variants. + */ template -std::pair, std::array> -initial_shape_and_size_funcs( - const ShapeMapOptions& shape_options, - double deformed_radius); +FunctionsOfTimeMap get_shape_and_size( + const std::variant, + FromVolumeFileShapeSize>& shape_map_options, + double initial_time, double shape_expiration_time, + double size_expiration_time, double deformed_radius); } // namespace domain::creators::time_dependent_options diff --git a/src/Domain/Creators/TimeDependentOptions/Sphere.cpp b/src/Domain/Creators/TimeDependentOptions/Sphere.cpp index 8916fd1d6785..3b9f25fcb445 100644 --- a/src/Domain/Creators/TimeDependentOptions/Sphere.cpp +++ b/src/Domain/Creators/TimeDependentOptions/Sphere.cpp @@ -31,34 +31,11 @@ namespace domain::creators::sphere { -TimeDependentMapOptions::RotationMapOptions::RotationMapOptions() = default; -TimeDependentMapOptions::RotationMapOptions::RotationMapOptions( - const std::array, 3> initial_values_in, - const double decay_timescale_in) - : initial_values(initial_values_in), decay_timescale(decay_timescale_in) {} - -TimeDependentMapOptions::ExpansionMapOptions::ExpansionMapOptions() = default; -TimeDependentMapOptions::ExpansionMapOptions::ExpansionMapOptions( - const std::array initial_values_in, - const double decay_timescale_in, - const std::array initial_values_outer_boundary_in, - const double decay_timescale_outer_boundary_in) - : initial_values(initial_values_in), - decay_timescale(decay_timescale_in), - initial_values_outer_boundary(initial_values_outer_boundary_in), - decay_timescale_outer_boundary(decay_timescale_outer_boundary_in) {} - -TimeDependentMapOptions::TranslationMapOptions::TranslationMapOptions() = - default; -TimeDependentMapOptions::TranslationMapOptions::TranslationMapOptions( - const std::array, 3> initial_values_in) - : initial_values(initial_values_in) {} - TimeDependentMapOptions::TimeDependentMapOptions( - const double initial_time, std::optional shape_map_options, - std::optional rotation_map_options, - std::optional expansion_map_options, - std::optional translation_map_options, + const double initial_time, ShapeMapOptionType shape_map_options, + RotationMapOptionType rotation_map_options, + ExpansionMapOptionType expansion_map_options, + TranslationMapOptionType translation_map_options, const bool transition_rot_scale_trans) : initial_time_(initial_time), shape_map_options_(std::move(shape_map_options)), @@ -81,6 +58,8 @@ TimeDependentMapOptions::create_functions_of_time( std::unordered_map expiration_times{ {size_name, std::numeric_limits::infinity()}, {shape_name, std::numeric_limits::infinity()}, + {rotation_name, std::numeric_limits::infinity()}, + {expansion_name, std::numeric_limits::infinity()}, {translation_name, std::numeric_limits::infinity()}}; // If we have control systems, overwrite these expiration times with the ones @@ -90,95 +69,35 @@ TimeDependentMapOptions::create_functions_of_time( } if (shape_map_options_.has_value()) { - auto [shape_funcs, size_funcs] = - time_dependent_options::initial_shape_and_size_funcs( - shape_map_options_.value(), deformed_radius_); - - // ShapeMap FunctionOfTime - result[shape_name] = - std::make_unique>( - initial_time_, std::move(shape_funcs), - expiration_times.at(shape_name)); - - // Size FunctionOfTime (used in ShapeMap) - result[size_name] = - std::make_unique>( - initial_time_, std::move(size_funcs), - expiration_times.at(size_name)); + // ShapeMap and Size FunctionOfTime (used in ShapeMap) + auto shape_and_size = time_dependent_options::get_shape_and_size( + shape_map_options_.value(), initial_time_, + expiration_times.at(shape_name), expiration_times.at(size_name), + deformed_radius_); + result.merge(shape_and_size); } // ExpansionMap FunctionOfTime if (expansion_map_options_.has_value()) { - result[expansion_name] = - std::make_unique( - std::array{ - {{gsl::at(expansion_map_options_.value().initial_values, 0)}, - {gsl::at(expansion_map_options_.value().initial_values, 1)}, - {gsl::at(expansion_map_options_.value().initial_values, 2)}}}, - initial_time_, expansion_map_options_.value().decay_timescale); + auto expansion_functions_of_time = time_dependent_options::get_expansion( + expansion_map_options_.value(), initial_time_, + expiration_times.at(expansion_name)); - // ExpansionMap in the Outer regionFunctionOfTime - result[expansion_outer_boundary_name] = std::make_unique< - FunctionsOfTime::SettleToConstant>( - std::array{ - {{gsl::at( - expansion_map_options_.value().initial_values_outer_boundary, - 0)}, - {gsl::at( - expansion_map_options_.value().initial_values_outer_boundary, - 1)}, - {gsl::at( - expansion_map_options_.value().initial_values_outer_boundary, - 2)}}}, - initial_time_, - expansion_map_options_.value().decay_timescale_outer_boundary); + result.merge(expansion_functions_of_time); } - DataVector initial_quaternion_value{4, 0.0}; - DataVector initial_quaternion_first_derivative_value{4, 0.0}; - DataVector initial_quaternion_second_derivative_value{4, 0.0}; - // RotationMap FunctionOfTime if (rotation_map_options_.has_value()) { - for (size_t i = 0; i < 4; i++) { - initial_quaternion_value[i] = - gsl::at(gsl::at(rotation_map_options_.value().initial_values, 0), i); - initial_quaternion_first_derivative_value[i] = - gsl::at(gsl::at(rotation_map_options_.value().initial_values, 1), i); - initial_quaternion_second_derivative_value[i] = - gsl::at(gsl::at(rotation_map_options_.value().initial_values, 2), i); - } - result[rotation_name] = - std::make_unique( - std::array{ - std::move(initial_quaternion_value), - std::move(initial_quaternion_first_derivative_value), - std::move(initial_quaternion_second_derivative_value)}, - initial_time_, rotation_map_options_.value().decay_timescale); + result[rotation_name] = time_dependent_options::get_rotation( + rotation_map_options_.value(), initial_time_, + expiration_times.at(rotation_name)); } - DataVector initial_translation_center{3, 0.0}; - DataVector initial_translation_velocity{3, 0.0}; - DataVector initial_translation_acceleration{3, 0.0}; - // Translation FunctionOfTime if (translation_map_options_.has_value()) { - for (size_t i = 0; i < 3; i++) { - initial_translation_center[i] = gsl::at( - gsl::at(translation_map_options_.value().initial_values, 0), i); - initial_translation_velocity[i] = gsl::at( - gsl::at(translation_map_options_.value().initial_values, 1), i); - initial_translation_acceleration[i] = gsl::at( - gsl::at(translation_map_options_.value().initial_values, 2), i); - } - result[translation_name] = - std::make_unique>( - initial_time_, - std::array{ - {std::move(initial_translation_center), - std::move(initial_translation_velocity), - std::move(initial_translation_acceleration)}}, - expiration_times.at(translation_name)); + result[translation_name] = time_dependent_options::get_translation( + translation_map_options_.value(), initial_time_, + expiration_times.at(translation_name)); } return result; @@ -190,6 +109,8 @@ void TimeDependentMapOptions::build_maps( const double outer_radius) { filled_ = filled; if (shape_map_options_.has_value()) { + const size_t l_max = time_dependent_options::l_max_from_shape_options( + shape_map_options_.value()); std::unique_ptr transition_func; @@ -230,12 +151,9 @@ void TimeDependentMapOptions::build_maps( /* outer_sphericity */ 1.0, static_cast(gsl::at(axes, j % 6))); } - gsl::at(shape_maps_, j) = ShapeMap{center, - shape_map_options_->l_max, - shape_map_options_->l_max, - std::move(transition_func), - shape_name, - size_name}; + gsl::at(shape_maps_, j) = + ShapeMap{center, l_max, l_max, std::move(transition_func), + shape_name, size_name}; } } else { // Shape map transitions from 1 to 0 from the inner radius to the first @@ -247,12 +165,9 @@ void TimeDependentMapOptions::build_maps( std::make_unique(inner_radius, shape_outer_radius); - shape_maps_[0] = ShapeMap{center, - shape_map_options_->l_max, - shape_map_options_->l_max, - std::move(transition_func), - shape_name, - size_name}; + shape_maps_[0] = + ShapeMap{center, l_max, l_max, std::move(transition_func), + shape_name, size_name}; } } diff --git a/src/Domain/Creators/TimeDependentOptions/Sphere.hpp b/src/Domain/Creators/TimeDependentOptions/Sphere.hpp index 6ae2aed27232..495fde452334 100644 --- a/src/Domain/Creators/TimeDependentOptions/Sphere.hpp +++ b/src/Domain/Creators/TimeDependentOptions/Sphere.hpp @@ -16,7 +16,10 @@ #include "Domain/CoordinateMaps/TimeDependent/RotScaleTrans.hpp" #include "Domain/CoordinateMaps/TimeDependent/Shape.hpp" #include "Domain/CoordinateMaps/TimeDependent/Translation.hpp" +#include "Domain/Creators/TimeDependentOptions/ExpansionMap.hpp" +#include "Domain/Creators/TimeDependentOptions/RotationMap.hpp" #include "Domain/Creators/TimeDependentOptions/ShapeMap.hpp" +#include "Domain/Creators/TimeDependentOptions/TranslationMap.hpp" #include "Domain/FunctionsOfTime/FunctionOfTime.hpp" #include "Domain/Structure/ObjectLabel.hpp" #include "Options/Auto.hpp" @@ -81,113 +84,18 @@ struct TimeDependentMapOptions { using ShapeMapOptions = time_dependent_options::ShapeMapOptions; + using ShapeMapOptionType = typename ShapeMapOptions::type::value_type; - struct RotationMapOptions { - using type = Options::Auto; - static std::string name() { return "RotationMap"; } - static constexpr Options::String help = { - "Options for a time-dependent rotation map about an arbitrary axis. " - "Specify 'None' to not use this map."}; - - struct InitialValues { - using type = std::array, 3>; - static constexpr Options::String help = { - "The initial values for the quaternion function and its first two " - "derivatives."}; - }; - - struct DecayTimescaleRotation { - using type = double; - static constexpr Options::String help = { - "The timescale for how fast the rotation approaches its asymptotic " - "value."}; - }; - - using options = tmpl::list; - RotationMapOptions(); - RotationMapOptions(std::array, 3> initial_values_in, - double decay_timescale_in); - - std::array, 3> initial_values{}; - double decay_timescale{std::numeric_limits::signaling_NaN()}; - }; - - struct ExpansionMapOptions { - using type = Options::Auto; - static std::string name() { return "ExpansionMap"; } - static constexpr Options::String help = { - "Options for the expansion map. Specify 'None' to not use this map."}; - - struct InitialValues { - using type = std::array; - static constexpr Options::String help = { - "Initial value and first two derivatives of expansion."}; - }; - - struct DecayTimescaleExpansion { - using type = double; - static constexpr Options::String help = { - "The timescale for how fast the expansion approaches " - "its asymptotic value."}; - }; + using RotationMapOptions = time_dependent_options::RotationMapOptions; + using RotationMapOptionType = typename RotationMapOptions::type::value_type; - struct InitialValuesOuterBoundary { - using type = std::array; - static constexpr Options::String help = { - "Initial value and first two derivatives of expansion outside the " - "transition region."}; - }; + using ExpansionMapOptions = time_dependent_options::ExpansionMapOptions; + using ExpansionMapOptionType = typename ExpansionMapOptions::type::value_type; - struct DecayTimescaleExpansionOuterBoundary { - using type = double; - static constexpr Options::String help = { - "The timescale for how fast the expansion approaches " - "its asymptotic value outside the transition region."}; - }; - - using options = tmpl::list; - ExpansionMapOptions(); - ExpansionMapOptions(std::array initial_values_in, - double decay_timescale_in, - std::array initial_values_outer_boundary_in, - double decay_timescale_outer_boundary_in); - - std::array initial_values{ - std::numeric_limits::signaling_NaN(), - std::numeric_limits::signaling_NaN(), - std::numeric_limits::signaling_NaN()}; - double decay_timescale{std::numeric_limits::signaling_NaN()}; - std::array initial_values_outer_boundary{ - std::numeric_limits::signaling_NaN(), - std::numeric_limits::signaling_NaN(), - std::numeric_limits::signaling_NaN()}; - double decay_timescale_outer_boundary{ - std::numeric_limits::signaling_NaN()}; - }; - - struct TranslationMapOptions { - using type = Options::Auto; - static std::string name() { return "TranslationMap"; } - static constexpr Options::String help = { - "Options for a time-dependent translation map in that keeps the " - "outer boundary fixed. Specify 'None' to not use this map."}; - - struct InitialValues { - using type = std::array, 3>; - static constexpr Options::String help = { - "Initial values for the translation map, its velocity and " - "acceleration."}; - }; - - using options = tmpl::list; - TranslationMapOptions(); - explicit TranslationMapOptions( - std::array, 3> initial_values_in); - - std::array, 3> initial_values{}; - }; + using TranslationMapOptions = + time_dependent_options::TranslationMapOptions<3>; + using TranslationMapOptionType = + typename TranslationMapOptions::type::value_type; struct TransitionRotScaleTrans { using type = bool; @@ -205,12 +113,12 @@ struct TimeDependentMapOptions { TimeDependentMapOptions() = default; - TimeDependentMapOptions( - double initial_time, std::optional shape_map_options, - std::optional rotation_map_options, - std::optional expansion_map_options, - std::optional translation_map_options, - bool transition_rot_scale_trans); + TimeDependentMapOptions(double initial_time, + ShapeMapOptionType shape_map_options, + RotationMapOptionType rotation_map_options, + ExpansionMapOptionType expansion_map_options, + TranslationMapOptionType translation_map_options, + bool transition_rot_scale_trans); /*! * \brief Create the function of time map using the options that were @@ -299,10 +207,10 @@ struct TimeDependentMapOptions { RotScaleTransMap inner_rot_scale_trans_map_{}; RotScaleTransMap transition_rot_scale_trans_map_{}; - std::optional shape_map_options_{}; - std::optional rotation_map_options_{}; - std::optional expansion_map_options_{}; - std::optional translation_map_options_{}; + ShapeMapOptionType shape_map_options_; + RotationMapOptionType rotation_map_options_; + ExpansionMapOptionType expansion_map_options_; + TranslationMapOptionType translation_map_options_; bool transition_rot_scale_trans_{false}; }; } // namespace domain::creators::sphere diff --git a/src/Domain/Creators/TimeDependentOptions/TranslationMap.cpp b/src/Domain/Creators/TimeDependentOptions/TranslationMap.cpp index 0c9528232939..5d3bdfc55abf 100644 --- a/src/Domain/Creators/TimeDependentOptions/TranslationMap.cpp +++ b/src/Domain/Creators/TimeDependentOptions/TranslationMap.cpp @@ -9,35 +9,75 @@ #include "DataStructures/DataVector.hpp" #include "Domain/Creators/TimeDependentOptions/FromVolumeFile.hpp" +#include "Domain/FunctionsOfTime/FunctionOfTime.hpp" +#include "Domain/FunctionsOfTime/PiecewisePolynomial.hpp" +#include "Options/Context.hpp" +#include "Options/ParseError.hpp" #include "Utilities/GenerateInstantiations.hpp" namespace domain::creators::time_dependent_options { template TranslationMapOptions::TranslationMapOptions( - std::variant, 3>, - FromVolumeFile> - values_from_options) { - if (std::holds_alternative, 3>>( - values_from_options)) { - auto& values = - std::get, 3>>(values_from_options); - for (size_t i = 0; i < initial_values.size(); i++) { - gsl::at(initial_values, i) = DataVector{Dim, 0.0}; - for (size_t j = 0; j < Dim; j++) { - gsl::at(initial_values, i)[j] = gsl::at(gsl::at(values, i), j); - } + const std::array, 3>& initial_values_in, + const Options::Context& context) { + if (initial_values_in.empty() or initial_values_in.size() > 3) { + PARSE_ERROR( + context, + "Must specify at least the value of the translation, and optionally " + "up to 2 time derivatives."); + } + + for (size_t i = 0; i < initial_values_in.size(); i++) { + gsl::at(initial_values, i) = DataVector{gsl::at(initial_values_in, i)}; + } +} + +template +std::unique_ptr get_translation( + const std::variant, FromVolumeFile>& + translation_map_options, + const double initial_time, const double expiration_time) { + const std::string name = "Translation"; + std::unique_ptr result{}; + + if (std::holds_alternative(translation_map_options)) { + const auto& from_vol_file = + std::get(translation_map_options); + const auto volume_fot = + from_vol_file.retrieve_function_of_time({name}, initial_time); + + // It must be a PiecewisePolynomial + if (UNLIKELY(dynamic_cast*>( + volume_fot.at(name).get()) == nullptr)) { + ERROR_NO_TRACE( + "Translation function of time read from volume data is not a " + "PiecewisePolynomial<2>. Cannot use it to initialize the translation " + "map."); } - } else if (std::holds_alternative>( - values_from_options)) { - auto& values_from_file = - std::get>(values_from_options); - initial_values = values_from_file.values; + + result = volume_fot.at(name)->create_at_time(initial_time, expiration_time); + } else if (std::holds_alternative>( + translation_map_options)) { + const auto& hard_coded_options = + std::get>(translation_map_options); + + result = std::make_unique>( + initial_time, hard_coded_options.initial_values, expiration_time); + } else { + ERROR("Unknown TranslationMap."); } + + return result; } #define DIM(data) BOOST_PP_TUPLE_ELEM(0, data) -#define INSTANTIATE(_, data) template class TranslationMapOptions; +#define INSTANTIATE(_, data) \ + template class TranslationMapOptions; \ + template std::unique_ptr \ + get_translation(const std::variant, \ + FromVolumeFile>& translation_map_options, \ + double initial_time, double expiration_time); GENERATE_INSTANTIATIONS(INSTANTIATE, (1, 2, 3)) diff --git a/src/Domain/Creators/TimeDependentOptions/TranslationMap.hpp b/src/Domain/Creators/TimeDependentOptions/TranslationMap.hpp index b12aa7379bac..c3dec791ff55 100644 --- a/src/Domain/Creators/TimeDependentOptions/TranslationMap.hpp +++ b/src/Domain/Creators/TimeDependentOptions/TranslationMap.hpp @@ -18,33 +18,51 @@ namespace domain::creators::time_dependent_options { /*! - * \brief Class to be used as an option for initializing translation map - * coefficients. + * \brief Class that holds hard coded translation map options from the input + * file. + * + * \details This class can also be used as an option tag with the \p type type + * alias, `name()` function, and \p help string. */ template struct TranslationMapOptions { - using type = Options::Auto; + using type = + Options::Auto, FromVolumeFile>, + Options::AutoLabel::None>; static std::string name() { return "TranslationMap"; } static constexpr Options::String help = { "Options for a time-dependent translation of the coordinates. Specify " "'None' to not use this map."}; struct InitialValues { - using type = std::variant, 3>, - FromVolumeFile>; + using type = std::array, 3>; static constexpr Options::String help = { - "Initial values for the translation map, its velocity and " - "acceleration."}; + "Initial values for the translation map. You can optionally specify " + "its first two time derivatives. If time derivatives aren't specified, " + "zero will be used."}; }; using options = tmpl::list; TranslationMapOptions() = default; // NOLINTNEXTLINE(google-explicit-constructor) - TranslationMapOptions(std::variant, 3>, - FromVolumeFile> - values_from_options); + TranslationMapOptions( + const std::array, 3>& initial_values_in, + const Options::Context& context = {}); std::array initial_values{}; }; + +/*! + * \brief Helper function that takes the variant of the translation map options, + * and returns the fully constructed translation function of time. + * + * \details Even if the function of time is read from a file, it will have a + * new \p initial_time and \p expiration_time. + */ +template +std::unique_ptr get_translation( + const std::variant, FromVolumeFile>& + translation_map_options, + double initial_time, double expiration_time); } // namespace domain::creators::time_dependent_options diff --git a/src/Domain/FunctionsOfTime/FunctionOfTime.hpp b/src/Domain/FunctionsOfTime/FunctionOfTime.hpp index 017fdef6ee51..4a58786d1bfc 100644 --- a/src/Domain/FunctionsOfTime/FunctionOfTime.hpp +++ b/src/Domain/FunctionsOfTime/FunctionOfTime.hpp @@ -49,6 +49,18 @@ class FunctionOfTime : public PUP::able { virtual auto get_clone() const -> std::unique_ptr = 0; + /// Create a FunctionOfTime at time \p t as if one had created a new + /// FunctionOfTime with \p t as its initial time, except the initial values of + /// new FunctionOfTime are the exact values of the old FunctionOfTime at time + /// \p t, and the new expriation is \p expiration_time. + /// + /// \details This defaults to just `get_clone()` since some functions can't be + /// updated/don't expire. + virtual std::unique_ptr create_at_time( + const double /*t*/, const double /*expiration_time*/) const { + return get_clone(); + } + /// Returns the domain of validity of the function. /// For FunctionsOfTime that allow a small amount of time extrapolation, /// `time_bounds` tells you the bounds including the allowed extrapolation diff --git a/src/Domain/FunctionsOfTime/IntegratedFunctionOfTime.cpp b/src/Domain/FunctionsOfTime/IntegratedFunctionOfTime.cpp index ea50ba15705e..19bfea9a5b4f 100644 --- a/src/Domain/FunctionsOfTime/IntegratedFunctionOfTime.cpp +++ b/src/Domain/FunctionsOfTime/IntegratedFunctionOfTime.cpp @@ -55,6 +55,13 @@ std::unique_ptr IntegratedFunctionOfTime::get_clone() const { return std::make_unique(*this); } +std::unique_ptr IntegratedFunctionOfTime::create_at_time( + const double t, const double expiration_time) const { + const auto initial_func_and_deriv = deriv_info_at_update_times_(t); + return std::make_unique( + t, initial_func_and_deriv.data, expiration_time, rotation_); +} + template std::array IntegratedFunctionOfTime::func_and_derivs(const double t) const { diff --git a/src/Domain/FunctionsOfTime/IntegratedFunctionOfTime.hpp b/src/Domain/FunctionsOfTime/IntegratedFunctionOfTime.hpp index f1a610a318f2..85a25a336a54 100644 --- a/src/Domain/FunctionsOfTime/IntegratedFunctionOfTime.hpp +++ b/src/Domain/FunctionsOfTime/IntegratedFunctionOfTime.hpp @@ -50,6 +50,9 @@ class IntegratedFunctionOfTime : public FunctionOfTime { auto get_clone() const -> std::unique_ptr override; + std::unique_ptr create_at_time( + double t, double expiration_time) const override; + std::array func(double t) const override; std::array func_and_deriv(double t) const override; [[noreturn]] std::array func_and_2_derivs( diff --git a/src/Domain/FunctionsOfTime/PiecewisePolynomial.cpp b/src/Domain/FunctionsOfTime/PiecewisePolynomial.cpp index 44ec4af1fe75..c68197ad897e 100644 --- a/src/Domain/FunctionsOfTime/PiecewisePolynomial.cpp +++ b/src/Domain/FunctionsOfTime/PiecewisePolynomial.cpp @@ -69,6 +69,13 @@ std::unique_ptr PiecewisePolynomial::get_clone() return std::make_unique(*this); } +template +std::unique_ptr PiecewisePolynomial::create_at_time( + const double t, const double expiration_time) const { + return std::make_unique(t, func_and_derivs(t), + expiration_time); +} + template template std::array diff --git a/src/Domain/FunctionsOfTime/PiecewisePolynomial.hpp b/src/Domain/FunctionsOfTime/PiecewisePolynomial.hpp index 93871e71ca16..f8c50e6874b0 100644 --- a/src/Domain/FunctionsOfTime/PiecewisePolynomial.hpp +++ b/src/Domain/FunctionsOfTime/PiecewisePolynomial.hpp @@ -44,6 +44,9 @@ class PiecewisePolynomial : public FunctionOfTime { auto get_clone() const -> std::unique_ptr override; + std::unique_ptr create_at_time( + double t, double expiration_time) const override; + // clang-tidy: google-runtime-references // clang-tidy: cppcoreguidelines-owning-memory,-warnings-as-errors WRAPPED_PUPable_decl_template(PiecewisePolynomial); // NOLINT @@ -71,6 +74,12 @@ class PiecewisePolynomial : public FunctionOfTime { /// an arbitrary time `t`. std::vector func_and_all_derivs(double t) const override; + /// Returns the function and `MaxDerivReturned` derivatives at + /// an arbitrary time `t`. + /// The function has multiple components. + template + std::array func_and_derivs(double t) const; + /// 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. @@ -103,12 +112,6 @@ class PiecewisePolynomial : public FunctionOfTime { void unpack_old_version(PUP::er& p, size_t version); - /// Returns the function and `MaxDerivReturned` derivatives at - /// an arbitrary time `t`. - /// The function has multiple components. - template - std::array func_and_derivs(double t) const; - void store_entry(double time_of_update, std::array func_and_derivs, double next_expiration_time); diff --git a/src/Domain/FunctionsOfTime/QuaternionFunctionOfTime.cpp b/src/Domain/FunctionsOfTime/QuaternionFunctionOfTime.cpp index a9a91ccd9ef6..1b7edd586143 100644 --- a/src/Domain/FunctionsOfTime/QuaternionFunctionOfTime.cpp +++ b/src/Domain/FunctionsOfTime/QuaternionFunctionOfTime.cpp @@ -32,6 +32,7 @@ #include "Utilities/ErrorHandling/Error.hpp" #include "Utilities/GenerateInstantiations.hpp" #include "Utilities/Gsl.hpp" +#include "Utilities/MakeArray.hpp" #include "Utilities/Serialization/PupBoost.hpp" #include "Utilities/StdHelpers.hpp" @@ -79,6 +80,14 @@ std::unique_ptr QuaternionFunctionOfTime::get_clone() return std::make_unique(*this); } +template +std::unique_ptr +QuaternionFunctionOfTime::create_at_time( + const double t, const double expiration_time) const { + return std::make_unique( + t, func(t), angle_f_of_t_.func_and_derivs(t), expiration_time); +} + template std::array QuaternionFunctionOfTime::time_bounds() const { return std::array{stored_quaternions_and_times_.initial_time(), diff --git a/src/Domain/FunctionsOfTime/QuaternionFunctionOfTime.hpp b/src/Domain/FunctionsOfTime/QuaternionFunctionOfTime.hpp index 97bdf0f99175..bd1890ad8170 100644 --- a/src/Domain/FunctionsOfTime/QuaternionFunctionOfTime.hpp +++ b/src/Domain/FunctionsOfTime/QuaternionFunctionOfTime.hpp @@ -42,6 +42,12 @@ namespace domain::FunctionsOfTime { * around the internal `PiecewisePolynomial::update` function with the addition * that it then updates the stored quaternions as well. * + * \note The initial rotation angles passed to the angle PiecewisePolynomial + * don't matter as we never actually use the angles themselves. We only use + * their derivatives (angular velocity) to determine map parameters. In theory + * we could determine each initial angle from the input axis-angle + * representation, but we don't need to. + * * The angle PiecewisePolynomial is accessible through the `angle_func`, * `angle_func_and_deriv`, and `angle_func_and_2_derivs` functions which * correspond to the function calls of a normal PiecewisePolynomial except @@ -79,6 +85,9 @@ class QuaternionFunctionOfTime : public FunctionOfTime { auto get_clone() const -> std::unique_ptr override; + std::unique_ptr create_at_time( + double t, double expiration_time) const override; + // clang-tidy: google-runtime-references // clang-tidy: cppcoreguidelines-owning-memory,-warnings-as-errors WRAPPED_PUPable_decl_template(QuaternionFunctionOfTime); // NOLINT diff --git a/src/Evolution/Ringdown/StrahlkorperCoefsInRingdownDistortedFrame.cpp b/src/Evolution/Ringdown/StrahlkorperCoefsInRingdownDistortedFrame.cpp index 58271989fbe3..bc25619cf100 100644 --- a/src/Evolution/Ringdown/StrahlkorperCoefsInRingdownDistortedFrame.cpp +++ b/src/Evolution/Ringdown/StrahlkorperCoefsInRingdownDistortedFrame.cpp @@ -10,6 +10,8 @@ #include "DataStructures/Matrix.hpp" #include "Domain/CoordinateMaps/Distribution.hpp" #include "Domain/Creators/Sphere.hpp" +#include "Domain/Creators/TimeDependentOptions/ExpansionMap.hpp" +#include "Domain/Creators/TimeDependentOptions/RotationMap.hpp" #include "Domain/Creators/TimeDependentOptions/Sphere.hpp" #include "Domain/StrahlkorperTransformations.hpp" #include "IO/H5/Dat.hpp" @@ -26,7 +28,7 @@ std::vector strahlkorper_coefs_in_ringdown_distorted_frame( const double settling_timescale, const std::array& exp_func_and_2_derivs, const std::array& exp_outer_bdry_func_and_2_derivs, - const std::array, 3>& rot_func_and_2_derivs) { + const std::vector>& rot_func_and_2_derivs) { // Read the AhC coefficients from the H5 file const std::vector>& ahc_inertial_h5 = ylm::read_surface_ylm( @@ -47,11 +49,11 @@ std::vector strahlkorper_coefs_in_ringdown_distorted_frame( // 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::ExpansionMapOptions + const domain::creators::time_dependent_options::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 + const domain::creators::time_dependent_options::RotationMapOptions rotation_map_options{rot_func_and_2_derivs, settling_timescale}; const domain::creators::sphere::TimeDependentMapOptions time_dependent_map_options{match_time, std::nullopt, diff --git a/src/Evolution/Ringdown/StrahlkorperCoefsInRingdownDistortedFrame.hpp b/src/Evolution/Ringdown/StrahlkorperCoefsInRingdownDistortedFrame.hpp index 835a1a7ee303..ea6c83b7d109 100644 --- a/src/Evolution/Ringdown/StrahlkorperCoefsInRingdownDistortedFrame.hpp +++ b/src/Evolution/Ringdown/StrahlkorperCoefsInRingdownDistortedFrame.hpp @@ -39,5 +39,5 @@ std::vector strahlkorper_coefs_in_ringdown_distorted_frame( double settling_timescale, const std::array& exp_func_and_2_derivs, const std::array& exp_outer_bdry_func_and_2_derivs, - const std::array, 3>& rot_func_and_2_derivs); + const std::vector>& rot_func_and_2_derivs); } // namespace evolution::Ringdown diff --git a/support/Pipelines/Bbh/Inspiral.yaml b/support/Pipelines/Bbh/Inspiral.yaml index 56d71b049c96..842e602340d4 100644 --- a/support/Pipelines/Bbh/Inspiral.yaml +++ b/support/Pipelines/Bbh/Inspiral.yaml @@ -118,9 +118,9 @@ DomainCreator: TimeDependentMaps: InitialTime: &InitialTime 0.0 ExpansionMap: - InitialValues: [1.0, {{ RadialExpansionVelocity }}] + InitialValues: [1.0, {{ RadialExpansionVelocity }}, 0.0] AsymptoticVelocityOuterBoundary: -1.0e-6 - DecayTimescaleOuterBoundaryVelocity: 50.0 + DecayTimescaleOuterBoundary: 50.0 RotationMap: InitialAngularVelocity: [0.0, 0.0, {{ InitialAngularVelocity }}] TranslationMap: diff --git a/support/Pipelines/Bbh/Ringdown.yaml b/support/Pipelines/Bbh/Ringdown.yaml index 2b062228ace8..11c8dcfe9ac1 100644 --- a/support/Pipelines/Bbh/Ringdown.yaml +++ b/support/Pipelines/Bbh/Ringdown.yaml @@ -65,13 +65,13 @@ DomainCreator: # and has not been tested for other configurations. SizeInitialValues: [0.0, -1.0, 0.0] RotationMap: - InitialValues: + InitialQuaternions: - [{{ Rotation0 }}, {{ Rotation1 }}, {{ Rotation2 }}, {{ Rotation3 }}] - [{{ dtRotation0 }}, {{ dtRotation1 }}, {{ dtRotation2 }}, {{ dtRotation3 }}] - [{{ dt2Rotation0 }}, {{ dt2Rotation1 }}, {{ dt2Rotation2}}, {{ dt2Rotation3 }}] - DecayTimescaleRotation: 1.0 + DecayTimescale: 1.0 ExpansionMap: # During the ringdown, only the outer boundary's expansion map is # applied. @@ -79,8 +79,8 @@ DomainCreator: InitialValuesOuterBoundary: [{{ ExpansionOuterBdry }}, {{ dtExpansionOuterBdry }}, {{ dt2ExpansionOuterBdry }}] - DecayTimescaleExpansion: 1.0 - DecayTimescaleExpansionOuterBoundary: 1.0 + DecayTimescale: 1.0 + DecayTimescaleOuterBoundary: 1.0 TranslationMap: InitialValues: [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]] TransitionRotScaleTrans: True diff --git a/tests/InputFiles/CurvedScalarWave/WorldtubeKerrSchild.yaml b/tests/InputFiles/CurvedScalarWave/WorldtubeKerrSchild.yaml index 3023e73d366c..9db35ba20b24 100644 --- a/tests/InputFiles/CurvedScalarWave/WorldtubeKerrSchild.yaml +++ b/tests/InputFiles/CurvedScalarWave/WorldtubeKerrSchild.yaml @@ -72,9 +72,9 @@ DomainCreator: TimeDependentMaps: InitialTime: 0. ExpansionMap: - InitialValues: [1., 0.] + InitialValues: [1., 0., 0.] AsymptoticVelocityOuterBoundary: 0. - DecayTimescaleOuterBoundaryVelocity: 1. + DecayTimescaleOuterBoundary: 1. RotationMap: InitialAngularVelocity: [0., 0., 0.08944271909999159] TranslationMap: diff --git a/tests/InputFiles/ExportCoordinates/InputTimeDependent3D.yaml b/tests/InputFiles/ExportCoordinates/InputTimeDependent3D.yaml index 6fef2193cac1..63038802510e 100644 --- a/tests/InputFiles/ExportCoordinates/InputTimeDependent3D.yaml +++ b/tests/InputFiles/ExportCoordinates/InputTimeDependent3D.yaml @@ -69,9 +69,9 @@ DomainCreator: TimeDependentMaps: InitialTime: 0.0 ExpansionMap: - InitialValues: [1.0, -4.6148457646200002e-05] + InitialValues: [1.0, -4.6148457646200002e-05, 0.0] AsymptoticVelocityOuterBoundary: -1.0e-6 - DecayTimescaleOuterBoundaryVelocity: 50.0 + DecayTimescaleOuterBoundary: 50.0 RotationMap: InitialAngularVelocity: [0.0, 0.0, 1.5264577062000000e-02] TranslationMap: diff --git a/tests/InputFiles/GeneralizedHarmonic/CylindricalBinaryBlackHole.yaml b/tests/InputFiles/GeneralizedHarmonic/CylindricalBinaryBlackHole.yaml index ea8f4f8ff7be..faf14f430912 100644 --- a/tests/InputFiles/GeneralizedHarmonic/CylindricalBinaryBlackHole.yaml +++ b/tests/InputFiles/GeneralizedHarmonic/CylindricalBinaryBlackHole.yaml @@ -58,9 +58,9 @@ DomainCreator: TimeDependentMaps: InitialTime: 0.0 ExpansionMap: - InitialValues: [1.0, -4.6148457646200002e-05] + InitialValues: [1.0, -4.6148457646200002e-05, 0.0] AsymptoticVelocityOuterBoundary: -1.0e-6 - DecayTimescaleOuterBoundaryVelocity: 50.0 + DecayTimescaleOuterBoundary: 50.0 RotationMap: InitialAngularVelocity: [0.0, 0.0, 1.5264577062000000e-02] TranslationMap: diff --git a/tests/InputFiles/GrMhd/GhValenciaDivClean/BinaryNeutronStar.yaml b/tests/InputFiles/GrMhd/GhValenciaDivClean/BinaryNeutronStar.yaml index effc521bc814..f3cf431fd8ba 100644 --- a/tests/InputFiles/GrMhd/GhValenciaDivClean/BinaryNeutronStar.yaml +++ b/tests/InputFiles/GrMhd/GhValenciaDivClean/BinaryNeutronStar.yaml @@ -64,9 +64,9 @@ DomainCreator: TimeDependentMaps: InitialTime: 0.0 ExpansionMap: - InitialValues: [1.0, 0.0] + InitialValues: [1.0, 0.0, 0.0] AsymptoticVelocityOuterBoundary: -1.0e-6 - DecayTimescaleOuterBoundaryVelocity: 50.0 + DecayTimescaleOuterBoundary: 50.0 RotationMap: InitialAngularVelocity: [0.0, 0.0, 0.00826225] TranslationMap: diff --git a/tests/Unit/ControlSystem/Test_Measurements.cpp b/tests/Unit/ControlSystem/Test_Measurements.cpp index 679c7966adfd..58441565a9b5 100644 --- a/tests/Unit/ControlSystem/Test_Measurements.cpp +++ b/tests/Unit/ControlSystem/Test_Measurements.cpp @@ -18,6 +18,7 @@ #include "Domain/Creators/Tags/Domain.hpp" #include "Domain/Creators/Tags/FunctionsOfTime.hpp" #include "Domain/Creators/TimeDependentOptions/BinaryCompactObject.hpp" +#include "Domain/Creators/TimeDependentOptions/RotationMap.hpp" #include "Domain/FunctionsOfTime/RegisterDerivedWithCharm.hpp" #include "Domain/Structure/ObjectLabel.hpp" #include "Framework/ActionTesting.hpp" @@ -323,7 +324,7 @@ SPECTRE_TEST_CASE("Unit.ControlSystem.FindTwoCenters", const domain::creators::bco::TimeDependentMapOptions time_dep_opts{ 0.0, std::nullopt, - domain::creators::bco::TimeDependentMapOptions::RotationMapOptions{ + domain::creators::time_dependent_options::RotationMapOptions{ std::array{0.0, 0.0, 0.1}}, std::nullopt, std::nullopt, diff --git a/tests/Unit/Domain/Creators/Test_BinaryCompactObject.cpp b/tests/Unit/Domain/Creators/Test_BinaryCompactObject.cpp index fef583f9dc30..835d5255f8f0 100644 --- a/tests/Unit/Domain/Creators/Test_BinaryCompactObject.cpp +++ b/tests/Unit/Domain/Creators/Test_BinaryCompactObject.cpp @@ -33,6 +33,7 @@ #include "Domain/FunctionsOfTime/PiecewisePolynomial.hpp" #include "Domain/FunctionsOfTime/QuaternionFunctionOfTime.hpp" #include "Domain/Structure/BlockNeighbor.hpp" +#include "Domain/Structure/ObjectLabel.hpp" #include "Framework/TestCreation.hpp" #include "Helpers/DataStructures/MakeWithRandomValues.hpp" #include "Helpers/Domain/BoundaryConditions/BoundaryCondition.hpp" @@ -404,9 +405,9 @@ std::string create_option_string( ? " TimeDependentMaps:\n" " InitialTime: 1.0\n" " ExpansionMap: \n" - " InitialValues: [1.0, -0.1]\n" + " InitialValues: [1.0, -0.1, 0.0]\n" " AsymptoticVelocityOuterBoundary: -0.1\n" - " DecayTimescaleOuterBoundaryVelocity: 5.0\n" + " DecayTimescaleOuterBoundary: 5.0\n" " RotationMap:\n" " InitialAngularVelocity: [0.0, 0.0, -0.2]\n" " TranslationMap:\n" @@ -839,6 +840,10 @@ void test_parse_errors() { // test_connectivity function. } +template +using HardcodedShape = + domain::creators::time_dependent_options::ShapeMapOptions; + void test_kerr_horizon_conforming() { INFO( "Check that inner radius is deformed to constant Boyer-Lindquist radius"); @@ -868,18 +873,13 @@ void test_kerr_horizon_conforming() { Distribution::Inverse, 120., domain::creators::bco::TimeDependentMapOptions{ - 0., - std::nullopt, - std::nullopt, - std::nullopt, - {{32_st, - domain::creators::time_dependent_options:: - KerrSchildFromBoyerLindquist{mass_A, spin_A}, - std::nullopt}}, - {{32_st, - domain::creators::time_dependent_options:: - KerrSchildFromBoyerLindquist{mass_B, spin_B}, - std::nullopt}}}}; + 0., std::nullopt, std::nullopt, std::nullopt, + HardcodedShape{ + 32_st, domain::creators::time_dependent_options:: + KerrSchildFromBoyerLindquist{mass_A, spin_A}}, + HardcodedShape{ + 32_st, domain::creators::time_dependent_options:: + KerrSchildFromBoyerLindquist{mass_B, spin_B}}}}; const auto domain = domain_creator.create_domain(); const auto functions_of_time = domain_creator.functions_of_time(); // Set up coordinates on an ellipsoid of constant Boyer-Lindquist radius diff --git a/tests/Unit/Domain/Creators/Test_CylindricalBinaryCompactObject.cpp b/tests/Unit/Domain/Creators/Test_CylindricalBinaryCompactObject.cpp index 7609cb76bf8a..7d97582dea65 100644 --- a/tests/Unit/Domain/Creators/Test_CylindricalBinaryCompactObject.cpp +++ b/tests/Unit/Domain/Creators/Test_CylindricalBinaryCompactObject.cpp @@ -27,12 +27,15 @@ #include "Domain/Creators/DomainCreator.hpp" #include "Domain/Creators/OptionTags.hpp" #include "Domain/Creators/TimeDependentOptions/BinaryCompactObject.hpp" +#include "Domain/Creators/TimeDependentOptions/RotationMap.hpp" +#include "Domain/Creators/TimeDependentOptions/ShapeMap.hpp" #include "Domain/Domain.hpp" #include "Domain/ExcisionSphere.hpp" #include "Domain/FunctionsOfTime/FixedSpeedCubic.hpp" #include "Domain/FunctionsOfTime/PiecewisePolynomial.hpp" #include "Domain/FunctionsOfTime/QuaternionFunctionOfTime.hpp" #include "Domain/Structure/BlockNeighbor.hpp" +#include "Domain/Structure/ObjectLabel.hpp" #include "Framework/TestCreation.hpp" #include "Helpers/Domain/BoundaryConditions/BoundaryCondition.hpp" #include "Helpers/Domain/Creators/TestHelpers.hpp" @@ -400,16 +403,18 @@ TimeDepOptions construct_time_dependent_options() { return TimeDepOptions{ expected_time, std::nullopt, - TimeDepOptions::RotationMapOptions{{initial_angular_velocity[0], - initial_angular_velocity[1], - initial_angular_velocity[2]}}, + domain::creators::time_dependent_options::RotationMapOptions{ + std::array{initial_angular_velocity[0], initial_angular_velocity[1], + initial_angular_velocity[2]}}, std::nullopt, - TimeDepOptions::ShapeMapOptions{ + domain::creators::time_dependent_options::ShapeMapOptions< + false, domain::ObjectLabel::A>{ 8_st, std::nullopt, {{initial_size_A_coefs[0][0], initial_size_A_coefs[1][0], initial_size_A_coefs[1][0]}}}, - TimeDepOptions::ShapeMapOptions{ + domain::creators::time_dependent_options::ShapeMapOptions< + false, domain::ObjectLabel::B>{ 8_st, std::nullopt, {{initial_size_B_coefs[0][0], initial_size_B_coefs[1][0], @@ -477,7 +482,8 @@ void test_parse_errors() { 25.0, false, 1_st, 3_st, TimeDepOptions{ 0.0, std::nullopt, - TimeDepOptions::RotationMapOptions{std::array{0.0, 0.0, 0.0}}, + domain::creators::time_dependent_options::RotationMapOptions< + false>{std::array{0.0, 0.0, 0.0}}, std::nullopt, std::nullopt, std::nullopt}, create_inner_boundary_condition(), create_outer_boundary_condition(), Options::Context{false, {}, 1, 1}), diff --git a/tests/Unit/Domain/Creators/Test_Sphere.cpp b/tests/Unit/Domain/Creators/Test_Sphere.cpp index c6c627fb0c52..4c6e3ac61060 100644 --- a/tests/Unit/Domain/Creators/Test_Sphere.cpp +++ b/tests/Unit/Domain/Creators/Test_Sphere.cpp @@ -38,12 +38,14 @@ #include "Domain/Creators/TimeDependence/RegisterDerivedWithCharm.hpp" #include "Domain/Creators/TimeDependence/UniformTranslation.hpp" #include "Domain/Creators/TimeDependentOptions/ShapeMap.hpp" +#include "Domain/Creators/TimeDependentOptions/TranslationMap.hpp" #include "Domain/Domain.hpp" #include "Domain/ElementMap.hpp" #include "Domain/FunctionsOfTime/PiecewisePolynomial.hpp" #include "Domain/Structure/BlockNeighbor.hpp" #include "Domain/Structure/Direction.hpp" #include "Domain/Structure/DirectionMap.hpp" +#include "Domain/Structure/ObjectLabel.hpp" #include "Domain/Structure/OrientationMap.hpp" #include "Framework/TestCreation.hpp" #include "Framework/TestHelpers.hpp" @@ -595,16 +597,18 @@ void test_sphere(const gsl::not_null gen) { if (time_dependent) { if (use_hard_coded_time_dep_options) { + using namespace domain::creators::time_dependent_options; // NOLINT translation_velocity = std::array{{0.001, -0.003, 0.005}}; - time_dependent_options = creators::sphere::TimeDependentMapOptions( + time_dependent_options = creators::sphere::TimeDependentMapOptions{ 1.0, - creators::sphere::TimeDependentMapOptions::ShapeMapOptions{ - l_max, std::nullopt}, - std::nullopt, std::nullopt, - creators::sphere::TimeDependentMapOptions::TranslationMapOptions{ - {std::array{0.0, 0.0, 0.0}, translation_velocity, - std::array{0.0, 0.0, 0.0}}}, - false); + ShapeMapOptions{l_max, + std::nullopt}, + std::nullopt, + std::nullopt, + TranslationMapOptions<3>{std::array{ + std::array{0.0, 0.0, 0.0}, translation_velocity, + std::array{0.0, 0.0, 0.0}}}, + false}; } else { time_dependent_options = std::make_unique< domain::creators::time_dependence::UniformTranslation<3, 0>>( @@ -717,8 +721,10 @@ void test_shape_distortion() { time_dependent_options = domain::creators::sphere::TimeDependentMapOptions{ time, - {{l_max, domain::creators::time_dependent_options:: - KerrSchildFromBoyerLindquist{mass, spin}}}, + domain::creators::time_dependent_options::ShapeMapOptions< + false, domain::ObjectLabel::None>{ + l_max, domain::creators::time_dependent_options:: + KerrSchildFromBoyerLindquist{mass, spin}}, std::nullopt, std::nullopt, std::nullopt, @@ -750,16 +756,18 @@ void test_shape_distortion() { time_dependent_options = domain::creators::sphere::TimeDependentMapOptions{ time, - {{l_max, - domain::creators::time_dependent_options::YlmsFromFile{ - h5_filename, std::vector{subfile_name}, time, std::nullopt, - false}, - // Constructing a strahlkorper from collocation radii will not exactly - // match the collocation points (see Strahlkorper constructor docs). - // For this reason, the 00 coef we calculate is not exact. This is why - // we just hard code the proper value here. If you change l_max, this - // value must also change. - std::array{-4.6442771561420703730e-01, 0.0, 0.0}}}, + domain::creators::time_dependent_options::ShapeMapOptions< + false, domain::ObjectLabel::None>{ + l_max, + domain::creators::time_dependent_options::YlmsFromFile{ + h5_filename, std::vector{subfile_name}, time, std::nullopt, + false}, + // Constructing a strahlkorper from collocation radii will not + // exactly match the collocation points (see Strahlkorper + // constructor docs). For this reason, the 00 coef we calculate is + // not exact. This is why we just hard code the proper value here. + // If you change l_max, this value must also change. + std::array{-4.6442771561420703730e-01, 0.0, 0.0}}, std::nullopt, std::nullopt, std::nullopt, diff --git a/tests/Unit/Domain/Creators/TimeDependentOptions/Test_BinaryCompactObject.cpp b/tests/Unit/Domain/Creators/TimeDependentOptions/Test_BinaryCompactObject.cpp index 409fedd0aaaa..29f6fc7c3b89 100644 --- a/tests/Unit/Domain/Creators/TimeDependentOptions/Test_BinaryCompactObject.cpp +++ b/tests/Unit/Domain/Creators/TimeDependentOptions/Test_BinaryCompactObject.cpp @@ -12,11 +12,15 @@ #include "DataStructures/DataVector.hpp" #include "Domain/Creators/BinaryCompactObject.hpp" #include "Domain/Creators/TimeDependentOptions/BinaryCompactObject.hpp" +#include "Domain/Creators/TimeDependentOptions/ExpansionMap.hpp" +#include "Domain/Creators/TimeDependentOptions/RotationMap.hpp" +#include "Domain/Creators/TimeDependentOptions/ShapeMap.hpp" #include "Domain/FunctionsOfTime/FixedSpeedCubic.hpp" #include "Domain/FunctionsOfTime/FunctionOfTime.hpp" #include "Domain/FunctionsOfTime/IntegratedFunctionOfTime.hpp" #include "Domain/FunctionsOfTime/PiecewisePolynomial.hpp" #include "Domain/FunctionsOfTime/QuaternionFunctionOfTime.hpp" +#include "Domain/Structure/ObjectLabel.hpp" #include "NumericalAlgorithms/SphericalHarmonics/Spherepack.hpp" #include "Utilities/CartesianProduct.hpp" #include "Utilities/Gsl.hpp" @@ -31,19 +35,19 @@ struct Inertial; namespace domain::creators::bco { template using ExpMapOptions = - typename TimeDependentMapOptions::ExpansionMapOptions; + typename TimeDependentMapOptions::ExpansionMapOptionType; template using RotMapOptions = - typename TimeDependentMapOptions::RotationMapOptions; + typename TimeDependentMapOptions::RotationMapOptionType; template using TransMapOptions = - typename TimeDependentMapOptions::TranslationMapOptions; + typename TimeDependentMapOptions::TranslationMapOptionType; template using ShapeMapAOptions = typename TimeDependentMapOptions< - IsCylindrical>::template ShapeMapOptions; + IsCylindrical>::template ShapeMapOptionType; template using ShapeMapBOptions = typename TimeDependentMapOptions< - IsCylindrical>::template ShapeMapOptions; + IsCylindrical>::template ShapeMapOptionType; namespace { // Test produce_all_maps for 1-4 maps using Expansion = domain::CoordinateMaps::TimeDependent::CubicScale<3>; @@ -120,23 +124,24 @@ void test(const bool include_expansion, const bool include_rotation, CAPTURE(include_translation); CAPTURE(include_shape_a); CAPTURE(include_shape_b); - std::optional> exp_map_options{}; - std::optional> rot_map_options{}; - std::optional> trans_map_options{}; - std::optional> shape_map_a_options{}; - std::optional> shape_map_b_options{}; + ExpMapOptions exp_map_options{}; + RotMapOptions rot_map_options{}; + TransMapOptions trans_map_options{}; + ShapeMapAOptions shape_map_a_options{}; + ShapeMapBOptions shape_map_b_options{}; - const std::array exp_values{1.0, 0.0}; + const std::array exp_values{1.0, 0.0, 0.0}; const double exp_outer_boundary_velocity = -0.01; const double exp_outer_boundary_timescale = 25.0; if (include_expansion) { - exp_map_options = ExpMapOptions{ - exp_values, exp_outer_boundary_velocity, exp_outer_boundary_timescale}; + exp_map_options = time_dependent_options::ExpansionMapOptions{ + exp_values, exp_outer_boundary_timescale, exp_outer_boundary_velocity}; } const std::array angular_velocity{0.2, -0.4, 0.6}; if (include_rotation) { - rot_map_options = RotMapOptions{angular_velocity}; + rot_map_options = + time_dependent_options::RotationMapOptions{angular_velocity}; } const std::array, 3> translation_values = { @@ -150,22 +155,30 @@ void test(const bool include_expansion, const bool include_rotation, const size_t l_max_A = 8; if (include_shape_a) { shape_map_a_options = - IsCylindrical ? ShapeMapAOptions{l_max_A, std::nullopt, - size_A_values} - : ShapeMapAOptions{ - l_max_A, std::nullopt, size_A_values, - transition_ends_at_cube_A}; + IsCylindrical + ? time_dependent_options::ShapeMapOptions< + not IsCylindrical, domain::ObjectLabel::A>{l_max_A, + std::nullopt, + size_A_values} + : time_dependent_options::ShapeMapOptions{ + l_max_A, std::nullopt, size_A_values, + transition_ends_at_cube_A}; } const std::array size_B_values{-0.001, -0.02, -0.3}; const size_t l_max_B = 10; if (include_shape_b) { shape_map_b_options = - IsCylindrical ? ShapeMapBOptions{l_max_B, std::nullopt, - size_B_values} - : ShapeMapBOptions{ - l_max_B, std::nullopt, size_B_values, - transition_ends_at_cube_B}; + IsCylindrical + ? time_dependent_options::ShapeMapOptions< + not IsCylindrical, domain::ObjectLabel::B>{l_max_B, + std::nullopt, + size_B_values} + : time_dependent_options::ShapeMapOptions{ + l_max_B, std::nullopt, size_B_values, + transition_ends_at_cube_B}; } const double initial_time = 1.5; @@ -221,27 +234,17 @@ void test(const bool include_expansion, const bool include_rotation, ExpBdryFoT expansion_outer_boundary{1.0, initial_time, exp_outer_boundary_velocity, exp_outer_boundary_timescale}; - RotFoT rotation{initial_time, - std::array{DataVector{1.0, 0.0, 0.0, 0.0}}, - std::array{{{3, 0.0}, - {gsl::at(angular_velocity, 0), - gsl::at(angular_velocity, 1), - gsl::at(angular_velocity, 2)}, - {3, 0.0}, - {3, 0.0}}}, - expiration_times.at( - TimeDependentMapOptions::rotation_name)}; + RotFoT rotation{ + initial_time, std::array{DataVector{1.0, 0.0, 0.0, 0.0}}, + std::array{ + {{3, 0.0}, DataVector{angular_velocity}, {3, 0.0}, {3, 0.0}}}, + expiration_times.at( + TimeDependentMapOptions::rotation_name)}; TransFoT translation{ initial_time, - std::array{{{gsl::at(translation_values, 0)[0], - gsl::at(translation_values, 0)[1], - gsl::at(translation_values, 0)[2]}, - {gsl::at(translation_values, 1)[0], - gsl::at(translation_values, 1)[1], - gsl::at(translation_values, 1)[2]}, - {gsl::at(translation_values, 2)[0], - gsl::at(translation_values, 2)[1], - gsl::at(translation_values, 2)[2]}}}, + std::array{DataVector{translation_values[0]}, + DataVector{translation_values[1]}, + DataVector{translation_values[2]}}, expiration_times.at( TimeDependentMapOptions::translation_name)}; SizeFoT size_A{initial_time, @@ -301,60 +304,59 @@ void test(const bool include_expansion, const bool include_rotation, RadiiType inner_outer_radii_A{}; RadiiType inner_outer_radii_B{}; - if constexpr (IsCylindrical) { - inner_outer_radii_A = std::array{0.8, 3.2}; - } else { - inner_outer_radii_A = std::array{0.8, 1.4, 3.2}; + if constexpr (IsCylindrical) { + inner_outer_radii_A = std::array{0.8, 3.2}; + } else { + inner_outer_radii_A = std::array{0.8, 1.4, 3.2}; + } + if constexpr (IsCylindrical) { + inner_outer_radii_B = std::array{0.5, 2.1}; + } else { + inner_outer_radii_B = std::array{0.5, 0.9, 2.1}; + } + const auto build_maps = [&time_dep_options, ¢ers, &inner_outer_radii_A, + &inner_outer_radii_B, &domain_envelope_radius, + &domain_outer_radius, &cube_a_center, + &cube_b_center, filled_A = not excise_A, + filled_B = not excise_B]() { + time_dep_options.build_maps(centers, cube_a_center, cube_b_center, + inner_outer_radii_A, inner_outer_radii_B, + filled_A, filled_B, domain_envelope_radius, + domain_outer_radius); + }; + + if constexpr (not IsCylindrical) { + if (include_shape_a and + (not excise_A and not transition_ends_at_cube_A)) { + CHECK_THROWS_WITH(build_maps(), + Catch::Matchers::ContainsSubstring( + "If the object is filled, the transition must " + "end at the cube.")); + continue; } - if constexpr (IsCylindrical) { - inner_outer_radii_B = std::array{0.5, 2.1}; - } else { - inner_outer_radii_B = std::array{0.5, 0.9, 2.1}; + if (include_shape_a and not excise_A) { + CHECK_THROWS_WITH( + build_maps(), + Catch::Matchers::ContainsSubstring( + "the centers of the inner and outer object are different")); + continue; + } + if (include_shape_b and + (not excise_B and not transition_ends_at_cube_B)) { + CHECK_THROWS_WITH(build_maps(), + Catch::Matchers::ContainsSubstring( + "If the object is filled, the transition must " + "end at the cube.")); + continue; } - const auto build_maps = [&time_dep_options, ¢ers, - &inner_outer_radii_A, &inner_outer_radii_B, - &domain_envelope_radius, &domain_outer_radius, - &cube_a_center, &cube_b_center, - filled_A = not excise_A, - filled_B = not excise_B]() { - time_dep_options.build_maps(centers, cube_a_center, cube_b_center, - inner_outer_radii_A, inner_outer_radii_B, - filled_A, filled_B, domain_envelope_radius, - domain_outer_radius); - }; - - if constexpr (not IsCylindrical) { - if (include_shape_a and - (not excise_A and not transition_ends_at_cube_A)) { - CHECK_THROWS_WITH(build_maps(), - Catch::Matchers::ContainsSubstring( - "If the object is filled, the transition must " - "end at the cube.")); - continue; - } - if (include_shape_a and not excise_A) { - CHECK_THROWS_WITH( - build_maps(), - Catch::Matchers::ContainsSubstring( - "the centers of the inner and outer object are different")); - continue; - } - if (include_shape_b and - (not excise_B and not transition_ends_at_cube_B)) { - CHECK_THROWS_WITH(build_maps(), - Catch::Matchers::ContainsSubstring( - "If the object is filled, the transition must " - "end at the cube.")); - continue; - } - if (include_shape_b and not excise_B) { - CHECK_THROWS_WITH( - build_maps(), - Catch::Matchers::ContainsSubstring( - "the centers of the inner and outer object are different")); - continue; - } + if (include_shape_b and not excise_B) { + CHECK_THROWS_WITH( + build_maps(), + Catch::Matchers::ContainsSubstring( + "the centers of the inner and outer object are different")); + continue; } + } // Now for each object, either we have included shape map options and // excision info, or we didn't include either. (i.e. include_shape_map_? @@ -577,14 +579,18 @@ void test_errors() { CHECK_THROWS_WITH( (TimeDependentMapOptions{ 1.0, std::nullopt, std::nullopt, std::nullopt, - ShapeMapAOptions{1, {}}, - ShapeMapBOptions{8, {}}}), + time_dependent_options::ShapeMapOptions< + not IsCylindrical, domain::ObjectLabel::A>{1, {}}, + time_dependent_options::ShapeMapOptions< + not IsCylindrical, domain::ObjectLabel::B>{8, {}}}), Catch::Matchers::ContainsSubstring("Initial LMax for object")); CHECK_THROWS_WITH( (TimeDependentMapOptions{ 1.0, std::nullopt, std::nullopt, std::nullopt, - ShapeMapAOptions{6, {}}, - ShapeMapBOptions{0, {}}}), + time_dependent_options::ShapeMapOptions< + not IsCylindrical, domain::ObjectLabel::A>{6, {}}, + time_dependent_options::ShapeMapOptions< + not IsCylindrical, domain::ObjectLabel::B>{0, {}}}), Catch::Matchers::ContainsSubstring("Initial LMax for object")); CHECK_THROWS_WITH((TimeDependentMapOptions{ 1.0, std::nullopt, std::nullopt, std::nullopt, @@ -604,10 +610,12 @@ void test_errors() { ([&radii]() { TimeDependentMapOptions time_dep_opts{ 1.0, - ExpMapOptions{{1.0, 0.0}, 0.0, 0.01}, - RotMapOptions{{0.0, 0.0, 0.0}}, std::nullopt, - ShapeMapAOptions{8, {}}, + std::nullopt, + std::nullopt, + time_dependent_options::ShapeMapOptions{ + 8, std::nullopt}, std::nullopt}; time_dep_opts.build_maps( std::array{std::array{5.0, 0.0, 0.0}, std::array{-5.0, 0.0, 0.0}}, @@ -637,12 +645,16 @@ void test_worldtube_fots() { const TimeDependentMapOptions worldtube_options{ initial_time, - ExpMapOptions{ - {initial_expansion, initial_expansion_deriv}, 1., 1.}, - RotMapOptions{{0.0, 0.0, 1.0}}, + time_dependent_options::ExpansionMapOptions{ + std::array{initial_expansion, initial_expansion_deriv, 0.0}, 1.0, + 1.0}, + time_dependent_options::RotationMapOptions{ + std::array{0.0, 0.0, 1.0}}, std::nullopt, - ShapeMapAOptions{2, {}, std::make_optional(size_a_opts)}, - ShapeMapBOptions{2, {}, std::make_optional(size_b_opts)}}; + time_dependent_options::ShapeMapOptions{ + 2, {}, std::make_optional(size_a_opts)}, + time_dependent_options::ShapeMapOptions{ + 2, {}, std::make_optional(size_b_opts)}}; const auto fots = worldtube_options.create_functions_of_time({}); CHECK(not fots.contains("Translation")); diff --git a/tests/Unit/Domain/Creators/TimeDependentOptions/Test_ExpansionMap.cpp b/tests/Unit/Domain/Creators/TimeDependentOptions/Test_ExpansionMap.cpp index 128294af4dcd..f64ff7c25f93 100644 --- a/tests/Unit/Domain/Creators/TimeDependentOptions/Test_ExpansionMap.cpp +++ b/tests/Unit/Domain/Creators/TimeDependentOptions/Test_ExpansionMap.cpp @@ -5,6 +5,7 @@ #include #include +#include #include #include #include @@ -16,7 +17,9 @@ #include "Domain/FunctionsOfTime/FunctionOfTime.hpp" #include "Domain/FunctionsOfTime/PiecewisePolynomial.hpp" #include "Domain/FunctionsOfTime/RegisterDerivedWithCharm.hpp" +#include "Domain/FunctionsOfTime/SettleToConstant.hpp" #include "Framework/TestCreation.hpp" +#include "Helpers/Domain/Creators/TimeDependent/TestHelpers.hpp" #include "IO/H5/AccessType.hpp" #include "IO/H5/File.hpp" #include "IO/H5/TensorData.hpp" @@ -28,171 +31,230 @@ #include "Utilities/MakeArray.hpp" #include "Utilities/Serialization/Serialize.hpp" +namespace domain::creators::time_dependent_options { namespace { +constexpr double infinity = std::numeric_limits::infinity(); void test_expansion_map_options() { { - const auto expansion_map_options = TestHelpers::test_creation< - domain::creators::time_dependent_options::ExpansionMapOptions>( - "InitialValues: [1.0, 2.0, 3.0]\n" - "InitialValuesOuterBoundary: [4.0, 5.0, 6.0]\n" - "DecayTimescaleOuterBoundary: 50\n" - "DecayTimescale: Auto\n" - "AsymptoticVelocityOuterBoundary: -1e-5"); - - CHECK(expansion_map_options.name() == "ExpansionMap"); - CHECK(expansion_map_options.initial_values == - std::array{DataVector{1.0}, DataVector{2.0}, DataVector{3.0}}); - CHECK(expansion_map_options.initial_values_outer_boundary == - std::array{DataVector{4.0}, DataVector{5.0}, DataVector{6.0}}); - CHECK(expansion_map_options.decay_timescale_outer_boundary == 50.0); - CHECK_FALSE(expansion_map_options.decay_timescale.has_value()); - CHECK(expansion_map_options.asymptotic_velocity_outer_boundary.has_value()); - CHECK(expansion_map_options.asymptotic_velocity_outer_boundary.value() == - -1e-5); + INFO("None"); + const auto expansion_map_options = + TestHelpers::test_option_tag>("None"); + + CHECK(not expansion_map_options.has_value()); } + + const auto non_settle_fots = []() { + INFO("Hardcoded AllowSettleFots = " + std::to_string(AllowSettleFoTs) + + ", Non-Settle Fots"); + const auto expansion_map_options = + TestHelpers::test_option_tag>( + "InitialValues: [1.0, 2.0, 3.0]\n" + "DecayTimescaleOuterBoundary: 50\n" + "AsymptoticVelocityOuterBoundary: -1e-5"); + + REQUIRE(expansion_map_options.has_value()); + CHECK(std::holds_alternative>( + expansion_map_options.value())); + + const auto& hardcoded_options = + std::get>( + expansion_map_options.value()); + + const std::array expected_values{DataVector{1.0}, DataVector{2.0}, + DataVector{3.0}}; + CHECK(hardcoded_options.initial_values == expected_values); + CHECK(hardcoded_options.initial_values_outer_boundary == + std::array{DataVector{1.0}, DataVector{0.0}, DataVector{0.0}}); + CHECK(hardcoded_options.decay_timescale_outer_boundary == 50.0); + CHECK_FALSE(hardcoded_options.decay_timescale.has_value()); + CHECK(hardcoded_options.asymptotic_velocity_outer_boundary == + std::optional{-1.e-5}); + + const auto expansion_fots = + get_expansion(expansion_map_options.value(), 0.3, 2.9); + const auto& expansion_ptr = expansion_fots.at("Expansion"); + const auto& expansion_outer_boundary_ptr = + expansion_fots.at("ExpansionOuterBoundary"); + + const auto* expansion = + dynamic_cast*>( + expansion_ptr.get()); + const auto* expansion_outer_boundary = + dynamic_cast( + expansion_outer_boundary_ptr.get()); + + REQUIRE(expansion != nullptr); + REQUIRE(expansion_outer_boundary != nullptr); + + CHECK(expansion->time_bounds() == std::array{0.3, 2.9}); + CHECK(expansion->func_and_2_derivs(0.3) == expected_values); + CHECK(expansion_outer_boundary->time_bounds() == std::array{0.3, infinity}); + CHECK(expansion_outer_boundary->decay_timescale() == 50.0); + CHECK(expansion_outer_boundary->velocity() == -1e-5); + }; + + non_settle_fots.template operator()(); + non_settle_fots.template operator()(); + { - const auto expansion_map_options = TestHelpers::test_creation< - domain::creators::time_dependent_options::ExpansionMapOptions>( - "InitialValues: [1.0, 2.0, 3.0]\n" - "InitialValuesOuterBoundary: [4.0, 5.0, 6.0]\n" - "DecayTimescaleOuterBoundary: 50\n" - "DecayTimescale: 40\n" - "AsymptoticVelocityOuterBoundary: Auto"); - - CHECK(expansion_map_options.name() == "ExpansionMap"); - CHECK(expansion_map_options.initial_values == - std::array{DataVector{1.0}, DataVector{2.0}, DataVector{3.0}}); - CHECK(expansion_map_options.initial_values_outer_boundary == + INFO("Hardcoded AllowSettleFots = true, Settle Fots"); + const auto expansion_map_options = + TestHelpers::test_option_tag>( + "InitialValues: [1.0, 2.0, 3.0]\n" + "InitialValuesOuterBoundary: [4.0, 5.0, 6.0]\n" + "DecayTimescaleOuterBoundary: 50\n" + "DecayTimescale: 40\n"); + + REQUIRE(expansion_map_options.has_value()); + CHECK(std::holds_alternative>( + expansion_map_options.value())); + + const auto& hardcoded_options = + std::get>(expansion_map_options.value()); + + const std::array expected_values{DataVector{1.0}, DataVector{2.0}, + DataVector{3.0}}; + CHECK(hardcoded_options.initial_values == expected_values); + CHECK(hardcoded_options.initial_values_outer_boundary == std::array{DataVector{4.0}, DataVector{5.0}, DataVector{6.0}}); - CHECK(expansion_map_options.decay_timescale_outer_boundary == 50.0); - CHECK(expansion_map_options.decay_timescale.has_value()); - CHECK(expansion_map_options.decay_timescale.value() == 40.0); + CHECK(hardcoded_options.decay_timescale_outer_boundary == 50.0); + CHECK(hardcoded_options.decay_timescale == std::optional{40.0}); CHECK_FALSE( - expansion_map_options.asymptotic_velocity_outer_boundary.has_value()); + hardcoded_options.asymptotic_velocity_outer_boundary.has_value()); + + const auto expansion_fots = + get_expansion(expansion_map_options.value(), 0.3, 2.9); + const auto& expansion_ptr = expansion_fots.at("Expansion"); + const auto& expansion_outer_boundary_ptr = + expansion_fots.at("ExpansionOuterBoundary"); + + const auto* expansion = + dynamic_cast( + expansion_ptr.get()); + const auto* expansion_outer_boundary = + dynamic_cast( + expansion_outer_boundary_ptr.get()); + + REQUIRE(expansion != nullptr); + REQUIRE(expansion_outer_boundary != nullptr); + + CHECK(expansion->time_bounds() == std::array{0.3, infinity}); + CHECK(expansion->func_and_2_derivs(0.3) == expected_values); + CHECK(expansion_outer_boundary->time_bounds() == std::array{0.3, infinity}); } - CHECK_THROWS_WITH( - (TestHelpers::test_creation< - domain::creators::time_dependent_options::ExpansionMapOptions>( - "InitialValues: [1.0, 2.0, 3.0]\n" - "InitialValuesOuterBoundary: [4.0, 5.0, 6.0]\n" - "DecayTimescaleOuterBoundary: 50\n" - "DecayTimescale: Auto\n" - "AsymptoticVelocityOuterBoundary: Auto")), - Catch::Matchers::ContainsSubstring( - "must specify one of DecayTimescale or " - "AsymptoticVelocityOuterBoundary, but not both.")); - CHECK_THROWS_WITH( - (TestHelpers::test_creation< - domain::creators::time_dependent_options::ExpansionMapOptions>( - "InitialValues: [1.0, 2.0, 3.0]\n" - "InitialValuesOuterBoundary: [4.0, 5.0, 6.0]\n" - "DecayTimescaleOuterBoundary: 50\n" - "DecayTimescale: 40\n" - "AsymptoticVelocityOuterBoundary: -1e-5")), - Catch::Matchers::ContainsSubstring( - "must specify one of DecayTimescale or " - "AsymptoticVelocityOuterBoundary, but not both.")); - CHECK_THROWS_WITH( - (TestHelpers::test_creation< - domain::creators::time_dependent_options::ExpansionMapOptions>( - "InitialValues: [1.0, 2.0, 3.0]\n" - "InitialValuesOuterBoundary: [4.0, 5.0, 6.0]\n" - "DecayTimescaleOuterBoundary: Auto\n" - "DecayTimescale: 40\n" - "AsymptoticVelocityOuterBoundary: Auto")), - Catch::Matchers::ContainsSubstring( - "When specifying the ExpansionMap initial outer " - "boundary values directly, you must also specify a " - "'DecayTimescaleOuterBoundary'.")); + const std::string filename{"Commencement.h5"}; + const std::string subfile_name{"VolumeData"}; + + domain::FunctionsOfTimeMap functions_of_time{}; + functions_of_time["Expansion"] = + std::make_unique>( + 0.0, std::array{DataVector{1.0}, DataVector{2.0}, DataVector{3.0}}, + 100.0); + functions_of_time["ExpansionOuterBoundary"] = + std::make_unique(0.0, 0.0, + -1e-5, 50.0); { - std::unordered_map> - functions_of_time{}; - functions_of_time["Expansion"] = - std::make_unique>( - 0.0, std::array{DataVector{1.0}, DataVector{2.0}, DataVector{3.0}}, - 100.0); - functions_of_time["ExpansionOuterBoundary"] = - std::make_unique(0.0, 0.0, - -1e-5, 50.0); - const std::string filename{"Commencement.h5"}; - const std::string subfile_name{"VolumeData"}; + INFO("FromVolumeFile non-Settle FoTs"); + if (file_system::check_if_file_exists(filename)) { file_system::rm(filename, true); } - { - h5::H5File h5_file{filename}; - 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, serialize(functions_of_time)); - } + TestHelpers::domain::creators::write_volume_data(filename, subfile_name, + functions_of_time); - { - const auto expansion_map_options = TestHelpers::test_creation< - domain::creators::time_dependent_options::ExpansionMapOptions>( - "DecayTimescaleOuterBoundary: 60\n" - "DecayTimescale: Auto\n" - "AsymptoticVelocityOuterBoundary: -2e-5\n" - "InitialValues:\n" - " H5Filename: " + - filename + "\n SubfileName: " + subfile_name + - "\n Time: 0.0\n" - "InitialValuesOuterBoundary:\n" - " H5Filename: " + - filename + "\n SubfileName: " + subfile_name + "\n Time: 0.0"); - CHECK(expansion_map_options.name() == "ExpansionMap"); - CHECK(expansion_map_options.initial_values == - std::array{DataVector{1.0}, DataVector{2.0}, DataVector{3.0}}); - CHECK(expansion_map_options.initial_values_outer_boundary == - std::array{DataVector{0.0}, DataVector{0.0}, DataVector{0.0}}); - CHECK(expansion_map_options.decay_timescale_outer_boundary == 60.0); - CHECK_FALSE(expansion_map_options.decay_timescale.has_value()); - CHECK( - expansion_map_options.asymptotic_velocity_outer_boundary.has_value()); - CHECK(expansion_map_options.asymptotic_velocity_outer_boundary.value() == - -2e-5); - } - { - const auto expansion_map_options = TestHelpers::test_creation< - domain::creators::time_dependent_options::ExpansionMapOptions>( - "DecayTimescaleOuterBoundary: Auto\n" - "DecayTimescale: Auto\n" - "AsymptoticVelocityOuterBoundary: Auto\n" - "InitialValues:\n" - " H5Filename: " + - filename + "\n SubfileName: " + subfile_name + - "\n Time: 0.0\n" - "InitialValuesOuterBoundary:\n" - " H5Filename: " + - filename + "\n SubfileName: " + subfile_name + "\n Time: 0.0"); - CHECK(expansion_map_options.name() == "ExpansionMap"); - CHECK(expansion_map_options.initial_values == - std::array{DataVector{1.0}, DataVector{2.0}, DataVector{3.0}}); - CHECK(expansion_map_options.initial_values_outer_boundary == - std::array{DataVector{0.0}, DataVector{0.0}, DataVector{0.0}}); - CHECK(expansion_map_options.decay_timescale_outer_boundary == 50.0); - CHECK_FALSE(expansion_map_options.decay_timescale.has_value()); - CHECK( - expansion_map_options.asymptotic_velocity_outer_boundary.has_value()); - CHECK(expansion_map_options.asymptotic_velocity_outer_boundary.value() == - -1e-5); - } + const auto expansion_map_options = + TestHelpers::test_option_tag>( + "H5Filename: Commencement.h5\n" + "SubfileName: VolumeData"); + + REQUIRE(expansion_map_options.has_value()); + CHECK( + std::holds_alternative(expansion_map_options.value())); + + const auto expansion_fots = + get_expansion(expansion_map_options.value(), 0.3, 100.0); + const auto& expansion_ptr = expansion_fots.at("Expansion"); + const auto& expansion_outer_boundary_ptr = + expansion_fots.at("ExpansionOuterBoundary"); + + const auto* expansion = + dynamic_cast*>( + expansion_ptr.get()); + const auto* expansion_outer_boundary = + dynamic_cast( + expansion_outer_boundary_ptr.get()); + + REQUIRE(expansion != nullptr); + REQUIRE(expansion_outer_boundary != nullptr); + + CHECK(expansion->time_bounds() == std::array{0.3, 100.0}); + CHECK(expansion->func_and_2_derivs(0.3) == + functions_of_time.at("Expansion")->func_and_2_derivs(0.3)); + CHECK( + expansion_outer_boundary->func_and_2_derivs(0.3) == + functions_of_time.at("ExpansionOuterBoundary")->func_and_2_derivs(0.3)); + CHECK(expansion_outer_boundary->time_bounds() == std::array{0.0, infinity}); + CHECK(expansion_outer_boundary->decay_timescale() == 50.0); + CHECK(expansion_outer_boundary->velocity() == -1e-5); + } + { + INFO("FromVolumeFile Settle FoTs"); if (file_system::check_if_file_exists(filename)) { file_system::rm(filename, true); } + + functions_of_time["Expansion"] = + std::make_unique( + std::array{DataVector{1.0}, DataVector{2.0}, DataVector{3.0}}, 0.0, + 100.0); + functions_of_time["ExpansionOuterBoundary"] = + std::make_unique( + std::array{DataVector{3.0}, DataVector{4.0}, DataVector{6.0}}, 0.0, + 100.0); + + TestHelpers::domain::creators::write_volume_data(filename, subfile_name, + functions_of_time); + + const auto expansion_map_options = + TestHelpers::test_option_tag>( + "H5Filename: Commencement.h5\n" + "SubfileName: VolumeData"); + + REQUIRE(expansion_map_options.has_value()); + CHECK( + std::holds_alternative(expansion_map_options.value())); + + const auto expansion_fots = + get_expansion(expansion_map_options.value(), 0.3, 100.0); + const auto& expansion_ptr = expansion_fots.at("Expansion"); + const auto& expansion_outer_boundary_ptr = + expansion_fots.at("ExpansionOuterBoundary"); + + const auto* expansion = + dynamic_cast( + expansion_ptr.get()); + const auto* expansion_outer_boundary = + dynamic_cast( + expansion_outer_boundary_ptr.get()); + + REQUIRE(expansion != nullptr); + REQUIRE(expansion_outer_boundary != nullptr); + + CHECK(expansion->time_bounds() == std::array{0.0, infinity}); + CHECK(expansion->func_and_2_derivs(0.3) == + functions_of_time.at("Expansion")->func_and_2_derivs(0.3)); + CHECK( + expansion_outer_boundary->func_and_2_derivs(0.3) == + functions_of_time.at("ExpansionOuterBoundary")->func_and_2_derivs(0.3)); + CHECK(expansion_outer_boundary->time_bounds() == std::array{0.0, infinity}); + } + + if (file_system::check_if_file_exists(filename)) { + file_system::rm(filename, true); } } } // namespace @@ -202,3 +264,4 @@ SPECTRE_TEST_CASE("Unit.Domain.Creators.TimeDependentOptions.ExpansionMap", domain::FunctionsOfTime::register_derived_with_charm(); test_expansion_map_options(); } +} // namespace domain::creators::time_dependent_options diff --git a/tests/Unit/Domain/Creators/TimeDependentOptions/Test_FromVolumeFile.cpp b/tests/Unit/Domain/Creators/TimeDependentOptions/Test_FromVolumeFile.cpp index 1a4bce0f99d7..730fdfe24524 100644 --- a/tests/Unit/Domain/Creators/TimeDependentOptions/Test_FromVolumeFile.cpp +++ b/tests/Unit/Domain/Creators/TimeDependentOptions/Test_FromVolumeFile.cpp @@ -1,8 +1,6 @@ // Distributed under the MIT License. // See LICENSE.txt for details. -#include "DataStructures/DataVector.hpp" -#include "Domain/Structure/ObjectLabel.hpp" #include "Framework/TestingFramework.hpp" #include @@ -11,13 +9,16 @@ #include #include +#include "DataStructures/DataVector.hpp" #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 "Domain/Structure/ObjectLabel.hpp" #include "Framework/TestCreation.hpp" +#include "Helpers/Domain/Creators/TimeDependent/TestHelpers.hpp" #include "IO/H5/AccessType.hpp" #include "IO/H5/File.hpp" #include "IO/H5/TensorData.hpp" @@ -30,35 +31,12 @@ #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() { +void test(const std::string& function_of_time_name) { 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> @@ -70,212 +48,25 @@ void test() { DataVector{3, 0.0}}, 100.0); - write_volume_data(filename, subfile_name, functions_of_time); + TestHelpers::domain::creators::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"); + domain::creators::time_dependent_options::FromVolumeFile>( + "H5Filename: " + filename + "\nSubfileName: " + subfile_name); + + const domain::FunctionsOfTimeMap fot_from_file = + from_volume_file.retrieve_function_of_time({function_of_time_name}, time); 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); - CHECK(from_volume_file.velocity_outer_boundary == velocity); - CHECK(from_volume_file.decay_timescale_outer_boundary == decay_timescale); - - 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); + CHECK(dynamic_cast&>( + *fot_from_file.at(function_of_time_name)) == + dynamic_cast&>( + *functions_of_time.at(function_of_time_name))); if (file_system::check_if_file_exists(filename)) { file_system::rm(filename, true); @@ -295,24 +86,31 @@ void test_errors() { } using FromVolumeFile = - domain::creators::time_dependent_options::FromVolumeFile< - domain::creators::time_dependent_options::names::Expansion>; + domain::creators::time_dependent_options::FromVolumeFile; + + auto from_volume_file = TestHelpers::test_creation( + "H5Filename: " + filename + "\nSubfileName: " + subfile_name); CHECK_THROWS_WITH( - (FromVolumeFile{filename, subfile_name, 0.0}), + (from_volume_file.retrieve_function_of_time({"Expansion"}, 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); + TestHelpers::domain::creators::write_volume_data(filename, subfile_name); + + from_volume_file = TestHelpers::test_creation( + "H5Filename: " + filename + "\nSubfileName: " + subfile_name); CHECK_THROWS_WITH( - (FromVolumeFile{filename, subfile_name, 0.0}), + (from_volume_file.retrieve_function_of_time({"Expansion"}, 0.0)), Catch::Matchers::ContainsSubstring( "Expansion: There are no functions of time in the subfile ")); subfile_name += "0"; + from_volume_file = TestHelpers::test_creation( + "H5Filename: " + filename + "\nSubfileName: " + subfile_name); std::unordered_map> functions_of_time{}; @@ -323,13 +121,17 @@ void test_errors() { DataVector{3, 0.0}}, 100.0); - write_volume_data(filename, subfile_name, functions_of_time); + TestHelpers::domain::creators::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 ")); + CHECK_THROWS_WITH( + (from_volume_file.retrieve_function_of_time({"Expansion"}, 0.0)), + Catch::Matchers::ContainsSubstring( + "No function of time named Expansion in the subfile ")); subfile_name += "0"; + from_volume_file = TestHelpers::test_creation( + "H5Filename: " + filename + "\nSubfileName: " + subfile_name); functions_of_time.clear(); functions_of_time["Expansion"] = std::make_unique>( @@ -338,14 +140,20 @@ void test_errors() { DataVector{3, 0.0}}, 1.0); - write_volume_data(filename, subfile_name, functions_of_time); + TestHelpers::domain::creators::write_volume_data(filename, subfile_name, + functions_of_time); CHECK_THROWS_WITH( - (FromVolumeFile{filename, subfile_name, 10.0}), + (from_volume_file.retrieve_function_of_time({"Expansion"}, 10.0)), Catch::Matchers::ContainsSubstring("Expansion: The requested time") and Catch::Matchers::ContainsSubstring( "is out of the range of the function of time")); + // Check that this is ok to call + const auto function_of_time = + from_volume_file.retrieve_function_of_time({"Expansion"}, std::nullopt); + (void)function_of_time; + if (file_system::check_if_file_exists(filename)) { file_system::rm(filename, true); } @@ -354,11 +162,7 @@ void test_errors() { SPECTRE_TEST_CASE("Unit.Domain.Creators.TimeDependentOptions.FromVolumeFile", "[Unit][Domain]") { domain::FunctionsOfTime::register_derived_with_charm(); - test(); - test(); - test(); - test>(); + test("Translation"); test_errors(); } } // namespace diff --git a/tests/Unit/Domain/Creators/TimeDependentOptions/Test_RotationMap.cpp b/tests/Unit/Domain/Creators/TimeDependentOptions/Test_RotationMap.cpp index 695ce464a617..1b66cb61b4bf 100644 --- a/tests/Unit/Domain/Creators/TimeDependentOptions/Test_RotationMap.cpp +++ b/tests/Unit/Domain/Creators/TimeDependentOptions/Test_RotationMap.cpp @@ -1,11 +1,11 @@ // Distributed under the MIT License. // See LICENSE.txt for details. -#include "Domain/FunctionsOfTime/QuaternionFunctionOfTime.hpp" #include "Framework/TestingFramework.hpp" #include #include +#include #include #include #include @@ -14,9 +14,11 @@ #include "Domain/Creators/TimeDependentOptions/FromVolumeFile.hpp" #include "Domain/Creators/TimeDependentOptions/RotationMap.hpp" #include "Domain/FunctionsOfTime/FunctionOfTime.hpp" -#include "Domain/FunctionsOfTime/PiecewisePolynomial.hpp" +#include "Domain/FunctionsOfTime/QuaternionFunctionOfTime.hpp" #include "Domain/FunctionsOfTime/RegisterDerivedWithCharm.hpp" +#include "Domain/FunctionsOfTime/SettleToConstantQuaternion.hpp" #include "Framework/TestCreation.hpp" +#include "Helpers/Domain/Creators/TimeDependent/TestHelpers.hpp" #include "IO/H5/AccessType.hpp" #include "IO/H5/File.hpp" #include "IO/H5/TensorData.hpp" @@ -28,56 +30,89 @@ #include "Utilities/MakeArray.hpp" #include "Utilities/Serialization/Serialize.hpp" +namespace domain::creators::time_dependent_options { namespace { -template -std::string make_array_str(const double value) { - std::stringstream ss{}; - ss << "[" << value; - if constexpr (Dim > 1) { - ss << ", " << value; - if constexpr (Dim > 2) { - ss << ", " << value; - } - } - ss << "]"; - - return ss.str(); -} - +constexpr double infinity = std::numeric_limits::infinity(); void test_rotation_map_options() { { - const auto rotation_map_options = TestHelpers::test_creation< - domain::creators::time_dependent_options::RotationMapOptions<2>>( - "InitialQuaternions: [[1.0, 0.0, 0.0, 0.0]]\n" - "InitialAngles: Auto\n" - "DecayTimescale: Auto\n"); - CHECK(rotation_map_options.name() == "RotationMap"); - std::array expected_quat_func = - make_array<3>(DataVector{4, 0.0}); - expected_quat_func[0][0] = 1.0; - const std::array expected_angle_func = - make_array<3>(DataVector{3, 0.0}); - CHECK(expected_quat_func == rotation_map_options.quaternions); - CHECK(expected_angle_func == rotation_map_options.angles); - CHECK_FALSE(rotation_map_options.decay_timescale.has_value()); + INFO("None"); + const auto rotation_map_options = + TestHelpers::test_option_tag>("None"); + + CHECK(not rotation_map_options.has_value()); } + const auto non_settle_fots = []() { + INFO("Hardcoded AllowSettleFots = " + std::to_string(AllowSettleFoTs) + + ", Non-Settle Fots"); + const auto rotation_map_options = + TestHelpers::test_option_tag>( + "InitialAngularVelocity: [0.0, 0.4, 0.1]\n"); + + REQUIRE(rotation_map_options.has_value()); + CHECK(std::holds_alternative>( + rotation_map_options.value())); + + const auto& hardcoded_options = + std::get>( + rotation_map_options.value()); + + CHECK(hardcoded_options.quaternions == + std::array{DataVector{1.0, 0.0, 0.0, 0.0}, DataVector{4, 0.0}, + DataVector{4, 0.0}}); + CHECK(hardcoded_options.angles == + std::array{DataVector{3, 0.0}, DataVector{0.0, 0.4, 0.1}, + DataVector{3, 0.0}, DataVector{3, 0.0}}); + CHECK_FALSE(hardcoded_options.decay_timescale.has_value()); + + const auto rotation_ptr = + get_rotation(rotation_map_options.value(), 0.3, 2.9); + + const auto* rotation = + dynamic_cast*>( + rotation_ptr.get()); + + REQUIRE(rotation != nullptr); + + CHECK(rotation->time_bounds() == std::array{0.3, 2.9}); + CHECK(rotation->func(0.3)[0] == DataVector{1.0, 0.0, 0.0, 0.0}); + }; + + non_settle_fots.template operator()(); + non_settle_fots.template operator()(); + { - const auto rotation_map_options = TestHelpers::test_creation< - domain::creators::time_dependent_options::RotationMapOptions<3>>( - "InitialQuaternions: [[1.0, 0.0, 0.0, 0.0]]\n" - "InitialAngles: [[0.1, 0.2, 0.3], [1.1, 1.2, 1.3]]\n" - "DecayTimescale: 50.0\n"); - CHECK(rotation_map_options.name() == "RotationMap"); - std::array expected_quat_func = - make_array<4>(DataVector{4, 0.0}); - expected_quat_func[0][0] = 1.0; - const std::array expected_angle_func{ - DataVector{0.1, 0.2, 0.3}, DataVector{1.1, 1.2, 1.3}, - DataVector{3, 0.0}, DataVector{3, 0.0}}; - CHECK(expected_quat_func == rotation_map_options.quaternions); - CHECK(expected_angle_func == rotation_map_options.angles); - CHECK(rotation_map_options.decay_timescale.has_value()); - CHECK(rotation_map_options.decay_timescale.value() == 50.0); + INFO("Hardcoded AllowSettleFots = true, Settle Fots"); + const auto rotation_map_options = + TestHelpers::test_option_tag>( + "InitialQuaternions: [[1.0, 0.0, 0.0, 0.0]]\n" + "DecayTimescale: 40\n"); + + REQUIRE(rotation_map_options.has_value()); + CHECK(std::holds_alternative>( + rotation_map_options.value())); + + const auto& hardcoded_options = + std::get>(rotation_map_options.value()); + + const std::array expected_values{DataVector{1.0}, DataVector{2.0}, + DataVector{3.0}}; + CHECK(hardcoded_options.quaternions == + std::array{DataVector{1.0, 0.0, 0.0, 0.0}, DataVector{4, 0.0}, + DataVector{4, 0.0}}); + CHECK(hardcoded_options.angles == make_array<4>(DataVector{3, 0.0})); + CHECK(hardcoded_options.decay_timescale == std::optional{40.0}); + + const auto rotation_ptr = + get_rotation(rotation_map_options.value(), 0.3, 2.9); + + const auto* rotation = + dynamic_cast( + rotation_ptr.get()); + + REQUIRE(rotation != nullptr); + + CHECK(rotation->time_bounds() == std::array{0.3, infinity}); + CHECK(rotation->func(0.3)[0] == DataVector{1.0, 0.0, 0.0, 0.0}); } std::unordered_map h5_file{filename}; - 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, serialize(functions_of_time)); - } { - const auto rotation_map_options = TestHelpers::test_creation< - domain::creators::time_dependent_options::RotationMapOptions<2>>( - "InitialQuaternions:\n" - " H5Filename: " + - filename + "\n SubfileName: " + subfile_name + - "\n Time: 0.0\n" - "InitialAngles: Auto\n" - "DecayTimescale: Auto\n"); - CHECK(rotation_map_options.name() == "RotationMap"); - // 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.1, 0.1, 0.1}, - DataVector{-0.03, 0.15, 0.15, 0.15}}; - std::array expected_angle{ - DataVector{3, 0.1}, DataVector{3, 0.2}, DataVector{3, 0.3}}; - - CHECK_ITERABLE_APPROX(expected_quaternion, - rotation_map_options.quaternions); - CHECK(expected_angle == rotation_map_options.angles); - CHECK_FALSE(rotation_map_options.decay_timescale.has_value()); + INFO("FromVolumeFile non-Settle FoTs"); + + if (file_system::check_if_file_exists(filename)) { + file_system::rm(filename, true); + } + + TestHelpers::domain::creators::write_volume_data(filename, subfile_name, + functions_of_time); + + const auto rotation_map_options = + TestHelpers::test_option_tag>( + "H5Filename: GoatCheese.h5\n" + "SubfileName: VolumeData"); + + REQUIRE(rotation_map_options.has_value()); + CHECK(std::holds_alternative(rotation_map_options.value())); + + const auto rotation_ptr = + get_rotation(rotation_map_options.value(), 0.3, 100.0); + + const auto* rotation = + dynamic_cast*>( + rotation_ptr.get()); + + REQUIRE(rotation != nullptr); + + CHECK(rotation->time_bounds() == std::array{0.3, 100.0}); + CHECK(rotation->func_and_2_derivs(0.3) == + functions_of_time.at("Rotation")->func_and_2_derivs(0.3)); } { - const auto rotation_map_options = TestHelpers::test_creation< - domain::creators::time_dependent_options::RotationMapOptions<3>>( - "InitialQuaternions:\n" - " H5Filename: " + - filename + "\n SubfileName: " + subfile_name + - "\n Time: 0.0\n" - "InitialAngles: [[0.11, 0.22, 0.33]]\n" - "DecayTimescale: Auto\n"); - CHECK(rotation_map_options.name() == "RotationMap"); - // 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.1, 0.1, 0.1}, - DataVector{-0.03, 0.15, 0.15, 0.15}, DataVector{4, 0.0}}; - std::array expected_angle{ - DataVector{0.11, 0.22, 0.33}, DataVector{3, 0.0}, DataVector{3, 0.0}, - DataVector{3, 0.0}}; - - CHECK_ITERABLE_APPROX(expected_quaternion, - rotation_map_options.quaternions); - CHECK(expected_angle == rotation_map_options.angles); - CHECK_FALSE(rotation_map_options.decay_timescale.has_value()); + INFO("FromVolumeFile Settle FoTs"); + + if (file_system::check_if_file_exists(filename)) { + file_system::rm(filename, true); + } + + functions_of_time["Rotation"] = + std::make_unique( + std::array{DataVector{1.0, 0.0, 0.0, 0.0}, DataVector{4, 2.0}, + DataVector{4, 3.0}}, + 0.0, 100.0); + + TestHelpers::domain::creators::write_volume_data(filename, subfile_name, + functions_of_time); + + const auto rotation_map_options = + TestHelpers::test_option_tag>( + "H5Filename: GoatCheese.h5\n" + "SubfileName: VolumeData"); + + REQUIRE(rotation_map_options.has_value()); + CHECK(std::holds_alternative(rotation_map_options.value())); + + const auto rotation_ptr = + get_rotation(rotation_map_options.value(), 0.3, 100.0); + + const auto* rotation = + dynamic_cast( + rotation_ptr.get()); + + REQUIRE(rotation != nullptr); + + CHECK(rotation->time_bounds() == std::array{0.0, infinity}); + CHECK(rotation->func_and_2_derivs(0.3) == + functions_of_time.at("Rotation")->func_and_2_derivs(0.3)); } if (file_system::check_if_file_exists(filename)) { @@ -171,3 +208,4 @@ SPECTRE_TEST_CASE("Unit.Domain.Creators.TimeDependentOptions.RotationMap", domain::FunctionsOfTime::register_derived_with_charm(); test_rotation_map_options(); } +} // namespace domain::creators::time_dependent_options diff --git a/tests/Unit/Domain/Creators/TimeDependentOptions/Test_ShapeMap.cpp b/tests/Unit/Domain/Creators/TimeDependentOptions/Test_ShapeMap.cpp index 5c4ec4c44868..c7b4118332e1 100644 --- a/tests/Unit/Domain/Creators/TimeDependentOptions/Test_ShapeMap.cpp +++ b/tests/Unit/Domain/Creators/TimeDependentOptions/Test_ShapeMap.cpp @@ -1,11 +1,9 @@ // Distributed under the MIT License. // See LICENSE.txt for details. -#include "DataStructures/DataVector.hpp" #include "Framework/TestingFramework.hpp" #include -#include #include #include #include @@ -13,8 +11,11 @@ #include #include #include +#include #include +#include "DataStructures/DataVector.hpp" +#include "Domain/Creators/TimeDependentOptions/FromVolumeFile.hpp" #include "Domain/Creators/TimeDependentOptions/ShapeMap.hpp" #include "Domain/FunctionsOfTime/FunctionOfTime.hpp" #include "Domain/FunctionsOfTime/PiecewisePolynomial.hpp" @@ -23,6 +24,7 @@ #include "Framework/TestCreation.hpp" #include "Framework/TestHelpers.hpp" #include "Helpers/DataStructures/MakeWithRandomValues.hpp" +#include "Helpers/Domain/Creators/TimeDependent/TestHelpers.hpp" #include "IO/H5/AccessType.hpp" #include "IO/H5/Dat.hpp" #include "IO/H5/File.hpp" @@ -35,6 +37,7 @@ #include "Utilities/FileSystem.hpp" #include "Utilities/Gsl.hpp" +namespace domain::creators::time_dependent_options { namespace { void test_kerr_schild_boyer_lindquist() { const auto kerr_schild_boyer_lindquist = TestHelpers::test_creation< @@ -65,148 +68,70 @@ void test_ylms_from_file() { CHECK_FALSE(ylms_from_file.set_l1_coefs_to_zero); } -void test_shape_map_options() { - { - constexpr bool include_transition_ends_at_cube = false; - constexpr domain::ObjectLabel object_label = domain::ObjectLabel::A; - CAPTURE(include_transition_ends_at_cube); - CAPTURE(object_label); - const auto shape_map_options = TestHelpers::test_creation< - domain::creators::time_dependent_options::ShapeMapOptions< - include_transition_ends_at_cube, object_label>>( - "LMax: 8\n" - "InitialValues: Spherical\n" - "SizeInitialValues: [0.5, 1.0, 2.4]"); - CHECK(shape_map_options.name() == "ShapeMapA"); - CHECK(shape_map_options.l_max == 8); - CHECK_FALSE(shape_map_options.initial_values.has_value()); - CHECK(shape_map_options.initial_size_values.has_value()); - CHECK(shape_map_options.initial_size_values.value() == - std::array{0.5, 1.0, 2.4}); - CHECK_FALSE(shape_map_options.transition_ends_at_cube); - } - { - constexpr bool include_transition_ends_at_cube = true; - constexpr domain::ObjectLabel object_label = domain::ObjectLabel::B; - CAPTURE(include_transition_ends_at_cube); - CAPTURE(object_label); - const auto shape_map_options = TestHelpers::test_creation< - domain::creators::time_dependent_options::ShapeMapOptions< - include_transition_ends_at_cube, object_label>>( - "LMax: 8\n" - "InitialValues:\n" - " Mass: 1.7\n" - " Spin: [0.45, 0.12, 0.34]\n" - "SizeInitialValues: Auto\n" - "TransitionEndsAtCube: True"); - CHECK(shape_map_options.name() == "ShapeMapB"); - CHECK(shape_map_options.l_max == 8); - CHECK(shape_map_options.initial_values.has_value()); - CHECK(std::holds_alternative( - shape_map_options.initial_values.value())); - CHECK_FALSE(shape_map_options.initial_size_values.has_value()); - CHECK(shape_map_options.transition_ends_at_cube); - } - { - constexpr bool include_transition_ends_at_cube = false; - constexpr domain::ObjectLabel object_label = domain::ObjectLabel::None; - CAPTURE(include_transition_ends_at_cube); - CAPTURE(object_label); - const auto shape_map_options = TestHelpers::test_creation< - domain::creators::time_dependent_options::ShapeMapOptions< - include_transition_ends_at_cube, object_label>>( - "LMax: 8\n" - "InitialValues:\n" - " H5Filename: TotalEclipseOfTheHeart.h5\n" - " SubfileNames:\n" - " - Ylm_coefs\n" - " MatchTime: 1.7\n" - " MatchTimeEpsilon: Auto\n" - " SetL1CoefsToZero: True\n" - " CheckFrame: True\n" - "SizeInitialValues: Auto"); - CHECK(shape_map_options.name() == "ShapeMap"); - CHECK(shape_map_options.l_max == 8); - CHECK(shape_map_options.initial_values.has_value()); - CHECK(std::holds_alternative< - domain::creators::time_dependent_options::YlmsFromFile>( - shape_map_options.initial_values.value())); - CHECK_FALSE(shape_map_options.initial_size_values.has_value()); - CHECK_FALSE(shape_map_options.transition_ends_at_cube); - } - { - constexpr bool include_transition_ends_at_cube = false; - constexpr domain::ObjectLabel object_label = domain::ObjectLabel::None; - CAPTURE(include_transition_ends_at_cube); - CAPTURE(object_label); - const auto shape_map_options = TestHelpers::test_creation< - domain::creators::time_dependent_options::ShapeMapOptions< - include_transition_ends_at_cube, object_label>>( - "LMax: 8\n" - "InitialValues:\n" - " DatFilename: PartialEclipseOfTheHeart.dat\n" - " MatchTime: 1.7\n" - " MatchTimeEpsilon: Auto\n" - " SetL1CoefsToZero: True\n" - "SizeInitialValues: Auto"); - CHECK(shape_map_options.name() == "ShapeMap"); - CHECK(shape_map_options.l_max == 8); - CHECK(shape_map_options.initial_values.has_value()); - CHECK(std::holds_alternative< - domain::creators::time_dependent_options::YlmsFromSpEC>( - shape_map_options.initial_values.value())); - CHECK_FALSE(shape_map_options.initial_size_values.has_value()); - CHECK_FALSE(shape_map_options.transition_ends_at_cube); - } -} - template void test_funcs(const gsl::not_null generator) { const double inner_radius = 0.5; const size_t l_max = 8; - // We choose a Schwarzschild BH so all coefs are zero and it's easy to check { - const auto shape_map_options = TestHelpers::test_creation< - domain::creators::time_dependent_options::ShapeMapOptions< - false, domain::ObjectLabel::None>>( + INFO("None"); + const auto shape_map_options = TestHelpers::test_option_tag< + ShapeMapOptions>("None"); + + CHECK(not shape_map_options.has_value()); + } + { + INFO("Spherical"); + const auto shape_map_options = TestHelpers::test_option_tag< + ShapeMapOptions>( "LMax: 8\n" - "InitialValues:\n" - " Mass: 1.0\n" - " Spin: [0.0, 0.0, 0.0]\n" - "SizeInitialValues: [0.5, 1.0, 2.4]"); + "InitialValues: Spherical\n" + "SizeInitialValues: Auto\n"); - const auto [shape_funcs, size_funcs] = - domain::creators::time_dependent_options::initial_shape_and_size_funcs( - shape_map_options, inner_radius); + REQUIRE(shape_map_options.has_value()); + REQUIRE(std::holds_alternative< + ShapeMapOptions>( + shape_map_options.value())); - for (size_t i = 0; i < shape_funcs.size(); i++) { - CHECK(gsl::at(shape_funcs, i) == - DataVector{ylm::Spherepack::spectral_size(l_max, l_max), 0.0}); - } - CHECK(size_funcs == std::array{DataVector{0.5}, DataVector{1.0}, - DataVector{2.4}, DataVector{0.0}}); + const FunctionsOfTimeMap shape_and_size = get_shape_and_size( + shape_map_options.value(), 0.1, 0.8, 0.9, inner_radius); + + CHECK(shape_and_size.at("Shape")->time_bounds() == std::array{0.1, 0.8}); + CHECK(shape_and_size.at("Size")->time_bounds() == std::array{0.1, 0.9}); + + CHECK(shape_and_size.at("Shape")->func_and_2_derivs(0.1) == + make_array<3>( + DataVector{ylm::Spherepack::spectral_size(l_max, l_max), 0.0})); + CHECK(shape_and_size.at("Size")->func_and_2_derivs(0.1) == + make_array<3>(DataVector{0.0})); } + // We choose a Schwarzschild BH so all coefs are zero and it's easy to check { - const auto shape_map_options = TestHelpers::test_creation< + INFO("KerrSchild"); + const auto shape_map_options = TestHelpers::test_option_tag< domain::creators::time_dependent_options::ShapeMapOptions< false, domain::ObjectLabel::None>>( "LMax: 8\n" "InitialValues:\n" " Mass: 1.0\n" " Spin: [0.0, 0.0, 0.0]\n" - "SizeInitialValues: Auto"); + "SizeInitialValues: [0.5, 1.0, 2.4]"); - const auto [shape_funcs, size_funcs] = - domain::creators::time_dependent_options::initial_shape_and_size_funcs( - shape_map_options, inner_radius); + REQUIRE(shape_map_options.has_value()); + REQUIRE(std::holds_alternative< + ShapeMapOptions>( + shape_map_options.value())); - for (size_t i = 0; i < shape_funcs.size(); i++) { - CHECK(gsl::at(shape_funcs, i) == - DataVector{ylm::Spherepack::spectral_size(l_max, l_max), 0.0}); - } - CHECK(size_funcs == std::array{DataVector{0.0}, DataVector{0.0}, - DataVector{0.0}, DataVector{0.0}}); + const FunctionsOfTimeMap shape_and_size = get_shape_and_size( + shape_map_options.value(), 0.1, 0.8, 0.9, inner_radius); + + CHECK(shape_and_size.at("Shape")->time_bounds() == std::array{0.1, 0.8}); + CHECK(shape_and_size.at("Size")->time_bounds() == std::array{0.1, 0.9}); + + CHECK(shape_and_size.at("Shape")->func_and_2_derivs(0.1) == + make_array<3>( + DataVector{ylm::Spherepack::spectral_size(l_max, l_max), 0.0})); + CHECK(shape_and_size.at("Size")->func_and_2_derivs(0.1) == + std::array{DataVector{0.5}, DataVector{1.0}, DataVector{2.4}}); } { const std::string test_filename{"TotalEclipseOfTheHeart.h5"}; @@ -261,9 +186,10 @@ void test_funcs(const gsl::not_null generator) { std::unique_ptr> functions_of_time{}; { - const auto shape_map_options = TestHelpers::test_creation< + INFO("YlmsFromFile specify size"); + const auto shape_map_options = TestHelpers::test_option_tag< domain::creators::time_dependent_options::ShapeMapOptions< - false, domain::ObjectLabel::None>>( + true, domain::ObjectLabel::B>>( "LMax: 8\n" "InitialValues:\n" " H5Filename: TotalEclipseOfTheHeart.h5\n" @@ -274,11 +200,26 @@ void test_funcs(const gsl::not_null generator) { " MatchTimeEpsilon: Auto\n" " SetL1CoefsToZero: True\n" " CheckFrame: False\n" - "SizeInitialValues: [1.1, 2.2, 3.3]"); + "SizeInitialValues: [1.1, 2.2, 3.3]\n" + "TransitionEndsAtCube: True"); - const auto [shape_funcs, size_funcs] = - domain::creators::time_dependent_options:: - initial_shape_and_size_funcs(shape_map_options, inner_radius); + REQUIRE(shape_map_options.has_value()); + REQUIRE( + std::holds_alternative>( + shape_map_options.value())); + + const FunctionsOfTimeMap shape_and_size = get_shape_and_size( + shape_map_options.value(), 0.1, 0.8, 0.9, inner_radius); + + CHECK(shape_and_size.at(shape_name)->time_bounds() == + std::array{0.1, 0.8}); + CHECK(shape_and_size.at(size_name)->time_bounds() == + std::array{0.1, 0.9}); + + const auto shape_funcs = + shape_and_size.at(shape_name)->func_and_2_derivs(0.1); + const auto size_funcs = + shape_and_size.at(size_name)->func_and_2_derivs(0.1); ylm::SpherepackIterator iter{l_max, l_max}; ylm::SpherepackIterator file_iter{file_l_max, file_l_max}; @@ -301,63 +242,25 @@ void test_funcs(const gsl::not_null generator) { } } } - CHECK(size_funcs == std::array{DataVector{1.1}, DataVector{2.2}, - DataVector{3.3}, DataVector{0.0}}); + CHECK(size_funcs == + std::array{DataVector{1.1}, DataVector{2.2}, DataVector{3.3}}); functions_of_time[shape_name] = - std::make_unique>( - time, shape_funcs, 100.0); - functions_of_time[size_name] = - std::make_unique>( - time, size_funcs, 100.0); - { - h5::H5File test_file(test_filename, true); - auto& vol_file = test_file.insert(volume_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, serialize(functions_of_time)); - } - } - - { - const auto shape_map_options = TestHelpers::test_creation< - domain::creators::time_dependent_options::ShapeMapOptions< - false, domain::ObjectLabel::B>>( - "LMax: 8\n" - "InitialValues:\n" - " H5Filename: TotalEclipseOfTheHeart.h5\n" - " SubfileName: VolumeData\n" - " Time: 1.7\n" - "SizeInitialValues: Auto"); + shape_and_size.at(shape_name)->get_clone(); + functions_of_time[size_name] = shape_and_size.at(size_name)->get_clone(); - const auto [shape_funcs, size_funcs] = - domain::creators::time_dependent_options:: - initial_shape_and_size_funcs(shape_map_options, inner_radius); - - CHECK(shape_funcs == - functions_of_time.at(shape_name)->func_and_2_derivs(time)); - auto all_size_funcs = - functions_of_time.at(size_name)->func_and_all_derivs(time); - CHECK(all_size_funcs.size() == size_funcs.size()); - for (size_t i = 0; i < size_funcs.size(); i++) { - CHECK(all_size_funcs[i] == gsl::at(size_funcs, i)); - } + // For the FromVolumeFile test later on + TestHelpers::domain::creators::write_volume_data( + test_filename, volume_subfile_name, functions_of_time); } // Already checked that shape funcs are correct. Here just check that size // funcs were automatically set to correct values { - const auto shape_map_options = TestHelpers::test_creation< + INFO("YlmsFromFile auto size"); + const auto shape_map_options = TestHelpers::test_option_tag< domain::creators::time_dependent_options::ShapeMapOptions< - false, domain::ObjectLabel::None>>( + true, domain::ObjectLabel::B>>( "LMax: 8\n" "InitialValues:\n" " H5Filename: TotalEclipseOfTheHeart.h5\n" @@ -368,17 +271,19 @@ void test_funcs(const gsl::not_null generator) { " MatchTimeEpsilon: Auto\n" " SetL1CoefsToZero: True\n" " CheckFrame: False\n" - "SizeInitialValues: Auto"); + "SizeInitialValues: Auto\n" + "TransitionEndsAtCube: true"); - const auto [shape_funcs, size_funcs] = - domain::creators::time_dependent_options:: - initial_shape_and_size_funcs(shape_map_options, inner_radius); + REQUIRE(shape_map_options.has_value()); + REQUIRE( + std::holds_alternative>( + shape_map_options.value())); + + const FunctionsOfTimeMap shape_and_size = get_shape_and_size( + shape_map_options.value(), 0.1, 0.8, 0.9, inner_radius); - for (size_t i = 0; i < shape_funcs.size(); i++) { - CHECK(gsl::at(shape_funcs, i)[0] == 0.0); - } CHECK_ITERABLE_APPROX( - size_funcs, + shape_and_size.at(size_name)->func_and_2_derivs(0.1), (std::array{ DataVector{(inner_radius - ylm::Spherepack::average( strahlkorpers[0].coefficients())) * @@ -386,7 +291,35 @@ void test_funcs(const gsl::not_null generator) { DataVector{(inner_radius - ylm::Spherepack::average( strahlkorpers[1].coefficients())) * 2.0 * sqrt(M_PI)}, - DataVector{0.0}, DataVector{0.0}})); + DataVector{0.0}})); + } + + { + INFO("FromVolumeFile"); + const auto shape_map_options = TestHelpers::test_option_tag< + domain::creators::time_dependent_options::ShapeMapOptions< + false, domain::ObjectLabel::B>>( + "TransitionEndsAtCube: false\n" + "H5Filename: TotalEclipseOfTheHeart.h5\n" + "SubfileName: VolumeData\n"); + + REQUIRE(shape_map_options.has_value()); + REQUIRE(std::holds_alternative< + FromVolumeFileShapeSize>( + shape_map_options.value())); + + const FunctionsOfTimeMap shape_and_size = get_shape_and_size( + shape_map_options.value(), 0.2, 1.1, 1.2, inner_radius); + + CHECK(shape_and_size.at(shape_name)->time_bounds() == + std::array{0.2, 1.1}); + CHECK(shape_and_size.at(size_name)->time_bounds() == + std::array{0.2, 1.2}); + + CHECK(shape_and_size.at(shape_name)->func_and_2_derivs(0.2) == + functions_of_time.at(shape_name)->func_and_2_derivs(0.2)); + CHECK(shape_and_size.at(size_name)->func_and_2_derivs(0.2) == + functions_of_time.at(size_name)->func_and_2_derivs(0.2)); } if (file_system::check_if_file_exists(test_filename)) { @@ -394,6 +327,7 @@ void test_funcs(const gsl::not_null generator) { } } { + INFO("YlmsFromSpec"); const std::string test_filename{"PartialEclipseOfTheHeart.dat"}; if (file_system::check_if_file_exists(test_filename)) { file_system::rm(test_filename, true); @@ -454,7 +388,7 @@ void test_funcs(const gsl::not_null generator) { } { - const auto shape_map_options = TestHelpers::test_creation< + const auto shape_map_options = TestHelpers::test_option_tag< domain::creators::time_dependent_options::ShapeMapOptions< false, domain::ObjectLabel::None>>( "LMax: 8\n" @@ -465,9 +399,17 @@ void test_funcs(const gsl::not_null generator) { " SetL1CoefsToZero: True\n" "SizeInitialValues: Auto"); - const auto [shape_funcs, size_funcs] = - domain::creators::time_dependent_options:: - initial_shape_and_size_funcs(shape_map_options, inner_radius); + REQUIRE(shape_map_options.has_value()); + REQUIRE(std::holds_alternative< + ShapeMapOptions>( + shape_map_options.value())); + + const FunctionsOfTimeMap shape_and_size = get_shape_and_size( + shape_map_options.value(), 0.2, 1.1, 1.2, inner_radius); + + const auto shape_funcs = + shape_and_size.at("Shape")->func_and_2_derivs(0.2); + const auto size_funcs = shape_and_size.at("Size")->func_and_2_derivs(0.2); ylm::SpherepackIterator iter{l_max, l_max}; ylm::SpherepackIterator file_iter{file_l_max, file_l_max}; @@ -498,11 +440,11 @@ void test_funcs(const gsl::not_null generator) { CHECK(size_funcs == std::array{DataVector{-1.0 * strahlkorper.coefficients()[0] * sqrt(0.5 * M_PI)}, - DataVector{0.0}, DataVector{0.0}, DataVector{0.0}}); + DataVector{0.0}, DataVector{0.0}}); } { - const auto shape_map_options = TestHelpers::test_creation< + const auto shape_map_options = TestHelpers::test_option_tag< domain::creators::time_dependent_options::ShapeMapOptions< false, domain::ObjectLabel::None>>( "LMax: 8\n" @@ -513,8 +455,8 @@ void test_funcs(const gsl::not_null generator) { " SetL1CoefsToZero: True\n" "SizeInitialValues: Auto"); CHECK_THROWS_WITH( - (domain::creators::time_dependent_options:: - initial_shape_and_size_funcs(shape_map_options, inner_radius)), + (get_shape_and_size(shape_map_options.value(), 0.2, 1.1, 1.2, + inner_radius)), Catch::Matchers::ContainsSubstring( "Unable to find requested time 1") and Catch::Matchers::ContainsSubstring( @@ -534,6 +476,6 @@ SPECTRE_TEST_CASE("Unit.Domain.Creators.TimeDependentOptions.ShapeMap", MAKE_GENERATOR(generator); test_kerr_schild_boyer_lindquist(); test_ylms_from_file(); - test_shape_map_options(); test_funcs(make_not_null(&generator)); } +} // namespace domain::creators::time_dependent_options diff --git a/tests/Unit/Domain/Creators/TimeDependentOptions/Test_TranslationMap.cpp b/tests/Unit/Domain/Creators/TimeDependentOptions/Test_TranslationMap.cpp index 090a8ff554da..92a4ee3e63ad 100644 --- a/tests/Unit/Domain/Creators/TimeDependentOptions/Test_TranslationMap.cpp +++ b/tests/Unit/Domain/Creators/TimeDependentOptions/Test_TranslationMap.cpp @@ -16,6 +16,7 @@ #include "Domain/FunctionsOfTime/PiecewisePolynomial.hpp" #include "Domain/FunctionsOfTime/RegisterDerivedWithCharm.hpp" #include "Framework/TestCreation.hpp" +#include "Helpers/Domain/Creators/TimeDependent/TestHelpers.hpp" #include "IO/H5/AccessType.hpp" #include "IO/H5/File.hpp" #include "IO/H5/TensorData.hpp" @@ -27,6 +28,7 @@ #include "Utilities/MakeArray.hpp" #include "Utilities/Serialization/Serialize.hpp" +namespace domain::creators::time_dependent_options { namespace { template std::string make_array_str(const double value) { @@ -46,16 +48,41 @@ std::string make_array_str(const double value) { template void test_translation_map_options() { { - const auto translation_map_options = TestHelpers::test_creation< + INFO("None"); + const auto translation_map_options = + TestHelpers::test_option_tag>("None"); + + CHECK(not translation_map_options.has_value()); + } + { + INFO("Hardcoded options"); + const auto translation_map_options = TestHelpers::test_option_tag< domain::creators::time_dependent_options::TranslationMapOptions>( "InitialValues: [" + make_array_str(1.0) + "," + make_array_str(2.0) + "," + make_array_str(3.0) + "]"); - CHECK(translation_map_options.name() == "TranslationMap"); - CHECK(translation_map_options.initial_values == - std::array{DataVector{Dim, 1.0}, DataVector{Dim, 2.0}, - DataVector{Dim, 3.0}}); + + REQUIRE(translation_map_options.has_value()); + CHECK(std::holds_alternative>( + translation_map_options.value())); + const std::array initial_values{DataVector{Dim, 1.0}, DataVector{Dim, 2.0}, + DataVector{Dim, 3.0}}; + CHECK(std::get>(translation_map_options.value()) + .initial_values == initial_values); + + const auto translation_ptr = + get_translation(translation_map_options.value(), 0.1, 65.8); + + const auto* translation = + dynamic_cast*>( + translation_ptr.get()); + + CHECK(translation != nullptr); + + CHECK(translation->time_bounds() == std::array{0.1, 65.8}); + CHECK(translation->func_and_2_derivs(0.1) == initial_values); } { + INFO("FromVolumeFile"); std::unordered_map> functions_of_time{}; @@ -71,31 +98,32 @@ void test_translation_map_options() { file_system::rm(filename, true); } - { - h5::H5File h5_file{filename}; - 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, serialize(functions_of_time)); - } + TestHelpers::domain::creators::write_volume_data(filename, subfile_name, + functions_of_time); - const auto translation_map_options = TestHelpers::test_creation< + const auto translation_map_options = TestHelpers::test_option_tag< domain::creators::time_dependent_options::TranslationMapOptions>( - "InitialValues:\n" - " H5Filename: " + - filename + "\n SubfileName: " + subfile_name + "\n Time: 0.0"); - CHECK(translation_map_options.name() == "TranslationMap"); - CHECK(translation_map_options.initial_values == - std::array{DataVector{Dim, 1.0}, DataVector{Dim, 2.0}, - DataVector{Dim, 3.0}}); + "H5Filename: " + filename + "\nSubfileName: " + subfile_name); + + REQUIRE(translation_map_options.has_value()); + CHECK(std::holds_alternative( + translation_map_options.value())); + const std::array initial_values{DataVector{Dim, 1.0}, DataVector{Dim, 2.0}, + DataVector{Dim, 3.0}}; + + const auto translation_ptr = + get_translation(translation_map_options.value(), 0.1, 65.8); + + const auto* translation = + dynamic_cast*>( + translation_ptr.get()); + + CHECK(translation != nullptr); + + CHECK(translation->time_bounds() == std::array{0.1, 65.8}); + CHECK_ITERABLE_APPROX( + translation->func_and_2_derivs(0.3), + functions_of_time.at("Translation")->func_and_2_derivs(0.3)); if (file_system::check_if_file_exists(filename)) { file_system::rm(filename, true); @@ -111,3 +139,4 @@ SPECTRE_TEST_CASE("Unit.Domain.Creators.TimeDependentOptions.TranslationMap", test_translation_map_options<2>(); test_translation_map_options<3>(); } +} // namespace domain::creators::time_dependent_options diff --git a/tests/Unit/Domain/FunctionsOfTime/Test_FixedSpeedCubic.cpp b/tests/Unit/Domain/FunctionsOfTime/Test_FixedSpeedCubic.cpp index 6b1c0bb996c1..7907978667ee 100644 --- a/tests/Unit/Domain/FunctionsOfTime/Test_FixedSpeedCubic.cpp +++ b/tests/Unit/Domain/FunctionsOfTime/Test_FixedSpeedCubic.cpp @@ -75,6 +75,11 @@ void test( CHECK(get_output(*fixed_speed_cubic) == "FixedSpeedCubic(t=10, f=1, v=0.4, tau^2=25)"); } + + CHECK(dynamic_cast( + *f_of_t->get_clone()) == + dynamic_cast( + *f_of_t->create_at_time(0.0, 0.0))); } void test_function(const double initial_function_value, diff --git a/tests/Unit/Domain/FunctionsOfTime/Test_IntegratedFunctionOfTime.cpp b/tests/Unit/Domain/FunctionsOfTime/Test_IntegratedFunctionOfTime.cpp index 5afff2c65230..65c38630bc6b 100644 --- a/tests/Unit/Domain/FunctionsOfTime/Test_IntegratedFunctionOfTime.cpp +++ b/tests/Unit/Domain/FunctionsOfTime/Test_IntegratedFunctionOfTime.cpp @@ -59,6 +59,20 @@ void test(const gsl::not_null f_of_t, CHECK(f_of_t->expiration_after(t) == approx(t + 0.5 * dt)); CHECK(*f_of_t_derived != f_of_t_derived_copy); } + + const size_t index = positions.size() - 8; + const auto new_time = static_cast(index) * dt; + const auto new_expiration = new_time + dt; + const auto copy_at_time = f_of_t->create_at_time(new_time, new_expiration); + + CAPTURE(rotation); + const double check_time = new_time; + CHECK_ITERABLE_APPROX(copy_at_time->func(check_time), + f_of_t_derived->func(check_time)); + CHECK_ITERABLE_APPROX(copy_at_time->func_and_deriv(check_time), + f_of_t_derived->func_and_deriv(check_time)); + const auto t_bounds = copy_at_time->time_bounds(); + CHECK(t_bounds == std::array{new_time, new_expiration}); } void test_out_of_order_update() { diff --git a/tests/Unit/Domain/FunctionsOfTime/Test_PiecewisePolynomial.cpp b/tests/Unit/Domain/FunctionsOfTime/Test_PiecewisePolynomial.cpp index d32315b0e979..bdaed9a5480d 100644 --- a/tests/Unit/Domain/FunctionsOfTime/Test_PiecewisePolynomial.cpp +++ b/tests/Unit/Domain/FunctionsOfTime/Test_PiecewisePolynomial.cpp @@ -72,6 +72,13 @@ void test(const gsl::not_null f_of_t, CHECK(t_bounds[0] == 0.0); CHECK(t_bounds[1] == 4.8); CHECK(f_of_t->expiration_after(0.0) == dt); + + const double new_time = 1.0; + const double new_expiration = 2.1; + const auto copy_at_time = f_of_t->create_at_time(new_time, new_expiration); + const auto new_t_bounds = copy_at_time->time_bounds(); + CHECK(new_t_bounds[0] == new_time); + CHECK(new_t_bounds[1] == new_expiration); } template @@ -173,21 +180,33 @@ void check_func_and_derivs( const FunctionsOfTime::PiecewisePolynomial& f_of_t, const quartic& q) { const double t = 0.5; - 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)); - } + + const auto& do_check = [&](const auto& f_of_t_to_check) { + CHECK_ITERABLE_APPROX(f_of_t_to_check.func(t), q.func_and_derivs<0>(t)); + CHECK_ITERABLE_APPROX(f_of_t_to_check.func_and_deriv(t), + q.func_and_derivs<1>(t)); + CHECK_ITERABLE_APPROX(f_of_t_to_check.func_and_2_derivs(t), + q.func_and_derivs<2>(t)); + auto func_and_all_derivs = f_of_t_to_check.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_ITERABLE_APPROX(func_and_all_derivs[i], + gsl::at(func_and_2_derivs, i)); + } + }; + + do_check(f_of_t); test_copy_semantics(f_of_t); auto f_of_t_copy = f_of_t; test_move_semantics(std::move(f_of_t_copy), f_of_t); + + const auto copy_at_time = f_of_t.create_at_time(0.4, 1.1); + + do_check(*copy_at_time); } void test_func_and_derivs() { diff --git a/tests/Unit/Domain/FunctionsOfTime/Test_QuaternionFunctionOfTime.cpp b/tests/Unit/Domain/FunctionsOfTime/Test_QuaternionFunctionOfTime.cpp index 9bbb9c5877fa..bd23b747b2b9 100644 --- a/tests/Unit/Domain/FunctionsOfTime/Test_QuaternionFunctionOfTime.cpp +++ b/tests/Unit/Domain/FunctionsOfTime/Test_QuaternionFunctionOfTime.cpp @@ -189,18 +189,25 @@ SPECTRE_TEST_CASE("Unit.Domain.FunctionsOfTime.QuaternionFunctionOfTime", DataVector{{cos(1.0), 0.0, 0.0, sin(1.0)}}, DataVector{{-0.5 * sin(1.0), 0.0, 0.0, 0.5 * cos(1.0)}}}; - CHECK_ITERABLE_APPROX(qfot.quat_func(2.0), expected_func); - CHECK_ITERABLE_APPROX(qfot.quat_func_and_deriv(2.0), - expected_func_and_deriv); - CHECK_ITERABLE_APPROX(qfot.func(2.0), expected_func); - CHECK_ITERABLE_APPROX(qfot.func_and_deriv(2.0), expected_func_and_deriv); - CHECK_ITERABLE_APPROX(qfot_ptr->func(2.0), qfot.func(2.0)); - CHECK_ITERABLE_APPROX(qfot_ptr->func_and_deriv(2.0), - qfot.func_and_deriv(2.0)); - CHECK_ITERABLE_APPROX(qfot_serialized_deserialized.func(2.0), - qfot.func(2.0)); - CHECK_ITERABLE_APPROX(qfot_serialized_deserialized.func_and_deriv(2.0), - qfot.func_and_deriv(2.0)); + const auto do_check = [&](const auto& qfot_to_check, + const auto& ptr_to_check, + const auto& serialized_to_check) { + CHECK_ITERABLE_APPROX(qfot_to_check.quat_func(2.0), expected_func); + CHECK_ITERABLE_APPROX(qfot_to_check.quat_func_and_deriv(2.0), + expected_func_and_deriv); + CHECK_ITERABLE_APPROX(qfot_to_check.func(2.0), expected_func); + CHECK_ITERABLE_APPROX(qfot_to_check.func_and_deriv(2.0), + expected_func_and_deriv); + CHECK_ITERABLE_APPROX(ptr_to_check->func(2.0), qfot_to_check.func(2.0)); + CHECK_ITERABLE_APPROX(ptr_to_check->func_and_deriv(2.0), + qfot_to_check.func_and_deriv(2.0)); + CHECK_ITERABLE_APPROX(serialized_to_check.func(2.0), + qfot_to_check.func(2.0)); + CHECK_ITERABLE_APPROX(serialized_to_check.func_and_deriv(2.0), + qfot_to_check.func_and_deriv(2.0)); + }; + + do_check(qfot, qfot_ptr, qfot_serialized_deserialized); CHECK_THROWS_WITH( qfot.quat_func_and_3_derivs(2.0), @@ -210,6 +217,18 @@ SPECTRE_TEST_CASE("Unit.Domain.FunctionsOfTime.QuaternionFunctionOfTime", test_copy_semantics(qfot); test_move_semantics(std::move(qfot_serialized_deserialized), qfot); + + const auto copy_at_time = qfot.create_at_time(0.7, 2.8); + const auto clone_copy_at_time = copy_at_time->get_clone(); + const auto copy_serialized = serialize_and_deserialize( + dynamic_cast< + const domain::FunctionsOfTime::QuaternionFunctionOfTime<2>&>( + *copy_at_time)); + + do_check(dynamic_cast< + const domain::FunctionsOfTime::QuaternionFunctionOfTime<2>&>( + *copy_at_time), + clone_copy_at_time, copy_serialized); } { diff --git a/tests/Unit/Domain/FunctionsOfTime/Test_SettleToConstant.cpp b/tests/Unit/Domain/FunctionsOfTime/Test_SettleToConstant.cpp index 0333a51700cb..8af0e73ee2a1 100644 --- a/tests/Unit/Domain/FunctionsOfTime/Test_SettleToConstant.cpp +++ b/tests/Unit/Domain/FunctionsOfTime/Test_SettleToConstant.cpp @@ -50,6 +50,11 @@ void test( CHECK(t_bounds[1] == std::numeric_limits::infinity()); CHECK(f_of_t->expiration_after(match_time) == std::numeric_limits::infinity()); + + CHECK(dynamic_cast( + *f_of_t->get_clone()) == + dynamic_cast( + *f_of_t->create_at_time(0.0, 0.0))); } } // namespace diff --git a/tests/Unit/Domain/FunctionsOfTime/Test_SettleToConstantQuaternion.cpp b/tests/Unit/Domain/FunctionsOfTime/Test_SettleToConstantQuaternion.cpp index 2fdafcc552c3..babac9156ea7 100644 --- a/tests/Unit/Domain/FunctionsOfTime/Test_SettleToConstantQuaternion.cpp +++ b/tests/Unit/Domain/FunctionsOfTime/Test_SettleToConstantQuaternion.cpp @@ -57,6 +57,12 @@ void test( CHECK(t_bounds[1] == std::numeric_limits::infinity()); CHECK(f_of_t->expiration_after(match_time) == std::numeric_limits::infinity()); + + CHECK( + dynamic_cast( + *f_of_t->get_clone()) == + dynamic_cast( + *f_of_t->create_at_time(0.0, 0.0))); } } // namespace diff --git a/tests/Unit/Evolution/Ringdown/Test_StrahlkorperCoefsInRingdownDistortedFrame.cpp b/tests/Unit/Evolution/Ringdown/Test_StrahlkorperCoefsInRingdownDistortedFrame.cpp index 23d7a5fa67ef..7404799d9a1c 100644 --- a/tests/Unit/Evolution/Ringdown/Test_StrahlkorperCoefsInRingdownDistortedFrame.cpp +++ b/tests/Unit/Evolution/Ringdown/Test_StrahlkorperCoefsInRingdownDistortedFrame.cpp @@ -15,8 +15,12 @@ #include "DataStructures/DataVector.hpp" #include "Domain/CoordinateMaps/Distribution.hpp" #include "Domain/Creators/Sphere.hpp" +#include "Domain/Creators/TimeDependentOptions/ExpansionMap.hpp" +#include "Domain/Creators/TimeDependentOptions/RotationMap.hpp" +#include "Domain/Creators/TimeDependentOptions/ShapeMap.hpp" #include "Domain/Creators/TimeDependentOptions/Sphere.hpp" #include "Domain/StrahlkorperTransformations.hpp" +#include "Domain/Structure/ObjectLabel.hpp" #include "Evolution/Ringdown/StrahlkorperCoefsInRingdownDistortedFrame.hpp" #include "Framework/TestHelpers.hpp" #include "Helpers/DataStructures/MakeWithRandomValues.hpp" @@ -86,7 +90,7 @@ SPECTRE_TEST_CASE( for (size_t i = 0; i < 4; ++i) { gsl::at(initial_unit_quaternion, i) /= initial_unit_quaternion_magnitude; } - const std::array, 3> rot_func_and_2_derivs{ + const std::vector> rot_func_and_2_derivs{ initial_unit_quaternion, make_with_random_values>(make_not_null(&generator), make_not_null(&fot_dist)), @@ -96,13 +100,14 @@ SPECTRE_TEST_CASE( std::uniform_real_distribution settling_dist{0.5, 1.5}; const double settling_timescale{settling_dist(generator)}; - const domain::creators::sphere::TimeDependentMapOptions::ShapeMapOptions + const domain::creators::time_dependent_options::ShapeMapOptions< + false, domain::ObjectLabel::None> shape_map_options{l_max, std::nullopt}; - const domain::creators::sphere::TimeDependentMapOptions::ExpansionMapOptions + const domain::creators::time_dependent_options::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 + const domain::creators::time_dependent_options::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, diff --git a/tests/Unit/Evolution/Systems/Cce/Test_OptionTags.cpp b/tests/Unit/Evolution/Systems/Cce/Test_OptionTags.cpp index d21998532637..e6a6c3a74dd3 100644 --- a/tests/Unit/Evolution/Systems/Cce/Test_OptionTags.cpp +++ b/tests/Unit/Evolution/Systems/Cce/Test_OptionTags.cpp @@ -108,13 +108,13 @@ SPECTRE_TEST_CASE("Unit.Evolution.Systems.Cce.OptionTags", "[Unit][Cce]") { "100.0") == 100.0); CHECK(TestHelpers::test_option_tag("4.0") == - Options::Auto{4.0}); + std::optional{4.0}); CHECK(TestHelpers::test_option_tag("Auto") == - Options::Auto{}); + std::optional{}); CHECK(TestHelpers::test_option_tag("2.0") == - Options::Auto{2.0}); + std::optional{2.0}); CHECK(TestHelpers::test_option_tag("Auto") == - Options::Auto{}); + std::optional{}); CHECK(TestHelpers::test_option_tag( "OptionTagsCceR0100.h5") == "OptionTagsCceR0100.h5"); CHECK(TestHelpers::test_option_tag< diff --git a/tests/Unit/Framework/TestCreation.hpp b/tests/Unit/Framework/TestCreation.hpp index e5226680b5d1..97f130310953 100644 --- a/tests/Unit/Framework/TestCreation.hpp +++ b/tests/Unit/Framework/TestCreation.hpp @@ -8,6 +8,7 @@ #include #include +#include "Options/Auto.hpp" #include "Options/Options.hpp" #include "Options/ParseOptions.hpp" #include "Options/Protocols/FactoryCreation.hpp" @@ -52,6 +53,23 @@ struct AddGroups> { } }; +template +struct option_return; + +template +struct option_return { + using type = typename OptionTag::type::value_type; +}; + +template +struct option_return { + using type = typename OptionTag::type; +}; + +template +using option_return_t = typename option_return< + OptionTag, tt::is_a_v>::type; + template struct SingleFactoryMetavariables { struct factory_creation @@ -106,7 +124,7 @@ T test_creation(const std::string& construction_string) { /// /// \see TestHelpers::test_creation() template -typename OptionTag::type test_option_tag( +TestCreation_detail::option_return_t test_option_tag( const std::string& construction_string) { Options::Parser> options(""); options.parse( diff --git a/tests/Unit/Helpers/Domain/CMakeLists.txt b/tests/Unit/Helpers/Domain/CMakeLists.txt index b3ff4bcef766..d05614b12589 100644 --- a/tests/Unit/Helpers/Domain/CMakeLists.txt +++ b/tests/Unit/Helpers/Domain/CMakeLists.txt @@ -5,6 +5,7 @@ set(LIBRARY "DomainHelpers") set(LIBRARY_SOURCES DomainTestHelpers.cpp + Creators/TimeDependent/TestHelpers.cpp ) add_spectre_library(${LIBRARY} ${SPECTRE_TEST_LIBS_TYPE} ${LIBRARY_SOURCES}) diff --git a/tests/Unit/Helpers/Domain/Creators/TimeDependent/TestHelpers.cpp b/tests/Unit/Helpers/Domain/Creators/TimeDependent/TestHelpers.cpp new file mode 100644 index 000000000000..d5997c77c06a --- /dev/null +++ b/tests/Unit/Helpers/Domain/Creators/TimeDependent/TestHelpers.cpp @@ -0,0 +1,45 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Helpers/Domain/Creators/TimeDependent/TestHelpers.hpp" + +#include +#include +#include +#include +#include +#include + +#include "DataStructures/DataVector.hpp" +#include "Domain/FunctionsOfTime/FunctionOfTime.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/Serialization/Serialize.hpp" + +namespace TestHelpers::domain::creators { +void write_volume_data( + const std::string& filename, const std::string& subfile_name, + const std::unordered_map< + std::string, + std::unique_ptr<::domain::FunctionsOfTime::FunctionOfTime>>& + 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)}); +} +} // namespace TestHelpers::domain::creators diff --git a/tests/Unit/Helpers/Domain/Creators/TimeDependent/TestHelpers.hpp b/tests/Unit/Helpers/Domain/Creators/TimeDependent/TestHelpers.hpp new file mode 100644 index 000000000000..aa7252ad176e --- /dev/null +++ b/tests/Unit/Helpers/Domain/Creators/TimeDependent/TestHelpers.hpp @@ -0,0 +1,19 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include +#include +#include + +#include "Domain/FunctionsOfTime/FunctionOfTime.hpp" + +namespace TestHelpers::domain::creators { +void write_volume_data( + const std::string& filename, const std::string& subfile_name, + const std::unordered_map< + std::string, + std::unique_ptr<::domain::FunctionsOfTime::FunctionOfTime>>& + functions_of_time = {}); +} // namespace TestHelpers::domain::creators diff --git a/tests/Unit/Helpers/Evolution/Systems/CurvedScalarWave/Worldtube/TestHelpers.cpp b/tests/Unit/Helpers/Evolution/Systems/CurvedScalarWave/Worldtube/TestHelpers.cpp index d496ab2fa529..2942a5f5641b 100644 --- a/tests/Unit/Helpers/Evolution/Systems/CurvedScalarWave/Worldtube/TestHelpers.cpp +++ b/tests/Unit/Helpers/Evolution/Systems/CurvedScalarWave/Worldtube/TestHelpers.cpp @@ -65,11 +65,11 @@ std::unique_ptr> worldtube_binary_compact_object( " TimeDependentMaps:\n" " InitialTime: 0.0\n" " ExpansionMap: \n" - " InitialValues: [1.0, 0.0]\n" + " InitialValues: [1.0, 0.0, 0.0]\n" " AsymptoticVelocityOuterBoundary: 0.0\n" - " DecayTimescaleOuterBoundaryVelocity: 1.0\n" + " DecayTimescaleOuterBoundary: 1.0\n" " RotationMap:\n" - " InitialAngularVelocity: [0.0, 0.0," + + " InitialAngularVelocity: [0.0, 0.0, " + angular_velocity_stream.str() + "]\n" " TranslationMap: None\n" diff --git a/tests/Unit/ParallelAlgorithms/Interpolation/Test_ComputeDestVars.cpp b/tests/Unit/ParallelAlgorithms/Interpolation/Test_ComputeDestVars.cpp index 89aca0198a26..350a530fd8b2 100644 --- a/tests/Unit/ParallelAlgorithms/Interpolation/Test_ComputeDestVars.cpp +++ b/tests/Unit/ParallelAlgorithms/Interpolation/Test_ComputeDestVars.cpp @@ -11,10 +11,12 @@ #include "DataStructures/Variables.hpp" #include "Domain/Creators/Sphere.hpp" #include "Domain/Creators/Tags/FunctionsOfTime.hpp" +#include "Domain/Creators/TimeDependentOptions/ShapeMap.hpp" #include "Domain/Creators/TimeDependentOptions/Sphere.hpp" #include "Domain/FunctionsOfTime/FunctionOfTime.hpp" #include "Domain/FunctionsOfTime/RegisterDerivedWithCharm.hpp" #include "Domain/Structure/ElementId.hpp" +#include "Domain/Structure/ObjectLabel.hpp" #include "Framework/ActionTesting.hpp" #include "NumericalAlgorithms/Spectral/Mesh.hpp" #include "Parallel/GlobalCache.hpp" @@ -108,9 +110,13 @@ void test() { domain::FunctionsOfTime::register_derived_with_charm(); using TDMO = domain::creators::sphere::TimeDependentMapOptions; - TDMO time_dep_opts{0.0, TDMO::ShapeMapOptions{2, std::nullopt}, - std::nullopt, std::nullopt, - std::nullopt, true}; + TDMO time_dep_opts{0.0, + domain::creators::time_dependent_options::ShapeMapOptions< + false, domain::ObjectLabel::None>{2, std::nullopt}, + std::nullopt, + std::nullopt, + std::nullopt, + true}; const auto domain_creator = domain::creators::Sphere( 0.9, 4.9, domain::creators::Sphere::Excision{}, 1_st, 7_st, false, {},