Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add Skew map #6446

Open
wants to merge 4 commits into
base: develop
Choose a base branch
from
+904 −0
Open

Add Skew map #6446

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
367ae71
Add Skew map
knelli2 Jan 17, 2025
e6725b2
fixup. Skew
knelli2 Jan 23, 2025
1ee436b
fixup. Compressed image
knelli2 Jan 23, 2025
15a6529
fixup. Skew
knelli2 Jan 23, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
59 changes: 30 additions & 29 deletions src/Domain/CoordinateMaps/TimeDependent/Skew.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,9 +44,7 @@ Skew::Skew(std::string function_of_time_name,
template <typename T>
std::array<tt::remove_cvref_wrap_t<T>, 3> Skew::operator()(
const std::array<T, 3>& source_coords, const double time,
const std::unordered_map<
std::string, std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
functions_of_time) const {
const domain::FunctionsOfTimeMap& functions_of_time) const {
return map_and_velocity_helper(source_coords, time, functions_of_time, false);
}

Expand All @@ -68,9 +66,12 @@ std::optional<std::array<double, 3>> Skew::inverse(

// Short ciruit if we're at the outer radius
const double target_radius = magnitude(target_coords);
if (target_radius >= outer_radius_ and
equal_within_roundoff(target_radius, outer_radius_)) {
return target_coords;
if (target_radius >= outer_radius_) {
if (equal_within_roundoff(target_radius, outer_radius_)) {
return target_coords;
} else {
return std::nullopt;
}
}

// Another short circuit
Expand Down Expand Up @@ -106,12 +107,14 @@ std::optional<std::array<double, 3>> Skew::inverse(
const double lower_func = root_func(lower_x);
const double upper_func = root_func(upper_x);

if (equal_within_roundoff(lower_func, 0.0)) {
// Since we adjusted the bounds by eps above, we have a slightly larger eps
// here to account for the arithmetic of the root func
if (equal_within_roundoff(lower_func, 0.0, 5 * eps)) {
temporary_source_coord[0] = lower_x;
return temporary_source_coord;
}

if (equal_within_roundoff(upper_func, 0.0)) {
if (equal_within_roundoff(upper_func, 0.0, 5 * eps)) {
temporary_source_coord[0] = upper_x;
return temporary_source_coord;
}
Expand All @@ -132,18 +135,14 @@ std::optional<std::array<double, 3>> Skew::inverse(
template <typename T>
std::array<tt::remove_cvref_wrap_t<T>, 3> Skew::frame_velocity(
const std::array<T, 3>& source_coords, const double time,
const std::unordered_map<
std::string, std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
functions_of_time) const {
const domain::FunctionsOfTimeMap& functions_of_time) const {
return map_and_velocity_helper(source_coords, time, functions_of_time, true);
}

template <typename T>
tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> Skew::jacobian(
const std::array<T, 3>& source_coords, const double time,
const std::unordered_map<
std::string, std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
functions_of_time) const {
const domain::FunctionsOfTimeMap& functions_of_time) const {
using ResultT = tt::remove_cvref_wrap_t<T>;

const ResultT width = get_width(source_coords);
Expand Down Expand Up @@ -171,9 +170,7 @@ tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> Skew::jacobian(
template <typename T>
tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> Skew::inv_jacobian(
const std::array<T, 3>& source_coords, const double time,
const std::unordered_map<
std::string, std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
functions_of_time) const {
const domain::FunctionsOfTimeMap& functions_of_time) const {
return determinant_and_inverse(
jacobian(source_coords, time, functions_of_time))
.second;
Expand Down Expand Up @@ -232,18 +229,22 @@ template <typename T>
std::array<tt::remove_cvref_wrap_t<T>, 3> Skew::get_width_deriv(
const std::array<T, 3>& source_coords_cvref) const {
using ResultT = tt::remove_cvref_wrap_t<T>;
std::array<ResultT, 3> source_coords{};
for (size_t i = 0; i < 3; i++) {
if constexpr (std::is_same_v<ResultT, DataVector>) {
// NOLINTBEGIN
gsl::at(source_coords, i)
.set_data_ref(make_not_null(&const_cast<DataVector&>(
dereference_wrapper(gsl::at(source_coords_cvref, i)))));
// NOLINTEND
} else {
gsl::at(source_coords, i) = gsl::at(source_coords_cvref, i);
// Can't do math with the cvref typed data so we convert it to the result
// type. If it's a DataVector then we use views. If it's doubles, just a copy
const std::array<ResultT, 3> source_coords = [&source_coords_cvref]() {
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you just do this:

const auto& x = dereference_wrapper(source_coords_cvref[0]);
const auto& y = dereference_wrapper(source_coords_cvref[1]);
const auto& z = dereference_wrapper(source_coords_cvref[2]);

Then you don't need any of these const-casts etc.

Copy link
Member

@nilsvu nilsvu Jan 23, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Edit: I think in other maps we just do math with cref(DataVector) and that seems to work, see Wedge.cpp. Haven't checked the details.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No I want to do math with the entire std::array, not the individual components. However, that math won't work between a std::array<DataVector> and std::array<std::reference_wrapper<const DataVector>>. That's why I const cast.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ok const-casts are pretty dangerous so maybe sprinkle some dereferencing in StdArrayHelpers.hpp instead?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We use const_cast in a lot of places for this exact purpose (to set data references). Quick VSCode/git search shows 136 occurrences in 64 files

std::array<ResultT, 3> result{};
for (size_t i = 0; i < 3; i++) {
if constexpr (std::is_same_v<ResultT, DataVector>) {
// NOLINTBEGIN
gsl::at(result, i).set_data_ref(make_not_null(&const_cast<DataVector&>(
dereference_wrapper(gsl::at(source_coords_cvref, i)))));
// NOLINTEND
} else {
gsl::at(result, i) = gsl::at(source_coords_cvref, i);
}
}
}
return result;
}();

const std::array<ResultT, 3> centered_source_coords = source_coords - center_;

Expand All @@ -259,7 +260,7 @@ std::array<tt::remove_cvref_wrap_t<T>, 3> Skew::get_width_deriv(
center_[1] * centered_source_coords_dot_center,
center_[2] * centered_source_coords_dot_center};
// We don't multiply grad_delta by 2 like in the dox because it will just
// canel with the factor of 1/2 in grad_width
// cancel with the factor of 1/2 in grad_width
grad_delta += centered_source_coords *
outer_radius_squared_minus_center_radius_squared_;

Expand Down
10 changes: 6 additions & 4 deletions src/Domain/CoordinateMaps/TimeDependent/Skew.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@ namespace domain::CoordinateMaps::TimeDependent {
* system.
*
* ### Mapped Coordinates
* The Skew coordinte map is given by the mapping
* The Skew coordinate map is given by the mapping
*
* \begin{align}
* \label{eq:map}
Expand All @@ -43,8 +43,10 @@ namespace domain::CoordinateMaps::TimeDependent {
* \end{align}
*
* where $F_y(t)$ and $F_z(t)$ are the angles within the $(x,y)$ plane between
* the undistored $x$-axis and the skewed $\bar{y}$-axis at the \p center,
* represented by `domain::FunctionsOfTime::FunctionOfTime`s.
* the undistorted $x$-axis and the skewed $\bar{y}$-axis at the \p center,
* represented by `domain::FunctionsOfTime::FunctionOfTime`s. The actual
* function of time should have two components; the first corresponds to $y$ and
* the second corresponds to $z$.
*
* $W(\vec{x})$ is a spatial function that should be 1 at the \p center (i.e.
* maximally skewed between the two objects) and fall off to 0 at the
Expand Down Expand Up @@ -141,7 +143,7 @@ namespace domain::CoordinateMaps::TimeDependent {
* \end{equation}
*
* We can bound the root by noticing that in $\ref{eq:map}$, if we substitute
* the extramal values of $W = 0$ and $W = 1$, we get bounds on $\bar{x}$ which
* the extremal values of $W = 0$ and $W = 1$, we get bounds on $\bar{x}$ which
* we can turn into bounds on $x$:
*
* \begin{align}
Expand Down
73 changes: 69 additions & 4 deletions tests/Unit/Domain/CoordinateMaps/TimeDependent/Test_Skew.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,9 @@
#include "Framework/TestingFramework.hpp"

#include <array>
#include <cmath>
#include <cstddef>
#include <limits>
#include <optional>
#include <random>
#include <string>
Expand All @@ -25,10 +27,12 @@

namespace domain {
namespace {
constexpr size_t deriv_order = 2;
using Polynomial = domain::FunctionsOfTime::PiecewisePolynomial<deriv_order>;
using FoftPtr = std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>;

template <typename Generator>
void test(const gsl::not_null<Generator*> generator) {
static constexpr size_t deriv_order = 2;

const double initial_time = 0.5;
double t = initial_time + 0.05;
const double dt = 0.6;
Expand All @@ -42,8 +46,6 @@ void test(const gsl::not_null<Generator*> generator) {
std::uniform_real_distribution<double> point_dist{-5.0, 5.0};
// NOLINTEND

using Polynomial = domain::FunctionsOfTime::PiecewisePolynomial<deriv_order>;
using FoftPtr = std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>;
std::unordered_map<std::string, FoftPtr> f_of_t_list{};
f_of_t_list[function_of_time_name] = std::make_unique<Polynomial>(
initial_time,
Expand Down Expand Up @@ -113,11 +115,74 @@ void test(const gsl::not_null<Generator*> generator) {
CHECK(skew_map == skew_map_deserialized);
CHECK_FALSE(skew_map != skew_map_deserialized);
}

void test_specific_points() {
const std::string function_of_time_name{"Skew"};
const double time = 0.0;
std::unordered_map<std::string, FoftPtr> f_of_t_list{};
// Use pi/4 so math is easy
f_of_t_list[function_of_time_name] = std::make_unique<Polynomial>(
time,
std::array{DataVector{M_PI_4, M_PI_4}, DataVector{2, 0.0},
DataVector{2, 0.0}},
std::numeric_limits<double>::infinity());

const std::array<double, 3> center{-1.0, 1.0, 0.0};
const double outer_radius = 100.0;

const CoordinateMaps::TimeDependent::Skew skew_map{function_of_time_name,
center, outer_radius};

std::array<double, 3> test_point{};
{
INFO("Center");
test_point = center;
CAPTURE(test_point);

const auto mapped_point = skew_map(test_point, time, f_of_t_list);
CHECK_ITERABLE_APPROX(mapped_point, (std::array{0.0, 1.0, 0.0}));
const auto inverse_point =
skew_map.inverse(mapped_point, time, f_of_t_list);
CHECK(inverse_point.has_value());
CHECK_ITERABLE_APPROX(inverse_point.value(), test_point);
// Jacobian not defined at the center
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

do we not need it?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Unclear. SpEC didn't handle this specially (and thus would get a division by zero error if the center was passed to the jacobian).

}
{
INFO("Outer radius");
test_point = std::array{0.0, 0.0, outer_radius};
CAPTURE(test_point);

const auto mapped_point = skew_map(test_point, time, f_of_t_list);
CHECK_ITERABLE_APPROX(mapped_point, test_point);
const auto inverse_point =
skew_map.inverse(mapped_point, time, f_of_t_list);
CHECK(inverse_point.has_value());
CHECK_ITERABLE_APPROX(inverse_point.value(), test_point);
const auto jacobian = skew_map.jacobian(test_point, time, f_of_t_list);
CHECK_ITERABLE_APPROX(jacobian, identity<3>(0.0));
}
{
INFO("Outer radius plus eps");
const double eps = std::numeric_limits<double>::epsilon();
test_point = std::array{0.0, 0.0, outer_radius + eps};
CAPTURE(test_point);

const auto mapped_point = skew_map(test_point, time, f_of_t_list);
CHECK_ITERABLE_APPROX(mapped_point, test_point);
const auto inverse_point =
skew_map.inverse(mapped_point, time, f_of_t_list);
CHECK(inverse_point.has_value());
CHECK_ITERABLE_APPROX(inverse_point.value(), test_point);
const auto jacobian = skew_map.jacobian(test_point, time, f_of_t_list);
CHECK_ITERABLE_APPROX(jacobian, identity<3>(0.0));
}
}
} // namespace

SPECTRE_TEST_CASE("Unit.Domain.CoordinateMaps.TimeDependent.Skew",
"[Domain][Unit]") {
MAKE_GENERATOR(generator);
test(make_not_null(&generator));
test_specific_points();
}
} // namespace domain
Loading