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

Jump to bottom

Add rotation, expansion and translation maps to Sphere domain #5758

Merged
merged 1 commit into from
Feb 26, 2024
Merged

Add rotation, expansion and translation maps to Sphere domain #5758

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
30 changes: 25 additions & 5 deletions src/Domain/Creators/Sphere.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -267,11 +267,25 @@ Sphere::Sphere(
// Build the maps. We only apply the maps in the inner-most shell. The
// inner radius is what's passed in, but the outer radius is the outer
// radius of the inner-most shell so we have to see how many shells we
// have
// have.
// The transition region for the rotation, expansion and translation maps
// occurs only in the outer-most shell. We require at least one radial
// partition, as the transition requires spherical shape.

if (radial_partitioning_.empty()) {
PARSE_ERROR(context,
"The hard-coded translation map requires at least two "
"shells. Use at least one radial partition.");
}

std::get<sphere::TimeDependentMapOptions>(time_dependent_options_.value())
.build_maps(std::array{0.0, 0.0, 0.0}, inner_radius_,
radial_partitioning_.empty() ? outer_radius_
: radial_partitioning_[0]);
.build_maps(
std::array{0.0, 0.0, 0.0},
// Inner shell radii for the Shape map
std::pair<double, double>{inner_radius_, radial_partitioning_[0]},
// Outer shell radii for the transition region
std::pair<double, double>{radial_partitioning_.back(),
outer_radius_});
}
}
}
Expand Down Expand Up @@ -374,9 +388,15 @@ Domain<3> Sphere::create_domain() const {
time_dependent_options_.value());

// First shell gets the distorted maps.
// Last shell gets the transition region for the rotation, expansion and
// translation maps
for (size_t block_id = 0; block_id < num_blocks_; block_id++) {
const bool include_distorted_map_in_first_shell =
block_id < num_blocks_per_shell_;
// False if block_id is in the last shell
const bool use_rigid =
block_id + num_blocks_per_shell_ + (fill_interior_ ? 1 : 0) <
num_blocks_;
block_maps_grid_to_distorted[block_id] =
hard_coded_options.grid_to_distorted_map(
include_distorted_map_in_first_shell);
Expand All @@ -385,7 +405,7 @@ Domain<3> Sphere::create_domain() const {
include_distorted_map_in_first_shell);
block_maps_grid_to_inertial[block_id] =
hard_coded_options.grid_to_inertial_map(
include_distorted_map_in_first_shell);
include_distorted_map_in_first_shell, use_rigid);
}
} else {
const auto& time_dependence = std::get<std::unique_ptr<
Expand Down
179 changes: 153 additions & 26 deletions src/Domain/Creators/SphereTimeDependentMaps.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,16 +21,24 @@
#include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
#include "Domain/FunctionsOfTime/PiecewisePolynomial.hpp"
#include "Domain/FunctionsOfTime/QuaternionFunctionOfTime.hpp"
#include "Domain/FunctionsOfTime/SettleToConstant.hpp"
#include "Domain/FunctionsOfTime/SettleToConstantQuaternion.hpp"
#include "NumericalAlgorithms/SphericalHarmonics/Spherepack.hpp"
#include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrHorizon.hpp"
#include "Utilities/ErrorHandling/Error.hpp"

namespace domain::creators::sphere {

TimeDependentMapOptions::TimeDependentMapOptions(
const double initial_time, const ShapeMapOptions& shape_map_options)
const double initial_time, const ShapeMapOptions& shape_map_options,
std::optional<RotationMapOptions> rotation_map_options,
std::optional<ExpansionMapOptions> expansion_map_options,
std::optional<TranslationMapOptions> translation_map_options)
: initial_time_(initial_time),
initial_l_max_(shape_map_options.l_max),
initial_shape_values_(shape_map_options.initial_values) {}
shape_map_options_(shape_map_options),
rotation_map_options_(rotation_map_options),
expansion_map_options_(expansion_map_options),
translation_map_options_(translation_map_options) {}

std::unordered_map<std::string,
std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>
Expand All @@ -46,7 +54,8 @@ TimeDependentMapOptions::create_functions_of_time(
// their initial expiration time to infinity (i.e. not expiring)
std::unordered_map<std::string, double> expiration_times{
{size_name, std::numeric_limits<double>::infinity()},
{shape_name, std::numeric_limits<double>::infinity()}};
{shape_name, std::numeric_limits<double>::infinity()},
{translation_name, std::numeric_limits<double>::infinity()}};

// If we have control systems, overwrite these expiration times with the ones
// supplied by the control system
Expand All @@ -55,16 +64,19 @@ TimeDependentMapOptions::create_functions_of_time(
}

DataVector shape_zeros{
ylm::Spherepack::spectral_size(initial_l_max_, initial_l_max_), 0.0};
ylm::Spherepack::spectral_size(shape_map_options_.l_max,
shape_map_options_.l_max),
0.0};
DataVector shape_func{};
DataVector size_func{1, 0.0};

