Skip to content

Commit

Permalink
Add translation map to sphere
Browse files Browse the repository at this point in the history 
Add rotation and expansion options

Add test for translation time dependence in Sphere

Fix RotScaleTrans initial values and other fixes

Fix

Fix
  • Loading branch information
@guilara
guilara committed Feb 26, 2024
1 parent 6dc9765 commit 399ec27
Show file tree
Hide file tree
Showing 7 changed files with 544 additions and 146 deletions.
30 changes: 25 additions & 5 deletions src/Domain/Creators/Sphere.cpp
View file
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
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

0 comments on commit 399ec27

Please sign in to comment.