From 6656db23c767158c828b5c3c8908c2c02debb88a Mon Sep 17 00:00:00 2001 From: Kyle Nelli Date: Thu, 16 Jan 2025 19:28:53 -0800 Subject: [PATCH] Add Skew map --- .../TimeDependent/CMakeLists.txt | 2 + .../CoordinateMaps/TimeDependent/Skew.cpp | 408 ++++++++++++++++++ .../CoordinateMaps/TimeDependent/Skew.hpp | 245 +++++++++++ .../TimeDependent/CMakeLists.txt | 1 + .../TimeDependent/Test_Skew.cpp | 245 +++++++++++ 5 files changed, 901 insertions(+) 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/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..1567a85ac7d2 --- /dev/null +++ b/src/Domain/CoordinateMaps/TimeDependent/Skew.cpp @@ -0,0 +1,408 @@ +// 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 + +#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_(square(outer_radius_)), + f_of_t_names_({f_of_t_name_}) { + if (magnitude(center_) >= outer_radius) { + ERROR("Center of Skew map " + << center_ << " with radius " << magnitude(center_) + << " must be within the outer radius " << outer_radius_); + } +} + +template +std::array, 3> Skew::operator()( + const std::array& source_coords, const double time, + const domain::FunctionsOfTimeMap& 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 domain::FunctionsOfTimeMap& 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 circuit if we're at the outer radius + const double target_radius = magnitude(target_coords); + if (equal_within_roundoff(target_radius, outer_radius_)) { + return target_coords; + } else if (target_radius >= outer_radius_) { + return std::nullopt; + } + + // Another short circuit. Target y & z are the same as source y & z + const double tan_sum = tan(func[0]) * (target_coords[1] - center_[1]) + + tan(func[1]) * (target_coords[2] - center_[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-y_C) + tan(f_z)*(z-z_C)) and 0<=w<=1, then + // x <= \bar{x} <= x + (tan(f_y)*(y-y_C) + tan(f_z)*(z-z_C)), ==> + // 0 <= \bar{x} - x <= (tan(f_y)*(y-y_C) + tan(f_z)*(z-z_C)), ==> + // -\bar{x} <= - x <= -\bar{x} + (tan(f_y)*(y-y_C) + tan(f_z)*(z-z_C)), ==> + // \bar{x} >= x >= \bar{x} - (tan(f_y)*(y-y_C) + tan(f_z)*(z-z_C)) + // 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); + + // 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, 5 * eps)) { + 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 domain::FunctionsOfTimeMap& 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 domain::FunctionsOfTimeMap& functions_of_time) const { + using ResultT = tt::remove_cvref_wrap_t; + + check_for_singular_map(source_coords, time, functions_of_time); + + 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] - center_[1]) + + tan_func[1] * (source_coords[2] - center_[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 domain::FunctionsOfTimeMap& functions_of_time) const { + return determinant_and_inverse( + jacobian(source_coords, time, functions_of_time)) + .second; +} + +template +tt::remove_cvref_wrap_t Skew::get_width( + const std::array& source_coords) const { + using ResultT = tt::remove_cvref_wrap_t; + // Will be reused for result + ResultT lambda = dot(source_coords, source_coords) / outer_radius_squared_; + + ResultT& result = lambda; + + // Go point by point because we have to check roundoff + for (size_t i = 0; i < get_size(result); i++) { + if (equal_within_roundoff(get_element(lambda, i), 0.0)) { + get_element(result, i) = 1.0; + } else if (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(dereference_wrapper(source_coords[0]), i), + get_element(dereference_wrapper(source_coords[1]), i), + get_element(dereference_wrapper(source_coords[2]), i)}; + ERROR("Skew map: 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) const { + using ResultT = tt::remove_cvref_wrap_t; + + const ResultT lambda = + dot(source_coords, source_coords) / outer_radius_squared_; + + std::array grad_width{}; + for (size_t i = 0; i < 3; i++) { + gsl::at(grad_width, i) = gsl::at(source_coords, i); + } + // Factors of two cancel from grad lambda and pi/2 + grad_width *= -M_PI / outer_radius_squared_ * sin(M_PI * lambda); + + for (size_t i = 0; i < get_size(lambda); i++) { + // sin(lambda * M_PI) is zero at both lambda=0 and lambda=1 + if (equal_within_roundoff(get_element(lambda, i), 0.0) or + 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 or get_element(lambda, i) > 1.0) { + using ::operator<<; + const std::vector bad_point{ + get_element(dereference_wrapper(source_coords[0]), i), + get_element(dereference_wrapper(source_coords[1]), i), + get_element(dereference_wrapper(source_coords[2]), i)}; + ERROR("Skew map: 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 FunctionsOfTimeMap& functions_of_time, + const bool return_velocity) const { + using ResultT = tt::remove_cvref_wrap_t; + + check_for_singular_map(source_coords, time, functions_of_time); + + 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] - center_[1]) + + deriv[1] * (1.0 + square(tan(func[1]))) * + (source_coords[2] - center_[2])); + } else { + result[0] += + width * (tan(func_and_deriv[0][0]) * (source_coords[1] - center_[1]) + + tan(func_and_deriv[0][1]) * (source_coords[2] - center_[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 these because they are uniquely determined + if (p.isUnpacking()) { + outer_radius_squared_ = square(outer_radius_); + f_of_t_names_.clear(); + f_of_t_names_.insert(f_of_t_name_); + } +} + +template +void Skew::check_for_singular_map( + const std::array& source_coords, const double time, + const FunctionsOfTimeMap& functions_of_time) const { + (void)source_coords; + (void)time; + (void)functions_of_time; + +#ifdef SPECTRE_DEBUG + using ResultT = tt::remove_cvref_wrap_t; + const ResultT lambda = + dot(source_coords, source_coords) / outer_radius_squared_; + + const auto& function_of_time = functions_of_time.at(f_of_t_name_); + const auto 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."); + const ResultT tan_sum = tan(func[0]) * (source_coords[1] - center_[1]) + + tan(func[1]) * (source_coords[2] - center_[2]); + + const ResultT x_deriv_of_mapped_x_coord = + 1.0 - tan_sum * M_PI * dereference_wrapper(source_coords[0]) / + outer_radius_squared_ * sin(M_PI * lambda); + + for (size_t i = 0; i < get_size(lambda); i++) { + if (get_element(x_deriv_of_mapped_x_coord, i) < 0.0) { + using ::operator<<; + const std::vector bad_point{ + get_element(dereference_wrapper(source_coords[0]), i), + get_element(dereference_wrapper(source_coords[1]), i), + get_element(dereference_wrapper(source_coords[2]), i)}; + ERROR("Skew map is singular. Found for point " + << bad_point << ". Important values are:\n Outer radius, R = " + << outer_radius_ << "\n |x|^2/R^2 = " << get_element(lambda, i) + << "\n tan(F_y)*y + tan(F_z)*z = " << tan_sum + << "\n F_y = " << func[0] << "\n F_z = " << func[1] + << "\n tan(F_y) = " << tan(func[0]) + << "\n tan(F_z) = " << tan(func[1]) << "\n"); + } + } +#endif +} + +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 domain::FunctionsOfTimeMap& functions_of_time) const; \ + template std::array, 3> \ + Skew::frame_velocity( \ + const std::array& source_coords, const double time, \ + const domain::FunctionsOfTimeMap& functions_of_time) const; \ + template tnsr::Ij, 3, Frame::NoFrame> \ + Skew::jacobian(const std::array& source_coords, double time, \ + const domain::FunctionsOfTimeMap& functions_of_time) const; \ + template tnsr::Ij, 3, Frame::NoFrame> \ + Skew::inv_jacobian( \ + const std::array& source_coords, double time, \ + const domain::FunctionsOfTimeMap& functions_of_time) 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) \ + 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; \ + template void Skew::check_for_singular_map( \ + const std::array& source_coords, double time, \ + const domain::FunctionsOfTimeMap& functions_of_time) 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..01f88fed0afe --- /dev/null +++ b/src/Domain/CoordinateMaps/TimeDependent/Skew.hpp @@ -0,0 +1,245 @@ +// 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 is + * intended to be used in the BinaryCompactObject domain. + * + * ### Mapped Coordinates + * The Skew coordinate map is given by the mapping + * + * \begin{align} + * \label{eq:map} + * \bar{x} &= x + W(\vec{x})\left(\tan(F_y(t))(y-y_C) + + * \tan(F_z(t))(z-z_C)\right) \\ + * \bar{y} &= y \\ + * \bar{z} &= z + * \end{align} + * + * where $\vec{x}_C = (x_C, y_C, z_C)$ is the \p center of the skew map (which + * is different than the origin of the coordinate system), and $F_y(t)$ and + * $F_z(t)$ are the angles within the $(x,y)$ plane between the undistorted + * $x$-axis and the skewed $\bar{y}$-axis at the origin, 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 $\vec{x}_C$ (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`. The reason + * that $W(\vec{x})$ *should* be 1 at $\vec{x}_C$, and doesn't *need* to be 1 at + * $\vec{x}_C$ is because having $W(\vec{x})$ centered at the origin is much + * better for the smoothness of higher derivatives of $W(\vec{x})$ than if it + * were centered at $\vec{x}_C$. The important part is that the map is left + * invariant along the $x=x_C$ line and that $W(\vec{x}_C)\approx 1$ for the + * skew control system, both of which are satisfied if $W(\vec{x})$ is centered + * at the origin and $|\vec{x}_C|\ll R$. + * + * With that in mind, $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 origin, to 1 at $R$. We choose $\lambda(\vec{x})$ to quadratically go + * from 0 at the origin to 1 at $R$ with the form + * + * \begin{equation} + * \label{eq:lambda} + * \lambda(\vec{x}) = \frac{|\vec{x}|^2}{R^2}. + * \end{equation} + * + * A quadratic form was chosen for $\lambda$ rather than higher powers of 2 to + * avoid needing too much resolution to resolve a steep falloff in the function. + * + * \note If the quantity $S = \tan(F_y(t))(y-y_C) + \tan(F_z(t))(z-z_C)$ becomes + * too large, then the map becomes singular because multiple $x$ will be mapped + * to the same $\bar{x}$. This is due to the fact we are adding a linear term to + * a cosine term with different weights. If $S$ is too large, the cosine term + * will overpower the linear and the map will become singular. + * + * ### 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-y_C) + \tan(F_z(t))(z-z_C)\right) - + * \bar{x} + * \end{equation} + * + * We can bound the root by noticing that in $\ref{eq:map}$, if we substitute + * 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} + * x &<=& \bar{x} &<=& x + (\tan(F_y)(y-y_C) + \tan(F_z)(z-z_C)) \\ + * 0 &<=& \bar{x} - x &<=& (\tan(F_y)(y-y_C) + \tan(F_z)(z-z_C)) \\ + * -\bar{x} &<=& - x &<=& -\bar{x} + (\tan(F_y)(y-y_C) + \tan(F_z)(z-z_C)) \\ + * \bar{x} &>=& x &>=& \bar{x} - (\tan(F_y)(y-y_C) + \tan(F_z)(z-z_C)) + * \end{align} + * + * where on each line we just made simple arithmetic operations. We pad each + * bound by $10^{-14}$ just to avoid roundoff issues. If either of the bounds is + * within roundoff of zero, the map is the identity 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-y_C) + + * \dot{F}_z(t)(1 + \tan^2(F_z(t))(z-z_C))\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-y_C) + + * \tan(F_z(t))(z-z_C)\right), \\ + * {W^0}_1 &= \frac{\partial(\bar{x}-x)}{\partial y} &= \frac{\partial + * W(\vec{x})}{\partial y}&\left(\tan(F_y(t))(y-y_C) + + * \tan(F_z(t))(z-z_C)\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-y_C) + + * \tan(F_z(t))(z-z_C)\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}$) is given by + * + * \begin{equation} + * \frac{\partial \lambda(\vec{x})}{\partial x^i} = \frac{2x^i}{R^2} + * \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_width( + const std::array& source_coords) const; + + template + std::array, 3> get_width_deriv( + const std::array& source_coords) 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; + + template + void check_for_singular_map( + const std::array& source_coords, double time, + const FunctionsOfTimeMap& functions_of_time) 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_{}; + 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..edc58222d56a --- /dev/null +++ b/tests/Unit/Domain/CoordinateMaps/TimeDependent/Test_Skew.cpp @@ -0,0 +1,245 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Framework/TestingFramework.hpp" + +#include +#include +#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 { +constexpr size_t deriv_order = 2; +using Polynomial = domain::FunctionsOfTime::PiecewisePolynomial; +using FoftPtr = std::unique_ptr; + +template +void test(const gsl::not_null generator) { + 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 + + 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 double outer_radius = outer_radius_dist(*generator); + const std::array center{ + point_dist(*generator), point_dist(*generator), point_dist(*generator)}; + 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); +} + +void test_specific_points() { + const std::string function_of_time_name{"Skew"}; + const double time = 0.0; + std::unordered_map f_of_t_list{}; + // Use pi/4 so math is easy + f_of_t_list[function_of_time_name] = std::make_unique( + time, + std::array{DataVector{2, M_PI_4}, DataVector{2, 0.0}, DataVector{2, 0.0}}, + std::numeric_limits::infinity()); + + // Make it unit so math is easy + const std::array center{0.0, -1.0 / sqrt(2.0), 1.0 / sqrt(2.0)}; + const double outer_radius = 100.0; + + const CoordinateMaps::TimeDependent::Skew skew_map{function_of_time_name, + center, outer_radius}; + + std::array test_point{}; + { + INFO("Center"); + test_point = center; + CAPTURE(test_point); + + const auto mapped_point = skew_map(test_point, time, f_of_t_list); + // Should be exact + CHECK(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); + auto expected_jacobian = identity<3>(0.0); + // Because the angles are pi/4 + get<0, 1>(expected_jacobian) = + 0.5 * (1.0 + cos(M_PI / square(outer_radius))); + get<0, 2>(expected_jacobian) = get<0, 1>(expected_jacobian); + CHECK_ITERABLE_APPROX(jacobian, expected_jacobian); + } + { + INFO("Half way along y-axis"); + test_point = std::array{0.0, 0.5 * outer_radius, 0.0}; + CAPTURE(test_point); + + const auto mapped_point = skew_map(test_point, time, f_of_t_list); + const double falloff = 0.25 * (2.0 + sqrt(2.0)); + const double tan_sum = test_point[1]; + // Should be exact + CHECK_ITERABLE_APPROX( + mapped_point, + (std::array{falloff * tan_sum, test_point[1], test_point[2]})); + 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); + auto expected_jacobian = identity<3>(0.0); + get<0, 2>(expected_jacobian) = falloff; + get<0, 1>(expected_jacobian) = + -0.5 * M_PI * tan_sum / (sqrt(2.0) * outer_radius) + + get<0, 2>(expected_jacobian); + CHECK_ITERABLE_APPROX(jacobian, expected_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::epsilon(); + test_point = std::array{0.0, outer_radius + eps, 0.0}; + 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)); + } +} + +void test_errors() { +#ifdef SPECTRE_DEBUG + const std::string function_of_time_name{"Skew"}; + const double time = 0.0; + std::unordered_map f_of_t_list{}; + // Use values close to PI/2 to trigger the error + f_of_t_list[function_of_time_name] = std::make_unique( + time, + std::array{DataVector{2, M_PI_2 - 1.e-3}, DataVector{2, 0.0}, + DataVector{2, 0.0}}, + std::numeric_limits::infinity()); + + const double outer_radius = 100.0; + + const CoordinateMaps::TimeDependent::Skew skew_map{ + function_of_time_name, std::array{0.0, 0.0, 0.0}, outer_radius}; + + const std::array point{50.0, 20.0, 30.0}; + + CHECK_THROWS_WITH((skew_map(point, time, f_of_t_list)), + Catch::Matchers::ContainsSubstring("Skew map is singular")); +#endif +} +} // namespace + +SPECTRE_TEST_CASE("Unit.Domain.CoordinateMaps.TimeDependent.Skew", + "[Domain][Unit]") { + MAKE_GENERATOR(generator); + test(make_not_null(&generator)); + test_specific_points(); + test_errors(); +} +} // namespace domain