diff --git a/src/Domain/Creators/BinaryCompactObject.cpp b/src/Domain/Creators/BinaryCompactObject.cpp index 690c5f2be20f..edbb52ff2eb1 100644 --- a/src/Domain/Creators/BinaryCompactObject.cpp +++ b/src/Domain/Creators/BinaryCompactObject.cpp @@ -564,12 +564,14 @@ Domain<3> BinaryCompactObject::create_domain() const { // Inside the excision sphere we add the grid to inertial map from the outer // shell. This allows the center of the excisions/horizons to be mapped // properly to the inertial frame. - if (is_excised_a_) { + if (is_excised_a_ and + grid_to_inertial_block_maps[number_of_blocks_ - 1] != nullptr) { domain.inject_time_dependent_map_for_excision_sphere( "ExcisionSphereA", grid_to_inertial_block_maps[number_of_blocks_ - 1]->get_clone()); } - if (is_excised_b_) { + if (is_excised_b_ and + grid_to_inertial_block_maps[number_of_blocks_ - 1] != nullptr) { domain.inject_time_dependent_map_for_excision_sphere( "ExcisionSphereB", grid_to_inertial_block_maps[number_of_blocks_ - 1]->get_clone()); @@ -620,7 +622,8 @@ Domain<3> BinaryCompactObject::create_domain() const { time_dependent_options_ ->distorted_to_inertial_map( block_for_distorted_frame); - } else { + } else if (grid_to_inertial_block_maps[number_of_blocks_ - 1] != + nullptr) { // No distorted frame grid_to_inertial_block_maps[block] = grid_to_inertial_block_maps[number_of_blocks_ - 1]->get_clone(); @@ -628,6 +631,9 @@ Domain<3> BinaryCompactObject::create_domain() const { } // Finally, inject the time dependent maps into the corresponding blocks for (size_t block = 0; block < number_of_blocks_; ++block) { + if (grid_to_inertial_block_maps[block] == nullptr) { + continue; + } domain.inject_time_dependent_map_for_block( block, std::move(grid_to_inertial_block_maps[block]), std::move(grid_to_distorted_block_maps[block]), diff --git a/src/Domain/Creators/BinaryCompactObjectHelpers.cpp b/src/Domain/Creators/BinaryCompactObjectHelpers.cpp index 34d06ed29c29..df254130d42a 100644 --- a/src/Domain/Creators/BinaryCompactObjectHelpers.cpp +++ b/src/Domain/Creators/BinaryCompactObjectHelpers.cpp @@ -22,6 +22,7 @@ #include "Domain/FunctionsOfTime/QuaternionFunctionOfTime.hpp" #include "NumericalAlgorithms/SphericalHarmonics/Spherepack.hpp" #include "Options/ParseError.hpp" +#include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrHorizon.hpp" #include "Utilities/ErrorHandling/Assert.hpp" #include "Utilities/ErrorHandling/Error.hpp" #include "Utilities/GenerateInstantiations.hpp" @@ -145,34 +146,75 @@ TimeDependentMapOptions::create_functions_of_time( } // Size and Shape FunctionOfTime for objects A and B - for (size_t i = 0; i < shape_names.size(); i++) { - if (i == 0 ? shape_options_A_.has_value() : shape_options_B_.has_value()) { - const auto make_initial_size_values = [](const auto& lambda_options) { - return std::array{ - {{gsl::at(lambda_options.value().initial_size_values, 0)}, - {gsl::at(lambda_options.value().initial_size_values, 1)}, - {gsl::at(lambda_options.value().initial_size_values, 2)}, - {0.0}}}; - }; + const auto build_shape_and_size_fot = [&result, &expiration_times, this]( + const auto& shape_options, + const double inner_radius, + const std::string& shape_name, + const std::string& size_name) { + const DataVector shape_zeros{ylm::Spherepack::spectral_size( + shape_options.l_max, shape_options.l_max), + 0.0}; + DataVector shape_func{}; + DataVector size_func{1, shape_options.initial_size_values[0]}; + if (shape_options.initial_values.has_value()) { + if (std::holds_alternative( + shape_options.initial_values.value())) { + const ylm::Spherepack ylm{shape_options.l_max, shape_options.l_max}; + const auto& mass_and_spin = + std::get( + shape_options.initial_values.value()); + const DataVector radial_distortion = + inner_radius - + get(gr::Solutions::kerr_schild_radius_from_boyer_lindquist( + inner_radius, ylm.theta_phi_points(), mass_and_spin.mass, + mass_and_spin.spin)); + shape_func = ylm.phys_to_spec(radial_distortion); + // Transform from SPHEREPACK to actual Ylm for size func + if (size_func[0] != 0.0) { + ERROR( + "Initial value for size map must be zero, because it is " + "overridden by the initial shape map values."); + } + size_func[0] = shape_func[0] * sqrt(0.5 * M_PI); + // Set l=0 for shape map to 0 because size is going to be used + shape_func[0] = 0.0; + } + } else { + shape_func = shape_zeros; + } - const size_t initial_l_max = i == 0 ? shape_options_A_.value().l_max - : shape_options_B_.value().l_max; - const std::array initial_size_values = - i == 0 ? make_initial_size_values(shape_options_A_) - : make_initial_size_values(shape_options_B_); - const DataVector shape_zeros{ - ylm::Spherepack::spectral_size(initial_l_max, initial_l_max), 0.0}; - - result[gsl::at(shape_names, i)] = - std::make_unique>( - initial_time_, - std::array{shape_zeros, shape_zeros, shape_zeros}, - expiration_times.at(gsl::at(shape_names, i))); - result[gsl::at(size_names, i)] = - std::make_unique>( - initial_time_, initial_size_values, - expiration_times.at(gsl::at(size_names, i))); + result[shape_name] = + std::make_unique>( + initial_time_, + std::array{std::move(shape_func), shape_zeros, + shape_zeros}, + expiration_times.at(shape_name)); + result[size_name] = + std::make_unique>( + initial_time_, + std::array{{std::move(size_func), + {shape_options.initial_size_values[1]}, + {shape_options.initial_size_values[2]}, + {0.0}}}, + expiration_times.at(size_name)); + }; + if (shape_options_A_.has_value()) { + if (not inner_radii_[0].has_value()) { + ERROR( + "A shape map was specified for object A, but no inner radius is " + "available. The object must be enclosed by a sphere."); } + build_shape_and_size_fot(shape_options_A_.value(), *inner_radii_[0], + shape_names[0], size_names[0]); + } + if (shape_options_B_.has_value()) { + if (not inner_radii_[1].has_value()) { + ERROR( + "A shape map was specified for object B, but no inner radius is " + "available. The object must be enclosed by a sphere."); + } + build_shape_and_size_fot(shape_options_B_.value(), *inner_radii_[1], + shape_names[1], size_names[1]); } return result; @@ -206,6 +248,8 @@ void TimeDependentMapOptions::build_maps( << ", but no time dependent map options were specified " "for that object."); } + // Store the inner radii for creating functions of time + gsl::at(inner_radii_, i) = radii[0]; const size_t initial_l_max = i == 0 ? shape_options_A_.value().l_max : shape_options_B_.value().l_max; @@ -420,13 +464,7 @@ TimeDependentMapOptions::grid_to_inertial_map( } else if (rotation_map_.has_value()) { return std::make_unique>(rotation_map_.value()); } else { - ERROR( - "Requesting grid to inertial map without a distorted frame and " - "without a Rotation or Expansion map for object " - << Object - << ". This means there are no time dependent maps. If you don't want " - "time dependent maps, specify 'None' for TimeDependentMapOptions. " - "Otherwise specify at least one time dependent map."); + return nullptr; } } } diff --git a/src/Domain/Creators/BinaryCompactObjectHelpers.hpp b/src/Domain/Creators/BinaryCompactObjectHelpers.hpp index 91d199d02ad7..19a8cedad41a 100644 --- a/src/Domain/Creators/BinaryCompactObjectHelpers.hpp +++ b/src/Domain/Creators/BinaryCompactObjectHelpers.hpp @@ -16,6 +16,7 @@ #include "Domain/CoordinateMaps/TimeDependent/Rotation.hpp" #include "Domain/CoordinateMaps/TimeDependent/Shape.hpp" #include "Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/ShapeMapTransitionFunction.hpp" +#include "Domain/Creators/SphereTimeDependentMaps.hpp" #include "Domain/FunctionsOfTime/FunctionOfTime.hpp" #include "Domain/Structure/ObjectLabel.hpp" #include "Options/Auto.hpp" @@ -68,6 +69,15 @@ struct ShapeMapOptions { "zero."}; }; + struct InitialValues { + using type = + Options::Auto, + sphere::Spherical>; + static constexpr Options::String help = { + "Initial Ylm coefficients for the shape map. Specify 'Spherical' for " + "all coefficients to be initialized to zero."}; + }; + struct SizeInitialValues { using type = std::array; static constexpr Options::String help = { @@ -82,13 +92,15 @@ struct ShapeMapOptions { "be 0 at the outer radius of the inner sphere around the object"}; }; - using common_options = tmpl::list; + using common_options = tmpl::list; using options = tmpl::conditional_t< IsCylindrical, common_options, tmpl::push_back>; size_t l_max{}; + std::optional> + initial_values{}; std::array initial_size_values{}; bool transition_ends_at_cube{false}; }; @@ -394,6 +406,7 @@ struct TimeDependentMapOptions { std::optional rotation_options_{}; std::optional> shape_options_A_{}; std::optional> shape_options_B_{}; + std::array, 2> inner_radii_{}; // Maps std::optional expansion_map_{}; diff --git a/src/Domain/Creators/SphereTimeDependentMaps.cpp b/src/Domain/Creators/SphereTimeDependentMaps.cpp index 2fbdc34c3c14..2e82c28352c0 100644 --- a/src/Domain/Creators/SphereTimeDependentMaps.cpp +++ b/src/Domain/Creators/SphereTimeDependentMaps.cpp @@ -78,10 +78,10 @@ TimeDependentMapOptions::create_functions_of_time( const auto& mass_and_spin = std::get( 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, - mass_and_spin.spin)) / - inner_radius; + inner_radius - + get(gr::Solutions::kerr_schild_radius_from_boyer_lindquist( + inner_radius, ylm.theta_phi_points(), mass_and_spin.mass, + mass_and_spin.spin)); shape_func = ylm.phys_to_spec(radial_distortion); // Transform from SPHEREPACK to actual Ylm for size func size_func[0] = shape_func[0] * sqrt(0.5 * M_PI); diff --git a/src/Domain/Creators/TimeDependence/Shape.cpp b/src/Domain/Creators/TimeDependence/Shape.cpp index 004e7a5c1cd6..8500149d436f 100644 --- a/src/Domain/Creators/TimeDependence/Shape.cpp +++ b/src/Domain/Creators/TimeDependence/Shape.cpp @@ -178,8 +178,8 @@ Shape