Skip to content

Commit

Permalink
fixup
Browse files Browse the repository at this point in the history
  • Loading branch information
nilsvu committed Jul 11, 2024
1 parent 4608e30 commit 3369b36
Showing 1 changed file with 22 additions and 18 deletions.
40 changes: 22 additions & 18 deletions src/IO/Exporter/Exporter.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -129,7 +129,7 @@ void interpolate_to_points(
}

// Data structure for extrapolation of tensor components into excisions
template <size_t ExtrapolationOrder>
template <size_t NumExtrapolationAnchors>
struct ExtrapolationInfo {
// Index of the target point in the result array (y target)
size_t target_index;
Expand All @@ -138,14 +138,15 @@ struct ExtrapolationInfo {
// Coordinate of the target point (x target)
double target_point;
// Coordinates of the anchor points (x data)
std::array<double, ExtrapolationOrder> anchors;
std::array<double, NumExtrapolationAnchors> anchors;
};

template <size_t Dim, size_t ExtrapolationOrder>
template <size_t Dim, size_t NumExtrapolationAnchors>
void add_extrapolation_anchors(
const gsl::not_null<std::vector<BlockLogicalCoords<Dim>>*>
block_logical_coords,
const gsl::not_null<std::vector<ExtrapolationInfo<ExtrapolationOrder>>*>
const gsl::not_null<
std::vector<ExtrapolationInfo<NumExtrapolationAnchors>>*>
extrapolation_info,
const ExcisionSphere<Dim>& excision_sphere, const Domain<Dim>& domain,
const std::array<std::vector<double>, Dim>& target_points,
Expand Down Expand Up @@ -173,16 +174,16 @@ void add_extrapolation_anchors(
// for all points are the same (`radial_anchors_grid_frame`) and we can
// avoid copying them for each point.
const double excision_radius = excision_sphere.radius();
std::array<double, ExtrapolationOrder> radial_anchors_grid_frame{};
for (size_t i = 0; i < ExtrapolationOrder; ++i) {
std::array<double, NumExtrapolationAnchors> radial_anchors_grid_frame{};
for (size_t i = 0; i < NumExtrapolationAnchors; ++i) {
gsl::at(radial_anchors_grid_frame, i) =
excision_radius +
static_cast<double>(i) * extrapolation_spacing * excision_radius;
}
#pragma omp parallel num_threads(num_threads)
{
std::vector<BlockLogicalCoords<Dim>> extra_block_logical_coords{};
std::vector<ExtrapolationInfo<ExtrapolationOrder>>
std::vector<ExtrapolationInfo<NumExtrapolationAnchors>>
extra_extrapolation_info{};
#pragma omp for nowait
for (size_t s = 0; s < block_logical_coords->size(); ++s) {
Expand Down Expand Up @@ -253,8 +254,8 @@ void add_extrapolation_anchors(
}
auto x_spherical_distorted = cartesian_to_spherical(*x_distorted);
tnsr::I<DataVector, Dim, Frame::Spherical<Frame::Grid>>
anchors_grid_spherical{ExtrapolationOrder};
for (size_t i = 0; i < ExtrapolationOrder; ++i) {
anchors_grid_spherical{NumExtrapolationAnchors};
for (size_t i = 0; i < NumExtrapolationAnchors; ++i) {
get<0>(anchors_grid_spherical)[i] =
gsl::at(radial_anchors_grid_frame, i);
}
Expand All @@ -272,7 +273,7 @@ void add_extrapolation_anchors(
// Add the anchor points to the the interpolation target points
auto anchors_block_logical = block_logical_coordinates(
domain, anchors_grid, time, functions_of_time);
for (size_t i = 0; i < ExtrapolationOrder; ++i) {
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.");
Expand Down Expand Up @@ -301,8 +302,9 @@ void add_extrapolation_anchors(
anchors_distorted.get(d) -= excision_sphere.center().get(d);
}
// Compute magnitude(anchors_distorted) and store in std::array
std::array<double, ExtrapolationOrder> radial_anchors_distorted_frame{};
for (size_t i = 0; i < ExtrapolationOrder; ++i) {
std::array<double, NumExtrapolationAnchors>
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]);
Expand All @@ -318,7 +320,7 @@ void add_extrapolation_anchors(
{// Target index is the point we're extrapolating to
s,
// Source index will be adjusted below
extra_extrapolation_info.size() * ExtrapolationOrder,
extra_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)});
Expand Down Expand Up @@ -488,9 +490,10 @@ std::vector<std::vector<double>> interpolate_to_points(
} // omp parallel

// Add anchor points for extrapolation into excisions
constexpr size_t extrapolation_order = 8; // Number of anchor points
constexpr size_t num_extrapolation_anchors = 8;
const double extrapolation_spacing = 0.3;
std::vector<ExtrapolationInfo<extrapolation_order>> extrapolation_info{};
std::vector<ExtrapolationInfo<num_extrapolation_anchors>>
extrapolation_info{};
if (extrapolate_into_excisions) {
for (const auto& [name, excision_sphere] : domain.excision_spheres()) {
add_extrapolation_anchors(make_not_null(&block_logical_coords),
Expand Down Expand Up @@ -530,13 +533,14 @@ std::vector<std::vector<double>> interpolate_to_points(
for (const auto& extrapolation : extrapolation_info) {
double extrapolation_error = 0.;
for (size_t i = 0; i < tensor_components.size(); ++i) {
intrp::polynomial_interpolation<extrapolation_order - 1>(
intrp::polynomial_interpolation<num_extrapolation_anchors - 1>(
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,
extrapolation_order),
gsl::make_span(extrapolation.anchors.data(), extrapolation_order));
num_extrapolation_anchors),
gsl::make_span(extrapolation.anchors.data(),
num_extrapolation_anchors));
}
}
// Clear the anchor points from the result
Expand Down

0 comments on commit 3369b36

Please sign in to comment.