Skip to content

Commit

Permalink
Use common time dep opts in BCO domains
Browse files Browse the repository at this point in the history
  • Loading branch information
knelli2 committed Sep 19, 2024
1 parent 32dd542 commit 9fbb81d
Show file tree
Hide file tree
Showing 11 changed files with 228 additions and 197 deletions.
108 changes: 55 additions & 53 deletions src/Domain/Creators/TimeDependentOptions/BinaryCompactObject.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -42,9 +42,9 @@ TimeDependentMapOptions<IsCylindrical>::TimeDependentMapOptions(
std::optional<ShapeMapOptions<domain::ObjectLabel::B>> shape_options_B,
const Options::Context& context)
: initial_time_(initial_time),
expansion_map_options_(expansion_map_options),
rotation_map_options_(rotation_map_options),
translation_map_options_(translation_map_options),
expansion_map_options_(std::move(expansion_map_options)),
rotation_map_options_(std::move(rotation_map_options)),
translation_map_options_(std::move(translation_map_options)),
shape_options_A_(shape_options_A),
shape_options_B_(shape_options_B) {
if (not(expansion_map_options_.has_value() or
Expand Down Expand Up @@ -82,37 +82,51 @@ TimeDependentMapOptions<IsCylindrical>::create_worldtube_functions_of_time()
std::unordered_map<std::string,
std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>
result{};

// The functions of time need to be valid only for the very first time step,
// after that they need to be updated by the worldtube singleton.
const double initial_expiration_time = initial_time_ + 1e-10;
if (not expansion_map_options_.has_value()) {
ERROR("Initial values for the expansion map need to be provided.");
}
if (not expansion_map_options_->asymptotic_velocity_outer_boundary
.has_value()) {
ERROR(
"The BinaryCompactObject domains don't support using a "
"SettleToConstant function of time for the expansion map. Set the "
"asymptotic velocity to a value and set the decay timescale to "
"'None'");
}
result[expansion_name] =
std::make_unique<FunctionsOfTime::IntegratedFunctionOfTime>(
initial_time_,
std::array<double, 2>{
{{gsl::at(expansion_map_options_.value().initial_values, 0)},
{gsl::at(expansion_map_options_.value().initial_values, 1)}}},
expansion_map_options_.value().initial_values[0][0],
expansion_map_options_.value().initial_values[1][0]},
initial_expiration_time, false);
result[expansion_outer_boundary_name] =
std::make_unique<FunctionsOfTime::FixedSpeedCubic>(
1.0, initial_time_,
expansion_map_options_.value().outer_boundary_velocity,
expansion_map_options_.value().outer_boundary_decay_time);
expansion_map_options_.value()
.asymptotic_velocity_outer_boundary.value(),
expansion_map_options_.value().decay_timescale_outer_boundary);

if (not rotation_map_options_.has_value()) {
ERROR(
"Initial values for the rotation map need to be provided when using "
"the worldtube.");
}
if (rotation_map_options_->decay_timescale.has_value()) {
ERROR(
"The BinaryCompactObject domains don't support using a "
"SettleToConstant function of time for the rotation map. Set the "
"decay timescale to 'None'");
}

result[rotation_name] =
std::make_unique<FunctionsOfTime::IntegratedFunctionOfTime>(
initial_time_,
std::array<double, 2>{
0.,
gsl::at(rotation_map_options_.value().initial_angular_velocity,
2)},
std::array<double, 2>{0., rotation_map_options_.value().angles[1][2]},
initial_expiration_time, true);

// Size and Shape FunctionOfTime for objects A and B. Only spherical excision
Expand Down Expand Up @@ -189,64 +203,52 @@ TimeDependentMapOptions<IsCylindrical>::create_functions_of_time(
// ExpansionMap FunctionOfTime for the function \f$a(t)\f$ in the
// domain::CoordinateMaps::TimeDependent::RotScaleTrans map
if (expansion_map_options_.has_value()) {
if (not expansion_map_options_->asymptotic_velocity_outer_boundary
.has_value()) {
ERROR(
"The BinaryCompactObject domains don't support using a "
"SettleToConstant function of time for the expansion map. Set the "
"asymptotic velocity to a value and set the decay timescale to "
"'None'");
}

result[expansion_name] =
std::make_unique<FunctionsOfTime::PiecewisePolynomial<2>>(
initial_time_,
std::array<DataVector, 3>{
{{gsl::at(expansion_map_options_.value().initial_values, 0)},
{gsl::at(expansion_map_options_.value().initial_values, 1)},
{0.0}}},
initial_time_, expansion_map_options_->initial_values,
expiration_times.at(expansion_name));

// ExpansionMap FunctionOfTime for the function \f$b(t)\f$ in the
// domain::CoordinateMaps::TimeDependent::RotScaleTrans map
result[expansion_outer_boundary_name] =
std::make_unique<FunctionsOfTime::FixedSpeedCubic>(
1.0, initial_time_,
expansion_map_options_.value().outer_boundary_velocity,
expansion_map_options_.value().outer_boundary_decay_time);
expansion_map_options_->initial_values_outer_boundary[0][0],
initial_time_,
expansion_map_options_->asymptotic_velocity_outer_boundary.value(),
expansion_map_options_->decay_timescale_outer_boundary);
}

// RotationMap FunctionOfTime for the rotation angles about each
// axis. The initial rotation angles don't matter as we never
// actually use the angles themselves. We only use their derivatives
// (omega) to determine map parameters. In theory we could determine
// each initial angle from the input axis-angle representation, but
// we don't need to.
// axis.
if (rotation_map_options_.has_value()) {
result[rotation_name] = std::make_unique<
FunctionsOfTime::QuaternionFunctionOfTime<3>>(
initial_time_,
std::array<DataVector, 1>{DataVector{1.0, 0.0, 0.0, 0.0}},
std::array<DataVector, 4>{
{{3, 0.0},
{gsl::at(rotation_map_options_.value().initial_angular_velocity,
0),
gsl::at(rotation_map_options_.value().initial_angular_velocity,
1),
gsl::at(rotation_map_options_.value().initial_angular_velocity,
2)},
{3, 0.0},
{3, 0.0}}},
expiration_times.at(rotation_name));
if (rotation_map_options_->decay_timescale.has_value()) {
ERROR(
"The BinaryCompactObject domains don't support using a "
"SettleToConstant function of time for the rotation map. Set the "
"decay timescale to 'None'");
}

result[rotation_name] =
std::make_unique<FunctionsOfTime::QuaternionFunctionOfTime<3>>(
initial_time_, std::array{rotation_map_options_->quaternions[0]},
rotation_map_options_->angles, expiration_times.at(rotation_name));
}

// TranslationMap FunctionOfTime
if (translation_map_options_.has_value()) {
result[translation_name] = std::make_unique<
FunctionsOfTime::PiecewisePolynomial<2>>(
initial_time_,
std::array<DataVector, 3>{
{{gsl::at(translation_map_options_.value().initial_values, 0)[0],
gsl::at(translation_map_options_.value().initial_values, 0)[1],
gsl::at(translation_map_options_.value().initial_values, 0)[2]},
{gsl::at(translation_map_options_.value().initial_values, 1)[0],
gsl::at(translation_map_options_.value().initial_values, 1)[1],
gsl::at(translation_map_options_.value().initial_values, 1)[2]},
{gsl::at(translation_map_options_.value().initial_values, 2)[0],
gsl::at(translation_map_options_.value().initial_values, 2)[1],
gsl::at(translation_map_options_.value().initial_values, 2)[2]}}},
expiration_times.at(translation_name));
result[translation_name] =
std::make_unique<FunctionsOfTime::PiecewisePolynomial<2>>(
initial_time_, translation_map_options_->initial_values,
expiration_times.at(translation_name));
}

// Size and Shape FunctionOfTime for objects A and B
Expand Down
101 changes: 11 additions & 90 deletions src/Domain/Creators/TimeDependentOptions/BinaryCompactObject.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -17,8 +17,10 @@
#include "Domain/CoordinateMaps/TimeDependent/Rotation.hpp"
#include "Domain/CoordinateMaps/TimeDependent/Shape.hpp"
#include "Domain/CoordinateMaps/TimeDependent/ShapeMapTransitionFunctions/ShapeMapTransitionFunction.hpp"
#include "Domain/Creators/TimeDependentOptions/ExpansionMap.hpp"
#include "Domain/Creators/TimeDependentOptions/RotationMap.hpp"
#include "Domain/Creators/TimeDependentOptions/ShapeMap.hpp"
#include "Domain/Creators/TimeDependentOptions/Sphere.hpp"
#include "Domain/Creators/TimeDependentOptions/TranslationMap.hpp"
#include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
#include "Domain/Structure/ObjectLabel.hpp"
#include "Options/Auto.hpp"
Expand Down Expand Up @@ -146,101 +148,20 @@ struct TimeDependentMapOptions {
/// The outer boundary radius of the map is always set to
/// the outer boundary of the Domain, so there is no option
/// here to set the outer boundary radius.
struct ExpansionMapOptions {
using type = Options::Auto<ExpansionMapOptions, Options::AutoLabel::None>;
static std::string name() { return "ExpansionMap"; }
static constexpr Options::String help = {
"Options for the expansion map. Specify 'None' to not use this map."};
struct InitialValues {
using type = std::array<double, 2>;
static constexpr Options::String help = {
"Initial value and deriv of expansion."};
};
struct AsymptoticVelocityOuterBoundary {
using type = double;
static constexpr Options::String help = {
"The asymptotic velocity of the outer boundary."};
};
struct DecayTimescaleOuterBoundaryVelocity {
using type = double;
static constexpr Options::String help = {
"The timescale for how fast the outer boundary velocity approaches "
"its asymptotic value."};
};
using options = tmpl::list<InitialValues, AsymptoticVelocityOuterBoundary,
DecayTimescaleOuterBoundaryVelocity>;
ExpansionMapOptions() = default;
ExpansionMapOptions(std::array<double, 2> initial_values_in,
double outer_boundary_velocity_in,
double outer_boundary_decay_time_in)
: initial_values(initial_values_in),
outer_boundary_velocity(outer_boundary_velocity_in),
outer_boundary_decay_time(outer_boundary_decay_time_in) {}

std::array<double, 2> initial_values{
std::numeric_limits<double>::signaling_NaN(),
std::numeric_limits<double>::signaling_NaN()};
double outer_boundary_velocity{
std::numeric_limits<double>::signaling_NaN()};
double outer_boundary_decay_time{
std::numeric_limits<double>::signaling_NaN()};
};

struct RotationMapOptions {
using type = Options::Auto<RotationMapOptions, Options::AutoLabel::None>;
static std::string name() { return "RotationMap"; }
static constexpr Options::String help = {
"Options for a time-dependent rotation map about an arbitrary axis. "
"Specify 'None' to not use this map."};

struct InitialAngularVelocity {
using type = std::array<double, 3>;
static constexpr Options::String help = {"The initial angular velocity."};
};

using options = tmpl::list<InitialAngularVelocity>;
using ExpansionMapOptions =
domain::creators::time_dependent_options::ExpansionMapOptions;

RotationMapOptions() = default;
explicit RotationMapOptions(
std::array<double, 3> initial_angular_velocity_in)
: initial_angular_velocity(initial_angular_velocity_in) {}

std::array<double, 3> initial_angular_velocity{};
};
/// \brief Options for the rotation map
using RotationMapOptions =
domain::creators::time_dependent_options::RotationMapOptions<3>;

/// \brief Options for the Translation Map, the outer radius is always set to
/// the outer boundary of the Domain, so there's no option needed for outer
/// boundary.
struct TranslationMapOptions {
using type = Options::Auto<TranslationMapOptions, Options::AutoLabel::None>;
static std::string name() { return "TranslationMap"; }
static constexpr Options::String help = {
"Options for a time-dependent translation map. Specify 'None' to not "
"use this map."};

struct InitialValues {
using type = std::array<std::array<double, 3>, 3>;
static constexpr Options::String help = {
"Initial position, velocity and acceleration."};
};

using options = tmpl::list<InitialValues>;
TranslationMapOptions() = default;
explicit TranslationMapOptions(
std::array<std::array<double, 3>, 3> initial_values_in)
: initial_values(initial_values_in) {}

std::array<std::array<double, 3>, 3> initial_values{};
};
using TranslationMapOptions =
domain::creators::time_dependent_options::TranslationMapOptions<3>;

// We use a type alias here instead of defining the ShapeMapOptions struct
// because there appears to be a bug in clang-10. If the definition of
// ShapeMapOptions is here inside TimeDependentMapOptions, on clang-10 there
// is a linking error that there is an undefined reference to
// Options::Option::parse_as<TimeDependentMapOptions<A>> (and B). This doesn't
// show up for GCC. If we put the definition of ShapeMapOptions outside of
// TimeDependentMapOptions and just use a type alias here, the linking error
// goes away.
/// \brief Options for the shape map
template <domain::ObjectLabel Object>
using ShapeMapOptions =
domain::creators::time_dependent_options::ShapeMapOptions<
Expand Down
12 changes: 9 additions & 3 deletions support/Pipelines/Bbh/Inspiral.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -101,11 +101,17 @@ DomainCreator:
TimeDependentMaps:
InitialTime: &InitialTime 0.0
ExpansionMap:
InitialValues: [1.0, {{ RadialExpansionVelocity }}]
InitialValues: [1.0, {{ RadialExpansionVelocity }}, 0.0]
InitialValuesOuterBoundary: [1.0, 0.0, 0.0]
DecayTimescale: Auto
DecayTimescaleOuterBoundary: 50.0
AsymptoticVelocityOuterBoundary: -1.0e-6
DecayTimescaleOuterBoundaryVelocity: 50.0
RotationMap:
InitialAngularVelocity: [0.0, 0.0, {{ InitialAngularVelocity }}]
InitialQuaternions: [[1.0, 0.0, 0.0, 0.0]]
InitialAngles:
- [0.0, 0.0, 0.0]
- [0.0, 0.0, {{ InitialAngularVelocity }}]
DecayTimescale: Auto
TranslationMap:
InitialValues: [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]
ShapeMapA:
Expand Down
12 changes: 9 additions & 3 deletions tests/InputFiles/CurvedScalarWave/WorldtubeKerrSchild.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -72,11 +72,17 @@ DomainCreator:
TimeDependentMaps:
InitialTime: 0.
ExpansionMap:
InitialValues: [1., 0.]
InitialValues: [1., 0., 0.]
InitialValuesOuterBoundary: [1., 0., 0.]
DecayTimescale: Auto
AsymptoticVelocityOuterBoundary: 0.
DecayTimescaleOuterBoundaryVelocity: 1.
DecayTimescaleOuterBoundary: 1.
RotationMap:
InitialAngularVelocity: [0., 0., 0.08944271909999159]
InitialQuaternions: [[1., 0., 0., 0.]]
InitialAngles:
- [0., 0., 0.]
- [0., 0., 0.08944271909999159]
DecayTimescale: Auto
TranslationMap:
InitialValues: [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]
ShapeMapA:
Expand Down
12 changes: 9 additions & 3 deletions tests/InputFiles/ExportCoordinates/InputTimeDependent3D.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -68,11 +68,17 @@ DomainCreator:
TimeDependentMaps:
InitialTime: 0.0
ExpansionMap:
InitialValues: [1.0, -4.6148457646200002e-05]
InitialValues: [1.0, -4.6148457646200002e-05, 0.0]
InitialValuesOuterBoundary: [1.0, 0.0, 0.0]
DecayTimescale: Auto
AsymptoticVelocityOuterBoundary: -1.0e-6
DecayTimescaleOuterBoundaryVelocity: 50.0
DecayTimescaleOuterBoundary: 50.0
RotationMap:
InitialAngularVelocity: [0.0, 0.0, 1.5264577062000000e-02]
InitialQuaternions: [[1.0, 0.0, 0.0, 0.0]]
InitialAngles:
- [0.0, 0.0, 0.0,]
- [0.0, 0.0, 1.5264577062000000e-02]
DecayTimescale: Auto
TranslationMap:
InitialValues: [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]
ShapeMapA:
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -58,11 +58,17 @@ DomainCreator:
TimeDependentMaps:
InitialTime: 0.0
ExpansionMap:
InitialValues: [1.0, -4.6148457646200002e-05]
InitialValues: [1.0, -4.6148457646200002e-05, 0.0]
InitialValuesOuterBoundary: [1.0, 0.0, 0.0]
DecayTimescale: Auto
DecayTimescaleOuterBoundary: 50.0
AsymptoticVelocityOuterBoundary: -1.0e-6
DecayTimescaleOuterBoundaryVelocity: 50.0
RotationMap:
InitialAngularVelocity: [0.0, 0.0, 1.5264577062000000e-02]
InitialQuaternions: [[1.0, 0.0, 0.0, 0.0]]
InitialAngles:
- [0.0, 0.0, 0.0]
- [0.0, 0.0, 1.5264577062000000e-02]
DecayTimescale: Auto
TranslationMap:
InitialValues: [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]
ShapeMapA:
Expand Down
12 changes: 9 additions & 3 deletions tests/InputFiles/GrMhd/GhValenciaDivClean/BinaryNeutronStar.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -64,11 +64,17 @@ DomainCreator:
TimeDependentMaps:
InitialTime: 0.0
ExpansionMap:
InitialValues: [1.0, 0.0]
InitialValues: [1.0, 0.0, 0.0]
InitialValuesOuterBoundary: [1.0, 0.0, 0.0]
DecayTimescale: Auto
AsymptoticVelocityOuterBoundary: -1.0e-6
DecayTimescaleOuterBoundaryVelocity: 50.0
DecayTimescaleOuterBoundary: 50.0
RotationMap:
InitialAngularVelocity: [0.0, 0.0, 0.00826225]
InitialQuaternions: [[1.0, 0.0, 0.0, 0.0]]
InitialAngles:
- [0.0, 0.0, 0.0]
- [0.0, 0.0, 0.00826225]
DecayTimescale: Auto
TranslationMap:
InitialValues: [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]
ShapeMapA: None
Expand Down
Loading

0 comments on commit 9fbb81d

Please sign in to comment.