From 367ae712ed19a1d599360489f7f49bd2ad17b0f6 Mon Sep 17 00:00:00 2001 From: Kyle Nelli Date: Thu, 16 Jan 2025 19:28:53 -0800 Subject: [PATCH] Add Skew map --- docs/Images/Domain/CoordinateMaps/Skew.svg | 164 +++++++ .../TimeDependent/CMakeLists.txt | 2 + .../CoordinateMaps/TimeDependent/Skew.cpp | 423 ++++++++++++++++++ .../CoordinateMaps/TimeDependent/Skew.hpp | 292 ++++++++++++ .../TimeDependent/CMakeLists.txt | 1 + .../TimeDependent/Test_Skew.cpp | 123 +++++ 6 files changed, 1005 insertions(+) create mode 100644 docs/Images/Domain/CoordinateMaps/Skew.svg create mode 100644 src/Domain/CoordinateMaps/TimeDependent/Skew.cpp create mode 100644 src/Domain/CoordinateMaps/TimeDependent/Skew.hpp create mode 100644 tests/Unit/Domain/CoordinateMaps/TimeDependent/Test_Skew.cpp diff --git a/docs/Images/Domain/CoordinateMaps/Skew.svg b/docs/Images/Domain/CoordinateMaps/Skew.svg new file mode 100644 index 000000000000..32e59d62f6a6 --- /dev/null +++ b/docs/Images/Domain/CoordinateMaps/Skew.svg @@ -0,0 +1,164 @@ + +image/svg+xml \ No newline at end of file diff --git a/src/Domain/CoordinateMaps/TimeDependent/CMakeLists.txt b/src/Domain/CoordinateMaps/TimeDependent/CMakeLists.txt index 013135709c1e..547951788fe1 100644 --- a/src/Domain/CoordinateMaps/TimeDependent/CMakeLists.txt +++ b/src/Domain/CoordinateMaps/TimeDependent/CMakeLists.txt @@ -9,6 +9,7 @@ spectre_target_sources( RotationMatrixHelpers.cpp RotScaleTrans.cpp Shape.cpp + Skew.cpp SphericalCompression.cpp Translation.cpp ) @@ -24,6 +25,7 @@ spectre_target_headers( RotationMatrixHelpers.hpp RotScaleTrans.hpp Shape.hpp + Skew.hpp SphericalCompression.hpp Translation.hpp ) diff --git a/src/Domain/CoordinateMaps/TimeDependent/Skew.cpp b/src/Domain/CoordinateMaps/TimeDependent/Skew.cpp new file mode 100644 index 000000000000..0e207fdd6d95 --- /dev/null +++ b/src/Domain/CoordinateMaps/TimeDependent/Skew.cpp @@ -0,0 +1,423 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Domain/CoordinateMaps/TimeDependent/Skew.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "DataStructures/DataVector.hpp" +#include "DataStructures/Tensor/EagerMath/DeterminantAndInverse.hpp" +#include "DataStructures/Tensor/Identity.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "Domain/FunctionsOfTime/FunctionOfTime.hpp" +#include "NumericalAlgorithms/RootFinding/TOMS748.hpp" +#include "Utilities/ContainerHelpers.hpp" +#include "Utilities/DereferenceWrapper.hpp" +#include "Utilities/EqualWithinRoundoff.hpp" +#include "Utilities/GenerateInstantiations.hpp" +#include "Utilities/Gsl.hpp" +#include "Utilities/MakeArray.hpp" +#include "Utilities/MakeWithValue.hpp" +#include "Utilities/StdArrayHelpers.hpp" +#include "Utilities/StdHelpers.hpp" +#include "Utilities/StlStreamDeclarations.hpp" +#include "Utilities/TypeTraits/RemoveReferenceWrapper.hpp" + +namespace domain::CoordinateMaps::TimeDependent { + +Skew::Skew(std::string function_of_time_name, + const std::array& center, const double outer_radius) + : f_of_t_name_(std::move(function_of_time_name)), + center_(center), + outer_radius_(outer_radius), + outer_radius_squared_minus_center_radius_squared_(square(outer_radius) - + dot(center_, center_)), + f_of_t_names_({f_of_t_name_}) {} + +template +std::array, 3> Skew::operator()( + const std::array& source_coords, const double time, + const std::unordered_map< + std::string, std::unique_ptr>& + functions_of_time) const { + return map_and_velocity_helper(source_coords, time, functions_of_time, false); +} + +std::optional> Skew::inverse( + const std::array& target_coords, const double time, + const std::unordered_map< + std::string, std::unique_ptr>& + functions_of_time) const { + const auto& function_of_time = functions_of_time.at(f_of_t_name_); + const DataVector func = function_of_time->func_and_deriv(time)[0]; + ASSERT(func.size() == 2, "Expected a function of time with size 2, not " + << func.size() << " in the Skew map."); + + // Short circuit if there's no function of time + if (equal_within_roundoff(func[0], 0.0) and + equal_within_roundoff(func[1], 0.0)) { + return target_coords; + } + + // 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; + } + + // Another short circuit + const double tan_sum = + tan(func[0]) * target_coords[1] + tan(func[1]) * target_coords[2]; + if (equal_within_roundoff(tan_sum, 0.0)) { + return target_coords; + } + + std::array temporary_source_coord = target_coords; + + const auto root_func = [&](const double source_coord_x) -> double { + temporary_source_coord[0] = source_coord_x; + const double width = get_width(temporary_source_coord); + return source_coord_x + width * tan_sum - target_coords[0]; + }; + + // \bar{x} -> x + w * (tan(f_y)*y + tan(f_z)*z) and 0<=w<=1 so then + // x <= \bar{x} <= x + (tan(f_y)*y + tan(f_z)*z) which imples + // 0 <= \bar{x} - x <= (tan(f_y)*y + tan(f_z)*z) which imples + // -\bar{x} <= - x <= -\bar{x} + (tan(f_y)*y + tan(f_z)*z) which imples + // \bar{x} >= x >= \bar{x} - (tan(f_y)*y + tan(f_z)*z) + // This last line is our bracket for the root + + // We give a small epsilon to give the root find a little buffer. It is highly + // unlikely there will be another root within these bounds + const double eps = std::numeric_limits::epsilon() * 100.0; + + const double x_0 = target_coords[0] - tan_sum - eps; + const double x_1 = target_coords[0] + eps; + const double lower_x = std::min(x_0, x_1); + const double upper_x = std::max(x_0, x_1); + const double lower_func = root_func(lower_x); + const double upper_func = root_func(upper_x); + + if (equal_within_roundoff(lower_func, 0.0)) { + temporary_source_coord[0] = lower_x; + return temporary_source_coord; + } + + if (equal_within_roundoff(upper_func, 0.0)) { + temporary_source_coord[0] = upper_x; + return temporary_source_coord; + } + + if (lower_func * upper_func > 0.0) { + ERROR("No root in Skew map inverse!. Target coord = " + << target_coords << ", lower_x = " << lower_x + << ", lower_func = " << lower_func << ", upper_x = " << upper_x + << ", upper_func = " << upper_func << ", tan_sum = " << tan_sum); + } + + temporary_source_coord[0] = RootFinder::toms748( + root_func, lower_x, upper_x, lower_func, upper_func, eps, eps); + + return temporary_source_coord; +} + +template +std::array, 3> Skew::frame_velocity( + const std::array& source_coords, const double time, + const std::unordered_map< + std::string, std::unique_ptr>& + functions_of_time) const { + return map_and_velocity_helper(source_coords, time, functions_of_time, true); +} + +template +tnsr::Ij, 3, Frame::NoFrame> Skew::jacobian( + const std::array& source_coords, const double time, + const std::unordered_map< + std::string, std::unique_ptr>& + functions_of_time) const { + using ResultT = tt::remove_cvref_wrap_t; + + const ResultT width = get_width(source_coords); + const std::array width_deriv = get_width_deriv(source_coords); + + const auto& function_of_time = functions_of_time.at(f_of_t_name_); + const DataVector func = function_of_time->func_and_deriv(time)[0]; + ASSERT(func.size() == 2, "Expected a function of time with size 2, not " + << func.size() << " in the Skew map."); + + auto result = identity<3>(dereference_wrapper(source_coords[0])); + + const DataVector tan_func = tan(func); + // Temporarily use component to avoid allocation + auto& tan_sum = get<0, 2>(result); + tan_sum = tan_func[0] * source_coords[1] + tan_func[1] * source_coords[2]; + + get<0, 0>(result) += width_deriv[0] * tan_sum; + get<0, 1>(result) = width_deriv[1] * tan_sum + width * tan_func[0]; + get<0, 2>(result) = width_deriv[2] * tan_sum + width * tan_func[1]; + + return result; +} + +template +tnsr::Ij, 3, Frame::NoFrame> Skew::inv_jacobian( + const std::array& source_coords, const double time, + const std::unordered_map< + std::string, std::unique_ptr>& + functions_of_time) const { + return determinant_and_inverse( + jacobian(source_coords, time, functions_of_time)) + .second; +} + +template +tt::remove_cvref_wrap_t Skew::get_delta( + const std::array& centered_source_coords) const { + return square(dot(centered_source_coords, center_)) + + dot(centered_source_coords, centered_source_coords) * + outer_radius_squared_minus_center_radius_squared_; +} + +template +tt::remove_cvref_wrap_t Skew::get_width( + const std::array& source_coords) const { + using ResultT = tt::remove_cvref_wrap_t; + std::array centered_source_coords{}; + for (size_t i = 0; i < 3; i++) { + gsl::at(centered_source_coords, i) = + gsl::at(source_coords, i) - gsl::at(center_, i); + } + // Will be reused for result + ResultT lambda = (dot(centered_source_coords, center_) + + sqrt(get_delta(centered_source_coords))) / + outer_radius_squared_minus_center_radius_squared_; + + ResultT& result = lambda; + + for (size_t i = 0; i < get_size(result); i++) { + if (get_element(lambda, i) <= 0.0 and + equal_within_roundoff(get_element(lambda, i), 0.0)) { + get_element(result, i) = 1.0; + } else if (get_element(lambda, i) >= 1.0 and + equal_within_roundoff(get_element(lambda, i), 1.0)) { + get_element(result, i) = 0.0; + } else if (get_element(lambda, i) > 0.0 and get_element(lambda, i) < 1.0) { + get_element(result, i) = 0.5 * (1.0 + cos(M_PI * get_element(lambda, i))); + } else { + using ::operator<<; + const std::vector bad_point{ + get_element(centered_source_coords[0], i), + get_element(centered_source_coords[1], i), + get_element(centered_source_coords[2], i)}; + ERROR("Skew map: Centered source coordinate " + << bad_point << " not in valid region. Lambda is " + << get_element(lambda, i) << ". Center is " << center_ + << ". Outer radius is " << outer_radius_); + } + } + + return result; +} + +template +std::array, 3> Skew::get_width_deriv( + const std::array& source_coords_cvref) const { + using ResultT = tt::remove_cvref_wrap_t; + std::array source_coords{}; + for (size_t i = 0; i < 3; i++) { + if constexpr (std::is_same_v) { + // NOLINTBEGIN + gsl::at(source_coords, i) + .set_data_ref(make_not_null(&const_cast( + dereference_wrapper(gsl::at(source_coords_cvref, i))))); + // NOLINTEND + } else { + gsl::at(source_coords, i) = gsl::at(source_coords_cvref, i); + } + } + + const std::array centered_source_coords = source_coords - center_; + + const ResultT delta = get_delta(centered_source_coords); + + const ResultT lambda = (dot(centered_source_coords, center_) + sqrt(delta)) / + outer_radius_squared_minus_center_radius_squared_; + + const ResultT centered_source_coords_dot_center = + dot(centered_source_coords, center_); + std::array grad_delta{ + center_[0] * centered_source_coords_dot_center, + 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 + grad_delta += centered_source_coords * + outer_radius_squared_minus_center_radius_squared_; + + // Start off with grad_width = grad_lambda. We ignore the factor of 1/2 on the + // second term because it will cancel with the factor of 2 from grad_delta + std::array grad_width = + (center_ + grad_delta / (sqrt(delta))) / + outer_radius_squared_minus_center_radius_squared_; + for (size_t i = 0; i < get_size(lambda); i++) { + // sin(lambda * M_PI) is zero at both lambda=0 and lambda=1 + if ((get_element(lambda, i) <= 0.0 and + equal_within_roundoff(get_element(lambda, i), 0.0)) or + (get_element(lambda, i) >= 1.0 and + equal_within_roundoff(get_element(lambda, i), 1.0))) { + for (size_t j = 0; j < 3; j++) { + get_element(gsl::at(grad_width, j), i) = 0.0; + } + } else if (get_element(lambda, i) > 0.0 and get_element(lambda, i) < 1.0) { + for (size_t j = 0; j < 3; j++) { + // grad_width already has the factor of grad_lambda + get_element(gsl::at(grad_width, j), i) *= + -0.5 * M_PI * sin(M_PI * get_element(lambda, i)); + } + } else { + using ::operator<<; + const std::vector bad_point{ + get_element(centered_source_coords[0], i), + get_element(centered_source_coords[1], i), + get_element(centered_source_coords[2], i)}; + ERROR("Skew map: Centered source coordinate " + << bad_point << " not in valid region. Lambda is " + << get_element(lambda, i) << ". Center is " << center_ + << ". Outer radius is " << outer_radius_); + } + } + + return grad_width; +} + +template +std::array, 3> Skew::map_and_velocity_helper( + const std::array& source_coords, const double time, + const std::unordered_map< + std::string, std::unique_ptr>& + functions_of_time, + const bool return_velocity) const { + using ResultT = tt::remove_cvref_wrap_t; + std::array result{}; + if (return_velocity) { + result = make_array<3>( + make_with_value(dereference_wrapper(source_coords[0]), 0.0)); + } else { + for (size_t i = 0; i < 3; i++) { + gsl::at(result, i) = gsl::at(source_coords, i); + } + } + + // Currently result is the source coords, but with the proper return type + const ResultT width = get_width(source_coords); + + const auto& function_of_time = functions_of_time.at(f_of_t_name_); + const auto func_and_deriv = function_of_time->func_and_deriv(time); + ASSERT(func_and_deriv[0].size() == 2, + "Expected a function of time with size 2, not " + << func_and_deriv[0].size() << " in the Skew map."); + + if (return_velocity) { + const auto& func = func_and_deriv[0]; + const auto& deriv = func_and_deriv[1]; + + result[0] = + width * (deriv[0] * (1.0 + square(tan(func[0]))) * source_coords[1] + + deriv[1] * (1.0 + square(tan(func[1]))) * source_coords[2]); + } else { + result[0] += width * (tan(func_and_deriv[0][0]) * source_coords[1] + + tan(func_and_deriv[0][1]) * source_coords[2]); + } + + return result; +} + +void Skew::pup(PUP::er& p) { + size_t version = 0; + p | version; + // Remember to increment the version number when making changes to this + // function. Retain support for unpacking data written by previous versions + // whenever possible. See `Domain` docs for details. + if (version >= 0) { + p | f_of_t_name_; + p | center_; + p | outer_radius_; + } + + // No need to pup this because it is uniquely determined by f_of_t_name_ + if (p.isUnpacking()) { + outer_radius_squared_minus_center_radius_squared_ = + square(outer_radius_) - dot(center_, center_); + f_of_t_names_.clear(); + f_of_t_names_.insert(f_of_t_name_); + } +} + +bool operator==(const Skew& lhs, const Skew& rhs) { + return lhs.f_of_t_name_ == rhs.f_of_t_name_ and lhs.center_ == rhs.center_ and + lhs.outer_radius_ == rhs.outer_radius_; +} + +bool operator!=(const Skew& lhs, const Skew& rhs) { return not(lhs == rhs); } + +#define DTYPE(data) BOOST_PP_TUPLE_ELEM(0, data) + +#define INSTANTIATE(_, data) \ + template std::array, 3> \ + Skew::operator()( \ + const std::array& source_coords, double time, \ + const std::unordered_map< \ + std::string, \ + std::unique_ptr>& \ + functions_of_time) const; \ + template std::array, 3> \ + Skew::frame_velocity( \ + const std::array& source_coords, const double time, \ + const std::unordered_map< \ + std::string, \ + std::unique_ptr>& \ + functions_of_time) const; \ + template tnsr::Ij, 3, Frame::NoFrame> \ + Skew::jacobian( \ + const std::array& source_coords, double time, \ + const std::unordered_map< \ + std::string, \ + std::unique_ptr>& \ + functions_of_time) const; \ + template tnsr::Ij, 3, Frame::NoFrame> \ + Skew::inv_jacobian( \ + const std::array& source_coords, double time, \ + const std::unordered_map< \ + std::string, \ + std::unique_ptr>& \ + functions_of_time) const; \ + template tt::remove_cvref_wrap_t Skew::get_delta( \ + const std::array& centered_source_coords) const; \ + template tt::remove_cvref_wrap_t Skew::get_width( \ + const std::array& source_coords) const; \ + template std::array, 3> \ + Skew::get_width_deriv(const std::array& source_coords_cvref) \ + const; \ + template std::array, 3> \ + Skew::map_and_velocity_helper( \ + const std::array& source_coords, double time, \ + const domain::FunctionsOfTimeMap& functions_of_time, \ + bool return_velocity) const; + +GENERATE_INSTANTIATIONS(INSTANTIATE, (double, DataVector, + std::reference_wrapper, + std::reference_wrapper)) + +#undef DTYPE +#undef INSTANTIATE + +} // namespace domain::CoordinateMaps::TimeDependent diff --git a/src/Domain/CoordinateMaps/TimeDependent/Skew.hpp b/src/Domain/CoordinateMaps/TimeDependent/Skew.hpp new file mode 100644 index 000000000000..77890eeb16ea --- /dev/null +++ b/src/Domain/CoordinateMaps/TimeDependent/Skew.hpp @@ -0,0 +1,292 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +#include "DataStructures/Tensor/TypeAliases.hpp" +#include "Domain/FunctionsOfTime/FunctionOfTime.hpp" +#include "Utilities/TypeTraits/RemoveReferenceWrapper.hpp" + +/// \cond +namespace PUP { +class er; +} // namespace PUP +/// \endcond + +namespace domain::CoordinateMaps::TimeDependent { +/*! + * \ingroup CoordMapsTimeDependentGroup + * \brief %Time dependent coordinate map that keeps the $y$ and $z$ coordinates + * unchanged, but will distort (or skew) the $x$ coordinate based on functions + * of time and radial distance to the origin. + * + * \details This coordinate map is only available in 3 dimensions and assumes + * all coordinates are given with respect to the origin of the coordinate + * system. + * + * ### Mapped Coordinates + * The Skew coordinte map is given by the mapping + * + * \begin{align} + * \label{eq:map} + * \bar{x} &= x + W(\vec{x})\left(\tan(F_y(t))y + \tan(F_z(t))z\right) \\ + * \bar{y} &= y \\ + * \bar{z} &= z + * \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. + * + * $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 + * \p outer_radius, $R$. Typically the \p outer_radius should be set to the + * envelope radius for the `domain::creators::BinaryCompactObject`. Here, $W$ is + * chosen to be + * + * \begin{equation} + * \label{eq:W} + * W(\vec{x}) = \frac{1}{2}\left(1 + \cos(\pi\lambda(\vec{x}))\right). + * \end{equation} + * + * When the $\cos$ term is $-1$, then $W = 0$, and when the $\cos$ term is 1, + * then $W = 1$. Therefore, the function $\lambda(\vec{x})$ must go from 0 at + * the \p center, to 1 at $R$. In order to derive + * $\lambda(\vec{x})$, we introduce three vectors, shown in the figure, all + * originating from the origin: $\vec{x}$ is our source coordinate, $\vec{x}_0$ + * is the \p center, and $\vec{x}_P$ is the vector to the intersection point + * between the \p outer_radius and the ray connecting $\vec{x}_0$ and $\vec{x}$. + * Therefore, $|\vec{x}_P| = R$. + * + * \image html Skew.svg "Vectors for Skew map" + * + * We choose $\lambda(\vec{x})$ to linearly go from 0 at $\vec{x}_0$ to 1 at + * $\vec{x}_P$ with the form + * + * \begin{equation} + * \label{eq:lambda} + * \lambda(\vec{x}) = \frac{|\vec{x} - \vec{x}_0|}{|\vec{x}_P - \vec{x}_0|} + * \end{equation} + * + * We must now find a functional form for $|\vec{x}_P - \vec{x}_0|$ that depends + * on $\vec{x}$. We define + * + * \begin{equation} + * \label{eq:x_p_x_0} + * \vec{x}_P - \vec{x}_0 = S\frac{\vec{x}-\vec{x}_0}{|\vec{x}-\vec{x}_0|}, + * \end{equation} + * + * since the ray from $\vec{x}_0$ to $\vec{x}_P$ is just a rescaling of the ray + * from $\vec{x}_0$ to $\vec{x}$. This implies that Eq. $\ref{eq:lambda}$ + * becomes + * + * \begin{equation} + * \label{eq:lambda_S} + * \lambda(\vec{x}) = \frac{|\vec{x} - \vec{x}_0|}{|S|}. + * \end{equation} + * + * Now using Eq. $\ref{eq:x_p_x_0}$, we take the magnitude squared of + * $\vec{x}_P$ and the fact that $|\vec{x}_P|^2 = R^2$ to get + * + * \begin{equation} + * \label{eq:roots} + * 0 = S^2 + 2\frac{(\vec{x}-\vec{x}_0)\cdot\vec{x}_0}{|\vec{x}-\vec{x}_0}S + * + (|\vec{x}_0|^2 - R^2). + * \end{equation} + * + * Using the [alternative + * form](https://en.wikipedia.org/wiki/Quadratic_formula#Equivalent_formulations) + * of the quadratic forumla, choosing the positive root, and simplifying a bit, + * we arrive at + * + * \begin{equation} + * \label{eq:S} + * S = \frac{(R^2 - |\vec{x}_0|^2)|\vec{x} - \vec{x}_0|}{(\vec{x} - + * \vec{x}_0)\cdot\vec{x}_0 + \sqrt{[(\vec{x} - \vec{x}_0)\cdot\vec{x}_0]^2 + + * |\vec{x} - \vec{x}_0|^2(R^2 - |\vec{x}_0|^2)}}. + * \end{equation} + * + * Substituting this into Eq. $\ref{eq:lambda_S}$, we get + * + * \begin{equation} + * \label{eq:lambda_S_2} + * \lambda(\vec{x}) = \frac{(\vec{x} - \vec{x}_0)\cdot\vec{x}_0 + + * \sqrt{\Delta}}{R^2 - |\vec{x}_0|^2}. + * \end{equation} + * + * where we have defined + * + * \begin{equation} + * \label{eq:delta} + * \Delta = [(\vec{x} - \vec{x}_0)\cdot\vec{x}_0]^2 + |\vec{x} - + * \vec{x}_0|^2(R^2 - |\vec{x}_0|^2). + * \end{equation} + * + * ### Inverse + * To find the inverse, we need to solve a 1D root find for the $x$ component of + * the coordinate. The inverse of the $\bar{y}$ and $\bar{z}$ coordinates are + * trivial because there was no mapping. The equation we need to find the root + * of is + * + * \begin{equation} + * 0 = x + W(\vec{x})\left(\tan(F_y(t))y + \tan(F_z(t))z\right) - \bar{x} + * \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 + * we can turn into bounds on $x$: + * + * \begin{align} + * x &<=& \bar{x} &<=& x + (\tan(F_y)y + \tan(F_z)z) \\ + * 0 &<=& \bar{x} - x &<=& (\tan(F_y)y + \tan(F_z)z) \\ + * -\bar{x} &<=& - x &<=& -\bar{x} + (\tan(F_y)y + \tan(F_z)z) \\ + * \bar{x} &>=& x &>=& \bar{x} - (\tan(F_y)y + \tan(F_z)z) + * \end{align} + * + * where on each line we just made simple arithmetic operations. We pad each + * bound by $1e-14$ just to avoid roundoff issues. If either of the bounds is + * within roundoff of zero, the map is the idenitiy at that point and we forgo + * the root find. The root that is found is the original $x$ coordinate. + * + * ### Frame Velocity + * Taking the time derivative of $\ref{eq:map}$, the frame velocity is + * + * \begin{align} + * \label{eq:frame_vel} + * \dot{\bar{x}} &= W(\vec{x})\left(\dot{F}_y(t)(1 + \tan^2(F_y(t)))y + + * \dot{F}_z(t)(1 + \tan^2(F_z(t))z)\right) \\ + * \dot{\bar{y}} &= 0 \\ + * \dot{\bar{z}} &= 0 + * \end{align} + * + * ### Jacobian and Inverse Jacobian + * Considering the first terms in each equation of $\ref{eq:map}$, part of the + * jacobian will be the identity matrix. The rest will come from only the $x$ + * equation. Therefore we can express the jacobian as + * + * \begin{equation} + * \frac{\partial\bar{x}^i}{\partial x^j} = \delta^i_j + {W^i}_j + * \end{equation} + * + * where all components of ${W^i}_j$ are zero except the following + * + * \begin{align} + * {W^0}_0 &= \frac{\partial(\bar{x}-x)}{\partial x} &= \frac{\partial + * W(\vec{x})}{\partial x}&\left(\tan(F_y(t))y + \tan(F_z(t))z\right), \\ + * {W^0}_1 &= \frac{\partial(\bar{x}-x)}{\partial y} &= \frac{\partial + * W(\vec{x})}{\partial y}&\left(\tan(F_y(t))y + \tan(F_z(t))z\right) + + * W\tan(F_y(t)), \\ + * {W^0}_2 &= \frac{\partial(\bar{x}-x)}{\partial z} &= \frac{\partial + * W(\vec{x})}{\partial z}&\left(\tan(F_y(t))y + \tan(F_z(t))z\right) + + * W\tan(F_z(t)). + * \end{align} + * + * The gradient of $W(\vec{x})$ (Eq. $\ref{eq:W}$) is given by + * + * \begin{equation} + * \frac{\partial W(\vec{x})}{\partial x^i} = + * -\frac{\pi}{2}\frac{\partial\lambda(\vec{x})}{\partial + * x^i}\sin(\pi\lambda(\vec{x})). + * \end{equation} + * + * The gradient of $\lambda(\vec{x})$ (Eq. $\ref{eq:lambda_S_2}$) is given by + * + * \begin{equation} + * \frac{\partial \lambda(\vec{x})}{\partial x^i} = + * \frac{1}{R^2 - |\vec{x}_0|^2}\left(x_0^i + + * \frac{1}{2\sqrt{\Delta}}\frac{\partial \Delta}{\partial x^i}\right). + * \end{equation} + * + * And finally, the gradient of $\Delta$ (Eq. $\ref{eq:delta}$) is given by + * + * \begin{equation} + * \frac{\partial \Delta}{\partial x^i} = 2x_0^i\left[\left(\vec{x} - + * \vec{x}_0\right)\cdot\vec{x}_0\right] + 2(R^2 - |\vec{x}_0|^2)(x^i - + * x_0^i). \end{equation} + * + * The inverse jacobian is computed by numerically inverting the jacobian. + */ +class Skew { + public: + static constexpr size_t dim = 3; + + Skew(std::string function_of_time_name, const std::array& center, + double outer_radius); + Skew() = default; + + template + std::array, 3> operator()( + const std::array& source_coords, double time, + const domain::FunctionsOfTimeMap& functions_of_time) const; + + /// The inverse function is only callable with doubles because the inverse + /// might fail if called for a point out of range, and it is unclear + /// what should happen if the inverse were to succeed for some points in a + /// DataVector but fail for other points. + std::optional> inverse( + const std::array& target_coords, double time, + const domain::FunctionsOfTimeMap& functions_of_time) const; + + template + std::array, 3> frame_velocity( + const std::array& source_coords, double time, + const domain::FunctionsOfTimeMap& functions_of_time) const; + + template + tnsr::Ij, 3, Frame::NoFrame> jacobian( + const std::array& source_coords, double time, + const domain::FunctionsOfTimeMap& functions_of_time) const; + + template + tnsr::Ij, 3, Frame::NoFrame> inv_jacobian( + const std::array& source_coords, double time, + const domain::FunctionsOfTimeMap& functions_of_time) const; + + // NOLINTNEXTLINE(google-runtime-references) + void pup(PUP::er& p); + + static bool is_identity() { return false; } + + const std::unordered_set& function_of_time_names() const { + return f_of_t_names_; + } + + private: + template + tt::remove_cvref_wrap_t get_delta( + const std::array& centered_source_coords) const; + + template + tt::remove_cvref_wrap_t get_width( + const std::array& source_coords) const; + + template + std::array, 3> get_width_deriv( + const std::array& source_coords_cvref) const; + + template + std::array, 3> map_and_velocity_helper( + const std::array& source_coords, double time, + const domain::FunctionsOfTimeMap& functions_of_time, + bool return_velocity) const; + + // NOLINTNEXTLINE(readability-redundant-declaration) + friend bool operator==(const Skew& lhs, const Skew& rhs); + std::string f_of_t_name_; + std::array center_{}; + double outer_radius_{}; + double outer_radius_squared_minus_center_radius_squared_{}; + std::unordered_set f_of_t_names_; +}; + +bool operator!=(const Skew& lhs, const Skew& rhs); + +} // namespace domain::CoordinateMaps::TimeDependent diff --git a/tests/Unit/Domain/CoordinateMaps/TimeDependent/CMakeLists.txt b/tests/Unit/Domain/CoordinateMaps/TimeDependent/CMakeLists.txt index 45a0be51bdbf..ef2605eb61bd 100644 --- a/tests/Unit/Domain/CoordinateMaps/TimeDependent/CMakeLists.txt +++ b/tests/Unit/Domain/CoordinateMaps/TimeDependent/CMakeLists.txt @@ -10,6 +10,7 @@ set(LIBRARY_SOURCES Test_RotationMatrixHelpers.cpp Test_RotScaleTrans.cpp Test_Shape.cpp + Test_Skew.cpp Test_SphericalCompression.cpp Test_Translation.cpp ) diff --git a/tests/Unit/Domain/CoordinateMaps/TimeDependent/Test_Skew.cpp b/tests/Unit/Domain/CoordinateMaps/TimeDependent/Test_Skew.cpp new file mode 100644 index 000000000000..2bd29064d047 --- /dev/null +++ b/tests/Unit/Domain/CoordinateMaps/TimeDependent/Test_Skew.cpp @@ -0,0 +1,123 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Framework/TestingFramework.hpp" + +#include +#include +#include +#include +#include +#include + +#include "DataStructures/DataVector.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "Domain/CoordinateMaps/TimeDependent/Skew.hpp" +#include "Domain/FunctionsOfTime/FunctionOfTime.hpp" +#include "Domain/FunctionsOfTime/PiecewisePolynomial.hpp" +#include "Framework/TestHelpers.hpp" +#include "Helpers/DataStructures/MakeWithRandomValues.hpp" +#include "Helpers/Domain/CoordinateMaps/TestMapHelpers.hpp" +#include "Utilities/ConstantExpressions.hpp" +#include "Utilities/Gsl.hpp" +#include "Utilities/StdArrayHelpers.hpp" +#include "Utilities/TypeTraits.hpp" + +namespace domain { +namespace { +template +void test(const gsl::not_null 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; + const double expiration_time = 20.0; + + const std::string function_of_time_name{"Skew"}; + + // NOLINTBEGIN + std::uniform_real_distribution fot_dist{-0.01, 0.01}; + std::uniform_real_distribution outer_radius_dist{50.0, 150.0}; + std::uniform_real_distribution point_dist{-5.0, 5.0}; + // NOLINTEND + + using Polynomial = domain::FunctionsOfTime::PiecewisePolynomial; + using FoftPtr = std::unique_ptr; + std::unordered_map f_of_t_list{}; + f_of_t_list[function_of_time_name] = std::make_unique( + initial_time, + std::array{make_with_random_values( + generator, make_not_null(&fot_dist), DataVector{2, 0.0}), + make_with_random_values( + generator, make_not_null(&fot_dist), DataVector{2, 0.0}), + DataVector{2, 0.0}}, + expiration_time); + + const std::array center{ + point_dist(*generator), point_dist(*generator), point_dist(*generator)}; + const double outer_radius = outer_radius_dist(*generator); + CAPTURE(center); + CAPTURE(outer_radius); + + const CoordinateMaps::TimeDependent::Skew skew_map{function_of_time_name, + center, outer_radius}; + + CHECK(skew_map.function_of_time_names().contains(function_of_time_name)); + + // test serialized/deserialized map + const auto skew_map_deserialized = serialize_and_deserialize(skew_map); + CHECK(skew_map_deserialized.function_of_time_names().contains( + function_of_time_name)); + + while (t < expiration_time) { + CAPTURE(t); + const std::array point_xi{ + point_dist(*generator), point_dist(*generator), point_dist(*generator)}; + + const std::array dv_point_xi{ + make_with_random_values( + generator, make_not_null(&point_dist), DataVector{10, 0.0}), + make_with_random_values( + generator, make_not_null(&point_dist), DataVector{10, 0.0}), + make_with_random_values( + generator, make_not_null(&point_dist), DataVector{10, 0.0})}; + + const auto run_checks = [&](const auto& points) { + CAPTURE(points); + test_jacobian(skew_map, points, t, f_of_t_list); + test_inv_jacobian(skew_map, points, t, f_of_t_list); + test_frame_velocity(skew_map, points, t, f_of_t_list); + + test_jacobian(skew_map_deserialized, points, t, f_of_t_list); + test_inv_jacobian(skew_map_deserialized, points, t, f_of_t_list); + test_frame_velocity(skew_map_deserialized, points, t, f_of_t_list); + }; + + run_checks(point_xi); + test_coordinate_map_argument_types(skew_map, point_xi, t, f_of_t_list); + test_coordinate_map_argument_types(skew_map_deserialized, point_xi, t, + f_of_t_list); + test_inverse_map(skew_map, point_xi, t, f_of_t_list); + test_inverse_map(skew_map_deserialized, point_xi, t, f_of_t_list); + run_checks(dv_point_xi); + + t += dt; + } + + // Check inequivalence operator + CHECK_FALSE(skew_map != skew_map); + CHECK_FALSE(skew_map_deserialized != skew_map_deserialized); + + // Check serialization + CHECK(skew_map == skew_map_deserialized); + CHECK_FALSE(skew_map != skew_map_deserialized); +} +} // namespace + +SPECTRE_TEST_CASE("Unit.Domain.CoordinateMaps.TimeDependent.Skew", + "[Domain][Unit]") { + MAKE_GENERATOR(generator); + test(make_not_null(&generator)); +} +} // namespace domain