From e4a5d54600d33f46088ebbf90bc160a094432ddb Mon Sep 17 00:00:00 2001 From: Guillermo Lara Date: Wed, 1 Nov 2023 10:09:17 +0000 Subject: [PATCH] Add boost to KerrSchild Add more Lorentz boost functions Add instantiation of time deriv of spatial metric Add boosted Kerr initial data Expand doc Fix initialization of lorentz boost matrix Remove unused reference Add test of boosted spacetime metric Throw error for unsupported tag Co-authored-by: Nils Vu --- .../GeneralRelativity/CMakeLists.txt | 1 + .../GeneralRelativity/KerrSchild.cpp | 295 ++++++++++++++---- .../GeneralRelativity/KerrSchild.hpp | 276 ++++++++++------ .../TimeDerivativeOfSpatialMetric.cpp | 2 +- .../SpecialRelativity/LorentzBoostMatrix.cpp | 102 +++++- .../SpecialRelativity/LorentzBoostMatrix.hpp | 43 ++- support/Pipelines/Bbh/InitialData.yaml | 2 + .../CurvedScalarWave/WorldtubeKerrSchild.yaml | 1 + .../GeneralizedHarmonic/KerrSchild.yaml | 1 + tests/InputFiles/Xcts/BinaryBlackHole.yaml | 2 + tests/InputFiles/Xcts/KerrSchild.yaml | 1 + .../Test_ApparentHorizon.cpp | 11 +- .../Test_BackgroundSpacetime.cpp | 3 +- .../Actions/Test_SetInitialData.cpp | 3 +- .../Test_AnalyticChristoffel.cpp | 3 +- .../GeneralRelativity/VerifyGrSolution.hpp | 50 +++ .../GeneralRelativity/Test_KerrSchild.cpp | 202 ++++++++++-- .../GeneralRelativity/Test_WrappedGr.cpp | 22 +- .../AnalyticSolutions/Xcts/Test_Kerr.cpp | 3 +- .../Test_LorentzBoostMatrix.cpp | 94 ++++++ 20 files changed, 923 insertions(+), 194 deletions(-) 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); }