From 043da72cb85220a486db30c18414674c0b36a9a9 Mon Sep 17 00:00:00 2001 From: Nils Vu Date: Tue, 19 Nov 2024 14:19:01 -0800 Subject: [PATCH 1/3] Fix a bug in SphereTransition --- .../ShapeMapTransitionFunctions/SphereTransition.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) 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; From 4b0e979ff6122b4cb2d41d45c382817a6a464fd5 Mon Sep 17 00:00:00 2001 From: Nils Vu Date: Tue, 19 Nov 2024 14:31:03 -0800 Subject: [PATCH 2/3] Add block_logical_coordinates_in_excision --- src/Domain/BlockLogicalCoordinates.cpp | 57 ++++++++++++- src/Domain/BlockLogicalCoordinates.hpp | 26 +++++- src/Domain/ElementLogicalCoordinates.cpp | 53 ++++++------ src/Domain/ElementLogicalCoordinates.hpp | 5 ++ ...Test_BlockAndElementLogicalCoordinates.cpp | 83 +++++++++++++++++++ 5 files changed, 191 insertions(+), 33 deletions(-) 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/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/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); } From 482207fa777138d4962a8d97e7ae070f6f55fd14 Mon Sep 17 00:00:00 2001 From: Nils Vu Date: Tue, 19 Nov 2024 15:09:30 -0800 Subject: [PATCH 3/3] Exporter: add mode to extrapolate from nearest element Uses the Lagrange basis in the nearest element to extrapolate into the excision region. This can be unstable when we want to extrapolate far into the excision because we're extrapolating over many logical element sizes. However, it 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. --- src/IO/Exporter/CMakeLists.txt | 1 + src/IO/Exporter/Exporter.cpp | 99 ++++++++++++++-- src/IO/Exporter/Exporter.hpp | 80 +++++++++++-- src/IO/Exporter/InterpolateToPoints.hpp | 5 +- src/IO/Exporter/Python/Bindings.cpp | 21 +++- src/IO/Exporter/Python/InterpolateToPoints.py | 20 +++- .../InitialDataUtilities/NumericData.cpp | 11 +- .../InitialDataUtilities/NumericData.hpp | 17 +-- src/Visualization/Python/PlotAlongLine.py | 24 ++-- src/Visualization/Python/PlotSlice.py | 24 ++-- tests/Unit/IO/Exporter/Test_Exporter.cpp | 111 +++++++++++++----- .../InitialDataUtilities/Test_NumericData.cpp | 4 +- 12 files changed, 323 insertions(+), 94 deletions(-) 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/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);