diff --git a/src/PointwiseFunctions/AnalyticSolutions/GeneralRelativity/CMakeLists.txt b/src/PointwiseFunctions/AnalyticSolutions/GeneralRelativity/CMakeLists.txt index 524b8d27a52c..282c7e567208 100644 --- a/src/PointwiseFunctions/AnalyticSolutions/GeneralRelativity/CMakeLists.txt +++ b/src/PointwiseFunctions/AnalyticSolutions/GeneralRelativity/CMakeLists.txt @@ -44,6 +44,7 @@ target_link_libraries( GeneralRelativity Interpolation Serialization + SpecialRelativity ) add_subdirectory(Python) diff --git a/src/PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrSchild.cpp b/src/PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrSchild.cpp index 854d009dc7b9..d84ff61ca0e4 100644 --- a/src/PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrSchild.cpp +++ b/src/PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrSchild.cpp @@ -20,8 +20,11 @@ #include "PointwiseFunctions/GeneralRelativity/ExtrinsicCurvature.hpp" #include "PointwiseFunctions/GeneralRelativity/InverseSpacetimeMetric.hpp" #include "PointwiseFunctions/GeneralRelativity/Tags.hpp" +#include "PointwiseFunctions/GeneralRelativity/TimeDerivativeOfSpatialMetric.hpp" +#include "PointwiseFunctions/SpecialRelativity/LorentzBoostMatrix.hpp" #include "Utilities/ConstantExpressions.hpp" #include "Utilities/ContainerHelpers.hpp" +#include "Utilities/ErrorHandling/Assert.hpp" #include "Utilities/GenerateInstantiations.hpp" #include "Utilities/Gsl.hpp" #include "Utilities/MakeWithValue.hpp" @@ -34,13 +37,14 @@ KerrSchild::KerrSchild(CkMigrateMessage* /*msg*/) {} KerrSchild::KerrSchild(const double mass, const std::array& dimensionless_spin, const std::array& center, + const std::array& boost_velocity, const Options::Context& context) : mass_(mass), - // clang-tidy: do not std::move trivial types. - dimensionless_spin_(dimensionless_spin), // NOLINT - // clang-tidy: do not std::move trivial types. + dimensionless_spin_(dimensionless_spin), center_(center), - zero_spin_(dimensionless_spin_ == std::array{{0., 0., 0.}}) { + boost_velocity_(boost_velocity), + zero_spin_(dimensionless_spin_ == std::array{{0., 0., 0.}}), + zero_velocity_(boost_velocity_ == std::array{{0., 0., 0.}}) { const double spin_magnitude = magnitude(dimensionless_spin_); if (spin_magnitude > 1.0) { PARSE_ERROR(context, "Spin magnitude must be < 1. Given spin: " @@ -56,7 +60,9 @@ void KerrSchild::pup(PUP::er& p) { p | mass_; p | dimensionless_spin_; p | center_; + p | boost_velocity_; p | zero_spin_; + p | zero_velocity_; } template @@ -68,13 +74,25 @@ template void KerrSchild::IntermediateComputer::operator()( const gsl::not_null*> x_minus_center, const gsl::not_null /*cache*/, - internal_tags::x_minus_center /*meta*/) const { + internal_tags::x_minus_center_unboosted /*meta*/) const { *x_minus_center = x_; for (size_t d = 0; d < 3; ++d) { x_minus_center->get(d) -= gsl::at(solution_.center(), d); } } +template +void KerrSchild::IntermediateComputer::operator()( + const gsl::not_null*> x_minus_center_boosted, + const gsl::not_null cache, + internal_tags::x_minus_center /*meta*/) const { + const auto& x_minus_center = cache->get_var( + *this, internal_tags::x_minus_center_unboosted{}); + // Boost the coordinates to the rest frame + sr::lorentz_boost(x_minus_center_boosted, x_minus_center, 0., + solution_.boost_velocity()); +} + template void KerrSchild::IntermediateComputer::operator()( const gsl::not_null*> a_dot_x, @@ -251,9 +269,9 @@ void KerrSchild::IntermediateComputer::operator()( template void KerrSchild::IntermediateComputer::operator()( - const gsl::not_null*> deriv_H, + const gsl::not_null*> deriv_H, const gsl::not_null cache, - internal_tags::deriv_H /*meta*/) const { + internal_tags::deriv_H_unboosted /*meta*/) const { const auto& deriv_log_r = cache->get_var(*this, internal_tags::deriv_log_r{}); const auto& deriv_H_temp1 = @@ -262,12 +280,23 @@ void KerrSchild::IntermediateComputer::operator()( get(cache->get_var(*this, internal_tags::deriv_H_temp2{})); const auto spin_a = solution_.dimensionless_spin() * solution_.mass(); + get<0>(*deriv_H) = 0.; for (size_t i = 0; i < 3; ++i) { - deriv_H->get(i) = + deriv_H->get(i + 1) = deriv_H_temp1 * deriv_log_r.get(i) - deriv_H_temp2 * gsl::at(spin_a, i); } } +template +void KerrSchild::IntermediateComputer::operator()( + const gsl::not_null*> deriv_H_boosted, + const gsl::not_null cache, + internal_tags::deriv_H /*meta*/) const { + const auto& deriv_H = cache->get_var( + *this, internal_tags::deriv_H_unboosted{}); + sr::lorentz_boost(deriv_H_boosted, deriv_H, -solution_.boost_velocity()); +} + template void KerrSchild::IntermediateComputer::operator()( const gsl::not_null*> denom, @@ -297,9 +326,9 @@ void KerrSchild::IntermediateComputer::operator()( template void KerrSchild::IntermediateComputer::operator()( - const gsl::not_null*> null_form, + const gsl::not_null*> null_form, const gsl::not_null cache, - internal_tags::null_form /*meta*/) const { + internal_tags::null_form_unboosted /*meta*/) const { const auto& a_dot_x_over_r = get(cache->get_var(*this, internal_tags::a_dot_x_over_r{})); const auto& r = get(cache->get_var(*this, internal_tags::r{})); @@ -309,71 +338,102 @@ void KerrSchild::IntermediateComputer::operator()( get(cache->get_var(*this, internal_tags::denom{})); const auto spin_a = solution_.dimensionless_spin() * solution_.mass(); + get<0>(*null_form) = 1.; for (size_t i = 0; i < 3; ++i) { const size_t cross_product_index_1 = (i + 1) % 3; const size_t cross_product_index_2 = (i + 2) % 3; - null_form->get(i) = denom * (r * x_minus_center.get(i) - - gsl::at(spin_a, cross_product_index_1) * - x_minus_center.get(cross_product_index_2) + - gsl::at(spin_a, cross_product_index_2) * - x_minus_center.get(cross_product_index_1) + - a_dot_x_over_r * gsl::at(spin_a, i)); + null_form->get(i + 1) = + denom * (r * x_minus_center.get(i) - + gsl::at(spin_a, cross_product_index_1) * + x_minus_center.get(cross_product_index_2) + + gsl::at(spin_a, cross_product_index_2) * + x_minus_center.get(cross_product_index_1) + + a_dot_x_over_r * gsl::at(spin_a, i)); } } template void KerrSchild::IntermediateComputer::operator()( - const gsl::not_null*> deriv_null_form, + const gsl::not_null*> null_form_boosted, const gsl::not_null cache, - internal_tags::deriv_null_form /*meta*/) const { + internal_tags::null_form /*meta*/) const { + const auto& null_form = cache->get_var( + *this, internal_tags::null_form_unboosted{}); + sr::lorentz_boost(null_form_boosted, null_form, -solution_.boost_velocity()); +} + +template +void KerrSchild::IntermediateComputer::operator()( + const gsl::not_null*> deriv_null_form, + const gsl::not_null cache, + internal_tags::deriv_null_form_unboosted /*meta*/) const { const auto spin_a = solution_.dimensionless_spin() * solution_.mass(); const auto& denom = get(cache->get_var(*this, internal_tags::denom{})); const auto& r = get(cache->get_var(*this, internal_tags::r{})); const auto& x_minus_center = cache->get_var(*this, internal_tags::x_minus_center{}); - const auto& null_form = - cache->get_var(*this, internal_tags::null_form{}); + const auto& null_form = cache->get_var( + *this, internal_tags::null_form_unboosted{}); const auto& a_dot_x_over_rsquared = get( cache->get_var(*this, internal_tags::a_dot_x_over_rsquared{})); const auto& deriv_log_r = cache->get_var(*this, internal_tags::deriv_log_r{}); + for (size_t d = 0; d < 4; ++d) { + deriv_null_form->get(d, 0) = 0.; + deriv_null_form->get(0, d) = 0.; + } for (size_t i = 0; i < 3; i++) { for (size_t j = 0; j < 3; j++) { - deriv_null_form->get(j, i) = + deriv_null_form->get(j + 1, i + 1) = denom * (gsl::at(spin_a, i) * gsl::at(spin_a, j) / r + - (x_minus_center.get(i) - 2.0 * r * null_form.get(i) - + (x_minus_center.get(i) - 2.0 * r * null_form.get(i + 1) - a_dot_x_over_rsquared * gsl::at(spin_a, i)) * deriv_log_r.get(j) * r); if (i == j) { - deriv_null_form->get(j, i) += denom * r; + deriv_null_form->get(j + 1, i + 1) += denom * r; } else { // add denom*epsilon^ijk a_k size_t k = (j + 1) % 3; if (k == i) { // j+1 = i (cyclic), so choose minus sign k++; k = k % 3; // and set k to be neither i nor j - deriv_null_form->get(j, i) -= denom * gsl::at(spin_a, k); + deriv_null_form->get(j + 1, i + 1) -= denom * gsl::at(spin_a, k); } else { // i+1 = j (cyclic), so choose plus sign - deriv_null_form->get(j, i) += denom * gsl::at(spin_a, k); + deriv_null_form->get(j + 1, i + 1) += denom * gsl::at(spin_a, k); } } } } } +template +void KerrSchild::IntermediateComputer::operator()( + const gsl::not_null*> deriv_null_form_boosted, + const gsl::not_null cache, + internal_tags::deriv_null_form /*meta*/) const { + const auto& deriv_null_form = cache->get_var( + *this, internal_tags::deriv_null_form_unboosted{}); + // Inverse-boost the first index, because it represents the derivative which + // is given in the boosted frame + sr::lorentz_boost(deriv_null_form_boosted, deriv_null_form, + -solution_.boost_velocity()); +} + template void KerrSchild::IntermediateComputer::operator()( const gsl::not_null*> lapse_squared, const gsl::not_null cache, internal_tags::lapse_squared /*meta*/) const { - if (solution_.zero_spin()) { + if (solution_.zero_spin() and solution_.zero_velocity()) { const auto& r = get(cache->get_var(*this, internal_tags::r{})); get(*lapse_squared) = 1.0 / (1.0 + 2.0 * solution_.mass() / r); return; } const auto& H = get(cache->get_var(*this, internal_tags::H{})); - get(*lapse_squared) = 1.0 / (1.0 + 2.0 * square(null_vector_0_) * H); + const auto& null_form = + cache->get_var(*this, internal_tags::null_form{}); + get(*lapse_squared) = 1.0 / (1.0 + 2.0 * square(get<0>(null_form)) * H); } template @@ -394,8 +454,7 @@ void KerrSchild::IntermediateComputer::operator()( const auto& lapse = get(cache->get_var(*this, gr::Tags::Lapse{})); const auto& lapse_squared = get(cache->get_var(*this, internal_tags::lapse_squared{})); - get(*deriv_lapse_multiplier) = - -square(null_vector_0_) * lapse * lapse_squared; + get(*deriv_lapse_multiplier) = -lapse * lapse_squared; } template @@ -404,10 +463,12 @@ void KerrSchild::IntermediateComputer::operator()( const gsl::not_null cache, internal_tags::shift_multiplier /*meta*/) const { const auto& H = get(cache->get_var(*this, internal_tags::H{})); + const auto& null_form = + cache->get_var(*this, internal_tags::null_form{}); const auto& lapse_squared = get(cache->get_var(*this, internal_tags::lapse_squared{})); - get(*shift_multiplier) = -2.0 * null_vector_0_ * H * lapse_squared; + get(*shift_multiplier) = 2.0 * get<0>(null_form) * H * lapse_squared; } template @@ -415,7 +476,7 @@ void KerrSchild::IntermediateComputer::operator()( const gsl::not_null*> shift, const gsl::not_null cache, gr::Tags::Shift /*meta*/) const { - if (solution_.zero_spin()) { + if (solution_.zero_spin() and solution_.zero_velocity()) { const auto& x_minus_center = cache->get_var(*this, internal_tags::x_minus_center{}); const auto& r_squared = @@ -434,7 +495,7 @@ void KerrSchild::IntermediateComputer::operator()( get(cache->get_var(*this, internal_tags::shift_multiplier{})); for (size_t i = 0; i < 3; ++i) { - shift->get(i) = shift_multiplier * null_form.get(i); + shift->get(i) = shift_multiplier * null_form.get(i + 1); } } @@ -443,7 +504,7 @@ void KerrSchild::IntermediateComputer::operator()( const gsl::not_null*> deriv_shift, const gsl::not_null cache, DerivShift /*meta*/) const { - if (solution_.zero_spin()) { + if (solution_.zero_spin() and solution_.zero_velocity()) { const auto& x_minus_center = cache->get_var(*this, internal_tags::x_minus_center{}); const auto& r = get(cache->get_var(*this, internal_tags::r{})); @@ -476,12 +537,15 @@ void KerrSchild::IntermediateComputer::operator()( for (size_t m = 0; m < 3; ++m) { for (size_t i = 0; i < 3; ++i) { - deriv_shift->get(m, i) = 4.0 * cube(null_vector_0_) * H * - null_form.get(i) * square(lapse_squared) * - deriv_H.get(m) - - 2.0 * null_vector_0_ * lapse_squared * - (null_form.get(i) * deriv_H.get(m) + - H * deriv_null_form.get(m, i)); + deriv_shift->get(m, i) = + 2.0 * lapse_squared * + (get<0>(null_form) * null_form.get(i + 1) * deriv_H.get(m + 1) + + H * get<0>(null_form) * deriv_null_form.get(m + 1, i + 1) + + H * null_form.get(i + 1) * deriv_null_form.get(m + 1, 0)) - + 4.0 * H * get<0>(null_form) * null_form.get(i + 1) * + square(lapse_squared) * + (square(get<0>(null_form)) * deriv_H.get(m + 1) + + 2. * H * get<0>(null_form) * deriv_null_form.get(m + 1, 0)); } } } @@ -491,7 +555,7 @@ void KerrSchild::IntermediateComputer::operator()( const gsl::not_null*> spatial_metric, const gsl::not_null cache, gr::Tags::SpatialMetric /*meta*/) const { - if (solution_.zero_spin()) { + if (solution_.zero_spin() and solution_.zero_velocity()) { const auto& x_minus_center = cache->get_var(*this, internal_tags::x_minus_center{}); const auto& r = get(cache->get_var(*this, internal_tags::r{})); @@ -516,7 +580,7 @@ void KerrSchild::IntermediateComputer::operator()( spatial_metric->get(i, i) = 1.; for (size_t j = i; j < 3; ++j) { // Symmetry spatial_metric->get(i, j) += - 2.0 * H * null_form.get(i) * null_form.get(j); + 2.0 * H * null_form.get(i + 1) * null_form.get(j + 1); } } } @@ -531,7 +595,7 @@ void KerrSchild::IntermediateComputer::operator()( get(cache->get_var(*this, internal_tags::lapse_squared{})); const auto& null_form = cache->get_var(*this, internal_tags::null_form{}); - if (solution_.zero_spin()) { + if (solution_.zero_spin() and solution_.zero_velocity()) { const auto& x_minus_center = cache->get_var(*this, internal_tags::x_minus_center{}); const auto& r = get(cache->get_var(*this, internal_tags::r{})); @@ -551,8 +615,9 @@ void KerrSchild::IntermediateComputer::operator()( for (size_t i = 0; i < 3; ++i) { for (size_t j = i; j < 3; ++j) { // Symmetry - inverse_spatial_metric->get(i, j) = - -2.0 * H * lapse_squared * null_form.get(i) * null_form.get(j); + inverse_spatial_metric->get(i, j) = -2.0 * H * lapse_squared * + null_form.get(i + 1) * + null_form.get(j + 1); } inverse_spatial_metric->get(i, i) += 1.; } @@ -575,10 +640,11 @@ void KerrSchild::IntermediateComputer::operator()( for (size_t j = i; j < 3; ++j) { // Symmetry for (size_t m = 0; m < 3; ++m) { deriv_spatial_metric->get(m, i, j) = - 2.0 * null_form.get(i) * null_form.get(j) * deriv_H.get(m) + + 2.0 * null_form.get(i + 1) * null_form.get(j + 1) * + deriv_H.get(m + 1) + 2.0 * H * - (null_form.get(i) * deriv_null_form.get(m, j) + - null_form.get(j) * deriv_null_form.get(m, i)); + (null_form.get(i + 1) * deriv_null_form.get(m + 1, j + 1) + + null_form.get(j + 1) * deriv_null_form.get(m + 1, i + 1)); } } } @@ -587,9 +653,62 @@ void KerrSchild::IntermediateComputer::operator()( template void KerrSchild::IntermediateComputer::operator()( const gsl::not_null*> dt_spatial_metric, - const gsl::not_null /*cache*/, + const gsl::not_null cache, ::Tags::dt> /*meta*/) const { - std::fill(dt_spatial_metric->begin(), dt_spatial_metric->end(), 0.); + if (solution_.zero_velocity()) { + std::fill(dt_spatial_metric->begin(), dt_spatial_metric->end(), 0.); + return; + } + const auto& ex_curvature = + cache->get_var(*this, gr::Tags::ExtrinsicCurvature{}); + const auto& lapse = cache->get_var(*this, gr::Tags::Lapse{}); + const auto& shift = + cache->get_var(*this, gr::Tags::Shift{}); + const auto& deriv_shift = + cache->get_var(*this, DerivShift{}); + const auto& spatial_metric = + cache->get_var(*this, gr::Tags::SpatialMetric{}); + const auto& deriv_spatial_metric = + cache->get_var(*this, DerivSpatialMetric{}); + // Reconstruct from the extrinsic curvature and shift + gr::time_derivative_of_spatial_metric(dt_spatial_metric, lapse, shift, + deriv_shift, spatial_metric, + deriv_spatial_metric, ex_curvature); +} + +template +void KerrSchild::IntermediateComputer::operator()( + const gsl::not_null*> null_form_dot_deriv_H, + const gsl::not_null cache, + internal_tags::null_form_dot_deriv_H /*meta*/) const { + const auto& deriv_H = + cache->get_var(*this, internal_tags::deriv_H{}); + const auto& null_form = + cache->get_var(*this, internal_tags::null_form{}); + get(*null_form_dot_deriv_H) = 0.0; + for (size_t i = 0; i < 3; ++i) { + get(*null_form_dot_deriv_H) += null_form.get(i + 1) * deriv_H.get(i + 1); + } +} + +template +void KerrSchild::IntermediateComputer::operator()( + const gsl::not_null*> + null_form_dot_deriv_null_form, + const gsl::not_null cache, + internal_tags::null_form_dot_deriv_null_form /*meta*/) + const { + const auto& deriv_null_form = + cache->get_var(*this, internal_tags::deriv_null_form{}); + const auto& null_form = + cache->get_var(*this, internal_tags::null_form{}); + for (size_t j = 0; j < 3; ++j) { + null_form_dot_deriv_null_form->get(j) = 0.0; + for (size_t i = 0; i < 3; ++i) { + null_form_dot_deriv_null_form->get(j) += + null_form.get(i + 1) * deriv_null_form.get(i + 1, j + 1); + } + } } template @@ -597,14 +716,52 @@ void KerrSchild::IntermediateComputer::operator()( const gsl::not_null*> extrinsic_curvature, const gsl::not_null cache, gr::Tags::ExtrinsicCurvature /*meta*/) const { - gr::extrinsic_curvature( - extrinsic_curvature, cache->get_var(*this, gr::Tags::Lapse{}), - cache->get_var(*this, gr::Tags::Shift{}), - cache->get_var(*this, DerivShift{}), - cache->get_var(*this, gr::Tags::SpatialMetric{}), - cache->get_var(*this, - ::Tags::dt>{}), - cache->get_var(*this, DerivSpatialMetric{})); + const auto& lapse = cache->get_var(*this, gr::Tags::Lapse{}); + const auto& H = cache->get_var(*this, internal_tags::H{}); + const auto& deriv_H = + cache->get_var(*this, internal_tags::deriv_H{}); + const auto& null_form = + cache->get_var(*this, internal_tags::null_form{}); + const auto& deriv_null_form = + cache->get_var(*this, internal_tags::deriv_null_form{}); + const auto& null_form_dot_deriv_H = + cache->get_var(*this, internal_tags::null_form_dot_deriv_H{}); + const auto& null_form_dot_deriv_null_form = cache->get_var( + *this, internal_tags::null_form_dot_deriv_null_form{}); + // We use a formula with l^0. However, since we actually use l_0, there is a + // sign difference. This can be seen by lowering the index with the + // Minkowski metric, which can be done for the null form given the form of the + // KerrSchild metric. + const double l0_relative_sign = -1.0; + for (size_t i = 0; i < 3; ++i) { + for (size_t j = i; j < 3; ++j) { + extrinsic_curvature->get(i, j) = + (-1.0 / get(lapse)) * + (null_form.get(i + 1) * null_form.get(j + 1) * get<0>(deriv_H) + + get(H) * + (null_form.get(i + 1) * deriv_null_form.get(0, j + 1) + + null_form.get(j + 1) * deriv_null_form.get(0, i + 1))) - + get(lapse) * + (get(H) * + (l0_relative_sign * null_form.get(i + 1) * + deriv_null_form.get(j + 1, 0) + + l0_relative_sign * null_form.get(j + 1) * + deriv_null_form.get(i + 1, 0) + + l0_relative_sign * get<0>(null_form) * + (deriv_null_form.get(i + 1, j + 1) + + deriv_null_form.get(j + 1, i + 1) + + 2.0 * get(H) * + (null_form.get(i + 1) * + null_form_dot_deriv_null_form.get(j) + + null_form.get(j + 1) * + null_form_dot_deriv_null_form.get(i)) + + 2.0 * null_form.get(i + 1) * null_form.get(j + 1) * + get(null_form_dot_deriv_H))) + + l0_relative_sign * get<0>(null_form) * + (null_form.get(i + 1) * deriv_H.get(j + 1) + + null_form.get(j + 1) * deriv_H.get(i + 1))); + } + } } template @@ -640,7 +797,7 @@ tnsr::i KerrSchild::IntermediateVars::get_var( const IntermediateComputer& computer, DerivLapse /*meta*/) { - if (computer.solution().zero_spin()) { + if (computer.solution().zero_spin() and computer.solution().zero_velocity()) { const auto& x_minus_center = get_var(computer, internal_tags::x_minus_center{}); const auto& r = get(get_var(computer, internal_tags::r{})); @@ -659,13 +816,21 @@ KerrSchild::IntermediateVars::get_var( return d_lapse; } tnsr::i result{}; + const auto& H = get_var(computer, internal_tags::H{}); const auto& deriv_H = get_var(computer, internal_tags::deriv_H{}); const auto& deriv_lapse_multiplier = get(get_var(computer, internal_tags::deriv_lapse_multiplier{})); + const auto& null_form = + get_var(computer, internal_tags::null_form{}); + const auto& deriv_null_form = + get_var(computer, internal_tags::deriv_null_form{}); for (size_t i = 0; i < 3; ++i) { - result.get(i) = deriv_lapse_multiplier * deriv_H.get(i); + result.get(i) = + deriv_lapse_multiplier * + (square(get<0>(null_form)) * deriv_H.get(i + 1) + + 2. * get(H) * get<0>(null_form) * deriv_null_form.get(i + 1, 0)); } return result; } @@ -702,11 +867,13 @@ KerrSchild::IntermediateVars::get_var( gr::Tags::DerivDetSpatialMetric /*meta*/) { const auto& deriv_H = get_var(computer, internal_tags::deriv_H{}); + const auto& null_form = + get_var(computer, internal_tags::null_form{}); auto result = make_with_value>(get<0>(deriv_H), 0.); for (size_t i = 0; i < 3; ++i) { - result.get(i) = 2.0 * square(null_vector_0_) * deriv_H.get(i); + result.get(i) = 2.0 * square(get<0>(null_form)) * deriv_H.get(i + 1); } return result; @@ -716,7 +883,7 @@ template Scalar KerrSchild::IntermediateVars::get_var( const IntermediateComputer& computer, gr::Tags::TraceExtrinsicCurvature /*meta*/) { - if (computer.solution().zero_spin()) { + if (computer.solution().zero_spin() and computer.solution().zero_velocity()) { const double m = computer.solution().mass(); const auto& r = get(get_var(computer, internal_tags::r{})); Scalar result(get_size(r)); @@ -735,7 +902,7 @@ tnsr::I KerrSchild::IntermediateVars::get_var( const IntermediateComputer& computer, gr::Tags::TraceSpatialChristoffelSecondKind /*meta*/) { - if (computer.solution().zero_spin()) { + if (computer.solution().zero_spin() and computer.solution().zero_velocity()) { const double m = computer.solution().mass(); const auto& r = get(get_var(computer, internal_tags::r{})); const auto& r_squared = @@ -763,6 +930,10 @@ tnsr::Abb KerrSchild::IntermediateVars::get_var( const IntermediateComputer& computer, gr::Tags::SpacetimeChristoffelSecondKind /*meta*/) { + ASSERT( + computer.solution().zero_velocity(), + "The tag gr::Tags::SpacetimeChristoffelSecondKind is not supported for a " + "boosted KerrSchild solution."); const auto& lapse = get_var(computer, gr::Tags::Lapse{}); const auto& shift = get_var(computer, gr::Tags::Shift{}); const auto& spatial_metric = diff --git a/src/PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrSchild.hpp b/src/PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrSchild.hpp index 0a3ce12e0def..6e14da7780d3 100644 --- a/src/PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrSchild.hpp +++ b/src/PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrSchild.hpp @@ -214,6 +214,30 @@ namespace Solutions { * * Right now we use (\f$\ref{eq:sphertocartsimple}\f$), but we may * wish to use the other transformation in the future. + * + * ## Boost of the Kerr-Schild solution + * + * We add initial momentum to the solution by applying a Lorentz boost to the + * metric. Since the Kerr-Schild metric can be expressed covariantly in terms of + * the Minkowski metric, a scalar function and a one form, we construct the + * metric in the rest frame of the black hole and then apply an inverse boost to + * each of the covariant objects individually. Notice that we also need to + * appropriately boost the coordinates to the to the rest frame before computing + * the metric. + * + * \warning While technically the boosted Kerr-Schild metric is dependent on + * both the time and space coordinates, we have implemented it only at $t = 0$ + * as in SpEC. Therefore it is technically not an analytic solution and should + * not be used to compute errors with respect to it. + * + * Moreover, since the boosted solution is intended for use as initial data, + * we do not compute the time derivatives of the lapse and shift in the boosted + * frame but set them to zero. + * + * Consequently, the gr::Tags::SpacetimeChristoffelSecondKind computed here, + * corresponds to the boosted Kerr-Schild for the gauge where lapse and shift + * have vanishing derivatives. + * */ class KerrSchild : public AnalyticSolution<3_st>, public MarkAsAnalyticSolution { @@ -233,12 +257,18 @@ class KerrSchild : public AnalyticSolution<3_st>, static constexpr Options::String help = { "The [x,y,z] center of the black hole"}; }; - using options = tmpl::list; + struct Velocity { + using type = std::array; + static constexpr Options::String help = { + "The [x,y,z] boost velocity of the black hole"}; + }; + using options = tmpl::list; static constexpr Options::String help{ "Black hole in Kerr-Schild coordinates"}; KerrSchild(double mass, const std::array& dimensionless_spin, const std::array& center, + const std::array& boost_velocity = {{0., 0., 0.}}, const Options::Context& context = {}); explicit KerrSchild(CkMigrateMessage* /*msg*/); @@ -282,55 +312,72 @@ class KerrSchild : public AnalyticSolution<3_st>, const std::array& dimensionless_spin() const { return dimensionless_spin_; } + const std::array& boost_velocity() const { + return boost_velocity_; + } bool zero_spin() const { return zero_spin_; } + bool zero_velocity() const { return zero_velocity_; } struct internal_tags { template - using x_minus_center = ::Tags::TempI<0, 3, Frame, DataType>; + using x_minus_center_unboosted = ::Tags::TempI<0, 3, Frame, DataType>; + template + using x_minus_center = ::Tags::TempI<1, 3, Frame, DataType>; template - using a_dot_x = ::Tags::TempScalar<1, DataType>; + using a_dot_x = ::Tags::TempScalar<2, DataType>; template - using a_dot_x_squared = ::Tags::TempScalar<2, DataType>; + using a_dot_x_squared = ::Tags::TempScalar<3, DataType>; template - using half_xsq_minus_asq = ::Tags::TempScalar<3, DataType>; + using half_xsq_minus_asq = ::Tags::TempScalar<4, DataType>; template - using r_squared = ::Tags::TempScalar<4, DataType>; + using r_squared = ::Tags::TempScalar<5, DataType>; template - using r = ::Tags::TempScalar<5, DataType>; + using r = ::Tags::TempScalar<6, DataType>; template - using a_dot_x_over_rsquared = ::Tags::TempScalar<6, DataType>; + using a_dot_x_over_rsquared = ::Tags::TempScalar<7, DataType>; template - using deriv_log_r_denom = ::Tags::TempScalar<7, DataType>; + using deriv_log_r_denom = ::Tags::TempScalar<8, DataType>; template - using deriv_log_r = ::Tags::Tempi<8, 3, Frame, DataType>; + using deriv_log_r = ::Tags::Tempi<9, 3, Frame, DataType>; template - using H_denom = ::Tags::TempScalar<9, DataType>; + using H_denom = ::Tags::TempScalar<10, DataType>; template - using H = ::Tags::TempScalar<10, DataType>; + using H = ::Tags::TempScalar<11, DataType>; template - using deriv_H_temp1 = ::Tags::TempScalar<11, DataType>; + using deriv_H_temp1 = ::Tags::TempScalar<12, DataType>; template - using deriv_H_temp2 = ::Tags::TempScalar<12, DataType>; + using deriv_H_temp2 = ::Tags::TempScalar<13, DataType>; + template + using deriv_H_unboosted = ::Tags::Tempa<14, 3, Frame, DataType>; template - using deriv_H = ::Tags::Tempi<13, 3, Frame, DataType>; + using deriv_H = ::Tags::Tempa<15, 3, Frame, DataType>; template - using denom = ::Tags::TempScalar<14, DataType>; + using denom = ::Tags::TempScalar<16, DataType>; template - using a_dot_x_over_r = ::Tags::TempScalar<15, DataType>; + using a_dot_x_over_r = ::Tags::TempScalar<17, DataType>; + template + using null_form_unboosted = ::Tags::Tempa<18, 3, Frame, DataType>; template - using null_form = ::Tags::Tempi<16, 3, Frame, DataType>; + using null_form = ::Tags::Tempa<19, 3, Frame, DataType>; template - using deriv_null_form = ::Tags::Tempij<17, 3, Frame, DataType>; + using deriv_null_form_unboosted = ::Tags::Tempab<20, 3, Frame, DataType>; + template + using deriv_null_form = ::Tags::Tempab<21, 3, Frame, DataType>; template - using lapse_squared = ::Tags::TempScalar<18, DataType>; + using null_form_dot_deriv_H = ::Tags::TempScalar<22, DataType>; + template + using null_form_dot_deriv_null_form = ::Tags::Tempi<23, 3, Frame, DataType>; template - using deriv_lapse_multiplier = ::Tags::TempScalar<19, DataType>; + using lapse_squared = ::Tags::TempScalar<24, DataType>; template - using shift_multiplier = ::Tags::TempScalar<20, DataType>; + using deriv_lapse_multiplier = ::Tags::TempScalar<25, DataType>; + template + using shift_multiplier = ::Tags::TempScalar<26, DataType>; }; template using CachedBuffer = CachedTempBuffer< + internal_tags::x_minus_center_unboosted, internal_tags::x_minus_center, internal_tags::a_dot_x, internal_tags::a_dot_x_squared, @@ -342,10 +389,15 @@ class KerrSchild : public AnalyticSolution<3_st>, internal_tags::H_denom, internal_tags::H, internal_tags::deriv_H_temp1, internal_tags::deriv_H_temp2, + internal_tags::deriv_H_unboosted, internal_tags::deriv_H, internal_tags::denom, internal_tags::a_dot_x_over_r, + internal_tags::null_form_unboosted, internal_tags::null_form, + internal_tags::deriv_null_form_unboosted, internal_tags::deriv_null_form, + internal_tags::null_form_dot_deriv_H, + internal_tags::null_form_dot_deriv_null_form, internal_tags::lapse_squared, gr::Tags::Lapse, internal_tags::deriv_lapse_multiplier, internal_tags::shift_multiplier, @@ -369,143 +421,175 @@ class KerrSchild : public AnalyticSolution<3_st>, const KerrSchild& solution() const { return solution_; } void operator()( - const gsl::not_null*> x_minus_center, - const gsl::not_null /*cache*/, + gsl::not_null*> x_minus_center, + gsl::not_null /*cache*/, + internal_tags::x_minus_center_unboosted /*meta*/) + const; + + void operator()( + gsl::not_null*> x_minus_center_boosted, + gsl::not_null /*cache*/, internal_tags::x_minus_center /*meta*/) const; - void operator()(const gsl::not_null*> a_dot_x, - const gsl::not_null cache, + void operator()(gsl::not_null*> a_dot_x, + gsl::not_null cache, internal_tags::a_dot_x /*meta*/) const; - void operator()(const gsl::not_null*> a_dot_x_squared, - const gsl::not_null cache, + void operator()(gsl::not_null*> a_dot_x_squared, + gsl::not_null cache, internal_tags::a_dot_x_squared /*meta*/) const; - void operator()(const gsl::not_null*> half_xsq_minus_asq, - const gsl::not_null cache, + void operator()(gsl::not_null*> half_xsq_minus_asq, + gsl::not_null cache, internal_tags::half_xsq_minus_asq /*meta*/) const; - void operator()(const gsl::not_null*> r_squared, - const gsl::not_null cache, + void operator()(gsl::not_null*> r_squared, + gsl::not_null cache, internal_tags::r_squared /*meta*/) const; - void operator()(const gsl::not_null*> r, - const gsl::not_null cache, + void operator()(gsl::not_null*> r, + gsl::not_null cache, internal_tags::r /*meta*/) const; void operator()( - const gsl::not_null*> a_dot_x_over_rsquared, - const gsl::not_null cache, + gsl::not_null*> a_dot_x_over_rsquared, + gsl::not_null cache, internal_tags::a_dot_x_over_rsquared /*meta*/) const; - void operator()(const gsl::not_null*> deriv_log_r_denom, - const gsl::not_null cache, + void operator()(gsl::not_null*> deriv_log_r_denom, + gsl::not_null cache, internal_tags::deriv_log_r_denom /*meta*/) const; - void operator()( - const gsl::not_null*> deriv_log_r, - const gsl::not_null cache, - internal_tags::deriv_log_r /*meta*/) const; + void operator()(gsl::not_null*> deriv_log_r, + gsl::not_null cache, + internal_tags::deriv_log_r /*meta*/) const; - void operator()(const gsl::not_null*> H_denom, - const gsl::not_null cache, + void operator()(gsl::not_null*> H_denom, + gsl::not_null cache, internal_tags::H_denom /*meta*/) const; - void operator()(const gsl::not_null*> H, - const gsl::not_null cache, + void operator()(gsl::not_null*> H, + gsl::not_null cache, internal_tags::H /*meta*/) const; - void operator()(const gsl::not_null*> deriv_H_temp1, - const gsl::not_null cache, + void operator()(gsl::not_null*> deriv_H_temp1, + gsl::not_null cache, internal_tags::deriv_H_temp1 /*meta*/) const; - void operator()(const gsl::not_null*> deriv_H_temp2, - const gsl::not_null cache, + void operator()(gsl::not_null*> deriv_H_temp2, + gsl::not_null cache, internal_tags::deriv_H_temp2 /*meta*/) const; - void operator()(const gsl::not_null*> deriv_H, - const gsl::not_null cache, + void operator()( + gsl::not_null*> deriv_H, + gsl::not_null cache, + internal_tags::deriv_H_unboosted /*meta*/) const; + + void operator()(gsl::not_null*> deriv_H_boosted, + gsl::not_null cache, internal_tags::deriv_H /*meta*/) const; - void operator()(const gsl::not_null*> denom, - const gsl::not_null cache, + void operator()(gsl::not_null*> denom, + gsl::not_null cache, internal_tags::denom /*meta*/) const; - void operator()(const gsl::not_null*> a_dot_x_over_r, - const gsl::not_null cache, + void operator()(gsl::not_null*> a_dot_x_over_r, + gsl::not_null cache, internal_tags::a_dot_x_over_r /*meta*/) const; - void operator()(const gsl::not_null*> null_form, - const gsl::not_null cache, - internal_tags::null_form /*meta*/) const; + void operator()( + gsl::not_null*> null_form, + gsl::not_null cache, + internal_tags::null_form_unboosted /*meta*/) const; + + void operator()( + gsl::not_null*> null_form_boosted, + gsl::not_null cache, + internal_tags::null_form /*meta*/) const; + + void operator()( + gsl::not_null*> deriv_null_form, + gsl::not_null cache, + internal_tags::deriv_null_form_unboosted /*meta*/) + const; void operator()( - const gsl::not_null*> deriv_null_form, - const gsl::not_null cache, + gsl::not_null*> deriv_null_form_boosted, + gsl::not_null cache, internal_tags::deriv_null_form /*meta*/) const; - void operator()(const gsl::not_null*> lapse_squared, - const gsl::not_null cache, + void operator()(gsl::not_null*> lapse_squared, + gsl::not_null cache, internal_tags::lapse_squared /*meta*/) const; - void operator()(const gsl::not_null*> lapse, - const gsl::not_null cache, + void operator()(gsl::not_null*> lapse, + gsl::not_null cache, gr::Tags::Lapse /*meta*/) const; void operator()( - const gsl::not_null*> deriv_lapse_multiplier, - const gsl::not_null cache, + gsl::not_null*> deriv_lapse_multiplier, + gsl::not_null cache, internal_tags::deriv_lapse_multiplier /*meta*/) const; - void operator()(const gsl::not_null*> shift_multiplier, - const gsl::not_null cache, + void operator()(gsl::not_null*> shift_multiplier, + gsl::not_null cache, internal_tags::shift_multiplier /*meta*/) const; - void operator()(const gsl::not_null*> shift, - const gsl::not_null cache, + void operator()(gsl::not_null*> shift, + gsl::not_null cache, gr::Tags::Shift /*meta*/) const; - void operator()( - const gsl::not_null*> deriv_shift, - const gsl::not_null cache, - DerivShift /*meta*/) const; + void operator()(gsl::not_null*> deriv_shift, + gsl::not_null cache, + DerivShift /*meta*/) const; - void operator()( - const gsl::not_null*> spatial_metric, - const gsl::not_null cache, - gr::Tags::SpatialMetric /*meta*/) const; + void operator()(gsl::not_null*> spatial_metric, + gsl::not_null cache, + gr::Tags::SpatialMetric /*meta*/) const; void operator()( - const gsl::not_null*> spatial_metric, - const gsl::not_null cache, + gsl::not_null*> spatial_metric, + gsl::not_null cache, gr::Tags::InverseSpatialMetric /*meta*/) const; - void operator()(const gsl::not_null*> - deriv_spatial_metric, - const gsl::not_null cache, - DerivSpatialMetric /*meta*/) const; + void operator()( + gsl::not_null*> deriv_spatial_metric, + gsl::not_null cache, + DerivSpatialMetric /*meta*/) const; void operator()( - const gsl::not_null*> dt_spatial_metric, - const gsl::not_null cache, + gsl::not_null*> dt_spatial_metric, + gsl::not_null cache, ::Tags::dt> /*meta*/) const; void operator()( - const gsl::not_null*> extrinsic_curvature, - const gsl::not_null cache, + gsl::not_null*> null_form_dot_deriv_H, + gsl::not_null cache, + internal_tags::null_form_dot_deriv_H /*meta*/) const; + + void operator()( + gsl::not_null*> + null_form_dot_deriv_null_form, + gsl::not_null cache, + internal_tags::null_form_dot_deriv_null_form /*meta*/) + const; + + void operator()( + gsl::not_null*> extrinsic_curvature, + gsl::not_null cache, gr::Tags::ExtrinsicCurvature /*meta*/) const; void operator()( - const gsl::not_null*> + gsl::not_null*> spatial_christoffel_first_kind, - const gsl::not_null cache, + gsl::not_null cache, gr::Tags::SpatialChristoffelFirstKind /*meta*/) const; void operator()( - const gsl::not_null*> + gsl::not_null*> spatial_christoffel_second_kind, - const gsl::not_null cache, + gsl::not_null cache, gr::Tags::SpatialChristoffelSecondKind /*meta*/) const; @@ -570,14 +654,18 @@ class KerrSchild : public AnalyticSolution<3_st>, make_array(std::numeric_limits::signaling_NaN()); std::array center_ = make_array(std::numeric_limits::signaling_NaN()); + std::array boost_velocity_ = + make_array(std::numeric_limits::signaling_NaN()); bool zero_spin_{}; + bool zero_velocity_{}; }; SPECTRE_ALWAYS_INLINE bool operator==(const KerrSchild& lhs, const KerrSchild& rhs) { return lhs.mass() == rhs.mass() and lhs.dimensionless_spin() == rhs.dimensionless_spin() and - lhs.center() == rhs.center(); + lhs.center() == rhs.center() and + lhs.boost_velocity() == rhs.boost_velocity(); } SPECTRE_ALWAYS_INLINE bool operator!=(const KerrSchild& lhs, diff --git a/src/PointwiseFunctions/GeneralRelativity/TimeDerivativeOfSpatialMetric.cpp b/src/PointwiseFunctions/GeneralRelativity/TimeDerivativeOfSpatialMetric.cpp index 36f344683f9d..d8f3459510a5 100644 --- a/src/PointwiseFunctions/GeneralRelativity/TimeDerivativeOfSpatialMetric.cpp +++ b/src/PointwiseFunctions/GeneralRelativity/TimeDerivativeOfSpatialMetric.cpp @@ -81,7 +81,7 @@ tnsr::ii time_derivative_of_spatial_metric( extrinsic_curvature); GENERATE_INSTANTIATIONS(INSTANTIATE, (1, 2, 3), (double, DataVector), - (Frame::Grid, Frame::Inertial)) + (Frame::Grid, Frame::Distorted, Frame::Inertial)) #undef DIM #undef DTYPE diff --git a/src/PointwiseFunctions/SpecialRelativity/LorentzBoostMatrix.cpp b/src/PointwiseFunctions/SpecialRelativity/LorentzBoostMatrix.cpp index c58eb504dab1..24e503e80b90 100644 --- a/src/PointwiseFunctions/SpecialRelativity/LorentzBoostMatrix.cpp +++ b/src/PointwiseFunctions/SpecialRelativity/LorentzBoostMatrix.cpp @@ -3,6 +3,7 @@ #include "PointwiseFunctions/SpecialRelativity/LorentzBoostMatrix.hpp" +#include #include #include #include @@ -58,6 +59,75 @@ void lorentz_boost_matrix( (*boost_matrix).get(i + 1, i + 1) += 1.0; } } + +template +tnsr::Ab lorentz_boost_matrix( + const std::array& velocity) { + tnsr::I velocity_as_tensor(velocity); + return lorentz_boost_matrix(velocity_as_tensor); +} + +template +void lorentz_boost( + const gsl::not_null*> result, + const tnsr::I& vector, + const double vector_component_0, + const std::array& velocity) { + if (velocity == make_array(0.)) { + *result = vector; + return; + } + // Inverse matrix with respect to the boost applied to one forms + const auto boost_matrix = lorentz_boost_matrix(-velocity); + for (size_t i = 0; i < SpatialDim; ++i) { + result->get(i) = boost_matrix.get(i + 1, 0) * vector_component_0; + for (size_t j = 0; j < SpatialDim; ++j) { + result->get(i) += boost_matrix.get(i + 1, j + 1) * vector.get(j); + } + } +} + +template +void lorentz_boost( + const gsl::not_null*> result, + const tnsr::a& one_form, + const std::array& velocity) { + if (velocity == make_array(0.)) { + *result = one_form; + return; + } + const auto boost_matrix = lorentz_boost_matrix(velocity); + for (size_t i = 0; i < SpatialDim + 1; ++i) { + result->get(i) = 0.; + for (size_t j = 0; j < SpatialDim + 1; ++j) { + result->get(i) += boost_matrix.get(i, j) * one_form.get(j); + } + } +} + +template +void lorentz_boost( + const gsl::not_null*> result, + const tnsr::ab& tensor, + const std::array& velocity) { + if (velocity == make_array(0.)) { + *result = tensor; + return; + } + const auto boost_matrix = lorentz_boost_matrix(velocity); + for (size_t i = 0; i < SpatialDim + 1; ++i) { + for (size_t k = 0; k < SpatialDim + 1; ++k) { + result->get(i, k) = 0.; + for (size_t j = 0; j < SpatialDim + 1; ++j) { + for (size_t l = 0; l < SpatialDim + 1; ++l) { + result->get(i, k) += boost_matrix.get(i, j) * boost_matrix.get(k, l) * + tensor.get(j, l); + } + } + } + } +} + } // namespace sr // Explicit Instantiations @@ -70,10 +140,40 @@ void lorentz_boost_matrix( template void sr::lorentz_boost_matrix( \ gsl::not_null*> \ boost_matrix, \ - const tnsr::I& velocity); + const tnsr::I& velocity); \ + template tnsr::Ab \ + sr::lorentz_boost_matrix(const std::array& velocity); GENERATE_INSTANTIATIONS(INSTANTIATE, (1, 2, 3)) +#undef DIM +#undef INSTANTIATE + +#define DIM(data) BOOST_PP_TUPLE_ELEM(0, data) +#define DTYPE(data) BOOST_PP_TUPLE_ELEM(1, data) +#define FRAME(data) BOOST_PP_TUPLE_ELEM(2, data) + +#define INSTANTIATE(_, data) \ + template void sr::lorentz_boost( \ + const gsl::not_null*> \ + result, \ + const tnsr::I& one_form, \ + const double one_form_component_0, \ + const std::array& velocity); \ + template void sr::lorentz_boost( \ + const gsl::not_null*> \ + result, \ + const tnsr::a& one_form, \ + const std::array& velocity); \ + template void sr::lorentz_boost( \ + const gsl::not_null*> \ + result, \ + const tnsr::ab& tensor, \ + const std::array& velocity); + +GENERATE_INSTANTIATIONS(INSTANTIATE, (1, 2, 3), (double, DataVector), + (Frame::Grid, Frame::Distorted, Frame::Inertial)) + #undef DIM #undef DTYPE #undef FRAME diff --git a/src/PointwiseFunctions/SpecialRelativity/LorentzBoostMatrix.hpp b/src/PointwiseFunctions/SpecialRelativity/LorentzBoostMatrix.hpp index 51f68bdc4fc4..2038a28f38f9 100644 --- a/src/PointwiseFunctions/SpecialRelativity/LorentzBoostMatrix.hpp +++ b/src/PointwiseFunctions/SpecialRelativity/LorentzBoostMatrix.hpp @@ -3,9 +3,12 @@ #pragma once +#include #include -#include "DataStructures/Tensor/TypeAliases.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "Utilities/Gsl.hpp" +#include "Utilities/MakeArray.hpp" /// \cond namespace gsl { @@ -53,5 +56,43 @@ template void lorentz_boost_matrix( gsl::not_null*> boost_matrix, const tnsr::I& velocity); + +template +tnsr::Ab lorentz_boost_matrix( + const std::array& velocity); +/// @} + +/// @{ +/*! + * \ingroup SpecialRelativityGroup + * \brief Apply a Lorentz boost to the spatial part of a vector. + * \details This requires passing the 0th component of the vector as an + * additional argument. + */ +template +void lorentz_boost(gsl::not_null*> result, + const tnsr::I& vector, + double vector_component_0, + const std::array& velocity); + +/*! + * \brief Apply a Lorentz boost to a one form. + */ +template +void lorentz_boost( + gsl::not_null*> result, + const tnsr::a& one_form, + const std::array& velocity); + +/*! + * \brief Apply a Lorentz boost to each component of a rank-2 tensor with + * lower or covariant indices. + * \note In the future we might want to write a single function capable to boost + * a tensor of arbitrary rank. + */ +template +void lorentz_boost(gsl::not_null*> result, + const tnsr::ab& tensor, + const std::array& velocity); /// @} } // namespace sr diff --git a/support/Pipelines/Bbh/InitialData.yaml b/support/Pipelines/Bbh/InitialData.yaml index 37cf84e04b6b..f9cea1ffd21a 100644 --- a/support/Pipelines/Bbh/InitialData.yaml +++ b/support/Pipelines/Bbh/InitialData.yaml @@ -34,6 +34,7 @@ Background: &background - {{ DimensionlessSpinLeft_y }} - {{ DimensionlessSpinLeft_z }} Center: [0., 0., 0.] + Velocity: [0., 0., 0.] ObjectRight: &kerr_right KerrSchild: Mass: &mass_right {{ MassRight }} @@ -42,6 +43,7 @@ Background: &background - {{ DimensionlessSpinRight_y }} - {{ DimensionlessSpinRight_z }} Center: [0., 0., 0.] + Velocity: [0., 0., 0.] AngularVelocity: {{ OrbitalAngularVelocity }} Expansion: {{ RadialExpansionVelocity }} LinearVelocity: [0., 0., 0.] diff --git a/tests/InputFiles/CurvedScalarWave/WorldtubeKerrSchild.yaml b/tests/InputFiles/CurvedScalarWave/WorldtubeKerrSchild.yaml index b0bb15108215..c4614c340028 100644 --- a/tests/InputFiles/CurvedScalarWave/WorldtubeKerrSchild.yaml +++ b/tests/InputFiles/CurvedScalarWave/WorldtubeKerrSchild.yaml @@ -26,6 +26,7 @@ BackgroundSpacetime: Mass: 1. Center: [0., 0., 0.] Spin: [0., 0., 0.] + Velocity: [0., 0., 0.] ResourceInfo: AvoidGlobalProc0: false diff --git a/tests/InputFiles/GeneralizedHarmonic/KerrSchild.yaml b/tests/InputFiles/GeneralizedHarmonic/KerrSchild.yaml index 1ff6103c4721..46b128c7f7dc 100644 --- a/tests/InputFiles/GeneralizedHarmonic/KerrSchild.yaml +++ b/tests/InputFiles/GeneralizedHarmonic/KerrSchild.yaml @@ -74,6 +74,7 @@ InitialData: &InitialData Mass: 1.0 Spin: [0.0, 0.0, 0.0] Center: &Center [0.0, 0.0, 0.0] + Velocity: [0.0, 0.0, 0.0] DomainCreator: Sphere: diff --git a/tests/InputFiles/Xcts/BinaryBlackHole.yaml b/tests/InputFiles/Xcts/BinaryBlackHole.yaml index 3427f2196ae1..9595723e40ff 100644 --- a/tests/InputFiles/Xcts/BinaryBlackHole.yaml +++ b/tests/InputFiles/Xcts/BinaryBlackHole.yaml @@ -84,11 +84,13 @@ Background: &background Mass: 0.4229 Spin: [0., 0., 0.] Center: [0., 0., 0.] + Velocity: [0., 0., 0.] ObjectRight: &kerr_right KerrSchild: Mass: 0.4229 Spin: [0., 0., 0.] Center: [0., 0., 0.] + Velocity: [0., 0., 0.] AngularVelocity: 0.0144 Expansion: 0. LinearVelocity: [0., 0., 0.] diff --git a/tests/InputFiles/Xcts/KerrSchild.yaml b/tests/InputFiles/Xcts/KerrSchild.yaml index e5d58ce020d9..f1d21eb32cae 100644 --- a/tests/InputFiles/Xcts/KerrSchild.yaml +++ b/tests/InputFiles/Xcts/KerrSchild.yaml @@ -71,6 +71,7 @@ Background: &solution Mass: &KerrMass 1. Spin: &KerrSpin [0., 0., 0.6] Center: &KerrCenter [0., 0., 0.] + Velocity: [0., 0., 0.] InitialGuess: Flatness diff --git a/tests/Unit/Elliptic/Systems/Xcts/BoundaryConditions/Test_ApparentHorizon.cpp b/tests/Unit/Elliptic/Systems/Xcts/BoundaryConditions/Test_ApparentHorizon.cpp index 2e505a41235c..20c2a6daf64f 100644 --- a/tests/Unit/Elliptic/Systems/Xcts/BoundaryConditions/Test_ApparentHorizon.cpp +++ b/tests/Unit/Elliptic/Systems/Xcts/BoundaryConditions/Test_ApparentHorizon.cpp @@ -411,7 +411,8 @@ void test_consistency_with_kerr(const bool compute_expansion) { std::array rotation = -0.5 * dimensionless_spin / horizon_kerrschild_radius; CAPTURE(rotation); - const KerrSchild solution{mass, dimensionless_spin, {{0., 0., 0.}}}; + const KerrSchild solution{ + mass, dimensionless_spin, {{0., 0., 0.}}, {{0., 0., 0.}}}; const ApparentHorizon kerr_horizon{ center, rotation, std::make_unique(solution), // Check with and without the negative-expansion condition. Either the @@ -579,8 +580,8 @@ SPECTRE_TEST_CASE("Unit.Xcts.BoundaryConditions.ApparentHorizon", " Lapse: Auto\n" " NegativeExpansion: None\n"); test_creation({{1., 2., 3.}}, {{0.1, 0.2, 0.3}}, - {{2.3, {{0.4, 0.5, 0.6}}, {{0., 0., 0.}}}}, - {{3.4, {{0.3, 0.2, 0.1}}, {{0., 0., 0.}}}}, + {{2.3, {{0.4, 0.5, 0.6}}, {{0., 0., 0.}}, {{0., 0., 0.}}}}, + {{3.4, {{0.3, 0.2, 0.1}}, {{0., 0., 0.}}, {{0., 0., 0.}}}}, "ApparentHorizon:\n" " Center: [1., 2., 3.]\n" " Rotation: [0.1, 0.2, 0.3]\n" @@ -589,11 +590,13 @@ SPECTRE_TEST_CASE("Unit.Xcts.BoundaryConditions.ApparentHorizon", " Mass: 2.3\n" " Spin: [0.4, 0.5, 0.6]\n" " Center: [0., 0., 0.]\n" + " Velocity: [0., 0., 0.]\n" " NegativeExpansion:\n" " KerrSchild:\n" " Mass: 3.4\n" " Spin: [0.3, 0.2, 0.1]\n" - " Center: [0., 0., 0.]\n"); + " Center: [0., 0., 0.]\n" + " Velocity: [0., 0., 0.]\n"); test_with_random_values(); test_consistency_with_kerr(false); test_consistency_with_kerr(true); diff --git a/tests/Unit/Evolution/Systems/CurvedScalarWave/Test_BackgroundSpacetime.cpp b/tests/Unit/Evolution/Systems/CurvedScalarWave/Test_BackgroundSpacetime.cpp index 4b5b451a0746..a7eb5a92adcd 100644 --- a/tests/Unit/Evolution/Systems/CurvedScalarWave/Test_BackgroundSpacetime.cpp +++ b/tests/Unit/Evolution/Systems/CurvedScalarWave/Test_BackgroundSpacetime.cpp @@ -37,7 +37,8 @@ SPECTRE_TEST_CASE("Unit.Evolution.Systems.CurvedScalarWave.BackgroundSpacetime", const auto ks_background = gr::Solutions::KerrSchild(0.5, {{0.1, 0.2, 0.3}}, {{1.0, 3.0, 2.0}}); const std::string& ks_option_string = - "Mass: 0.5\nSpin: [0.1,0.2,0.3]\nCenter: [1.0,3.0,2.0]"; + "Mass: 0.5\nSpin: [0.1,0.2,0.3]\nCenter: [1.0,3.0,2.0]\nVelocity: " + "[0.0,0.0,0.0]"; test_option_parsing(ks_background, ks_option_string); const auto spherical_ks_background = gr::Solutions::SphericalKerrSchild( diff --git a/tests/Unit/Evolution/Systems/GeneralizedHarmonic/Actions/Test_SetInitialData.cpp b/tests/Unit/Evolution/Systems/GeneralizedHarmonic/Actions/Test_SetInitialData.cpp index 997b71338bd9..e2064e10aaac 100644 --- a/tests/Unit/Evolution/Systems/GeneralizedHarmonic/Actions/Test_SetInitialData.cpp +++ b/tests/Unit/Evolution/Systems/GeneralizedHarmonic/Actions/Test_SetInitialData.cpp @@ -319,7 +319,8 @@ SPECTRE_TEST_CASE("Unit.Evolution.Systems.Gh.NumericInitialData", "GeneralizedHarmonic(KerrSchild):\n" " Mass: 1.\n" " Spin: [0, 0, 0]\n" - " Center: [0, 0, 0]", + " Center: [0, 0, 0]\n" + " Velocity: [0, 0, 0]", false); } diff --git a/tests/Unit/Evolution/Systems/GeneralizedHarmonic/GaugeSourceFunctions/Test_AnalyticChristoffel.cpp b/tests/Unit/Evolution/Systems/GeneralizedHarmonic/GaugeSourceFunctions/Test_AnalyticChristoffel.cpp index eb66caf5e16e..c89d1784a956 100644 --- a/tests/Unit/Evolution/Systems/GeneralizedHarmonic/GaugeSourceFunctions/Test_AnalyticChristoffel.cpp +++ b/tests/Unit/Evolution/Systems/GeneralizedHarmonic/GaugeSourceFunctions/Test_AnalyticChristoffel.cpp @@ -122,7 +122,8 @@ void test_ks(const Mesh<3>& mesh) { " GeneralizedHarmonic(KerrSchild):\n" " Mass: 1.2\n" " Spin: [0.1, 0.2, 0.3]\n" - " Center: [-0.1, -0.2, -0.4]\n") + " Center: [-0.1, -0.2, -0.4]\n" + " Velocity: [0.0, 0.0, 0.0]\n") ->get_clone()); CHECK_FALSE(gauge_condition->is_harmonic()); diff --git a/tests/Unit/Helpers/PointwiseFunctions/AnalyticSolutions/GeneralRelativity/VerifyGrSolution.hpp b/tests/Unit/Helpers/PointwiseFunctions/AnalyticSolutions/GeneralRelativity/VerifyGrSolution.hpp index 7ac0b7c7afee..1e3320dabc23 100644 --- a/tests/Unit/Helpers/PointwiseFunctions/AnalyticSolutions/GeneralRelativity/VerifyGrSolution.hpp +++ b/tests/Unit/Helpers/PointwiseFunctions/AnalyticSolutions/GeneralRelativity/VerifyGrSolution.hpp @@ -426,5 +426,55 @@ void verify_consistency(const Solution& solution, const double time, derivative_delta), derivative_approx); } + +/// Check the consistency of quantities dependent on other metric quantities +/// returned by a solution. This includes checking pointwise relations such as +/// consistency of the metric and inverse and comparing returned and +/// numerical derivatives. +template +void verify_spatial_consistency(const Solution& solution, const double time, + const tnsr::I& position, + const double derivative_delta, + const double derivative_tolerance) { + using Lapse = gr::Tags::Lapse; + using Shift = gr::Tags::Shift; + using SpatialMetric = gr::Tags::SpatialMetric; + using SqrtDetSpatialMetric = gr::Tags::SqrtDetSpatialMetric; + using InverseSpatialMetric = gr::Tags::InverseSpatialMetric; + using ExtrinsicCurvature = gr::Tags::ExtrinsicCurvature; + using tags = + tmpl::list, + Tags::dt, detail::deriv, + Tags::dt, Tags::dt, detail::deriv>; + + auto derivative_approx = approx.epsilon(derivative_tolerance); + + const auto vars = solution.variables(position, time, tags{}); + + const auto numerical_metric_det_and_inverse = + determinant_and_inverse(get(vars)); + CHECK_ITERABLE_APPROX(get(get(vars)), + sqrt(get(numerical_metric_det_and_inverse.first))); + CHECK_ITERABLE_APPROX(get(vars), + numerical_metric_det_and_inverse.second); + + CHECK_ITERABLE_CUSTOM_APPROX( + SINGLE_ARG(get>(vars)), + detail::space_derivative(solution, position, time, + derivative_delta), + derivative_approx); + CHECK_ITERABLE_CUSTOM_APPROX( + SINGLE_ARG(get>(vars)), + detail::space_derivative(solution, position, time, + derivative_delta), + derivative_approx); + CHECK_ITERABLE_CUSTOM_APPROX( + SINGLE_ARG(get>(vars)), + detail::space_derivative(solution, position, time, + derivative_delta), + derivative_approx); +} + } // namespace VerifyGrSolution } // namespace TestHelpers diff --git a/tests/Unit/PointwiseFunctions/AnalyticSolutions/GeneralRelativity/Test_KerrSchild.cpp b/tests/Unit/PointwiseFunctions/AnalyticSolutions/GeneralRelativity/Test_KerrSchild.cpp index 8f5ccc9dd93b..0410a7e40904 100644 --- a/tests/Unit/PointwiseFunctions/AnalyticSolutions/GeneralRelativity/Test_KerrSchild.cpp +++ b/tests/Unit/PointwiseFunctions/AnalyticSolutions/GeneralRelativity/Test_KerrSchild.cpp @@ -33,6 +33,7 @@ #include "NumericalAlgorithms/Spectral/Quadrature.hpp" #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrSchild.hpp" #include "PointwiseFunctions/GeneralRelativity/Tags.hpp" +#include "PointwiseFunctions/SpecialRelativity/LorentzBoostMatrix.hpp" #include "Utilities/ConstantExpressions.hpp" #include "Utilities/MakeWithValue.hpp" #include "Utilities/TMPL.hpp" @@ -307,39 +308,64 @@ void test_tag_retrieval(const DataType& used_for_size) { const double mass = 1.234; const std::array spin{{0.1, -0.2, 0.3}}; const std::array center{{1.0, 2.0, 3.0}}; + const std::array boost_velocity{{0.2, -0.3, 0.1}}; const auto x = spatial_coords(used_for_size); const double t = 1.3; + using all_tags = typename gr::Solutions::KerrSchild::tags; + + // We test gr::Tags::SpacetimeChristoffelSecondKind + // separately as it is not supported for the boosted solution. + using tags_for_zero_boost = + tmpl::list>; + using regular_tags = tmpl::remove< + all_tags, gr::Tags::SpacetimeChristoffelSecondKind>; + // Evaluate solution - const gr::Solutions::KerrSchild solution(mass, spin, center); - TestHelpers::AnalyticSolutions::test_tag_retrieval( - solution, x, t, - typename gr::Solutions::KerrSchild::template tags{}); + const gr::Solutions::KerrSchild solution(mass, spin, center, boost_velocity); + TestHelpers::AnalyticSolutions::test_tag_retrieval(solution, x, t, + regular_tags{}); + + const gr::Solutions::KerrSchild solution_zero_boost(mass, spin, center); + TestHelpers::AnalyticSolutions::test_tag_retrieval(solution_zero_boost, x, t, + tags_for_zero_boost{}); } template -void test_einstein_solution() { +void test_einstein_solution(const bool use_non_zero_velocity) { + INFO("Verify KerrSchild solution satisfies Einstein equations"); // Parameters // ...for KerrSchild solution const double mass = 1.7; const std::array spin{{0.1, 0.2, 0.3}}; const std::array center{{0.3, 0.2, 0.4}}; + const std::array boost_velocity{{0.4, -0.51, -0.5}}; // ...for grid const std::array lower_bound{{0.8, 1.22, 1.30}}; const double time = -2.8; - gr::Solutions::KerrSchild solution(mass, spin, center); - TestHelpers::VerifyGrSolution::verify_consistency( - solution, time, tnsr::I{lower_bound}, 0.01, 1.0e-10); - if constexpr (std::is_same_v) { - // Don't look at time-independent solution in other than the inertial - // frame. - const size_t grid_size = 8; - const std::array upper_bound{{0.82, 1.24, 1.32}}; - TestHelpers::VerifyGrSolution::verify_time_independent_einstein_solution( - solution, grid_size, lower_bound, upper_bound, - std::numeric_limits::epsilon() * 1.e5); + if (use_non_zero_velocity) { + gr::Solutions::KerrSchild solution(mass, spin, center, boost_velocity); + // Time derivatives are non zero. However we do not code up the time + // dependence of the boosted solution for simplicity. Hence we only check + // consistency in the spatial derivatives. + TestHelpers::VerifyGrSolution::verify_spatial_consistency( + solution, time, tnsr::I{lower_bound}, 0.01, 1.0e-10); + } else { + gr::Solutions::KerrSchild solution(mass, spin, center); + TestHelpers::VerifyGrSolution::verify_consistency( + solution, time, tnsr::I{lower_bound}, 0.01, 1.0e-10); + if constexpr (std::is_same_v) { + // Don't look at time-independent solution in other than the inertial + // frame. + const size_t grid_size = 8; + const std::array upper_bound{{0.82, 1.24, 1.32}}; + TestHelpers::VerifyGrSolution::verify_time_independent_einstein_solution( + solution, grid_size, lower_bound, upper_bound, + std::numeric_limits::epsilon() * 1.e5); + } } + } template @@ -366,13 +392,127 @@ void test_zero_spin_optimization(const DataType& used_for_size) { }); } +template +void test_boosted_for_zero_velocity(const DataType& used_for_size) { + const double epsilon = std::numeric_limits::epsilon(); + const std::array boost_velocity{{0.0, 0.0, epsilon}}; + + gr::Solutions::KerrSchild solution_zero_velocity( + 3.0, {{0., 0., 0.}}, {{0.2, 0.3, 0.2}}); + gr::Solutions::KerrSchild solution_tiny_velocity( + 3.0, {{0., 0., 1e-50}}, {{0.2, 0.3, 0.2}}, boost_velocity); + CHECK(solution_zero_velocity.zero_velocity()); + CHECK(not solution_tiny_velocity.zero_velocity()); + const auto x = spatial_coords(used_for_size); + + // We do not test gr::Tags::SpacetimeChristoffelSecondKind + // as it is not supported for the boosted solution. + using all_tags = typename gr::Solutions::KerrSchild::tags; + using tags_to_test = tmpl::remove< + all_tags, gr::Tags::SpacetimeChristoffelSecondKind>; + + const auto all_tags_zero_velocity = + solution_zero_velocity.variables(x, 0., tags_to_test{}); + const auto all_tags_tiny_velocity = + solution_tiny_velocity.variables(x, 0., tags_to_test{}); + tmpl::for_each( + [&all_tags_zero_velocity, &all_tags_tiny_velocity](auto tag_v) { + using tag = tmpl::type_from; + CHECK_ITERABLE_APPROX(get(all_tags_zero_velocity), + get(all_tags_tiny_velocity)); + }); +} + +template +void test_boosted_spacetime_metric(const DataType& used_for_size) { + // Test that the spacetime metric of the boosted solution indeed + // corresponds to the boosted spacetime metric tensor + + const auto x = spatial_coords(used_for_size); + const double t = 0.0; + + // KerrSchild solution + const double mass = 1.0; + const std::array spin{{0.0, 0.0, 0.0}}; + const std::array center{{0.0, 0.0, 0.0}}; + gr::Solutions::KerrSchild solution(mass, spin, center); + + // Boosted KerrSchild solution + const std::array velocity{{0.5, 0.0, 0.0}}; + gr::Solutions::KerrSchild boosted_solution(mass, spin, center, velocity); + + // Need to evaluate at the boosted coordinates + // Boost coordinates + tnsr::I x_boosted; + sr::lorentz_boost(make_not_null(&x_boosted), x, 0., -velocity); + + // We do not test gr::Tags::SpacetimeChristoffelSecondKind + // as it is not supported for the boosted solution. + using all_tags = typename gr::Solutions::KerrSchild::tags; + using tags_to_test = tmpl::remove< + all_tags, gr::Tags::SpacetimeChristoffelSecondKind>; + + const auto vars = solution.variables(x_boosted, t, tags_to_test{}); + const auto& lapse = get>(vars); + const auto& shift = get>(vars); + const auto& spatial_metric = + get>(vars); + + // Compute spacetime metric + tnsr::aa spacetime_metric; + gr::spacetime_metric(make_not_null(&spacetime_metric), lapse, shift, + spatial_metric); + + auto spacetime_metric_as_asymmetric_tensor = + make_with_value>(x, 0.0); + for (size_t a = 0; a < 4; ++a) { + for (size_t b = 0; b < 4; ++b) { + spacetime_metric_as_asymmetric_tensor.get(a, b) = + spacetime_metric.get(a, b); + } + } + + tnsr::ab expected_boosted_spacetime_metric; + sr::lorentz_boost(make_not_null(&expected_boosted_spacetime_metric), + spacetime_metric_as_asymmetric_tensor, -velocity); + + // Boosted metric quantities + const auto boosted_vars = boosted_solution.variables(x, t, tags_to_test{}); + const auto& lapse_from_boosted_solution = + get>(boosted_vars); + const auto& shift_from_boosted_solution = + get>(boosted_vars); + const auto& spatial_metric_from_boosted_solution = + get>(boosted_vars); + + tnsr::aa spacetime_metric_from_boosted_solution; + gr::spacetime_metric(make_not_null(&spacetime_metric_from_boosted_solution), + lapse_from_boosted_solution, shift_from_boosted_solution, + spatial_metric_from_boosted_solution); + + auto spacetime_metric_from_boosted_solution_as_asymmetric_tensor = + make_with_value>(x, 0.0); + for (size_t a = 0; a < 4; ++a) { + for (size_t b = 0; b < 4; ++b) { + spacetime_metric_from_boosted_solution_as_asymmetric_tensor.get(a, b) = + spacetime_metric_from_boosted_solution.get(a, b); + } + } + + CHECK_ITERABLE_APPROX( + expected_boosted_spacetime_metric, + spacetime_metric_from_boosted_solution_as_asymmetric_tensor); +} + void test_serialize() { - gr::Solutions::KerrSchild solution(3.0, {{0.2, 0.3, 0.2}}, {{0.0, 3.0, 4.0}}); + gr::Solutions::KerrSchild solution(3.0, {{0.2, 0.3, 0.2}}, {{0.0, 3.0, 4.0}}, + {{0.1, 0.2, 0.3}}); test_serialization(solution); } void test_copy_and_move() { - gr::Solutions::KerrSchild solution(3.0, {{0.2, 0.3, 0.2}}, {{0.0, 3.0, 4.0}}); + gr::Solutions::KerrSchild solution(3.0, {{0.2, 0.3, 0.2}}, {{0.0, 3.0, 4.0}}, + {{0.1, 0.2, 0.3}}); test_copy_semantics(solution); auto solution_copy = solution; // clang-tidy: std::move of trivially copyable type @@ -383,9 +523,11 @@ void test_construct_from_options() { const auto created = TestHelpers::test_creation( "Mass: 0.5\n" "Spin: [0.1,0.2,0.3]\n" - "Center: [1.0,3.0,2.0]"); - CHECK(created == - gr::Solutions::KerrSchild(0.5, {{0.1, 0.2, 0.3}}, {{1.0, 3.0, 2.0}})); + "Center: [1.0,3.0,2.0]\n" + "Velocity: [0.5,0.4,0.3]"); + CHECK(created == gr::Solutions::KerrSchild(0.5, {{0.1, 0.2, 0.3}}, + {{1.0, 3.0, 2.0}}, + {{0.5, 0.4, 0.3}})); } } // namespace @@ -401,18 +543,26 @@ SPECTRE_TEST_CASE("Unit.PointwiseFunctions.AnalyticSolutions.Gr.KerrSchild", test_numerical_deriv_det_spatial_metric(DataVector(5)); test_tag_retrieval(DataVector(5)); test_tag_retrieval(0.0); - test_einstein_solution(); + test_einstein_solution(true); + test_einstein_solution(false); test_zero_spin_optimization(DataVector(5)); test_zero_spin_optimization(0.0); + test_boosted_for_zero_velocity(DataVector(5)); + test_boosted_for_zero_velocity(0.0); + test_boosted_spacetime_metric(0.0); test_schwarzschild(DataVector(5)); test_schwarzschild(0.0); test_numerical_deriv_det_spatial_metric(DataVector(5)); test_tag_retrieval(DataVector(5)); test_tag_retrieval(0.0); - test_einstein_solution(); + test_einstein_solution(true); + test_einstein_solution(false); test_zero_spin_optimization(DataVector(5)); test_zero_spin_optimization(0.0); + test_boosted_for_zero_velocity(DataVector(5)); + test_boosted_for_zero_velocity(0.0); + test_boosted_spacetime_metric(0.0); CHECK_THROWS_WITH( []() { @@ -429,13 +579,15 @@ SPECTRE_TEST_CASE("Unit.PointwiseFunctions.AnalyticSolutions.Gr.KerrSchild", CHECK_THROWS_WITH(TestHelpers::test_creation( "Mass: -0.5\n" "Spin: [0.1,0.2,0.3]\n" - "Center: [1.0,3.0,2.0]"), + "Center: [1.0,3.0,2.0]\n" + "Velocity: [0,0,0]"), Catch::Matchers::ContainsSubstring( "Value -0.5 is below the lower bound of 0")); CHECK_THROWS_WITH( TestHelpers::test_creation( "Mass: 0.5\n" "Spin: [1.1,0.9,0.3]\n" - "Center: [1.0,3.0,2.0]"), + "Center: [1.0,3.0,2.0]\n" + "Velocity: [0,0,0]"), Catch::Matchers::ContainsSubstring("Spin magnitude must be < 1")); } diff --git a/tests/Unit/PointwiseFunctions/AnalyticSolutions/GeneralRelativity/Test_WrappedGr.cpp b/tests/Unit/PointwiseFunctions/AnalyticSolutions/GeneralRelativity/Test_WrappedGr.cpp index 6c9050e47654..2aaf7728803a 100644 --- a/tests/Unit/PointwiseFunctions/AnalyticSolutions/GeneralRelativity/Test_WrappedGr.cpp +++ b/tests/Unit/PointwiseFunctions/AnalyticSolutions/GeneralRelativity/Test_WrappedGr.cpp @@ -44,10 +44,9 @@ void compare_different_wrapped_solutions(const double mass, template void test_copy_and_move(const SolutionType& solution) { test_copy_semantics(solution); - auto solution_copy = solution; auto solution_copy2 = solution; // clang-tidy: std::move of trivially copyable type - test_move_semantics(std::move(solution_copy2), solution_copy); // NOLINT + test_move_semantics(std::move(solution_copy2), solution); // NOLINT } template @@ -136,8 +135,27 @@ void test_construct_from_options() { const double mass = 0.5; const std::array spin{{0.1, 0.2, 0.3}}; const std::array center{{1.0, 3.0, 2.0}}; + CHECK(created == gh::Solutions::WrappedGr(mass, spin, center)); } + +template <> +void test_construct_from_options() { + const auto created = TestHelpers::test_creation< + gh::Solutions::WrappedGr>( + "Mass: 0.5\n" + "Spin: [0.1,0.2,0.3]\n" + "Center: [1.0,3.0,2.0]\n" + "Velocity: [0.1,-0.3,0.5]"); + const double mass = 0.5; + const std::array spin{{0.1, 0.2, 0.3}}; + const std::array center{{1.0, 3.0, 2.0}}; + const std::array velocity{{0.1, -0.3, 0.5}}; + + CHECK(created == gh::Solutions::WrappedGr( + mass, spin, center, velocity)); +} + } // namespace SPECTRE_TEST_CASE("Unit.PointwiseFunctions.AnalyticSolutions.Gr.WrappedGr", diff --git a/tests/Unit/PointwiseFunctions/AnalyticSolutions/Xcts/Test_Kerr.cpp b/tests/Unit/PointwiseFunctions/AnalyticSolutions/Xcts/Test_Kerr.cpp index 1e70c4d4e416..06acde4ae11a 100644 --- a/tests/Unit/PointwiseFunctions/AnalyticSolutions/Xcts/Test_Kerr.cpp +++ b/tests/Unit/PointwiseFunctions/AnalyticSolutions/Xcts/Test_Kerr.cpp @@ -59,7 +59,8 @@ SPECTRE_TEST_CASE("Unit.PointwiseFunctions.AnalyticSolutions.Xcts.Kerr", "KerrSchild:\n" " Mass: 0.43\n" " Spin: [0.1, 0.2, 0.3]\n" - " Center: [1., 2., 3.]"); + " Center: [1., 2., 3.]\n" + " Velocity: [0., 0., 0.]"); } } // namespace Xcts::Solutions diff --git a/tests/Unit/PointwiseFunctions/SpecialRelativity/Test_LorentzBoostMatrix.cpp b/tests/Unit/PointwiseFunctions/SpecialRelativity/Test_LorentzBoostMatrix.cpp index 5cf56b76f3d4..3a2f005da375 100644 --- a/tests/Unit/PointwiseFunctions/SpecialRelativity/Test_LorentzBoostMatrix.cpp +++ b/tests/Unit/PointwiseFunctions/SpecialRelativity/Test_LorentzBoostMatrix.cpp @@ -60,6 +60,97 @@ void test_lorentz_boost_matrix_analytic(const double& velocity_squared) { } CHECK_ITERABLE_APPROX(inverse_check, identity_matrix); } + +template +void test_lorentz_boost(const std::array velocity) { + const DataVector used_for_size{3., 4., 5.}; + + MAKE_GENERATOR(generator); + std::uniform_real_distribution<> distribution(-0.2, 0.2); + + auto covariant_vector = + make_with_random_values>( + make_not_null(&generator), make_not_null(&distribution), + used_for_size); + + tnsr::a boosted_covariant_vector; + sr::lorentz_boost( + make_not_null(&boosted_covariant_vector), covariant_vector, velocity); + + // Check that applying the inverse boost returns the original vector + tnsr::a unboosted_covariant_vector; + sr::lorentz_boost( + make_not_null(&unboosted_covariant_vector), boosted_covariant_vector, + -velocity); + CHECK_ITERABLE_APPROX(unboosted_covariant_vector, covariant_vector); + + // Check boost of the spatial vector + auto contravariant_vector = + make_with_random_values>( + make_not_null(&generator), make_not_null(&distribution), + used_for_size); + tnsr::I spatial_vector = + make_with_value>(used_for_size, 0.0); + for (size_t i = 0; i < SpatialDim; ++i) { + spatial_vector.get(i) = contravariant_vector.get(i + 1); + } + // Lower index, i.e. v_0 + const double vector_component_0 = get<0>(contravariant_vector); + + tnsr::I boosted_spatial_vector = + make_with_value>(used_for_size, 0.0); + sr::lorentz_boost( + make_not_null(&boosted_spatial_vector), spatial_vector, + vector_component_0, velocity); + + // We don't have an overload of the Lorentz boost for 4d vectors (just for one + // forms). But we can just change the sign of the velocity to boost it with + // the inverse Lorentz matrix for testing + tnsr::a boosted_contravariant_vector; + sr::lorentz_boost( + make_not_null(&boosted_contravariant_vector), contravariant_vector, + -velocity); + tnsr::I expected_spatial_vector = + make_with_value>(used_for_size, 0.0); + for (size_t i = 0; i < SpatialDim; ++i) { + expected_spatial_vector.get(i) = boosted_contravariant_vector.get(i + 1); + } + CHECK_ITERABLE_APPROX(expected_spatial_vector, boosted_spatial_vector); + + // Check boost on rank-2 matrix by constructing one with an outer tensor + // product, i.e. T_{ab} = v_{a} u_{b} + tnsr::ab tensor = + make_with_value>(used_for_size, + 0.0); + tnsr::ab boosted_tensor = + make_with_value>(used_for_size, + 0.0); + tnsr::ab expected_tensor = + make_with_value>(used_for_size, + 0.0); + // We need a second covariant vector + auto covariant_vector_2 = + make_with_random_values>( + make_not_null(&generator), make_not_null(&distribution), + used_for_size); + // We boost it as well + tnsr::a boosted_covariant_vector_2; + sr::lorentz_boost( + make_not_null(&boosted_covariant_vector_2), covariant_vector_2, velocity); + + for (size_t i = 0; i < SpatialDim + 1; ++i) { + for (size_t j = 0; j < SpatialDim + 1; ++j) { + tensor.get(i, j) = covariant_vector.get(i) * covariant_vector_2.get(j); + expected_tensor.get(i, j) = + boosted_covariant_vector.get(i) * boosted_covariant_vector_2.get(j); + } + } + + sr::lorentz_boost(make_not_null(&boosted_tensor), + tensor, velocity); + CHECK_ITERABLE_APPROX(expected_tensor, boosted_tensor); +} + } // namespace SPECTRE_TEST_CASE( @@ -85,4 +176,7 @@ SPECTRE_TEST_CASE( d = large_velocity_squared; CHECK_FOR_DOUBLES(test_lorentz_boost_matrix_analytic, (1, 2, 3)); + + const std::array velocity{{0.1, -0.4, 0.3}}; + test_lorentz_boost(velocity); }