diff --git a/src/Domain/BlockLogicalCoordinates.cpp b/src/Domain/BlockLogicalCoordinates.cpp index 091b665fe2e8..058295dbb7e3 100644 --- a/src/Domain/BlockLogicalCoordinates.cpp +++ b/src/Domain/BlockLogicalCoordinates.cpp @@ -21,7 +21,8 @@ template std::optional> block_logical_coordinates_single_point( const tnsr::I& input_point, const Block& block, - const double time, const domain::FunctionsOfTimeMap& functions_of_time) { + const double time, const domain::FunctionsOfTimeMap& functions_of_time, + const bool allow_extrapolation) { std::optional> logical_point{}; if (block.is_time_dependent()) { if constexpr (std::is_same_v) { @@ -125,7 +126,7 @@ block_logical_coordinates_single_point( logical_point->get(d) = -1.0; continue; } - if (abs(logical_point->get(d)) > 1.0) { + if (abs(logical_point->get(d)) > 1.0 and not allow_extrapolation) { return std::nullopt; } } @@ -133,6 +134,49 @@ block_logical_coordinates_single_point( return logical_point; } +template +BlockLogicalCoords block_logical_coordinates_in_excision( + const tnsr::I& input_point, + const ExcisionSphere& excision_sphere, + const std::vector>& blocks, const double time, + const domain::FunctionsOfTimeMap& functions_of_time) { + // Comment by NV (Apr 2024): This function can be made more robust by first + // checking if the point is inside the excision sphere at all, but the + // excision sphere doesn't currently have this information. It needs the + // grid-to-inertial map including the deformation of the excision sphere to + // determine this. + for (const auto& [block_id, direction] : + excision_sphere.abutting_directions()) { + auto x_logical = block_logical_coordinates_single_point( + input_point, blocks[block_id], time, functions_of_time, true); + if (not x_logical.has_value()) { + continue; + } + // Discard block if the point has angular logical coordinates outside the + // range [-1, 1] + for (size_t d = 0; d < Dim; ++d) { + if (d != direction.dimension() and std::abs(x_logical->get(d)) > 1.) { + x_logical = std::nullopt; + break; + } + } + if (not x_logical.has_value()) { + continue; + } + // Discard block if the point is radially inside the block or on the other + // side of the excision sphere + const double radial_distance_this_block = + x_logical->get(direction.dimension()) * direction.sign(); + if (radial_distance_this_block < 1.) { + continue; + } + // The checks above should leave only 1 valid block, so return that + return make_id_pair(domain::BlockId(block_id), + std::move(x_logical.value())); + } + return std::nullopt; +} + template std::vector> block_logical_coordinates( const Domain& domain, const tnsr::I& x, @@ -174,11 +218,18 @@ std::vector> block_logical_coordinates( block_logical_coordinates_single_point( \ const tnsr::I& input_point, \ const Block& block, const double time, \ - const domain::FunctionsOfTimeMap& functions_of_time); \ + const domain::FunctionsOfTimeMap& functions_of_time, \ + bool allow_extrapolation); \ template std::vector> \ block_logical_coordinates( \ const Domain& domain, \ const tnsr::I& x, const double time, \ + const domain::FunctionsOfTimeMap& functions_of_time); \ + template BlockLogicalCoords \ + block_logical_coordinates_in_excision( \ + const tnsr::I& input_point, \ + const ExcisionSphere& excision_sphere, \ + const std::vector>& blocks, const double time, \ const domain::FunctionsOfTimeMap& functions_of_time); GENERATE_INSTANTIATIONS(INSTANTIATE, (1, 2, 3), diff --git a/src/Domain/BlockLogicalCoordinates.hpp b/src/Domain/BlockLogicalCoordinates.hpp index b322f218ae2c..e2b8ee60b629 100644 --- a/src/Domain/BlockLogicalCoordinates.hpp +++ b/src/Domain/BlockLogicalCoordinates.hpp @@ -21,6 +21,8 @@ template class Domain; template class Block; +template +class ExcisionSphere; /// \endcond template @@ -69,10 +71,32 @@ auto block_logical_coordinates( const domain::FunctionsOfTimeMap& functions_of_time = {}) -> std::vector>; +/// \par Extrapolation +/// If `allow_extrapolation` is `true` (default is `false`), then the function +/// can return block-logical coordinates outside of [-1, 1]. This is useful for +/// extrapolating into excision regions. template std::optional> block_logical_coordinates_single_point( const tnsr::I& input_point, const Block& block, double time = std::numeric_limits::signaling_NaN(), - const domain::FunctionsOfTimeMap& functions_of_time = {}); + const domain::FunctionsOfTimeMap& functions_of_time = {}, + bool allow_extrapolation = false); /// @} + +/*! + * \brief Finds the block-logical coordinates of a point in an excision sphere. + * + * The point will be reported to belong to the nearest block that abuts the + * excision sphere and will have a logical coordinate outside the range [-1, 1] + * in the radial direction. This is useful for extrapolating into excision + * regions using the Lagrange polynomial basis of the nearest element. If the + * point is not within the excision sphere, the return value is std::nullopt. + */ +template +BlockLogicalCoords block_logical_coordinates_in_excision( + const tnsr::I& input_point, + const ExcisionSphere& excision_sphere, + const std::vector>& blocks, + double time = std::numeric_limits::signaling_NaN(), + const domain::FunctionsOfTimeMap& functions_of_time = {}); diff --git a/src/Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/SphereTransition.cpp b/src/Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/SphereTransition.cpp index 484f5941aef5..ce096373b815 100644 --- a/src/Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/SphereTransition.cpp +++ b/src/Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/SphereTransition.cpp @@ -65,10 +65,12 @@ std::optional SphereTransition::original_radius_over_radius( if ((original_radius + eps_) >= r_min_ and (original_radius - eps_) <= r_max_) { return std::optional{original_radius / mag}; - } else if ((a_ > 0.0 and mag > r_max_) or (a_ < 0.0 and mag < r_min_)) { + } else if ((a_ > 0.0 and original_radius > r_max_) or + (a_ < 0.0 and original_radius < r_min_)) { // a_ being positive is a sentinel for reverse (see constructor) return {1.0 + radial_distortion / mag}; - } else if ((a_ < 0.0 and mag > r_max_) or (a_ > 0.0 and mag < r_min_)) { + } else if ((a_ < 0.0 and original_radius > r_max_) or + (a_ > 0.0 and original_radius < r_min_)) { return {1.0}; } else { return std::nullopt; diff --git a/src/Domain/ElementLogicalCoordinates.cpp b/src/Domain/ElementLogicalCoordinates.cpp index d039556ec573..c61ce0472024 100644 --- a/src/Domain/ElementLogicalCoordinates.cpp +++ b/src/Domain/ElementLogicalCoordinates.cpp @@ -27,36 +27,29 @@ element_logical_coordinates( const ElementId& element_id) { tnsr::I x_element_logical{}; for (size_t d = 0; d < Dim; ++d) { - // Check if the point is outside the element - const double up = element_id.segment_id(d).endpoint(Side::Upper); - const double lo = element_id.segment_id(d).endpoint(Side::Lower); - if (x_block_logical.get(d) < lo or x_block_logical.get(d) > up) { + // Check if the point is inside the element in this dimension. + // The following conditions are valid: + // - The block-logical coordinates are within the element's bounds. + // - The block-logical coordinate is outside the block and we're in the + // outermost element in that direction. This means we want to extrapolate + // out of the block. + const auto& segment_id = element_id.segment_id(d); + const double up = segment_id.endpoint(Side::Upper); + const double lo = segment_id.endpoint(Side::Lower); + if ((x_block_logical.get(d) >= lo and x_block_logical.get(d) <= up) or + (x_block_logical.get(d) < -1. and segment_id.index() == 0) or + (x_block_logical.get(d) > 1. and + segment_id.index() == two_to_the(segment_id.refinement_level()) - 1)) { + // Map to element logical coords + x_element_logical.get(d) = + (2.0 * x_block_logical.get(d) - up - lo) / (up - lo); + } else { return std::nullopt; } - // Map to element logical coords - x_element_logical.get(d) = - (2.0 * x_block_logical.get(d) - up - lo) / (up - lo); } return x_element_logical; } -namespace { -// The segments bounds are binary fractions (i.e. the numerator is an -// integer and the denominator is a power of 2) so these floating point -// comparisons should be safe from roundoff problems -// Need to return true if on upper face of block -// Otherwise return true if point is within the segment or on the lower bound -bool segment_contains(const double x_block_logical, - const double lower_bound_block_logical, - const double upper_bound_block_logical) { - if (UNLIKELY(x_block_logical == upper_bound_block_logical)) { - return (upper_bound_block_logical == 1.0); - } - return (x_block_logical >= lower_bound_block_logical and - x_block_logical < upper_bound_block_logical); -} -} // namespace - template std::unordered_map, ElementLogicalCoordHolder> element_logical_coordinates( @@ -90,13 +83,15 @@ element_logical_coordinates( if (not x_elem.has_value()) { continue; } - // Disambiguate points on shared element boundaries + // Disambiguate points on shared element boundaries. To do this, + // associate boundary points with the element with the lower index. bool is_contained = true; for (size_t d = 0; d < Dim; ++d) { - const double up = element_id.segment_id(d).endpoint(Side::Upper); - const double lo = element_id.segment_id(d).endpoint(Side::Lower); - const double x_block_log = x_block_logical.get(d); - if (not segment_contains(x_block_log, lo, up)) { + const auto& segment_id = element_id.segment_id(d); + const double up = segment_id.endpoint(Side::Upper); + if (x_block_logical.get(d) == up and + segment_id.index() != + two_to_the(segment_id.refinement_level()) - 1) { is_contained = false; break; } diff --git a/src/Domain/ElementLogicalCoordinates.hpp b/src/Domain/ElementLogicalCoordinates.hpp index 4d6f8d811249..aa98a3dd345d 100644 --- a/src/Domain/ElementLogicalCoordinates.hpp +++ b/src/Domain/ElementLogicalCoordinates.hpp @@ -32,6 +32,11 @@ class IdPair; * have element logical coordinates of -1 or 1. See the other function overload * for handling multiple points and disambiguating points on shared element * boundaries. + * + * Points with block-logical coordinates outside [-1, 1] are considered to + * belong to the outermost element of the block in that direction, and will + * therefore return an element logical coordinate outside of [-1, 1] as well. + * This is useful for extrapolation. */ template std::optional> diff --git a/src/IO/Exporter/CMakeLists.txt b/src/IO/Exporter/CMakeLists.txt index 41b790d73dbc..3c02911b4b07 100644 --- a/src/IO/Exporter/CMakeLists.txt +++ b/src/IO/Exporter/CMakeLists.txt @@ -34,6 +34,7 @@ target_link_libraries( FunctionsOfTime H5 Interpolation + Options Serialization Utilities ) diff --git a/src/IO/Exporter/Exporter.cpp b/src/IO/Exporter/Exporter.cpp index 2cad6bfa4f6d..6fcf0a1a7b0f 100644 --- a/src/IO/Exporter/Exporter.cpp +++ b/src/IO/Exporter/Exporter.cpp @@ -21,6 +21,8 @@ #include "IO/H5/VolumeData.hpp" #include "NumericalAlgorithms/Interpolation/IrregularInterpolant.hpp" #include "NumericalAlgorithms/Interpolation/PolynomialInterpolation.hpp" +#include "Options/Context.hpp" +#include "Options/Options.hpp" #include "Utilities/FileSystem.hpp" #include "Utilities/GenerateInstantiations.hpp" #include "Utilities/GetOutput.hpp" @@ -28,6 +30,59 @@ #include "Utilities/Overloader.hpp" #include "Utilities/Serialization/Serialize.hpp" +namespace spectre::Exporter { + +std::ostream& operator<<(std::ostream& os, + const ExcisionExtrapolationMode& mode) { + switch (mode) { + case ExcisionExtrapolationMode::NoExtrapolation: + os << "NoExtrapolation"; + break; + case ExcisionExtrapolationMode::NearestElement: + os << "NearestElement"; + break; + case ExcisionExtrapolationMode::RadialAnchors: + os << "RadialAnchors"; + break; + default: + os << "Unknown"; + break; + } + return os; +} + +} // namespace spectre::Exporter + +template <> +spectre::Exporter::ExcisionExtrapolationMode +Options::create_from_yaml::create< + void>(const Options::Option& options) { + const auto type_read = options.parse_as(); + for (const auto value : + {spectre::Exporter::ExcisionExtrapolationMode::NoExtrapolation, + spectre::Exporter::ExcisionExtrapolationMode::NearestElement, + spectre::Exporter::ExcisionExtrapolationMode::RadialAnchors}) { + if (type_read == get_output(value)) { + return value; + } + } + PARSE_ERROR( + options.context(), + "Failed to convert \"" + << type_read + << "\" to spectre::Exporter::ExcisionExtrapolationMode.\nMust be one " + "of " + << get_output( + spectre::Exporter::ExcisionExtrapolationMode::NoExtrapolation) + << ", " + << get_output( + spectre::Exporter::ExcisionExtrapolationMode::NearestElement) + << ", or " + << get_output( + spectre::Exporter::ExcisionExtrapolationMode::RadialAnchors) + << "."); +} + // Ignore OpenMP pragmas when OpenMP is not enabled #pragma GCC diagnostic push #pragma GCC diagnostic ignored "-Wunknown-pragmas" @@ -329,7 +384,7 @@ std::vector> interpolate_to_points( const std::variant& observation, const std::vector& tensor_components, const std::array, Dim>& target_points, - const bool extrapolate_into_excisions, + const ExcisionExtrapolationMode excision_extrapolation_mode, const std::optional num_threads) { domain::creators::register_derived_with_charm(); domain::creators::time_dependence::register_derived_with_charm(); @@ -446,19 +501,37 @@ std::vector> interpolate_to_points( } } // for blocks if (block_logical_coords[s].has_value() or - not extrapolate_into_excisions) { + excision_extrapolation_mode == + ExcisionExtrapolationMode::NoExtrapolation) { continue; } - // The point wasn't found in any block. Check if it's in an excision and - // set up extrapolation if requested. + // The point wasn't found in any block. Extrapolate if requested. for (const auto& [name, excision_sphere] : domain.excision_spheres()) { - if (add_extrapolation_anchors( - make_not_null(&extra_block_logical_coords), - make_not_null(&extra_extrapolation_info), excision_sphere, - domain, target_point, time, functions_of_time, - extrapolation_spacing)) { - extra_extrapolation_info.back().target_index = s; - break; + if (excision_extrapolation_mode == + ExcisionExtrapolationMode::NearestElement) { + // Get the block-logical coordinates for the nearest element (the + // coordinates will be outside [-1, 1], so Lagrange extrapolation will + // happen) + block_logical_coords[s] = block_logical_coordinates_in_excision( + target_point, excision_sphere, domain.blocks(), time, + functions_of_time); + if (block_logical_coords[s].has_value()) { + break; + } + } else if (excision_extrapolation_mode == + ExcisionExtrapolationMode::RadialAnchors) { + // Set up extrapolation anchors for this point + if (add_extrapolation_anchors( + make_not_null(&extra_block_logical_coords), + make_not_null(&extra_extrapolation_info), excision_sphere, + domain, target_point, time, functions_of_time, + extrapolation_spacing)) { + extra_extrapolation_info.back().target_index = s; + break; + } + } else { + ERROR("Invalid excision extrapolation mode: " + << excision_extrapolation_mode); } } } // omp for target points @@ -503,7 +576,7 @@ std::vector> interpolate_to_points( } } - if (extrapolate_into_excisions) { + if (excision_extrapolation_mode == ExcisionExtrapolationMode::RadialAnchors) { // Extrapolate into excisions from the anchor points #pragma omp parallel for num_threads(resolved_num_threads) for (const auto& extrapolation : extrapolation_info) { @@ -540,7 +613,7 @@ std::vector> interpolate_to_points( const std::variant& observation, \ const std::vector& tensor_components, \ const std::array, DIM(data)>& target_points, \ - bool extrapolate_into_excisions, \ + ExcisionExtrapolationMode excision_extrapolation_mode, \ const std::optional num_threads); GENERATE_INSTANTIATIONS(INSTANTIATE, (1, 2, 3)) diff --git a/src/IO/Exporter/Exporter.hpp b/src/IO/Exporter/Exporter.hpp index a1a86c38425e..5f73f052ee3f 100644 --- a/src/IO/Exporter/Exporter.hpp +++ b/src/IO/Exporter/Exporter.hpp @@ -15,6 +15,14 @@ #include #include +/// \cond +namespace Options { +class Option; +template +struct create_from_yaml; +} // namespace Options +/// \endcond + /// Functions that are intended to be used by external programs, e.g. to /// interpolate data in volume files to target points. namespace spectre::Exporter { @@ -32,6 +40,46 @@ struct ObservationStep { int value; }; +enum class ExcisionExtrapolationMode { + /*! + * \brief No extrapolation into the excisions is performed + */ + NoExtrapolation, + /*! + * \brief Extrapolate into excisions from the nearest element + * + * Uses the Lagrange basis in the nearest element to extrapolate into the + * excision region. This can be unstable when the element is too small (e.g. + * when it is h-refined a lot) and/or has log or inverse radial map. In those + * cases the logical coordinate that we're extrapolating to can be many + * logical element sizes away. Also, this relies on inverting the neighboring + * block's map outside the block, which is not guaranteed to work. + * + * This mode can be useful for extrapolating only very slightly into the + * excision, e.g. if the excision in the initial data is the apparent horizon + * and the evolution needs a little space around it to find and track the + * horizon. + */ + NearestElement, + /*! + * \brief Extrapolate into excisions from radial anchor points + * + * Uses the strategy from \cite Etienne2008re (adjusted for distorted + * excisions): we choose uniformly spaced radial anchor points spaced as + * $\Delta r = 0.3 r_\mathrm{AH}$ in the grid frame (where the excision is + * spherical), then map the anchor points to the distorted frame (where we + * have the target point) and do a 7th order polynomial extrapolation into the + * excision region. + * + * This mode is the most stable and accurate to fill the entire excision for + * moving puncture codes. + */ + RadialAnchors, +}; + +std::ostream& operator<<(std::ostream& os, + const ExcisionExtrapolationMode& mode); + /*! * \brief Interpolate data in volume files to target points * @@ -47,15 +95,12 @@ struct ObservationStep { * "Lapse", "Shift_x", "Shift_y", "Shift_z", "SpatialMetric_xx", etc. * Look into the H5 file to see what components are available. * \param target_points The points to interpolate to, in inertial coordinates. - * \param extrapolate_into_excisions Enables extrapolation into excision regions - * of the domain (default is `false`). This can be useful to fill the excision - * region with (constraint-violating but smooth) data so it can be imported into - * moving puncture codes. Specifically, we implement the strategy used in - * \cite Etienne2008re adjusted for distorted excisions: we choose uniformly - * spaced radial anchor points spaced as $\Delta r = 0.3 r_\mathrm{AH}$ in the - * grid frame (where the excision is spherical), then map the anchor points to - * the distorted frame (where we have the target point) and do a 7th order - * polynomial extrapolation into the excision region. + * \param excision_extrapolation_mode Enables extrapolation into excision + * regions of the domain (default is `NoExtrapolation`). This can be useful to + * fill the excision region with (constraint-violating but smooth) data so it + * can be imported into moving puncture codes, or to extrapolate a little bit + * into the excision to help the evolution find and track the apparent horizon. + * See `ExcisionExtrapolationMode` for the available modes. * \param num_threads The number of threads to use if OpenMP is linked in. If * not specified, OpenMP will determine the number of threads automatically. * It's also possible to set the number of threads using the environment @@ -73,7 +118,22 @@ std::vector> interpolate_to_points( const std::variant& observation, const std::vector& tensor_components, const std::array, Dim>& target_points, - bool extrapolate_into_excisions = false, + ExcisionExtrapolationMode excision_extrapolation_mode = + ExcisionExtrapolationMode::NoExtrapolation, std::optional num_threads = std::nullopt); } // namespace spectre::Exporter + +template <> +struct Options::create_from_yaml { + template + static spectre::Exporter::ExcisionExtrapolationMode create( + const Options::Option& options) { + return create(options); + } +}; + +template <> +spectre::Exporter::ExcisionExtrapolationMode +Options::create_from_yaml::create< + void>(const Options::Option& options); diff --git a/src/IO/Exporter/InterpolateToPoints.hpp b/src/IO/Exporter/InterpolateToPoints.hpp index e4290cb5a1b7..3afdce2d779c 100644 --- a/src/IO/Exporter/InterpolateToPoints.hpp +++ b/src/IO/Exporter/InterpolateToPoints.hpp @@ -37,7 +37,8 @@ tuples::tagged_tuple_from_typelist interpolate_to_points( const std::string& subfile_name, const std::variant& observation, const tnsr::I& target_points, - bool extrapolate_into_excisions = false, + ExcisionExtrapolationMode excision_extrapolation_mode = + ExcisionExtrapolationMode::NoExtrapolation, std::optional num_threads = std::nullopt) { // Convert target_points to an array of vectors // Possible performance optimization: avoid copying the data by allowing @@ -65,7 +66,7 @@ tuples::tagged_tuple_from_typelist interpolate_to_points( // Interpolate! const auto interpolated_data = interpolate_to_points( volume_files_or_glob, subfile_name, observation, tensor_components, - target_points_array, extrapolate_into_excisions, num_threads); + target_points_array, excision_extrapolation_mode, num_threads); // Convert the interpolated data to a tagged_tuple tuples::tagged_tuple_from_typelist result{}; size_t component_index = 0; diff --git a/src/IO/Exporter/Python/Bindings.cpp b/src/IO/Exporter/Python/Bindings.cpp index 166f8d41cd9c..55a81f5c40b5 100644 --- a/src/IO/Exporter/Python/Bindings.cpp +++ b/src/IO/Exporter/Python/Bindings.cpp @@ -14,6 +14,14 @@ namespace py = pybind11; PYBIND11_MODULE(_Pybindings, m) { // NOLINT enable_segfault_handler(); + py::enum_( + m, "ExcisionExtrapolationMode") + .value("NoExtrapolation", + spectre::Exporter::ExcisionExtrapolationMode::NoExtrapolation) + .value("NearestElement", + spectre::Exporter::ExcisionExtrapolationMode::NearestElement) + .value("RadialAnchors", + spectre::Exporter::ExcisionExtrapolationMode::RadialAnchors); m.def( "interpolate_to_points", [](const std::variant, std::string>& @@ -21,7 +29,8 @@ PYBIND11_MODULE(_Pybindings, m) { // NOLINT const std::string& subfile_name, const size_t observation_id, const std::vector& tensor_components, std::vector> target_points, - bool extrapolate_into_excisions, + spectre::Exporter::ExcisionExtrapolationMode + excision_extrapolation_mode, const std::optional& num_threads) { const size_t dim = target_points.size(); const spectre::Exporter::ObservationId obs_id{observation_id}; @@ -29,17 +38,17 @@ PYBIND11_MODULE(_Pybindings, m) { // NOLINT return spectre::Exporter::interpolate_to_points( volume_files_or_glob, subfile_name, obs_id, tensor_components, make_array, 1>(std::move(target_points)), - extrapolate_into_excisions, num_threads); + excision_extrapolation_mode, num_threads); } else if (dim == 2) { return spectre::Exporter::interpolate_to_points( volume_files_or_glob, subfile_name, obs_id, tensor_components, make_array, 2>(std::move(target_points)), - extrapolate_into_excisions, num_threads); + excision_extrapolation_mode, num_threads); } else if (dim == 3) { return spectre::Exporter::interpolate_to_points( volume_files_or_glob, subfile_name, obs_id, tensor_components, make_array, 3>(std::move(target_points)), - extrapolate_into_excisions, num_threads); + excision_extrapolation_mode, num_threads); } else { ERROR("Invalid dimension of target points: " << dim @@ -50,6 +59,8 @@ PYBIND11_MODULE(_Pybindings, m) { // NOLINT }, py::arg("volume_files_or_glob"), py::arg("subfile_name"), py::arg("observation_id"), py::arg("tensor_components"), - py::arg("target_points"), py::arg("extrapolate_into_excisions") = false, + py::arg("target_points"), + py::arg("excision_extrapolation_mode") = + spectre::Exporter::ExcisionExtrapolationMode::NoExtrapolation, py::arg("num_threads") = std::nullopt); } diff --git a/src/IO/Exporter/Python/InterpolateToPoints.py b/src/IO/Exporter/Python/InterpolateToPoints.py index bb2cd9768abf..5d1b017e591d 100644 --- a/src/IO/Exporter/Python/InterpolateToPoints.py +++ b/src/IO/Exporter/Python/InterpolateToPoints.py @@ -8,7 +8,7 @@ import rich import spectre.IO.H5 as spectre_h5 -from spectre.IO.Exporter import interpolate_to_points +from spectre.IO.Exporter import ExcisionExtrapolationMode, interpolate_to_points from spectre.Visualization.OpenVolfiles import ( open_volfiles_command, parse_points, @@ -43,12 +43,20 @@ ) @click.option( "--extrapolate-into-excisions", - is_flag=True, + "excision_extrapolation_mode", + type=click.Choice(ExcisionExtrapolationMode.__members__), + default=ExcisionExtrapolationMode.NoExtrapolation.name, + show_default=True, + callback=lambda ctx, param, value: ExcisionExtrapolationMode.__members__[ + value + ], help=( - "Enables extrapolation into excision regions of the domain. " - "This can be useful to fill the excision region with " - "(constraint-violating but smooth) data so it can be imported into " - "moving puncture codes." + "Enables extrapolation into excision regions of the domain. This can be" + " useful to fill the excision region with (constraint-violating but" + " smooth) data so it can be imported into moving puncture codes (set to" + " 'AnchorPoints'), or to extrapolate a little bit into the excision to" + " help the evolution find and track the apparent horizon (set to" + " 'NearestElement')." ), ) @click.option( diff --git a/src/PointwiseFunctions/InitialDataUtilities/NumericData.cpp b/src/PointwiseFunctions/InitialDataUtilities/NumericData.cpp index 96fc3eeed230..5c29b1de0651 100644 --- a/src/PointwiseFunctions/InitialDataUtilities/NumericData.cpp +++ b/src/PointwiseFunctions/InitialDataUtilities/NumericData.cpp @@ -9,18 +9,19 @@ #include #include -NumericData::NumericData(std::string file_glob, std::string subgroup, - int observation_step, bool extrapolate_into_excisions) +NumericData::NumericData( + std::string file_glob, std::string subgroup, int observation_step, + spectre::Exporter::ExcisionExtrapolationMode excision_extrapolation_mode) : file_glob_(std::move(file_glob)), subgroup_(std::move(subgroup)), observation_step_(observation_step), - extrapolate_into_excisions_(extrapolate_into_excisions) {} + excision_extrapolation_mode_(excision_extrapolation_mode) {} bool operator==(const NumericData& lhs, const NumericData& rhs) { return lhs.file_glob() == rhs.file_glob() and lhs.subgroup() == rhs.subgroup() and lhs.observation_step() == rhs.observation_step() and - lhs.extrapolate_into_excisions() == rhs.extrapolate_into_excisions(); + lhs.excision_extrapolation_mode() == rhs.excision_extrapolation_mode(); } bool operator!=(const NumericData& lhs, const NumericData& rhs) { @@ -31,7 +32,7 @@ void NumericData::pup(PUP::er& p) { p | file_glob_; p | subgroup_; p | observation_step_; - p | extrapolate_into_excisions_; + p | excision_extrapolation_mode_; } namespace elliptic::analytic_data { diff --git a/src/PointwiseFunctions/InitialDataUtilities/NumericData.hpp b/src/PointwiseFunctions/InitialDataUtilities/NumericData.hpp index 3bed3389a83d..c672a0dcf98a 100644 --- a/src/PointwiseFunctions/InitialDataUtilities/NumericData.hpp +++ b/src/PointwiseFunctions/InitialDataUtilities/NumericData.hpp @@ -84,7 +84,7 @@ class NumericData { struct ExtrapolateIntoExcisions { static constexpr Options::String help = { "Whether to extrapolate data into excised regions"}; - using type = bool; + using type = spectre::Exporter::ExcisionExtrapolationMode; }; using options = tmpl::list; @@ -98,8 +98,10 @@ class NumericData { NumericData& operator=(NumericData&&) = default; ~NumericData() = default; - NumericData(std::string file_glob, std::string subgroup, int observation_step, - bool extrapolate_into_excisions); + NumericData( + std::string file_glob, std::string subgroup, int observation_step, + spectre::Exporter::ExcisionExtrapolationMode excision_extrapolation_mode = + spectre::Exporter::ExcisionExtrapolationMode::NoExtrapolation); const std::string& file_glob() const { return file_glob_; } @@ -107,8 +109,9 @@ class NumericData { int observation_step() const { return observation_step_; } - bool extrapolate_into_excisions() const { - return extrapolate_into_excisions_; + spectre::Exporter::ExcisionExtrapolationMode excision_extrapolation_mode() + const { + return excision_extrapolation_mode_; } template @@ -119,7 +122,7 @@ class NumericData { tmpl::list>( file_glob_, subgroup_, spectre::Exporter::ObservationStep{observation_step_}, x, - extrapolate_into_excisions_); + excision_extrapolation_mode_); } template @@ -138,7 +141,7 @@ class NumericData { std::string file_glob_{}; std::string subgroup_{}; int observation_step_{}; - bool extrapolate_into_excisions_{}; + spectre::Exporter::ExcisionExtrapolationMode excision_extrapolation_mode_{}; }; bool operator==(const NumericData& lhs, const NumericData& rhs); diff --git a/src/Visualization/Python/PlotAlongLine.py b/src/Visualization/Python/PlotAlongLine.py index 6e824cc7b3d5..c98ad39c51e5 100644 --- a/src/Visualization/Python/PlotAlongLine.py +++ b/src/Visualization/Python/PlotAlongLine.py @@ -11,7 +11,7 @@ import matplotlib.pyplot as plt import numpy as np -from spectre.IO.Exporter import interpolate_to_points +from spectre.IO.Exporter import ExcisionExtrapolationMode, interpolate_to_points from spectre.Visualization.OpenVolfiles import ( open_volfiles, open_volfiles_command, @@ -54,7 +54,7 @@ def plot_along_line( vars, line_start, line_end, - extrapolate_into_excisions=False, + excision_extrapolation_mode=ExcisionExtrapolationMode.NoExtrapolation, num_samples=200, num_threads=None, x_logscale=False, @@ -148,7 +148,7 @@ def update_plot(obs_id, obs_time): observation_id=obs_id, tensor_components=vars, target_points=target_coords, - extrapolate_into_excisions=extrapolate_into_excisions, + excision_extrapolation_mode=excision_extrapolation_mode, num_threads=num_threads, ) for y, var_name, line in zip(vars_on_line, vars, lines): @@ -197,12 +197,20 @@ def update_plot(obs_id, obs_time): ) @click.option( "--extrapolate-into-excisions", - is_flag=True, + "excision_extrapolation_mode", + type=click.Choice(ExcisionExtrapolationMode.__members__), + default=ExcisionExtrapolationMode.NoExtrapolation.name, + show_default=True, + callback=lambda ctx, param, value: getattr( + ExcisionExtrapolationMode, value + ), help=( - "Enables extrapolation into excision regions of the domain. " - "This can be useful to fill the excision region with " - "(constraint-violating but smooth) data so it can be imported into " - "moving puncture codes." + "Enables extrapolation into excision regions of the domain. This can be" + " useful to fill the excision region with (constraint-violating but" + " smooth) data so it can be imported into moving puncture codes (set to" + " 'AnchorPoints'), or to extrapolate a little bit into the excision to" + " help the evolution find and track the apparent horizon (set to" + " 'NearestElement')." ), ) @click.option( diff --git a/src/Visualization/Python/PlotSlice.py b/src/Visualization/Python/PlotSlice.py index 9394d6f7d9c1..564e986b8831 100644 --- a/src/Visualization/Python/PlotSlice.py +++ b/src/Visualization/Python/PlotSlice.py @@ -14,7 +14,7 @@ import rich import spectre.IO.H5 as spectre_h5 -from spectre.IO.Exporter import interpolate_to_points +from spectre.IO.Exporter import ExcisionExtrapolationMode, interpolate_to_points from spectre.Visualization.OpenVolfiles import ( open_volfiles, open_volfiles_command, @@ -85,7 +85,7 @@ def plot_slice( slice_extent, slice_normal, slice_up, - extrapolate_into_excisions=False, + excision_extrapolation_mode=ExcisionExtrapolationMode.NoExtrapolation, num_samples=[200, 200], num_threads=None, title=None, @@ -201,7 +201,7 @@ def plot_slice(obs_id, time): observation_id=obs_id, tensor_components=[var_name], target_points=target_coords.reshape(3, np.prod(num_samples)), - extrapolate_into_excisions=extrapolate_into_excisions, + excision_extrapolation_mode=excision_extrapolation_mode, num_threads=num_threads, )[0] ).reshape(num_samples) @@ -291,12 +291,20 @@ def update_plot(frame): ) @click.option( "--extrapolate-into-excisions", - is_flag=True, + "excision_extrapolation_mode", + type=click.Choice(ExcisionExtrapolationMode.__members__), + default=ExcisionExtrapolationMode.NoExtrapolation.name, + show_default=True, + callback=lambda ctx, param, value: getattr( + ExcisionExtrapolationMode, value + ), help=( - "Enables extrapolation into excision regions of the domain. " - "This can be useful to fill the excision region with " - "(constraint-violating but smooth) data so it can be imported into " - "moving puncture codes." + "Enables extrapolation into excision regions of the domain. This can be" + " useful to fill the excision region with (constraint-violating but" + " smooth) data so it can be imported into moving puncture codes (set to" + " 'AnchorPoints'), or to extrapolate a little bit into the excision to" + " help the evolution find and track the apparent horizon (set to" + " 'NearestElement')." ), ) @click.option( diff --git a/tests/Unit/Domain/Test_BlockAndElementLogicalCoordinates.cpp b/tests/Unit/Domain/Test_BlockAndElementLogicalCoordinates.cpp index f701015dbeee..52904d013330 100644 --- a/tests/Unit/Domain/Test_BlockAndElementLogicalCoordinates.cpp +++ b/tests/Unit/Domain/Test_BlockAndElementLogicalCoordinates.cpp @@ -33,6 +33,7 @@ #include "Domain/Structure/ElementId.hpp" #include "Domain/Structure/InitialElementIds.hpp" #include "Framework/TestHelpers.hpp" +#include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrHorizon.hpp" #include "Utilities/Gsl.hpp" #include "Utilities/MakeArray.hpp" #include "Utilities/StdHelpers.hpp" @@ -817,6 +818,86 @@ void test_block_logical_coordinates_with_roundoff_error() { CHECK(get<1>(block_logical_coords[5]->data) < 1.0); CHECK(get<0>(block_logical_coords[6]->data) < 1.0); } + +void test_excision_spheres(const bool time_dependent) { + CAPTURE(time_dependent); + const auto shell = domain::creators::Sphere( + 1., 3., domain::creators::Sphere::Excision{}, 0_st, 3_st, true, + std::nullopt, {}, domain::CoordinateMaps::Distribution::Linear, + ShellWedges::All, + time_dependent + ? std::make_optional(domain::creators::Sphere::TimeDepOptionType{ + domain::creators::sphere::TimeDependentMapOptions{ + 0.0, + {{12, + domain::creators::time_dependent_options:: + KerrSchildFromBoyerLindquist{1.0, {{0.0, 0.0, 0.9}}}, + std::nullopt}}, + std::nullopt, + std::nullopt, + std::nullopt, + false}}) + : std::nullopt); + const auto domain = shell.create_domain(); + const auto& excision_sphere = domain.excision_spheres().at("ExcisionSphere"); + const double time = 0.0; + const auto functions_of_time = shell.functions_of_time(); + + std::vector, BlockLogicalCoords<3>>> test_data{}; + // Outside the excision sphere + test_data.push_back({tnsr::I{{2., 0., 0.}}, std::nullopt}); + if (time_dependent) { + // On the boundary of the excision sphere + test_data.push_back( + {tnsr::I{ + {get(gr::Solutions::kerr_schild_radius_from_boyer_lindquist< + double>(1.0, {{M_PI_2, 0.0}}, 1.0, {{0.0, 0.0, 0.9}})), + 0., 0.}}, + {{domain::BlockId{4}, + tnsr::I{{0., 0., -1.}}}}}); + // Inside the excision sphere + test_data.push_back( + {tnsr::I{ + {get(gr::Solutions::kerr_schild_radius_from_boyer_lindquist< + double>(0.5, {{M_PI_2, 0.0}}, 1.0, {{0.0, 0.0, 0.9}})), + 0., 0.}}, + {{domain::BlockId{4}, + // If the shape map preserved the BL-KS transformation within the + // excision, the result should be -1.5. However, apparently it + // doesn't. + tnsr::I{ + {0., 0., -1.31579939060867179}}}}}); + } else { + // On the boundary of the excision sphere + test_data.push_back( + {tnsr::I{{1., 0., 0.}}, + {{domain::BlockId{4}, + tnsr::I{{0., 0., -1.}}}}}); + // Inside the excision sphere + test_data.push_back( + {tnsr::I{{0.5, 0., 0.}}, + {{domain::BlockId{4}, + tnsr::I{{0., 0., -1.5}}}}}); + } + + for (const auto& [input_point, expected] : test_data) { + CAPTURE(input_point); + const auto result = block_logical_coordinates_in_excision( + input_point, excision_sphere, domain.blocks(), time, functions_of_time); + CHECK(result.has_value() == expected.has_value()); + if (result.has_value()) { + CHECK(result->id.get_index() == expected->id.get_index()); + CHECK_ITERABLE_APPROX(result->data, expected->data); + const auto element_logical_coords = element_logical_coordinates( + result->data, ElementId<3>{result->id.get_index()}); + REQUIRE(element_logical_coords.has_value()); + const tnsr::I expected_element_logical{ + {get<0>(expected->data), get<1>(expected->data), + get<2>(expected->data)}}; + CHECK_ITERABLE_APPROX(*element_logical_coords, expected_element_logical); + } + } +} } // namespace SPECTRE_TEST_CASE("Unit.Domain.BlockAndElementLogicalCoords", @@ -834,4 +915,6 @@ SPECTRE_TEST_CASE("Unit.Domain.BlockAndElementLogicalCoords", test_block_logical_coordinates1fail(); test_element_ids_are_uniquely_determined(); test_block_logical_coordinates_with_roundoff_error(); + test_excision_spheres(false); + test_excision_spheres(true); } diff --git a/tests/Unit/IO/Exporter/Test_Exporter.cpp b/tests/Unit/IO/Exporter/Test_Exporter.cpp index d61db6939ccc..d182426e71d2 100644 --- a/tests/Unit/IO/Exporter/Test_Exporter.cpp +++ b/tests/Unit/IO/Exporter/Test_Exporter.cpp @@ -16,9 +16,13 @@ #include "DataStructures/Tensor/Tensor.hpp" #include "Domain/Creators/BinaryCompactObject.hpp" #include "Domain/Creators/Rectilinear.hpp" +#include "Domain/Creators/RegisterDerivedWithCharm.hpp" +#include "Domain/Creators/TimeDependence/RegisterDerivedWithCharm.hpp" #include "Domain/ElementMap.hpp" +#include "Domain/FunctionsOfTime/RegisterDerivedWithCharm.hpp" #include "Domain/Structure/InitialElementIds.hpp" #include "Evolution/Systems/ScalarWave/Tags.hpp" +#include "Framework/TestCreation.hpp" #include "IO/Exporter/Exporter.hpp" #include "IO/Exporter/InterpolateToPoints.hpp" #include "IO/H5/File.hpp" @@ -26,16 +30,38 @@ #include "IO/H5/VolumeData.hpp" #include "Informer/InfoFromBuild.hpp" #include "NumericalAlgorithms/Spectral/LogicalCoordinates.hpp" +#include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrHorizon.hpp" #include "Utilities/FileSystem.hpp" #include "Utilities/Serialization/Serialize.hpp" namespace spectre::Exporter { SPECTRE_TEST_CASE("Unit.IO.Exporter", "[Unit]") { + { + INFO("ExcisionExtrapolationMode enum"); + CHECK(get_output(ExcisionExtrapolationMode::NoExtrapolation) == + "NoExtrapolation"); + CHECK(get_output(ExcisionExtrapolationMode::NearestElement) == + "NearestElement"); + CHECK(get_output(ExcisionExtrapolationMode::RadialAnchors) == + "RadialAnchors"); + for (const auto mode : {ExcisionExtrapolationMode::NoExtrapolation, + ExcisionExtrapolationMode::NearestElement, + ExcisionExtrapolationMode::RadialAnchors}) { + const auto created = + TestHelpers::test_creation( + get_output(mode)); + CHECK(created == mode); + } + } + #ifdef _OPENMP // Disable OpenMP multithreading since multiple unit tests may run in parallel omp_set_num_threads(1); #endif + domain::creators::register_derived_with_charm(); + domain::creators::time_dependence::register_derived_with_charm(); + domain::FunctionsOfTime::register_derived_with_charm(); { INFO("Bundled volume data files"); const auto interpolated_data = interpolate_to_points<3>( @@ -119,6 +145,7 @@ SPECTRE_TEST_CASE("Unit.IO.Exporter", "[Unit]") { { INFO("Extrapolation into BBH excisions"); using Object = domain::creators::BinaryCompactObject::Object; + const std::array spin{{0.0, 0.0, 0.3}}; const domain::creators::BinaryCompactObject domain_creator{ Object{1., 4., 8., true, true}, Object{0.8, 2.5, -6., true, true}, @@ -131,7 +158,17 @@ SPECTRE_TEST_CASE("Unit.IO.Exporter", "[Unit]") { true, domain::CoordinateMaps::Distribution::Projective, domain::CoordinateMaps::Distribution::Inverse, - 120.}; + 120., + domain::creators::bco::TimeDependentMapOptions{ + 0.0, + std::nullopt, + std::nullopt, + std::nullopt, + {{12, + domain::creators::time_dependent_options:: + KerrSchildFromBoyerLindquist{1.0, spin}, + std::nullopt}}, + std::nullopt}}; const auto domain = domain_creator.create_domain(); const auto functions_of_time = domain_creator.functions_of_time(); const double time = 1.0; @@ -171,9 +208,13 @@ SPECTRE_TEST_CASE("Unit.IO.Exporter", "[Unit]") { serialize(functions_of_time)); } // Interpolate - tnsr::I target_points{{{{8.5, 8.6, 8.7, 8.9, 9., 10.}, - {0., 0., 0., 0., 0., 0.}, - {0., 0., 0., 0., 0., 0.}}}}; + const double deformed_radius = + get(gr::Solutions::kerr_schild_radius_from_boyer_lindquist( + 1.0, {{M_PI_2, 0.0}}, 1.0, spin)); + tnsr::I target_points{ + {{{8.5, 8.6, 8.7, 8.9, deformed_radius, 10.}, + {0., 0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 0., 0.}}}}; const size_t num_target_points = get<0>(target_points).size(); std::array, 3> target_points_array{}; for (size_t d = 0; d < 3; ++d) { @@ -182,31 +223,45 @@ SPECTRE_TEST_CASE("Unit.IO.Exporter", "[Unit]") { gsl::at(target_points_array, d)[i] = target_points.get(d)[i]; } } - const auto interpolated_data = interpolate_to_points<3>( - h5_file_name, "/VolumeData", ObservationId{123}, {"Psi"}, - target_points_array, true); - CHECK(interpolated_data.size() == 1); - CHECK(interpolated_data[0].size() == num_target_points); - // Check result - DataVector psi_interpolated(num_target_points); - DataVector psi_expected(num_target_points); - for (size_t i = 0; i < num_target_points; ++i) { - psi_interpolated[i] = interpolated_data[0][i]; - tnsr::I x_target{ - {{get<0>(target_points)[i], get<1>(target_points)[i], - get<2>(target_points)[i]}}}; - psi_expected[i] = func(x_target); - } - // These points are extrapolated and therefore less precise - Approx approx_extrapolated = Approx::custom().epsilon(1.e-1).scale(1.0); - for (size_t i = 0; i < 4; ++i) { - CHECK(psi_interpolated[i] == approx_extrapolated(psi_expected[i])); + for (const auto extrapolation_mode : + {ExcisionExtrapolationMode::NearestElement, + ExcisionExtrapolationMode::RadialAnchors}) { + CAPTURE(extrapolation_mode); + const auto interpolated_data = interpolate_to_points<3>( + h5_file_name, "/VolumeData", ObservationId{123}, {"Psi"}, + target_points_array, extrapolation_mode); + CHECK(interpolated_data.size() == 1); + CHECK(interpolated_data[0].size() == num_target_points); + // Check result + DataVector psi_interpolated(num_target_points); + DataVector psi_expected(num_target_points); + for (size_t i = 0; i < num_target_points; ++i) { + psi_interpolated[i] = interpolated_data[0][i]; + const tnsr::I x_target{ + {{get<0>(target_points)[i], get<1>(target_points)[i], + get<2>(target_points)[i]}}}; + psi_expected[i] = func(x_target); + } + // These points are extrapolated and therefore less precise + const Approx approx_extrapolated = + Approx::custom().epsilon(1.e-1).scale(1.0); + for (size_t i = 0; i < 4; ++i) { + if (i < 2 and + extrapolation_mode == ExcisionExtrapolationMode::NearestElement) { + // These points are too far away from the excision surface to be + // extrapolated precisely enough with the nearest element method + continue; + } + CAPTURE(i); + CHECK(psi_interpolated[i] == approx_extrapolated(psi_expected[i])); + } + // This point is exact + CHECK(psi_interpolated[4] == approx(psi_expected[4])); + // This point is interpolated + const Approx approx_interpolated = + Approx::custom().epsilon(1.e-6).scale(1.0); + CHECK(psi_interpolated[5] == approx_interpolated(psi_expected[5])); } - // This point is exact - CHECK(psi_interpolated[4] == approx(psi_expected[4])); - // This point is interpolated - Approx approx_interpolated = Approx::custom().epsilon(1.e-6).scale(1.0); - CHECK(psi_interpolated[5] == approx_interpolated(psi_expected[5])); // Delete the test file if (file_system::check_if_file_exists(h5_file_name)) { file_system::rm(h5_file_name, true); diff --git a/tests/Unit/PointwiseFunctions/InitialDataUtilities/Test_NumericData.cpp b/tests/Unit/PointwiseFunctions/InitialDataUtilities/Test_NumericData.cpp index ddd651c226c1..2683115b3ece 100644 --- a/tests/Unit/PointwiseFunctions/InitialDataUtilities/Test_NumericData.cpp +++ b/tests/Unit/PointwiseFunctions/InitialDataUtilities/Test_NumericData.cpp @@ -21,8 +21,8 @@ void test_numeric_data() { "ObservationStep: -1\n" "ExtrapolateIntoExcisions: False\n"; const auto created = TestHelpers::test_creation(option_string); - const elliptic::analytic_data::NumericData numeric_data{ - file_name, "element_data", -1, false}; + const elliptic::analytic_data::NumericData numeric_data{file_name, + "element_data", -1}; CHECK(created == numeric_data); test_serialization(numeric_data); test_copy_semantics(numeric_data);