if (initial_shape_values_.has_value()) {
if (shape_map_options_.initial_values.has_value()) {
if (std::holds_alternative<KerrSchildFromBoyerLindquist>(
initial_shape_values_.value())) {
const ylm::Spherepack ylm{initial_l_max_, initial_l_max_};
const auto& mass_and_spin =
std::get<KerrSchildFromBoyerLindquist>(initial_shape_values_.value());
shape_map_options_.initial_values.value())) {
const ylm::Spherepack ylm{shape_map_options_.l_max,
shape_map_options_.l_max};
const auto& mass_and_spin = std::get<KerrSchildFromBoyerLindquist>(
shape_map_options_.initial_values.value());
const DataVector radial_distortion =
1.0 - get(gr::Solutions::kerr_schild_radius_from_boyer_lindquist(
inner_radius, ylm.theta_phi_points(), mass_and_spin.mass,
Expand Down Expand Up @@ -101,31 +113,143 @@ TimeDependentMapOptions::create_functions_of_time(
{0.0}}},
expiration_times.at(size_name));

// ExpansionMap FunctionOfTime
if (expansion_map_options_.has_value()) {
result[expansion_name] =
std::make_unique<FunctionsOfTime::SettleToConstant>(
std::array<DataVector, 3>{
{{gsl::at(expansion_map_options_.value().initial_values, 0)},
{gsl::at(expansion_map_options_.value().initial_values, 1)},
{gsl::at(expansion_map_options_.value().initial_values, 2)}}},
initial_time_, expansion_map_options_.value().decay_timescale);

// ExpansionMap in the Outer regionFunctionOfTime
result[expansion_outer_boundary_name] = std::make_unique<
FunctionsOfTime::SettleToConstant>(
std::array<DataVector, 3>{
{{gsl::at(
expansion_map_options_.value().initial_values_outer_boundary,
0)},
{gsl::at(
expansion_map_options_.value().initial_values_outer_boundary,
1)},
{gsl::at(
expansion_map_options_.value().initial_values_outer_boundary,
2)}}},
initial_time_,
expansion_map_options_.value().decay_timescale_outer_boundary);
}

DataVector initial_quaternion_value{4, 0.0};
DataVector initial_quaternion_first_derivative_value{4, 0.0};
DataVector initial_quaternion_second_derivative_value{4, 0.0};

// RotationMap FunctionOfTime
if (rotation_map_options_.has_value()) {
for (size_t i = 0; i < 3; i++) {
initial_quaternion_value[i] =
gsl::at(gsl::at(rotation_map_options_.value().initial_values, 0), i);
initial_quaternion_first_derivative_value[i] =
gsl::at(gsl::at(rotation_map_options_.value().initial_values, 1), i);
initial_quaternion_second_derivative_value[i] =
gsl::at(gsl::at(rotation_map_options_.value().initial_values, 2), i);
}
result[rotation_name] =
std::make_unique<FunctionsOfTime::SettleToConstantQuaternion>(
std::array<DataVector, 3>{
std::move(initial_quaternion_value),
std::move(initial_quaternion_first_derivative_value),
std::move(initial_quaternion_second_derivative_value)},
initial_time_, rotation_map_options_.value().decay_timescale);
}

DataVector initial_translation_center{3, 0.0};
DataVector initial_translation_velocity{3, 0.0};
DataVector initial_translation_acceleration{3, 0.0};

// Translation FunctionOfTime
if (translation_map_options_.has_value()) {
for (size_t i = 0; i < 3; i++) {
initial_translation_center[i] = gsl::at(
gsl::at(translation_map_options_.value().initial_values, 0), i);
initial_translation_velocity[i] = gsl::at(
gsl::at(translation_map_options_.value().initial_values, 1), i);
initial_translation_acceleration[i] = gsl::at(
gsl::at(translation_map_options_.value().initial_values, 2), i);
}
result[translation_name] =
std::make_unique<FunctionsOfTime::PiecewisePolynomial<2>>(
initial_time_,
std::array<DataVector, 3>{
{std::move(initial_translation_center),
std::move(initial_translation_velocity),
std::move(initial_translation_acceleration)}},
expiration_times.at(translation_name));
}

return result;
}

void TimeDependentMapOptions::build_maps(const std::array<double, 3>& center,
const double inner_radius,
const double outer_radius) {
void TimeDependentMapOptions::build_maps(
const std::array<double, 3>& center,
std::pair<double, double> inner_shell_radii,
std::pair<double, double> outer_shell_radii) {
std::unique_ptr<domain::CoordinateMaps::ShapeMapTransitionFunctions::
ShapeMapTransitionFunction>
transition_func =
std::make_unique<domain::CoordinateMaps::ShapeMapTransitionFunctions::
SphereTransition>(inner_radius, outer_radius);
shape_map_ = ShapeMap{center, initial_l_max_,
initial_l_max_, std::move(transition_func),
shape_name, size_name};
SphereTransition>(inner_shell_radii.first,
inner_shell_radii.second);
shape_map_ = ShapeMap{center,
shape_map_options_.l_max,
shape_map_options_.l_max,
std::move(transition_func),
shape_name,
size_name};

inner_rot_scale_trans_map_ = RotScaleTransMap{
expansion_map_options_.has_value()
? std::optional<std::pair<
std::string, std::string>>{{expansion_name,
expansion_outer_boundary_name}}
: std::nullopt,
rotation_map_options_.has_value()
? std::optional<std::string>{rotation_name}
: std::nullopt,
translation_map_options_.has_value()
? std::optional<std::string>{translation_name}
: std::nullopt,
outer_shell_radii.first,
outer_shell_radii.second,
domain::CoordinateMaps::TimeDependent::RotScaleTrans<
3>::BlockRegion::Inner};

transition_rot_scale_trans_map_ = RotScaleTransMap{
expansion_map_options_.has_value()
? std::optional<std::pair<
std::string, std::string>>{{expansion_name,
expansion_outer_boundary_name}}
: std::nullopt,
rotation_map_options_.has_value()
? std::optional<std::string>{rotation_name}
: std::nullopt,
translation_map_options_.has_value()
? std::optional<std::string>{translation_name}
: std::nullopt,
outer_shell_radii.first,
outer_shell_radii.second,
domain::CoordinateMaps::TimeDependent::RotScaleTrans<
3>::BlockRegion::Transition};
}

// If you edit any of the functions below, be sure to update the documentation
// in the Sphere domain creator as well as this class' documentation.
TimeDependentMapOptions::MapType<Frame::Distorted, Frame::Inertial>
TimeDependentMapOptions::distorted_to_inertial_map(
const bool include_distorted_map) {
const bool include_distorted_map) const {
if (include_distorted_map) {
return std::make_unique<
IdentityForComposition<Frame::Distorted, Frame::Inertial>>(
IdentityMap{});
return std::make_unique<DistortedToInertialComposition>(
inner_rot_scale_trans_map_);
} else {
return nullptr;
}
Expand All @@ -142,13 +266,16 @@ TimeDependentMapOptions::grid_to_distorted_map(
}

TimeDependentMapOptions::MapType<Frame::Grid, Frame::Inertial>
TimeDependentMapOptions::grid_to_inertial_map(
const bool include_distorted_map) const {
TimeDependentMapOptions::grid_to_inertial_map(const bool include_distorted_map,
const bool use_rigid) const {
if (include_distorted_map) {
return std::make_unique<GridToInertialComposition>(shape_map_);
return std::make_unique<GridToInertialComposition>(
shape_map_, inner_rot_scale_trans_map_);
} else if (use_rigid) {
return std::make_unique<GridToInertialSimple>(inner_rot_scale_trans_map_);
} else {
return std::make_unique<
IdentityForComposition<Frame::Grid, Frame::Inertial>>(IdentityMap{});
return std::make_unique<GridToInertialSimple>(
transition_rot_scale_trans_map_);
}
}
} // namespace domain::creators::sphere
Loading
Loading