diff --git a/src/Domain/CoordinateMaps/TimeDependent/Translation.cpp b/src/Domain/CoordinateMaps/TimeDependent/Translation.cpp index de3f73ed2d86..cf64bb003f87 100644 --- a/src/Domain/CoordinateMaps/TimeDependent/Translation.cpp +++ b/src/Domain/CoordinateMaps/TimeDependent/Translation.cpp @@ -359,7 +359,8 @@ std::array, Dim> Translation::piecewise_helper( for (size_t i = 0; i < Dim; i++) { get_element(gsl::at(result, i), k) += gsl::at(func_or_deriv_of_time, i); } - ASSERT(max(magnitude(result)) <= outer_radius_.value(), + const double eps = std::numeric_limits::epsilon() * 100.0; + ASSERT(max(magnitude(result)) - eps <= outer_radius_.value(), "Coordinates translated outside outer radius, This should not " "happen"); } else if (get_element(radius, k) > inner_radius_.value() and @@ -373,7 +374,8 @@ std::array, Dim> Translation::piecewise_helper( get_element(gsl::at(result, i), k) += gsl::at(func_or_deriv_of_time, i) * radial_falloff_factor; } - ASSERT(max(magnitude(result)) <= outer_radius_.value(), + const double eps = std::numeric_limits::epsilon() * 100.0; + ASSERT(max(magnitude(result)) - eps <= outer_radius_.value(), "Coordinates translated outside outer radius, This should not " "happen"); } diff --git a/src/Domain/Creators/Sphere.cpp b/src/Domain/Creators/Sphere.cpp index b244b51f240c..5ded4bdcd91f 100644 --- a/src/Domain/Creators/Sphere.cpp +++ b/src/Domain/Creators/Sphere.cpp @@ -267,11 +267,23 @@ Sphere::Sphere( // Build the maps. We only apply the maps in the inner-most shell. The // inner radius is what's passed in, but the outer radius is the outer // radius of the inner-most shell so we have to see how many shells we - // have + // have. + // The translation map linear transition occurs only in the outer-most + // shell, such that the outer boundary is not translated. Require at least + // one radial partition, as the transition requires spherical shape. + + if (radial_partitioning_.empty()) { + PARSE_ERROR(context, + "The hard-coded translation map requires at least two " + "shells. Use at least one radial partition."); + } + std::get(time_dependent_options_.value()) .build_maps(std::array{0.0, 0.0, 0.0}, inner_radius_, - radial_partitioning_.empty() ? outer_radius_ - : radial_partitioning_[0]); + radial_partitioning_[0], + // Translation inner radius + std::pair{radial_partitioning_.back(), + outer_radius_}); } } } @@ -374,9 +386,14 @@ Domain<3> Sphere::create_domain() const { time_dependent_options_.value()); // First shell gets the distorted maps. + // Last shell gets translation with transition region for (size_t block_id = 0; block_id < num_blocks_; block_id++) { const bool include_distorted_map_in_first_shell = block_id < num_blocks_per_shell_; + // False if block_id is in the last shell + const bool use_rigid_translation = + block_id + num_blocks_per_shell_ + (fill_interior_ ? 1 : 0) < + num_blocks_; block_maps_grid_to_distorted[block_id] = hard_coded_options.grid_to_distorted_map( include_distorted_map_in_first_shell); @@ -385,7 +402,7 @@ Domain<3> Sphere::create_domain() const { include_distorted_map_in_first_shell); block_maps_grid_to_inertial[block_id] = hard_coded_options.grid_to_inertial_map( - include_distorted_map_in_first_shell); + include_distorted_map_in_first_shell, use_rigid_translation); } } else { const auto& time_dependence = std::get> @@ -46,7 +49,8 @@ TimeDependentMapOptions::create_functions_of_time( // their initial expiration time to infinity (i.e. not expiring) std::unordered_map expiration_times{ {size_name, std::numeric_limits::infinity()}, - {shape_name, std::numeric_limits::infinity()}}; + {shape_name, std::numeric_limits::infinity()}, + {translation_name, std::numeric_limits::infinity()}}; // If we have control systems, overwrite these expiration times with the ones // supplied by the control system @@ -101,12 +105,32 @@ TimeDependentMapOptions::create_functions_of_time( {0.0}}}, expiration_times.at(size_name)); + DataVector initial_translation_center_temp{3, 0.0}; + DataVector initial_translation_velocity_temp{3, 0.0}; + for (size_t i = 0; i < 3; i++) { + initial_translation_center_temp[i] = + gsl::at(initial_translation_values_.front(), i); + initial_translation_velocity_temp[i] = + gsl::at(initial_translation_values_.back(), i); + } + + // TranslationMap FunctionOfTime + result[translation_name] = + std::make_unique>( + initial_time_, + std::array{ + {std::move(initial_translation_center_temp), + std::move(initial_translation_velocity_temp), + {3, 0.0}}}, + expiration_times.at(translation_name)); + return result; } -void TimeDependentMapOptions::build_maps(const std::array& center, - const double inner_radius, - const double outer_radius) { +void TimeDependentMapOptions::build_maps( + const std::array& center, const double inner_radius, + const double outer_radius, + std::pair translation_transition_radii) { std::unique_ptr transition_func = @@ -115,17 +139,22 @@ void TimeDependentMapOptions::build_maps(const std::array& center, shape_map_ = ShapeMap{center, initial_l_max_, initial_l_max_, std::move(transition_func), shape_name, size_name}; + + rigid_translation_map_ = TranslationMap{translation_name}; + + linear_falloff_translation_map_ = + TranslationMap{translation_name, translation_transition_radii.first, + translation_transition_radii.second}; } // If you edit any of the functions below, be sure to update the documentation // in the Sphere domain creator as well as this class' documentation. TimeDependentMapOptions::MapType TimeDependentMapOptions::distorted_to_inertial_map( - const bool include_distorted_map) { + const bool include_distorted_map) const { if (include_distorted_map) { - return std::make_unique< - IdentityForComposition>( - IdentityMap{}); + return std::make_unique( + rigid_translation_map_); } else { return nullptr; } @@ -143,12 +172,15 @@ TimeDependentMapOptions::grid_to_distorted_map( TimeDependentMapOptions::MapType TimeDependentMapOptions::grid_to_inertial_map( - const bool include_distorted_map) const { + const bool include_distorted_map, const bool use_rigid_translation) const { if (include_distorted_map) { - return std::make_unique(shape_map_); + return std::make_unique(shape_map_, + rigid_translation_map_); + } else if (use_rigid_translation) { + return std::make_unique(rigid_translation_map_); } else { - return std::make_unique< - IdentityForComposition>(IdentityMap{}); + return std::make_unique( + linear_falloff_translation_map_); } } } // namespace domain::creators::sphere diff --git a/src/Domain/Creators/SphereTimeDependentMaps.hpp b/src/Domain/Creators/SphereTimeDependentMaps.hpp index e333d9acfad9..614933ed2713 100644 --- a/src/Domain/Creators/SphereTimeDependentMaps.hpp +++ b/src/Domain/Creators/SphereTimeDependentMaps.hpp @@ -14,6 +14,7 @@ #include "Domain/CoordinateMaps/CoordinateMap.hpp" #include "Domain/CoordinateMaps/Identity.hpp" #include "Domain/CoordinateMaps/TimeDependent/Shape.hpp" +#include "Domain/CoordinateMaps/TimeDependent/Translation.hpp" #include "Domain/FunctionsOfTime/FunctionOfTime.hpp" #include "Options/Auto.hpp" #include "Options/String.hpp" @@ -81,6 +82,7 @@ struct TimeDependentMapOptions { using IdentityMap = domain::CoordinateMaps::Identity<3>; // Time-dependent maps using ShapeMap = domain::CoordinateMaps::TimeDependent::Shape; + using TranslationMap = domain::CoordinateMaps::TimeDependent::Translation<3>; template using IdentityForComposition = @@ -88,14 +90,20 @@ struct TimeDependentMapOptions { using GridToDistortedComposition = domain::CoordinateMap; using GridToInertialComposition = - domain::CoordinateMap; + domain::CoordinateMap; + using GridToInertialSimple = + domain::CoordinateMap; + using DistortedToInertialComposition = + domain::CoordinateMap; public: using maps_list = tmpl::list, IdentityForComposition, IdentityForComposition, - GridToDistortedComposition, GridToInertialComposition>; + GridToDistortedComposition, GridToInertialSimple, + GridToInertialComposition, DistortedToInertialComposition>; /// \brief The initial time of the functions of time. struct InitialTime { @@ -131,7 +139,26 @@ struct TimeDependentMapOptions { std::optional> initial_values{}; }; - using options = tmpl::list; + struct TranslationMapOptions { + using type = TranslationMapOptions; + 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."}; + + struct InitialValues { + using type = std::array, 2>; + static constexpr Options::String help = { + "Initial values for the translation map."}; + }; + + using options = tmpl::list; + + std::array, 2> initial_values{}; + }; + + using options = + tmpl::list; static constexpr Options::String help{ "The options for all the hard-coded time dependent maps in the Sphere " "domain."}; @@ -139,7 +166,8 @@ struct TimeDependentMapOptions { TimeDependentMapOptions() = default; TimeDependentMapOptions(double initial_time, - const ShapeMapOptions& shape_map_options); + const ShapeMapOptions& shape_map_options, + const TranslationMapOptions& translation_map_options); /*! * \brief Create the function of time map using the options that were @@ -149,6 +177,7 @@ struct TimeDependentMapOptions { * * - Size: `PiecewisePolynomial<3>` * - Shape: `PiecewisePolynomial<2>` + * - Translation: `PiecewisePolynomial<2>` */ std::unordered_map> @@ -162,19 +191,22 @@ struct TimeDependentMapOptions { * Currently, this constructs a: * * - Shape: `Shape` (with a size function of time) + * - Translation: `Translation` */ - void build_maps(const std::array& center, double inner_radius, - double outer_radius); + void build_maps( + const std::array& center, double inner_radius, + double outer_radius, + std::pair translation_transition_radii); /*! * \brief This will construct the map from `Frame::Distorted` to * `Frame::Inertial`. * - * If the argument `include_distorted_map` is true, then this will be an - * identity map. If it is false, then this returns `nullptr`. + * If the argument `include_distorted_map` is true, then this will be a + * translation map. If it is false, then this returns `nullptr`. */ - static MapType distorted_to_inertial_map( - bool include_distorted_map); + MapType distorted_to_inertial_map( + bool include_distorted_map) const; /*! * \brief This will construct the map from `Frame::Grid` to @@ -195,16 +227,20 @@ struct TimeDependentMapOptions { * only be an identity map. */ MapType grid_to_inertial_map( - bool include_distorted_map) const; + bool include_distorted_map, bool use_rigid_translation) const; inline static const std::string size_name{"Size"}; inline static const std::string shape_name{"Shape"}; + inline static const std::string translation_name{"Translation"}; private: double initial_time_{std::numeric_limits::signaling_NaN()}; size_t initial_l_max_{}; ShapeMap shape_map_{}; + TranslationMap rigid_translation_map_{}; + TranslationMap linear_falloff_translation_map_{}; std::optional> initial_shape_values_{}; + std::array, 2> initial_translation_values_{}; }; } // namespace domain::creators::sphere diff --git a/tests/InputFiles/Xcts/KerrSchild.yaml b/tests/InputFiles/Xcts/KerrSchild.yaml index f99510c7602f..cb93995ba3f1 100644 --- a/tests/InputFiles/Xcts/KerrSchild.yaml +++ b/tests/InputFiles/Xcts/KerrSchild.yaml @@ -4,7 +4,7 @@ Executable: SolveXcts Testing: Check: parse;execute_check_output - Timeout: 40 + Timeout: 50 Priority: High ExpectedOutput: - KerrSchildReductions.h5 @@ -26,7 +26,7 @@ Background: &solution KerrSchild: Mass: &KerrMass 1. Spin: &KerrSpin [0., 0., 0.6] - Center: [0., 0., 0.] + Center: &Center [0., 0., 0.] InitialGuess: Flatness @@ -48,8 +48,8 @@ DomainCreator: UseEquiangularMap: True EquatorialCompression: None WhichWedges: All - RadialPartitioning: [] - RadialDistribution: [Logarithmic] + RadialPartitioning: [3.0] + RadialDistribution: [Linear, Logarithmic] TimeDependentMaps: InitialTime: 0. ShapeMap: @@ -57,6 +57,8 @@ DomainCreator: InitialValues: Mass: *KerrMass Spin: *KerrSpin + TranslationMap: + InitialValues: [*Center, [0., 0., 0.]] OuterBoundaryCondition: AnalyticSolution: Solution: *solution diff --git a/tests/Unit/Domain/Creators/Test_Sphere.cpp b/tests/Unit/Domain/Creators/Test_Sphere.cpp index 59ac389d7cd1..2989a66ed07f 100644 --- a/tests/Unit/Domain/Creators/Test_Sphere.cpp +++ b/tests/Unit/Domain/Creators/Test_Sphere.cpp @@ -139,6 +139,9 @@ std::string option_string( " ShapeMap:\n" " LMax: 10\n" " InitialValues: Spherical\n" + " TranslationMap:\n" + " InitialValues: [[0.0, 0.0, 0.0]," + " [0.001, -0.003, 0.005]]\n" : " TimeDependentMaps:\n" " UniformTranslation:\n" " InitialTime: 1.0\n" @@ -235,7 +238,8 @@ void test_sphere_construction( const bool expect_boundary_conditions = true, const std::vector& times = {0.}, const std::array& velocity = {{0., 0., 0.}}, - const std::array& size_values = {{0., 0., 0.}}) { + const std::array& size_values = {{0., 0., 0.}}, + const bool use_hard_coded_maps = true) { // check consistency of domain const auto domain = TestHelpers::domain::creators::test_domain_creator( sphere, expect_boundary_conditions, false, times); @@ -306,8 +310,8 @@ void test_sphere_construction( const double corner_distance_from_origin = [&block, &coords_on_spherical_partition, ¤t_time, &block_id, &num_blocks_per_shell, &inner_radius, &outer_radius, - &radial_partitioning, &size_values, &functions_of_time, - &velocity]() -> double { + &radial_partitioning, &size_values, &functions_of_time, &velocity, + &is_inner_cube, &num_blocks, &use_hard_coded_maps]() -> double { // use stationary map if independent of time if (not block.is_time_dependent()) { return get(magnitude( @@ -324,32 +328,76 @@ void test_sphere_construction( get(magnitude(inertial_location_time_dep)); const double delta_t = current_time - 1.0; - // Only when using hard coded time dependent maps do we have - // distorted maps in the inner shell. - if (block.has_distorted_frame()) { - CHECK(block_id < num_blocks_per_shell); - // This is the calculation of the inverse of the - // SphericalCompression map since the shape map is the identity - // currently - const double min_radius = inner_radius; - const double max_radius = radial_partitioning.size() > 0 - ? radial_partitioning[0] - : outer_radius; - const double func = - gsl::at(size_values, 0) + gsl::at(size_values, 1) * delta_t + - 0.5 * gsl::at(size_values, 2) * square(delta_t); - const double func_Y00 = func / sqrt(4.0 * M_PI); - const double max_minus_min = max_radius - min_radius; - const double scale = - (max_minus_min + func_Y00 * max_radius / target_radius) / - (max_minus_min + func_Y00); - - return target_radius * scale; - } else if (block.moving_mesh_grid_to_inertial_map().is_identity()) { - // This happens in our test if we have hard coded time dependent - // maps, but we aren't in the first shell. Then it's just the - // identity map - return target_radius; + double temp_radius = target_radius; + if (use_hard_coded_maps) { + // Undo the translation + if (block_id + num_blocks_per_shell < num_blocks and + not is_inner_cube) { + // For inner shells we are using a translation time dependence. + // origin in inertial frame (need to shift inertial + // coord by velocity * (final_time - initial_time)) + + temp_radius = sqrt(square(get<0>(inertial_location_time_dep) - + velocity[0] * (delta_t)) + + square(get<1>(inertial_location_time_dep) - + velocity[1] * (delta_t)) + + square(get<2>(inertial_location_time_dep) - + velocity[2] * (delta_t))); + } else { + // For the outer shell there is a linear transition that leaves + // the outer boundary fixed + const double min_radius = radial_partitioning.back(); + const double max_radius = outer_radius; + // The radial falloff factor is determined with respect to the + // grid frame + const auto grid_location_time_dep = + block.moving_mesh_logical_to_grid_map()( + coords_on_spherical_partition, current_time, + functions_of_time); + const double grid_target_radius = + get(magnitude(grid_location_time_dep)); + const double radial_falloff_factor = + (max_radius - grid_target_radius) / + (max_radius - min_radius); + + temp_radius = sqrt( + square(get<0>(inertial_location_time_dep) - + radial_falloff_factor * velocity[0] * (delta_t)) + + square(get<1>(inertial_location_time_dep) - + radial_falloff_factor * velocity[1] * (delta_t)) + + square(get<2>(inertial_location_time_dep) - + radial_falloff_factor * velocity[2] * (delta_t))); + } + + // Only when using hard coded time dependent maps do we have + // distorted maps in the inner shell. + if (block.has_distorted_frame()) { + CHECK(block_id < num_blocks_per_shell); + // This is the calculation of the inverse of the + // SphericalCompression map since the shape map is the identity + // currently + const double min_radius = inner_radius; + const double max_radius = not radial_partitioning.empty() + ? radial_partitioning[0] + : outer_radius; + const double func = + gsl::at(size_values, 0) + + gsl::at(size_values, 1) * delta_t + + 0.5 * gsl::at(size_values, 2) * square(delta_t); + const double func_Y00 = func / sqrt(4.0 * M_PI); + const double max_minus_min = max_radius - min_radius; + const double scale = + (max_minus_min + func_Y00 * max_radius / target_radius) / + (max_minus_min + func_Y00); + + return temp_radius * scale; + } else { + // This happens in our test if we have hard coded time dependent + // maps, but we aren't in the first shell. Then it's just the + // identity map + return temp_radius; + } + } else { // This means we are using a translation time dependence. // origin in inertial frame (need to shift inertial @@ -535,11 +583,15 @@ void test_parse_errors() { "an outflow-type boundary condition, you must use that.")); using TDMO = domain::creators::sphere::TimeDependentMapOptions; CHECK_THROWS_WITH( - creators::Sphere( - inner_radius, outer_radius, inner_cube, refinement, initial_extents, - use_equiangular_map, equatorial_compression, radial_partitioning, - radial_distribution, which_wedges, - TDMO{1.0, TDMO::ShapeMapOptions{5, std::nullopt}}, nullptr), + creators::Sphere(inner_radius, outer_radius, inner_cube, refinement, + initial_extents, use_equiangular_map, + equatorial_compression, radial_partitioning, + radial_distribution, which_wedges, + TDMO{1.0, TDMO::ShapeMapOptions{5, std::nullopt}, + TDMO::TranslationMapOptions{ + std::array{0.0, 0.0, 0.0}, + std::array{0.001, -0.003, 0.005}}}, + nullptr), Catch::Matchers::ContainsSubstring( "Currently cannot use hard-coded time dependent maps with an inner " "cube. Use a TimeDependence instead.")); @@ -580,7 +632,7 @@ void test_sphere() { const size_t l_max = 10; const std::vector times{1., 10.}; - for (auto [interior_index, equiangular, equatorial_compression, array_index, + for (auto [interior_index, equiangular, equatorial_compression, index, which_wedges, time_dependent, use_hard_coded_time_dep_options, with_boundary_conditions] : random_sample<5>( @@ -601,8 +653,6 @@ void test_sphere() { CAPTURE(fill_interior); CAPTURE(equiangular); CAPTURE(equatorial_compression.has_value()); // Whether map is present. - CAPTURE(radial_partitioning[array_index]); - CAPTURE(radial_distribution[array_index]); CAPTURE(which_wedges); CAPTURE(time_dependent); CAPTURE(with_boundary_conditions); @@ -614,28 +664,46 @@ void test_sphere() { } CAPTURE(use_hard_coded_time_dep_options); + // If we are using hard coded maps, we need at least two shells (or one + // radial partition) for the translation map. + const auto array_index = + (use_hard_coded_time_dep_options and + gsl::at(radial_partitioning, index).empty() + ? index + 1 + : index); + CAPTURE(gsl::at(radial_partitioning, array_index)); + CAPTURE(gsl::at(radial_distribution, array_index)); + if (which_wedges != ShellWedges::All and with_boundary_conditions) { continue; } creators::Sphere::RadialDistribution::type radial_distribution_variant; - if (radial_distribution[array_index].size() == 1) { - radial_distribution_variant = radial_distribution[array_index][0]; + if (gsl::at(radial_distribution, array_index).size() == 1) { + radial_distribution_variant = + gsl::at(radial_distribution, array_index)[0]; } else { - radial_distribution_variant = radial_distribution[array_index]; + radial_distribution_variant = gsl::at(radial_distribution, array_index); } std::optional time_dependent_options{}; + auto translation_velocity = std::array{{0., 0., 0.}}; + if (time_dependent) { if (use_hard_coded_time_dep_options) { + translation_velocity = std::array{{0.001, -0.003, 0.005}}; time_dependent_options = creators::sphere::TimeDependentMapOptions( - 1.0, creators::sphere::TimeDependentMapOptions::ShapeMapOptions{ - l_max, std::nullopt}); + 1.0, + creators::sphere::TimeDependentMapOptions::ShapeMapOptions{ + l_max, std::nullopt}, + creators::sphere::TimeDependentMapOptions::TranslationMapOptions{ + std::array{0.0, 0.0, 0.0}, translation_velocity}); } else { time_dependent_options = std::make_unique< domain::creators::time_dependence::UniformTranslation<3, 0>>( 1.0, velocity); + translation_velocity = velocity; } } @@ -647,23 +715,22 @@ void test_sphere() { initial_extents, equiangular, equatorial_compression, - radial_partitioning[array_index], + gsl::at(radial_partitioning, array_index), radial_distribution_variant, which_wedges, std::move(time_dependent_options), with_boundary_conditions ? create_boundary_condition(true) : nullptr}; - test_sphere_construction( - sphere, inner_radius, outer_radius, fill_interior, - radial_partitioning[array_index], which_wedges, - with_boundary_conditions, - time_dependent ? times : std::vector{0.}, - time_dependent ? velocity : std::array{{0., 0., 0.}}, - size_values); + test_sphere_construction(sphere, inner_radius, outer_radius, fill_interior, + gsl::at(radial_partitioning, array_index), + which_wedges, with_boundary_conditions, + time_dependent ? times : std::vector{0.}, + translation_velocity, size_values, + use_hard_coded_time_dep_options); TestHelpers::domain::creators::test_creation( option_string(inner_radius, outer_radius, interior, initial_refinement, initial_extents, equiangular, equatorial_compression, - radial_partitioning[array_index], - radial_distribution[array_index], which_wedges, + gsl::at(radial_partitioning, array_index), + gsl::at(radial_distribution, array_index), which_wedges, time_dependent, use_hard_coded_time_dep_options, with_boundary_conditions), sphere, with_boundary_conditions); diff --git a/tests/Unit/Domain/Creators/Test_SphereTimeDependentMaps.cpp b/tests/Unit/Domain/Creators/Test_SphereTimeDependentMaps.cpp index 239e292772f3..ceebb1030027 100644 --- a/tests/Unit/Domain/Creators/Test_SphereTimeDependentMaps.cpp +++ b/tests/Unit/Domain/Creators/Test_SphereTimeDependentMaps.cpp @@ -34,25 +34,29 @@ std::string create_option_string(const bool use_non_zero_shape) { (use_non_zero_shape ? " InitialValues:\n" " Mass: 1.0\n" " Spin: [0.0, 0.0, 0.0]\n"s - : " InitialValues: Spherical\n"s); + : " InitialValues: Spherical\n"s) + + "TranslationMap:\n" + " InitialValues: [[0.1, -3.2, 1.1], [-0.3, 0.5, -0.7]]\n"s; } void test(const bool use_non_zero_shape) { CAPTURE(use_non_zero_shape); const double initial_time = 1.5; const size_t l_max = 8; - TimeDependentMapOptions time_dep_options = - TestHelpers::test_creation( - create_option_string(use_non_zero_shape)); + auto time_dep_options = TestHelpers::test_creation( + create_option_string(use_non_zero_shape)); std::unordered_map expiration_times{ {TimeDependentMapOptions::shape_name, 15.5}, {TimeDependentMapOptions::size_name, + std::numeric_limits::infinity()}, + {TimeDependentMapOptions::translation_name, std::numeric_limits::infinity()}}; // These are hard-coded so this is just a regression test CHECK(TimeDependentMapOptions::size_name == "Size"s); CHECK(TimeDependentMapOptions::shape_name == "Shape"s); + CHECK(TimeDependentMapOptions::translation_name == "Translation"s); using PP2 = domain::FunctionsOfTime::PiecewisePolynomial<2>; using PP3 = domain::FunctionsOfTime::PiecewisePolynomial<3>; @@ -70,9 +74,19 @@ void test(const bool use_non_zero_shape) { std::array{shape_zeros, shape_zeros, shape_zeros}, expiration_times.at(TimeDependentMapOptions::shape_name)}; + DataVector initial_translation_center{{0.1, -3.2, 1.1}}; + DataVector initial_velocity{{-0.3, 0.5, -0.7}}; + PP2 translation_non_zero{ + initial_time, + std::array{ + {initial_translation_center, initial_velocity, {3, 0.0}}}, + expiration_times.at(TimeDependentMapOptions::translation_name)}; + const std::array center{-5.0, -0.01, -0.02}; const double inner_radius = 0.5; const double outer_radius = 2.1; + const double translation_inner_radius = 3.1; + const double translation_outer_radius = 3.9; const auto functions_of_time = time_dep_options.create_functions_of_time(inner_radius, expiration_times); @@ -86,35 +100,48 @@ void test(const bool use_non_zero_shape) { *functions_of_time.at(TimeDependentMapOptions::shape_name).get()) == shape_all_zero); + CHECK(dynamic_cast( + *functions_of_time.at(TimeDependentMapOptions::translation_name) + .get()) == translation_non_zero); + for (const bool include_distorted : make_array(true, false)) { - time_dep_options.build_maps(center, inner_radius, outer_radius); - - const auto grid_to_distorted_map = - time_dep_options.grid_to_distorted_map(include_distorted); - const auto grid_to_inertial_map = - time_dep_options.grid_to_inertial_map(include_distorted); - const auto distorted_to_inertial_map = - time_dep_options.distorted_to_inertial_map(include_distorted); - - // All of these maps are tested individually. Rather than going through the - // effort of coming up with a source coordinate and calculating analytically - // what we would get after it's mapped, we just check that whether it's - // supposed to be a nullptr and if it's not, if it's the identity and that - // the jacobians are time dependent. - const auto check_map = [](const auto& map, const bool is_null, - const bool is_identity) { - if (is_null) { - CHECK(map == nullptr); - } else { - CHECK(map->is_identity() == is_identity); - CHECK(map->inv_jacobian_is_time_dependent() == not is_identity); - CHECK(map->jacobian_is_time_dependent() == not is_identity); - } - }; - - check_map(grid_to_inertial_map, false, not include_distorted); - check_map(grid_to_distorted_map, not include_distorted, false); - check_map(distorted_to_inertial_map, not include_distorted, true); + for (const bool use_rigid_translation : make_array(true, false)) { + time_dep_options.build_maps( + center, inner_radius, outer_radius, + std::pair{translation_inner_radius, + translation_outer_radius}); + + const auto grid_to_distorted_map = + time_dep_options.grid_to_distorted_map(include_distorted); + const auto grid_to_inertial_map = time_dep_options.grid_to_inertial_map( + include_distorted, use_rigid_translation); + const auto distorted_to_inertial_map = + time_dep_options.distorted_to_inertial_map(include_distorted); + + // All of these maps are tested individually. Rather than going through + // the effort of coming up with a source coordinate and calculating + // analytically what we would get after it's mapped, we just check that + // whether it's supposed to be a nullptr and if it's not, if it's the + // identity and that the jacobians are time dependent. + const auto check_map = [](const auto& map, const bool is_null, + const bool is_identity) { + if (is_null) { + CHECK(map == nullptr); + } else { + CHECK(map->is_identity() == is_identity); + CHECK(map->inv_jacobian_is_time_dependent() == not is_identity); + CHECK(map->jacobian_is_time_dependent() == not is_identity); + } + }; + + // There is no null pointer in the grid to inertial map + check_map(grid_to_inertial_map, false, false); + + check_map(grid_to_distorted_map, not include_distorted, false); + + // If no shape distortion, there is only the translation + check_map(distorted_to_inertial_map, not include_distorted, false); + } } } diff --git a/tests/Unit/ParallelAlgorithms/Interpolation/Test_ComputeDestVars.cpp b/tests/Unit/ParallelAlgorithms/Interpolation/Test_ComputeDestVars.cpp index 1382722833c4..ec0d40c5e460 100644 --- a/tests/Unit/ParallelAlgorithms/Interpolation/Test_ComputeDestVars.cpp +++ b/tests/Unit/ParallelAlgorithms/Interpolation/Test_ComputeDestVars.cpp @@ -108,11 +108,18 @@ 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}}; + TDMO time_dep_opts{ + 0.0, TDMO::ShapeMapOptions{2, std::nullopt}, + TDMO::TranslationMapOptions{std::array{0.0, 0.0, 0.0}, + std::array{0.0, 0.0, 0.0}}}; const auto domain_creator = domain::creators::Sphere( - 0.9, 4.9, domain::creators::Sphere::Excision{}, 1_st, 7_st, false, {}, {}, - {}, {}, time_dep_opts); + 0.9, 4.9, domain::creators::Sphere::Excision{}, 1_st, 7_st, false, {}, + std::vector{2.0}, + std::vector{ + {domain::CoordinateMaps::Distribution::Linear, + domain::CoordinateMaps::Distribution::Linear}}, + {}, time_dep_opts); const Domain<3> domain = domain_creator.create_domain(); Parallel::GlobalCache cache{{domain_creator.functions_of_time()}};