From 35788b3c886e5e437df7b5c697ccc6969136b9c0 Mon Sep 17 00:00:00 2001 From: Mark Scheel Date: Mon, 30 Jan 2023 10:27:19 -0800 Subject: [PATCH] StrahlkorperCoords src can be Frame::Distorted. Formerly only Frame::Grid was allowed. --- .../StrahlkorperCoordsInDifferentFrame.cpp | 70 +++++++++---------- ...est_StrahlkorperCoordsInDifferentFrame.cpp | 53 +++++++++----- 2 files changed, 70 insertions(+), 53 deletions(-) diff --git a/src/ApparentHorizons/StrahlkorperCoordsInDifferentFrame.cpp b/src/ApparentHorizons/StrahlkorperCoordsInDifferentFrame.cpp index 2fbf39fdc92ba..52a6ee24926c8 100644 --- a/src/ApparentHorizons/StrahlkorperCoordsInDifferentFrame.cpp +++ b/src/ApparentHorizons/StrahlkorperCoordsInDifferentFrame.cpp @@ -8,6 +8,7 @@ #include "DataStructures/Tensor/Tensor.hpp" #include "DataStructures/Variables.hpp" #include "Domain/Block.hpp" +#include "Domain/BlockLogicalCoordinates.hpp" #include "Domain/CoordinateMaps/CoordinateMap.hpp" #include "Domain/Domain.hpp" #include "Domain/FunctionsOfTime/FunctionOfTime.hpp" @@ -15,6 +16,7 @@ #include "NumericalAlgorithms/SphericalHarmonics/Tags.hpp" #include "NumericalAlgorithms/SphericalHarmonics/YlmSpherepack.hpp" #include "Utilities/ContainerHelpers.hpp" +#include "Utilities/ErrorHandling/Assert.hpp" #include "Utilities/GenerateInstantiations.hpp" #include "Utilities/Gsl.hpp" @@ -27,6 +29,8 @@ void strahlkorper_coords_in_different_frame( std::string, std::unique_ptr>& functions_of_time, const double time) { + static_assert(std::is_same_v, + "Destination frame must currently be Inertial frame"); destructive_resize_components( dest_cartesian_coords, src_strahlkorper.ylm_spherepack().physical_size()); @@ -53,41 +57,34 @@ void strahlkorper_coords_in_different_frame( // We now wish to map src_cartesian_coords to the destination frame. // Each Block will have a different map, so the mapping must be done - // Block by Block. - for (const auto& block : domain.blocks()) { - // Once there are more possible source frames than the grid frame - // (i.e. the distorted frame), then this static_assert will change, - // and there will be an `if constexpr` below to treat the different - // possible source frames. - static_assert(std::is_same_v, - "Source frame must currently be Grid frame"); - static_assert(std::is_same_v, - "Destination frame must currently be Inertial frame"); - const auto& grid_to_inertial_map = block.moving_mesh_grid_to_inertial_map(); - const auto& logical_to_grid_map = block.moving_mesh_logical_to_grid_map(); - // Fill only those dest_cartesian_coords that are in this Block. - // Determine which coords are in this Block by checking if the - // inverse grid-to-logical map yields logical coords between -1 and 1. - for (size_t s = 0; s < get<0>(src_cartesian_coords).size(); ++s) { - const tnsr::I x_src{ - {get<0>(src_cartesian_coords)[s], get<1>(src_cartesian_coords)[s], - get<2>(src_cartesian_coords)[s]}}; - const auto x_logical = logical_to_grid_map.inverse(x_src); - // x_logical might be an empty std::optional. - if (x_logical.has_value() and get<0>(x_logical.value()) <= 1.0 and - get<0>(x_logical.value()) >= -1.0 and - get<1>(x_logical.value()) <= 1.0 and - get<1>(x_logical.value()) >= -1.0 and - get<2>(x_logical.value()) <= 1.0 and - get<2>(x_logical.value()) >= -1.0) { - const auto x_dest = - grid_to_inertial_map(x_src, time, functions_of_time); - get<0>(*dest_cartesian_coords)[s] = get<0>(x_dest); - get<1>(*dest_cartesian_coords)[s] = get<1>(x_dest); - get<2>(*dest_cartesian_coords)[s] = get<2>(x_dest); - // Note that if a point is on the boundary of two or more - // Blocks, it might get filled twice, but that's ok. - } + // Block by Block. Since each point in general has a different + // block (given by block_logical_coordinates), this means we must + // map point by point. + const auto block_logical_coords = block_logical_coordinates( + domain, src_cartesian_coords, time, functions_of_time); + for (size_t s = 0; s < block_logical_coords.size(); ++s) { + const tnsr::I x_src{{get<0>(src_cartesian_coords)[s], + get<1>(src_cartesian_coords)[s], + get<2>(src_cartesian_coords)[s]}}; + ASSERT(block_logical_coords[s].has_value(), + "Found a point (source coords " << x_src + << ") that is not in any Block."); + const auto& block = + domain.blocks()[block_logical_coords[s].value().id.get_index()]; + if constexpr (std::is_same_v) { + const auto x_dest = block.moving_mesh_grid_to_inertial_map()( + x_src, time, functions_of_time); + get<0>(*dest_cartesian_coords)[s] = get<0>(x_dest); + get<1>(*dest_cartesian_coords)[s] = get<1>(x_dest); + get<2>(*dest_cartesian_coords)[s] = get<2>(x_dest); + } else { + static_assert(std::is_same_v, + "Src frame must be Distorted if it is not Grid"); + const auto x_dest = block.moving_mesh_distorted_to_inertial_map()( + x_src, time, functions_of_time); + get<0>(*dest_cartesian_coords)[s] = get<0>(x_dest); + get<1>(*dest_cartesian_coords)[s] = get<1>(x_dest); + get<2>(*dest_cartesian_coords)[s] = get<2>(x_dest); } } } @@ -107,7 +104,8 @@ void strahlkorper_coords_in_different_frame( functions_of_time, \ const double time); -GENERATE_INSTANTIATIONS(INSTANTIATE, (::Frame::Grid), (::Frame::Inertial)) +GENERATE_INSTANTIATIONS(INSTANTIATE, (::Frame::Grid, ::Frame::Distorted), + (::Frame::Inertial)) #undef INSTANTIATE #undef DESTFRAME diff --git a/tests/Unit/ApparentHorizons/Test_StrahlkorperCoordsInDifferentFrame.cpp b/tests/Unit/ApparentHorizons/Test_StrahlkorperCoordsInDifferentFrame.cpp index 97e36dfc45573..a37d37ee3bf26 100644 --- a/tests/Unit/ApparentHorizons/Test_StrahlkorperCoordsInDifferentFrame.cpp +++ b/tests/Unit/ApparentHorizons/Test_StrahlkorperCoordsInDifferentFrame.cpp @@ -23,17 +23,18 @@ namespace { +template void test_strahlkorper_coords_in_different_frame() { const size_t grid_points_each_dimension = 5; // Set up a Strahlkorper corresponding to a Schwarzschild hole of - // mass 1, in the grid frame. + // mass 1, in the source frame. // Center the Strahlkorper at (0.03,0.02,0.01) so that we test a // nonzero center. - const std::array strahlkorper_grid_center = {0.03, 0.02, 0.01}; + const std::array strahlkorper_src_center = {0.03, 0.02, 0.01}; const size_t l_max = 8; - const Strahlkorper strahlkorper_grid(l_max, 2.0, - strahlkorper_grid_center); + const Strahlkorper strahlkorper_src(l_max, 2.0, + strahlkorper_src_center); // Create a Domain. // We choose a spherical shell domain extending from radius 1.9M to @@ -43,15 +44,32 @@ void test_strahlkorper_coords_in_different_frame() { std::vector radial_partitioning{}; std::vector radial_distribution{ domain::CoordinateMaps::Distribution::Linear}; + + // In computing the time_dependence, make sure that src-to-inertial + // velocity is (0.01,0.02,0.03) to agree with the analytic checks + // below. If src is distorted frame, then grid-to-distorted + // velocity doesn't matter for the value of the check (but does matter + // in terms of which points are in which blocks). + std::unique_ptr> + time_dependence; + if constexpr (std::is_same_v) { + time_dependence = std::make_unique< + domain::creators::time_dependence::UniformTranslation<3>>( + 0.0, std::array({{0.01, 0.02, 0.03}})); + } else { + static_assert(std::is_same_v, + "Src frame must be Distorted if it is not Grid"); + time_dependence = std::make_unique< + domain::creators::time_dependence::UniformTranslation<3>>( + 0.0, std::array({{-0.02, -0.01, -0.01}}), + std::array({{0.01, 0.02, 0.03}})); + } domain::creators::Shell domain_creator( 1.9, 2.9, 1, std::array{grid_points_each_dimension, grid_points_each_dimension}, false, {{1.0, 2}}, radial_partitioning, radial_distribution, - ShellWedges::All, - std::make_unique< - domain::creators::time_dependence::UniformTranslation<3>>( - 0.0, std::array({{0.01, 0.02, 0.03}}))); + ShellWedges::All, std::move(time_dependence)); Domain<3> domain = domain_creator.create_domain(); const auto functions_of_time = domain_creator.functions_of_time(); @@ -60,18 +78,18 @@ void test_strahlkorper_coords_in_different_frame() { tnsr::I inertial_coords{}; strahlkorper_coords_in_different_frame(make_not_null(&inertial_coords), - strahlkorper_grid, domain, + strahlkorper_src, domain, functions_of_time, time); - // Now compare with expected result, which is the grid-frame coords of + // Now compare with expected result, which is the src-frame coords of // the Strahlkorper translated by (0.005,0.01,0.015). - const auto grid_coords = StrahlkorperFunctions::cartesian_coords( - strahlkorper_grid, StrahlkorperFunctions::radius(strahlkorper_grid), + const auto src_coords = StrahlkorperFunctions::cartesian_coords( + strahlkorper_src, StrahlkorperFunctions::radius(strahlkorper_src), StrahlkorperFunctions::rhat( - StrahlkorperFunctions::theta_phi(strahlkorper_grid))); - CHECK_ITERABLE_APPROX(get<0>(grid_coords) + 0.005, get<0>(inertial_coords)); - CHECK_ITERABLE_APPROX(get<1>(grid_coords) + 0.01, get<1>(inertial_coords)); - CHECK_ITERABLE_APPROX(get<2>(grid_coords) + 0.015, get<2>(inertial_coords)); + StrahlkorperFunctions::theta_phi(strahlkorper_src))); + CHECK_ITERABLE_APPROX(get<0>(src_coords) + 0.005, get<0>(inertial_coords)); + CHECK_ITERABLE_APPROX(get<1>(src_coords) + 0.01, get<1>(inertial_coords)); + CHECK_ITERABLE_APPROX(get<2>(src_coords) + 0.015, get<2>(inertial_coords)); } SPECTRE_TEST_CASE("Unit.ApparentHorizons.StrahlkorperCoordsInDifferentFrame", @@ -79,7 +97,8 @@ SPECTRE_TEST_CASE("Unit.ApparentHorizons.StrahlkorperCoordsInDifferentFrame", domain::creators::register_derived_with_charm(); domain::creators::time_dependence::register_derived_with_charm(); domain::FunctionsOfTime::register_derived_with_charm(); - test_strahlkorper_coords_in_different_frame(); + test_strahlkorper_coords_in_different_frame(); + test_strahlkorper_coords_in_different_frame(); } } // namespace