diff --git a/docs/References.bib b/docs/References.bib index 9fb1554bf5b1..21a330d27ed8 100644 --- a/docs/References.bib +++ b/docs/References.bib @@ -843,6 +843,21 @@ @article{Dumbser2017okk year = "2018" } +@article{Etienne2008re, + author = {Etienne, Zachariah B. and Liu, Yuk Tung and Shapiro, Stuart + L. and Baumgarte, Thomas W.}, + title = {General relativistic simulations of black-hole-neutron-star + mergers: Effects of black-hole spin}, + eprint = {0812.2245}, + archiveprefix = {arXiv}, + primaryclass = {astro-ph}, + doi = {10.1103/PhysRevD.79.044024}, + journal = {Phys. Rev. D}, + volume = {79}, + pages = {044024}, + year = {2009} +} + @article{Etienne2010ui, author = "Etienne, Zachariah B. and Liu, Yuk Tung and Shapiro, Stuart L.", diff --git a/src/DataStructures/Tensor/EagerMath/CartesianToSpherical.cpp b/src/DataStructures/Tensor/EagerMath/CartesianToSpherical.cpp index bf10167a5fd8..67e789e9a507 100644 --- a/src/DataStructures/Tensor/EagerMath/CartesianToSpherical.cpp +++ b/src/DataStructures/Tensor/EagerMath/CartesianToSpherical.cpp @@ -85,7 +85,7 @@ tnsr::I spherical_to_cartesian( Frame::Spherical>& x); GENERATE_INSTANTIATIONS(INSTANTIATE, (double, DataVector), (1, 2, 3), - (Frame::Grid, Frame::Inertial)) + (Frame::Grid, Frame::Distorted, Frame::Inertial)) #undef DTYPE #undef DIM diff --git a/src/IO/Exporter/Exporter.cpp b/src/IO/Exporter/Exporter.cpp index 9a1389b5793b..794f670312b4 100644 --- a/src/IO/Exporter/Exporter.cpp +++ b/src/IO/Exporter/Exporter.cpp @@ -8,6 +8,7 @@ #include #endif // _OPENMP +#include "DataStructures/Tensor/EagerMath/CartesianToSpherical.hpp" #include "Domain/BlockLogicalCoordinates.hpp" #include "Domain/Creators/RegisterDerivedWithCharm.hpp" #include "Domain/Creators/TimeDependence/RegisterDerivedWithCharm.hpp" @@ -19,6 +20,7 @@ #include "IO/H5/TensorData.hpp" #include "IO/H5/VolumeData.hpp" #include "NumericalAlgorithms/Interpolation/IrregularInterpolant.hpp" +#include "NumericalAlgorithms/Interpolation/PolynomialInterpolation.hpp" #include "Utilities/FileSystem.hpp" #include "Utilities/GenerateInstantiations.hpp" #include "Utilities/GetOutput.hpp" @@ -123,7 +125,176 @@ void interpolate_to_points( } } } // omp for - } // omp parallel + } // omp parallel +} + +// Data structure for extrapolation of tensor components into excisions +template +struct ExtrapolationInfo { + // Index of the target point in the result array (y target) + size_t target_index; + // Index of the first anchor point in the source array (y data) + size_t source_index; + // Coordinate of the target point (x target) + double target_point; + // Coordinates of the anchor points (x data) + std::array anchors; +}; + +// Construct anchor points for extrapolation into excisions. +// +// Anchor points are placed radially in the grid frame around the excision +// sphere, which is spherical in the grid frame. Then the anchor points are +// added to the `block_logical_coords` so data is interpolated to them. Also an +// entry is added to `extrapolation_info` so the interpolated data on the anchor +// points can be extrapolated to the target point in the excision. +// +// This function returns `true` if the target point is inside the excision and +// the anchor points were added. Otherwise, returns `false` and the next +// excision should be tried. +// +// Alternatives and notes: +// - We could extrapolate directly from the nearest element, using the Lagrange +// polynomial basis that the element already has. However, this is unstable +// when the element is too small (possibly h-refined) and/or has log spacing. +// 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. +// - We could construct anchor points in different frames, e.g. in the inertial +// frame. These choices probably don't make a big difference. +// - We could do more performance optimizations for the case where the excision +// sphere is spherical. E.g. in that case the radial anchor points for all +// points are the same. +template +bool add_extrapolation_anchors( + const gsl::not_null>*> + block_logical_coords, + const gsl::not_null< + std::vector>*> + extrapolation_info, + const ExcisionSphere& excision_sphere, const Domain& domain, + const tnsr::I& target_point, + const double time, const domain::FunctionsOfTimeMap& functions_of_time, + const double extrapolation_spacing) { + // Get spherical coordinates around the excision in distorted frame. + // Note that the excision sphere doesn't have a distorted map, but its + // grid-to-inertial map is just the neighboring block's + // distorted-to-inertial map, so we can use either of those. + auto x_distorted = + [&excision_sphere, &target_point, &time, + &functions_of_time]() -> tnsr::I { + if (excision_sphere.is_time_dependent()) { + const auto x_grid = + excision_sphere.moving_mesh_grid_to_inertial_map().inverse( + target_point, time, functions_of_time); + ASSERT(x_grid.has_value(), + "Failed to invert grid-to-inertial map of excision sphere at " + << excision_sphere.center() << " for point " << target_point + << "."); + tnsr::I result{}; + for (size_t d = 0; d < Dim; ++d) { + result.get(d) = x_grid->get(d); + } + return result; + } else { + tnsr::I result{}; + for (size_t d = 0; d < Dim; ++d) { + result.get(d) = target_point.get(d); + } + return result; + } + }(); + for (size_t d = 0; d < Dim; ++d) { + x_distorted.get(d) -= excision_sphere.center().get(d); + } + auto x_spherical_distorted = cartesian_to_spherical(x_distorted); + // Construct anchor points in grid frame + tnsr::I> + anchors_grid_spherical{NumExtrapolationAnchors}; + for (size_t i = 0; i < NumExtrapolationAnchors; ++i) { + get<0>(anchors_grid_spherical)[i] = + 1. + static_cast(i) * extrapolation_spacing; + } + get<0>(anchors_grid_spherical) *= excision_sphere.radius(); + // The grid-distorted map preserves angles + if constexpr (Dim > 1) { + get<1>(anchors_grid_spherical) = get<1>(x_spherical_distorted); + } + if constexpr (Dim > 2) { + get<2>(anchors_grid_spherical) = get<2>(x_spherical_distorted); + } + auto anchors_grid = spherical_to_cartesian(anchors_grid_spherical); + for (size_t d = 0; d < Dim; ++d) { + anchors_grid.get(d) += excision_sphere.center().get(d); + } + // Map anchor points to block logical coordinates. These are needed for + // interpolation. + auto anchors_block_logical = + block_logical_coordinates(domain, anchors_grid, time, functions_of_time); + const size_t block_id = anchors_block_logical[0]->id.get_index(); + const auto& block = domain.blocks()[block_id]; + // Map anchor points to the distorted frame. We will extrapolate in + // the distorted frame because we have the target point in the distorted + // frame. The target point in the grid frame is undefined because there's + // no grid-distorted map in the excision sphere. We could extrapolate in + // the inertial frame, but that's just an unnecessary transformation. + auto anchors_distorted = + [&block, &anchors_grid, &time, + &functions_of_time]() -> tnsr::I { + if (block.has_distorted_frame()) { + return block.moving_mesh_grid_to_distorted_map()(anchors_grid, time, + functions_of_time); + } else { + tnsr::I result{}; + for (size_t d = 0; d < Dim; ++d) { + result.get(d) = anchors_grid.get(d); + } + return result; + } + }(); + for (size_t d = 0; d < Dim; ++d) { + anchors_distorted.get(d) -= excision_sphere.center().get(d); + } + // Compute magnitude(anchors_distorted) and store in std::array + std::array radial_anchors_distorted_frame{}; + for (size_t i = 0; i < NumExtrapolationAnchors; ++i) { + auto& radial_anchor_distorted_frame = + gsl::at(radial_anchors_distorted_frame, i); + radial_anchor_distorted_frame = square(get<0>(anchors_distorted)[i]); + if constexpr (Dim > 1) { + radial_anchor_distorted_frame += square(get<1>(anchors_distorted)[i]); + } + if constexpr (Dim > 2) { + radial_anchor_distorted_frame += square(get<2>(anchors_distorted)[i]); + } + radial_anchor_distorted_frame = sqrt(radial_anchor_distorted_frame); + } + // Return false if the target point is not inside the excision. It would be + // nice to do this earlier to avoid unnecessary work. Here we use the first + // anchor point in the distorted frame, which is at the excision boundary in + // the angular direction of the target point. + const double excision_radius_distorted = radial_anchors_distorted_frame[0]; + if (get<0>(x_spherical_distorted) > excision_radius_distorted) { + return false; + } + // Add the anchor points to the the interpolation target points + for (size_t i = 0; i < NumExtrapolationAnchors; ++i) { + ASSERT(anchors_block_logical[i].has_value(), + "Extrapolation anchor point is not in any block. This should " + "not happen."); + block_logical_coords->push_back(std::move(anchors_block_logical[i])); + } + // Add the anchor points to the extrapolation info + extrapolation_info->push_back( + {// Target index is the point we're extrapolating to. It will be set + // outside this function. + std::numeric_limits::max(), + // Source index will be adjusted outside this function as well. + extrapolation_info->size() * NumExtrapolationAnchors, + // Target point and anchors are the radii in the distorted frame + get<0>(x_spherical_distorted), + std::move(radial_anchors_distorted_frame)}); + return true; } // Determines the selected observation ID in the volume data file, given either @@ -158,6 +329,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 std::optional num_threads) { domain::creators::register_derived_with_charm(); domain::creators::time_dependence::register_derived_with_charm(); @@ -235,10 +407,21 @@ std::vector> interpolate_to_points( // through the domain. This is the most expensive part of the function, so we // parallelize the loop. std::vector> block_logical_coords(num_target_points); + // We also set up the extrapolation into excisions here. Anchor points are + // added to the `block_logical_coords` and additional information is collected + // in `extrapolation_info` for later extrapolation. + constexpr size_t num_extrapolation_anchors = 8; + const double extrapolation_spacing = 0.3; + std::vector> + extrapolation_info{}; #pragma omp parallel num_threads(resolved_num_threads) { + // Set up thread-local variables tnsr::I target_point{}; -#pragma omp for + std::vector> extra_block_logical_coords{}; + std::vector> + extra_extrapolation_info{}; +#pragma omp for nowait for (size_t s = 0; s < num_target_points; ++s) { for (size_t d = 0; d < Dim; ++d) { target_point.get(d) = gsl::at(target_points, d)[s]; @@ -252,17 +435,49 @@ std::vector> interpolate_to_points( break; } } // for blocks - } // omp for target points - } // omp parallel + if (block_logical_coords[s].has_value() or + not extrapolate_into_excisions) { + continue; + } + // The point wasn't found in any block. Check if it's in an excision and + // set up extrapolation 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; + } + } + } // omp for target points +#pragma omp critical + { + // Append the extra block logical coordinates and extrapolation info from + // this thread to the global vectors. Also set the source index to the + // index in `block_logical_coords` where we're going to insert the new + // coordinates. + for (auto& info : extra_extrapolation_info) { + info.source_index += block_logical_coords.size(); + } + block_logical_coords.insert(block_logical_coords.end(), + extra_block_logical_coords.begin(), + extra_block_logical_coords.end()); + extrapolation_info.insert(extrapolation_info.end(), + extra_extrapolation_info.begin(), + extra_extrapolation_info.end()); + } // omp critical + } // omp parallel // Allocate memory for result std::vector> result{}; result.reserve(tensor_components.size()); for (size_t i = 0; i < tensor_components.size(); ++i) { - result.emplace_back(num_target_points, + result.emplace_back(block_logical_coords.size(), std::numeric_limits::signaling_NaN()); } - std::vector filled_data(num_target_points, false); + std::vector filled_data(block_logical_coords.size(), false); // Process all volume files in serial, because loading data with H5 must be // done in serial anyway. Instead, the loop over elements within each file is @@ -277,6 +492,29 @@ std::vector> interpolate_to_points( break; } } + + if (extrapolate_into_excisions) { + // Extrapolate into excisions from the anchor points +#pragma omp parallel for num_threads(resolved_num_threads) + for (const auto& extrapolation : extrapolation_info) { + double extrapolation_error = 0.; + for (size_t i = 0; i < tensor_components.size(); ++i) { + intrp::polynomial_interpolation( + make_not_null(&result[i][extrapolation.target_index]), + make_not_null(&extrapolation_error), extrapolation.target_point, + // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-pointer-arithmetic) + gsl::make_span(result[i].data() + extrapolation.source_index, + num_extrapolation_anchors), + gsl::make_span(extrapolation.anchors.data(), + num_extrapolation_anchors)); + } + } + // Clear the anchor points from the result + for (size_t i = 0; i < tensor_components.size(); ++i) { + result[i].resize(num_target_points); + } + } + return result; } @@ -292,6 +530,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, \ 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 55778a856cf7..fdc784182bd8 100644 --- a/src/IO/Exporter/Exporter.hpp +++ b/src/IO/Exporter/Exporter.hpp @@ -45,6 +45,15 @@ 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 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 @@ -62,6 +71,7 @@ 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, std::optional num_threads = std::nullopt); } // namespace spectre::Exporter diff --git a/src/IO/Exporter/Python/Bindings.cpp b/src/IO/Exporter/Python/Bindings.cpp index 582eb2a7c13f..166f8d41cd9c 100644 --- a/src/IO/Exporter/Python/Bindings.cpp +++ b/src/IO/Exporter/Python/Bindings.cpp @@ -21,6 +21,7 @@ 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, const std::optional& num_threads) { const size_t dim = target_points.size(); const spectre::Exporter::ObservationId obs_id{observation_id}; @@ -28,17 +29,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)), - num_threads); + extrapolate_into_excisions, 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)), - num_threads); + extrapolate_into_excisions, 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)), - num_threads); + extrapolate_into_excisions, num_threads); } else { ERROR("Invalid dimension of target points: " << dim @@ -49,5 +50,6 @@ 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("num_threads") = std::nullopt); + py::arg("target_points"), py::arg("extrapolate_into_excisions") = false, + py::arg("num_threads") = std::nullopt); } diff --git a/src/IO/Exporter/Python/InterpolateToPoints.py b/src/IO/Exporter/Python/InterpolateToPoints.py index 74b94c4c2a96..bb2cd9768abf 100644 --- a/src/IO/Exporter/Python/InterpolateToPoints.py +++ b/src/IO/Exporter/Python/InterpolateToPoints.py @@ -41,6 +41,16 @@ "to a couple of target points." ), ) +@click.option( + "--extrapolate-into-excisions", + is_flag=True, + 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." + ), +) @click.option( "--output", "-o", type=click.Path(writable=True), help="Output text file" ) @@ -74,7 +84,7 @@ def interpolate_to_points_command( target_coords_file, output, delimiter, - num_threads, + **kwargs, ): """Interpolate volume data to target coordinates.""" # Load target coords from file @@ -97,7 +107,7 @@ def interpolate_to_points_command( observation_id=obs_id, tensor_components=vars, target_points=target_coords.T, - num_threads=num_threads, + **kwargs, ) ) diff --git a/src/Visualization/Python/PlotAlongLine.py b/src/Visualization/Python/PlotAlongLine.py index cea3be5bc6de..35aca3a8c5d1 100644 --- a/src/Visualization/Python/PlotAlongLine.py +++ b/src/Visualization/Python/PlotAlongLine.py @@ -69,6 +69,16 @@ def points_on_line( "Specify as comma-separated list, e.g. '1,0,0'. [required]" ), ) +@click.option( + "--extrapolate-into-excisions", + is_flag=True, + 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." + ), +) @click.option( "--num-samples", "-N", @@ -117,6 +127,7 @@ def plot_along_line_command( vars, line_start, line_end, + extrapolate_into_excisions, num_samples, num_threads, x_logscale, @@ -210,6 +221,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, num_threads=num_threads, ) for y, var_name, line in zip(vars_on_line, vars, lines): diff --git a/tests/Unit/IO/Exporter/Test_Exporter.cpp b/tests/Unit/IO/Exporter/Test_Exporter.cpp index 98b2a8dbd414..a61c1b9c556a 100644 --- a/tests/Unit/IO/Exporter/Test_Exporter.cpp +++ b/tests/Unit/IO/Exporter/Test_Exporter.cpp @@ -12,7 +12,12 @@ #endif // _OPENMP #include "DataStructures/DataVector.hpp" +#include "DataStructures/Tensor/EagerMath/Magnitude.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "Domain/Creators/BinaryCompactObject.hpp" #include "Domain/Creators/Rectangle.hpp" +#include "Domain/ElementMap.hpp" +#include "Domain/Structure/InitialElementIds.hpp" #include "IO/Exporter/Exporter.hpp" #include "IO/H5/File.hpp" #include "IO/H5/TensorData.hpp" @@ -89,6 +94,100 @@ SPECTRE_TEST_CASE("Unit.IO.Exporter", "[Unit]") { file_system::rm(h5_file_name, true); } } + { + INFO("Extrapolation into BBH excisions"); + using Object = domain::creators::BinaryCompactObject::Object; + const domain::creators::BinaryCompactObject domain_creator{ + Object{1., 4., 8., true, true}, + Object{0.8, 2.5, -6., true, true}, + 60., + 300., + 0_st, + 6_st, + true, + domain::CoordinateMaps::Distribution::Projective, + domain::CoordinateMaps::Distribution::Inverse, + 120.}; + const auto domain = domain_creator.create_domain(); + const auto functions_of_time = domain_creator.functions_of_time(); + const double time = 1.0; + const auto element_ids = + initial_element_ids(domain_creator.initial_refinement_levels()); + const Mesh<3> mesh{18, Spectral::Basis::Legendre, + Spectral::Quadrature::GaussLobatto}; + const auto xi = logical_coordinates(mesh); + std::vector element_volume_data{}; + const auto func = [](const auto& x) { + auto x1 = x; + get<0>(x1) -= 8; + auto x2 = x; + get<0>(x2) += 6; + const auto r1 = get(magnitude(x1)); + const auto r2 = get(magnitude(x2)); + return blaze::evaluate(exp(-square(r1)) + exp(-square(r2))); + }; + for (const auto& element_id : element_ids) { + ElementMap<3, Frame::Inertial> element_map{ + element_id, domain.blocks()[element_id.block_id()]}; + const auto x = element_map(xi, time, functions_of_time); + DataVector psi = func(x); + element_volume_data.push_back(ElementVolumeData{ + element_id, {TensorComponent{"Psi", std::move(psi)}}, mesh}); + } + // Write to file + const std::string h5_file_name{"Unit.IO.Exporter.Excision.h5"}; + if (file_system::check_if_file_exists(h5_file_name)) { + file_system::rm(h5_file_name, true); + } + { // scope to open and close file + h5::H5File h5_file(h5_file_name); + auto& volfile = h5_file.insert("/VolumeData", 0); + volfile.write_volume_data(123, time, element_volume_data, + serialize(domain), + 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 size_t num_target_points = get<0>(target_points).size(); + std::array, 3> target_points_array{}; + for (size_t d = 0; d < 3; ++d) { + gsl::at(target_points_array, d).resize(num_target_points); + for (size_t i = 0; i < num_target_points; ++i) { + 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])); + } + // 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); + } + } } } // namespace spectre::Exporter