Skip to content

Commit

Permalink
Merge pull request #6324 from nilsdeppe/bns_improvements_4
Browse files Browse the repository at this point in the history 
Support equiangular grid points for BNS domain
  • Loading branch information
@nilsdeppe
nilsdeppe authored Nov 27, 2024
2 parents 423ac2b + b5e1fef commit 08a4104
Show file tree
Hide file tree
Showing 9 changed files with 178 additions and 102 deletions.
155 changes: 98 additions & 57 deletions src/Domain/CoordinateMaps/Frustum.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,23 +30,25 @@ namespace domain::CoordinateMaps {
Frustum::Frustum(const std::array<std::array<double, 2>, 4>& face_vertices,
const double lower_bound, const double upper_bound,
OrientationMap<3> orientation_of_frustum,
const bool with_equiangular_map,
bool equiangular_map_at_outer, bool equiangular_map_at_inner,
const Distribution zeta_distribution,
const std::optional<double> distribution_value,
const double sphericity, const double transition_phi,
const double opening_angle)
// clang-tidy: trivially copyable
: orientation_of_frustum_(std::move(orientation_of_frustum)), // NOLINT
with_equiangular_map_(with_equiangular_map),
is_identity_(face_vertices ==
std::array<std::array<double, 2>, 4>{{{{-1.0, -1.0}},
{{1.0, 1.0}},
{{-1.0, -1.0}},
{{1.0, 1.0}}}} and
lower_bound == -1.0 and upper_bound == 1.0 and
orientation_of_frustum_.is_aligned() and
not with_equiangular_map and
zeta_distribution == Distribution::Linear),
equiangular_map_at_outer_{equiangular_map_at_outer},
equiangular_map_at_inner_{equiangular_map_at_inner},
is_identity_(
face_vertices ==
std::array<std::array<double, 2>, 4>{{{{-1.0, -1.0}},
{{1.0, 1.0}},
{{-1.0, -1.0}},
{{1.0, 1.0}}}} and
lower_bound == -1.0 and upper_bound == 1.0 and
orientation_of_frustum_.is_aligned() and
not(equiangular_map_at_inner_ or equiangular_map_at_outer_) and
zeta_distribution == Distribution::Linear),
zeta_distribution_(zeta_distribution),
sphericity_(sphericity) {
ASSERT(sphericity_ >= 0.0 and sphericity_ <= 1.0,
Expand Down Expand Up @@ -169,33 +171,41 @@ std::array<tt::remove_cvref_wrap_t<T>, 3> Frustum::operator()(
"Only the distributions Linear, Projective, and Logarithmic are "
"supported.");
}
const ReturnType cap_xi_zero = with_equiangular_map_ ? tan(M_PI_4 * xi) : xi;
const ReturnType& cap_xi_zero =
equiangular_map_at_inner_ ? tan(M_PI_4 * xi) : xi;
const double one_plus_phi_square = 1.0 + phi_ * phi_;
const ReturnType cap_xi_upper =
with_equiangular_map_
const ReturnType& cap_xi_upper =
equiangular_map_at_outer_
? one_plus_phi_square * one_over_tan_half_opening_angle_ *
tan(half_opening_angle_ * (xi + phi_) /
one_plus_phi_square) -
phi_
: xi;
const ReturnType cap_xi_transition = 0.5 * (1.0 + cap_zeta) * cap_xi_upper +
0.5 * (1.0 - cap_zeta) * cap_xi_zero;
const ReturnType cap_eta = with_equiangular_map_ ? tan(M_PI_4 * eta) : eta;
const ReturnType& cap_xi_transition = 0.5 * (1.0 + cap_zeta) * cap_xi_upper +
0.5 * (1.0 - cap_zeta) * cap_xi_zero;

const ReturnType& cap_eta_zero =
equiangular_map_at_inner_ ? tan(M_PI_4 * eta) : eta;
const ReturnType& cap_eta_upper =
equiangular_map_at_outer_ ? tan(M_PI_4 * eta) : eta;
const ReturnType& cap_eta_transition =
0.5 * (1.0 + cap_zeta) * cap_eta_upper +
0.5 * (1.0 - cap_zeta) * cap_eta_zero;

ReturnType physical_x =
sigma_x_ + delta_x_xi_ * cap_xi_transition +
(delta_x_zeta_ + delta_x_xi_zeta_ * cap_xi_transition) * cap_zeta;
ReturnType physical_y =
sigma_y_ + delta_y_eta_ * cap_eta +
(delta_y_zeta_ + delta_y_eta_zeta_ * cap_eta) * cap_zeta;
sigma_y_ + delta_y_eta_ * cap_eta_transition +
(delta_y_zeta_ + delta_y_eta_zeta_ * cap_eta_transition) * cap_zeta;
ReturnType physical_z = sigma_z_ + delta_z_zeta_ * cap_zeta;
if (sphericity_ > 0.0) {
const ReturnType upper_surface_x =
sigma_x_ + delta_x_xi_ * cap_xi_upper +
(delta_x_zeta_ + delta_x_xi_zeta_ * cap_xi_upper);
const ReturnType upper_surface_y =
sigma_y_ + delta_y_eta_ * cap_eta +
(delta_y_zeta_ + delta_y_eta_zeta_ * cap_eta);
sigma_y_ + delta_y_eta_ * cap_eta_upper +
(delta_y_zeta_ + delta_y_eta_zeta_ * cap_eta_upper);
const double upper_surface_z = sigma_z_ + delta_z_zeta_;
const ReturnType upper_surface_r =
sqrt(square(upper_surface_x) + square(upper_surface_y) +
Expand Down Expand Up @@ -236,7 +246,7 @@ std::optional<std::array<double, 3>> Frustum::inverse(
logical_coords[1] =
(physical_coords[1] - sigma_y_ - delta_y_zeta_ * logical_coords[2]) /
denom1;
if (with_equiangular_map_) {
if (equiangular_map_at_inner_ or equiangular_map_at_outer_) {
logical_coords[0] = atan(logical_coords[0]) / M_PI_4;
logical_coords[1] = atan(logical_coords[1]) / M_PI_4;
}
Expand Down Expand Up @@ -369,39 +379,57 @@ tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> Frustum::jacobian(
"Only the distributions Linear, Projective, and Logarithmic are "
"supported.");
}
const ReturnType& cap_xi_zero = with_equiangular_map_ ? tan(M_PI_4 * xi) : xi;
const ReturnType& cap_xi_zero =
equiangular_map_at_inner_ ? tan(M_PI_4 * xi) : xi;
const double one_plus_phi_square = 1.0 + phi_ * phi_;
const ReturnType cap_xi_upper =
with_equiangular_map_
const ReturnType& cap_xi_upper =
equiangular_map_at_outer_
? one_plus_phi_square * one_over_tan_half_opening_angle_ *
tan(half_opening_angle_ * (xi + phi_) /
one_plus_phi_square) -
phi_
: xi;
const ReturnType cap_xi_transition =
with_equiangular_map_ ? 0.5 * (1.0 + cap_zeta) * cap_xi_upper +
0.5 * (1.0 - cap_zeta) * cap_xi_zero
: xi;

const ReturnType& cap_eta = with_equiangular_map_ ? tan(M_PI_4 * eta) : eta;
const ReturnType cap_xi_zero_deriv =
with_equiangular_map_ ? M_PI_4 * (1.0 + square(cap_xi_zero))
: make_with_value<ReturnType>(xi, 1.0);
const ReturnType cap_xi_upper_deriv =
with_equiangular_map_
const ReturnType& cap_xi_transition =
(equiangular_map_at_outer_ or equiangular_map_at_inner_)
? 0.5 * (1.0 + cap_zeta) * cap_xi_upper +
0.5 * (1.0 - cap_zeta) * cap_xi_zero
: xi;

const ReturnType& cap_eta_zero =
equiangular_map_at_inner_ ? tan(M_PI_4 * eta) : eta;
const ReturnType& cap_eta_upper =
equiangular_map_at_outer_ ? tan(M_PI_4 * eta) : eta;
const ReturnType& cap_eta_transition =
0.5 * (1.0 + cap_zeta) * cap_eta_upper +
0.5 * (1.0 - cap_zeta) * cap_eta_zero;

const ReturnType& cap_xi_zero_deriv =
equiangular_map_at_inner_ ? M_PI_4 * (1.0 + square(cap_xi_zero))
: make_with_value<ReturnType>(xi, 1.0);
const ReturnType& cap_xi_upper_deriv =
equiangular_map_at_outer_
? one_over_tan_half_opening_angle_ * half_opening_angle_ *
(1.0 + square(tan(half_opening_angle_ * (xi + phi_) /
one_plus_phi_square)))
: make_with_value<ReturnType>(xi, 1.0);

const ReturnType cap_xi_transition_deriv =
with_equiangular_map_ ? 0.5 * (1.0 + cap_zeta) * cap_xi_upper_deriv +
0.5 * (1.0 - cap_zeta) * cap_xi_zero_deriv
: make_with_value<ReturnType>(xi, 1.0);
const ReturnType& cap_xi_transition_deriv =
(equiangular_map_at_inner_ or equiangular_map_at_outer_)
? 0.5 * (1.0 + cap_zeta) * cap_xi_upper_deriv +
0.5 * (1.0 - cap_zeta) * cap_xi_zero_deriv
: make_with_value<ReturnType>(xi, 1.0);

const ReturnType cap_eta_deriv = with_equiangular_map_
? M_PI_4 * (1.0 + square(cap_eta))
: make_with_value<ReturnType>(eta, 1.0);
const ReturnType& cap_eta_zero_deriv =
equiangular_map_at_inner_ ? M_PI_4 * (1.0 + square(cap_eta_zero))
: make_with_value<ReturnType>(eta, 1.0);
const ReturnType& cap_eta_upper_deriv =
equiangular_map_at_outer_ ? M_PI_4 * (1.0 + square(cap_eta_upper))
: make_with_value<ReturnType>(eta, 1.0);
const ReturnType& cap_eta_transition_deriv =
(equiangular_map_at_inner_ or equiangular_map_at_outer_)
? 0.5 * (1.0 + cap_zeta) * cap_eta_upper_deriv +
0.5 * (1.0 - cap_zeta) * cap_eta_zero_deriv
: make_with_value<ReturnType>(eta, 1.0);

auto jacobian_matrix =
make_with_value<tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame>>(
Expand All @@ -413,7 +441,7 @@ tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> Frustum::jacobian(
const size_t mapped_dim_0 = orientation_of_frustum_.inverse_map()(0);
jacobian_matrix.get(mapped_dim_0, 0) =
delta_x_xi_ + delta_x_xi_zeta_ * cap_zeta;
if (with_equiangular_map_) {
if (equiangular_map_at_inner_ or equiangular_map_at_outer_) {
jacobian_matrix.get(mapped_dim_0, 0) *= cap_xi_transition_deriv;
}
if (mapped_xi.side() == Side::Lower) {
Expand All @@ -426,8 +454,8 @@ tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> Frustum::jacobian(
const size_t mapped_dim_1 = orientation_of_frustum_.inverse_map()(1);
jacobian_matrix.get(mapped_dim_1, 1) =
delta_y_eta_ + delta_y_eta_zeta_ * cap_zeta;
if (with_equiangular_map_) {
jacobian_matrix.get(mapped_dim_1, 1) *= cap_eta_deriv;
if (equiangular_map_at_inner_ or equiangular_map_at_outer_) {
jacobian_matrix.get(mapped_dim_1, 1) *= cap_eta_transition_deriv;
}
if (mapped_eta.side() == Side::Lower) {
jacobian_matrix.get(mapped_dim_1, 1) *= -1.0;
Expand All @@ -440,7 +468,9 @@ tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> Frustum::jacobian(
{delta_x_zeta_ + delta_x_xi_zeta_ * cap_xi_transition +
(delta_x_xi_ + delta_x_xi_zeta_ * cap_zeta) * 0.5 *
(cap_xi_upper - cap_xi_zero),
delta_y_zeta_ + delta_y_eta_zeta_ * cap_eta,
delta_y_zeta_ + delta_y_eta_zeta_ * cap_eta_transition +
(delta_y_eta_ + delta_y_eta_zeta_ * cap_zeta) * 0.5 *
(cap_eta_upper - cap_eta_zero),
make_with_value<ReturnType>(zeta, delta_z_zeta_)}});

if (zeta_distribution_ != Distribution::Linear) {
Expand All @@ -458,9 +488,9 @@ tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> Frustum::jacobian(
delta_x_zeta_ +
delta_x_xi_zeta_ * cap_xi_upper;

const ReturnType flat_frustum_y = sigma_y_ + delta_y_eta_ * cap_eta +
const ReturnType flat_frustum_y = sigma_y_ + delta_y_eta_ * cap_eta_upper +
delta_y_zeta_ +
delta_y_eta_zeta_ * cap_eta;
delta_y_eta_zeta_ * cap_eta_upper;

const double flat_frustum_z = sigma_z_ + delta_z_zeta_;

Expand Down Expand Up @@ -489,7 +519,7 @@ tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> Frustum::jacobian(
-1.0 * r_over_mag_flat * flat_frustum_z_hat *
flat_frustum_x_hat}});

if (with_equiangular_map_) {
if (equiangular_map_at_inner_ or equiangular_map_at_outer_) {
delta_dX_dxi[0] *= cap_xi_upper_deriv;
delta_dX_dxi[1] *= cap_xi_upper_deriv;
delta_dX_dxi[2] *= cap_xi_upper_deriv;
Expand All @@ -510,10 +540,10 @@ tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> Frustum::jacobian(
-1.0 * r_over_mag_flat * flat_frustum_z_hat *
flat_frustum_y_hat}});

if (with_equiangular_map_) {
delta_dX_deta[0] *= cap_eta_deriv;
delta_dX_deta[1] *= cap_eta_deriv;
delta_dX_deta[2] *= cap_eta_deriv;
if (equiangular_map_at_inner_ or equiangular_map_at_outer_) {
delta_dX_deta[0] *= cap_eta_upper_deriv;
delta_dX_deta[1] *= cap_eta_upper_deriv;
delta_dX_deta[2] *= cap_eta_upper_deriv;
}

get<0, 1>(jacobian_matrix) += delta_dX_deta[0];
Expand Down Expand Up @@ -549,14 +579,24 @@ tnsr::Ij<tt::remove_cvref_wrap_t<T>, 3, Frame::NoFrame> Frustum::inv_jacobian(
}

void Frustum::pup(PUP::er& p) {
size_t version = 4;
size_t version = 5;
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 | orientation_of_frustum_;
p | with_equiangular_map_;
if (version < 5) {
bool with_equiangular_map{false};
p | with_equiangular_map;
if (p.isUnpacking()) {
equiangular_map_at_inner_ = with_equiangular_map;
equiangular_map_at_outer_ = with_equiangular_map;
}
} else {
p | equiangular_map_at_inner_;
p | equiangular_map_at_outer_;
}
p | is_identity_;
if (version < 2 /*is unpacking*/) {
bool with_projective_map = false; // unused
Expand Down Expand Up @@ -628,7 +668,8 @@ bool operator==(const Frustum& lhs, const Frustum& rhs) {
// above.
// `radius` is also not compared for similar reasons.
return lhs.orientation_of_frustum_ == rhs.orientation_of_frustum_ and
lhs.with_equiangular_map_ == rhs.with_equiangular_map_ and
lhs.equiangular_map_at_inner_ == rhs.equiangular_map_at_inner_ and
lhs.equiangular_map_at_outer_ == rhs.equiangular_map_at_outer_ and
lhs.is_identity_ == rhs.is_identity_ and
lhs.zeta_distribution_ == rhs.zeta_distribution_ and
lhs.sigma_x_ == rhs.sigma_x_ and
Expand Down
16 changes: 11 additions & 5 deletions src/Domain/CoordinateMaps/Frustum.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -253,7 +253,7 @@ namespace CoordinateMaps {
* \phi^2)\tan(\frac{\pi}{4}\frac{\xi+\phi}{1+\phi^2})-\phi\f]
*
* For \f$\phi=0\f$ this generalized map simplifies to the original equiangular
* map \f$Xi_0 = \tan(\frac{\pi}{4}\xi)\f$.
* map \f$\Xi_0 = \tan(\frac{\pi}{4}\xi)\f$.
*
* This function is compatible with the ability to bulge out the Frustum; the
* map and Jacobian remain unchanged, modulo this substitution.
Expand All @@ -270,8 +270,12 @@ class Frustum {
* \param upper_bound z distance from the origin to the upper base of the
* frustum.
* \param orientation_of_frustum The orientation of the frustum in 3D space.
* \param with_equiangular_map Determines whether to apply a tangent function
* mapping to the logical coordinates (true) or not (false).
* \param equiangular_map_at_outer Determines whether to apply a tangent
* function mapping to the logical coordinates (true) or not (false) at the
* larger side of the frustum.
* \param equiangular_map_at_inner Determines whether to apply a tangent
* function mapping to the logical coordinates (true) or not (false) at the
* smaller side of the frustum.
* \param zeta_distribution Whether to apply a linear, logarithmic, or
* projective mapping to the logical zeta coordinate.
* \param distribution_value In the case of a projective map, the projective
Expand All @@ -288,7 +292,8 @@ class Frustum {
Frustum(const std::array<std::array<double, 2>, 4>& face_vertices,
double lower_bound, double upper_bound,
OrientationMap<3> orientation_of_frustum,
bool with_equiangular_map = false,
bool equiangular_map_at_outer = false,
bool equiangular_map_at_inner = false,
Distribution zeta_distribution = Distribution::Linear,
std::optional<double> distribution_value = std::nullopt,
double sphericity = 0.0, double transition_phi = 0.0,
Expand Down Expand Up @@ -333,7 +338,8 @@ class Frustum {

OrientationMap<3> orientation_of_frustum_ =
OrientationMap<3>::create_aligned();
bool with_equiangular_map_{false};
bool equiangular_map_at_outer_{false};
bool equiangular_map_at_inner_{false};
bool is_identity_{false};
Distribution zeta_distribution_ = Distribution::Linear;
std::optional<double> zeta_distribution_value_{std::nullopt};
Expand Down
2 changes: 2 additions & 0 deletions src/Domain/Creators/BinaryCompactObject.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -530,6 +530,8 @@ Domain<3> BinaryCompactObject<UseWorldtube>::create_domain() const {
Maps maps_frustums = domain::make_vector_coordinate_map_base<
Frame::BlockLogical, Frame::Inertial, 3>(frustum_coordinate_maps(
length_inner_cube_, length_outer_cube_, use_equiangular_map_,
use_equiangular_map_ and not use_single_block_a_ and
not use_single_block_b_,
{{-translation_, -center_of_mass_offset_[0], -center_of_mass_offset_[1]}},
radial_distribution_envelope_,
radial_distribution_envelope_ ==
Expand Down
12 changes: 7 additions & 5 deletions src/Domain/Creators/FrustalCloak.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,11 +74,13 @@ FrustalCloak::FrustalCloak(
Domain<3> FrustalCloak::create_domain() const {
std::vector<std::unique_ptr<
CoordinateMapBase<Frame::BlockLogical, Frame::Inertial, 3>>>
coord_maps = domain::make_vector_coordinate_map_base<
Frame::BlockLogical, Frame::Inertial, 3>(frustum_coordinate_maps(
length_inner_cube_, length_outer_cube_, use_equiangular_map_,
origin_preimage_, domain::CoordinateMaps::Distribution::Projective,
projection_factor_));
coord_maps = domain::make_vector_coordinate_map_base<Frame::BlockLogical,
Frame::Inertial, 3>(
frustum_coordinate_maps(
length_inner_cube_, length_outer_cube_, use_equiangular_map_,
use_equiangular_map_, origin_preimage_,
domain::CoordinateMaps::Distribution::Projective,
projection_factor_));
return Domain<3>{std::move(coord_maps),
corners_for_biradially_layered_domains(
0, 1, false, false, {{1, 2, 3, 4, 5, 6, 7, 8}})};
Expand Down
Loading

0 comments on commit 08a4104

Please sign in to comment.