From 5f85c719db1d90204b78ac12c77f90e528eea629 Mon Sep 17 00:00:00 2001 From: Megum Rebelo Date: Sat, 4 Jan 2025 18:52:46 +0000 Subject: [PATCH 1/3] Start of BinaryWithGravitationalWaves class for XCTS - class to implement background data - superposition of two boosted black holes - no radiative content is added for now Additional BoundaryCondition class (Dirichlet): SuperposedBoostedBinary --- docs/References.bib | 76 ++ src/Elliptic/Executables/Xcts/SolveXcts.hpp | 4 + .../Xcts/BoundaryConditions/CMakeLists.txt | 2 + .../Xcts/BoundaryConditions/Factory.hpp | 14 +- .../SuperposedBoostedBinary.cpp | 300 +++++++ .../SuperposedBoostedBinary.hpp | 239 ++++++ .../Xcts/BinaryWithGravitationalWaves.cpp | 757 ++++++++++++++++++ .../Xcts/BinaryWithGravitationalWaves.hpp | 692 ++++++++++++++++ .../AnalyticData/Xcts/CMakeLists.txt | 3 + .../Xcts/BoundaryConditions/CMakeLists.txt | 1 + .../SuperposedBoostedBinary.py | 161 ++++ .../Test_SuperposedBoostedBinary.cpp | 181 +++++ .../Xcts/BinaryWithGravitationalWaves.py | 201 +++++ .../AnalyticData/Xcts/CMakeLists.txt | 1 + .../Test_BinaryWithGravitationalWaves.cpp | 152 ++++ 15 files changed, 2779 insertions(+), 5 deletions(-) create mode 100644 src/Elliptic/Systems/Xcts/BoundaryConditions/SuperposedBoostedBinary.cpp create mode 100644 src/Elliptic/Systems/Xcts/BoundaryConditions/SuperposedBoostedBinary.hpp create mode 100644 src/PointwiseFunctions/AnalyticData/Xcts/BinaryWithGravitationalWaves.cpp create mode 100644 src/PointwiseFunctions/AnalyticData/Xcts/BinaryWithGravitationalWaves.hpp create mode 100644 tests/Unit/Elliptic/Systems/Xcts/BoundaryConditions/SuperposedBoostedBinary.py create mode 100644 tests/Unit/Elliptic/Systems/Xcts/BoundaryConditions/Test_SuperposedBoostedBinary.cpp create mode 100644 tests/Unit/PointwiseFunctions/AnalyticData/Xcts/BinaryWithGravitationalWaves.py create mode 100644 tests/Unit/PointwiseFunctions/AnalyticData/Xcts/Test_BinaryWithGravitationalWaves.cpp diff --git a/docs/References.bib b/docs/References.bib index 660f9140838b..9a68a4766643 100644 --- a/docs/References.bib +++ b/docs/References.bib @@ -457,6 +457,20 @@ @article{Buchman:2012dw year = "2012" } +@article{Buonanno2005xu, + author = "Buonanno, Alessandra and Chen, Yanbei and Damour, Thibault", + title = "{Transition from inspiral to plunge in precessing binaries of + spinning black holes}", + eprint = "gr-qc/0508067", + archivePrefix = "arXiv", + doi = "10.1103/PhysRevD.74.104005", + journal = "Phys. Rev. D", + volume = "74", + pages = "104005", + year = "2006", + url = "https://doi.org/10.1103/PhysRevD.74.104005" +} + @article{Casoni2012, author = {Casoni, E. and Peraire, J. and Huerta, A.}, title = {One-dimensional shock-capturing for high-order discontinuous @@ -1345,6 +1359,21 @@ @Article{Hou2007 url = "https://doi.org/10.1007/s10915-006-9105-9" } +@article{Jaranowski1997ky, + author = "Jaranowski, Piotr and Schaefer, Gerhard", + title = "{Third postNewtonian higher order ADM Hamilton dynamics for + two-body point mass systems}", + eprint = "gr-qc/9712075", + archivePrefix = "arXiv", + doi = "10.1103/PhysRevD.57.7274", + url = "https://doi.org/10.1103/PhysRevD.57.7274", + journal = "Phys. Rev. D", + volume = "57", + pages = "7274--7291", + year = "1998", + note = "[Erratum: Phys.Rev.D 63, 029902 (2001)]" +} + @article{Jiang1996, author = {Guang-Shan Jiang and Chi-Wang Shu}, doi = {https://doi.org/10.1006/jcph.1996.0130}, @@ -1375,6 +1404,22 @@ @article{Kastaun2020uxr year = 2021 } +@article{Kelly2007uc, + author = "Kelly, Bernard J. and Tichy, Wolfgang and Campanelli, + Manuela and Whiting, Bernard F.", + title = "{Black hole puncture initial data with realistic gravitational + wave content}", + eprint = "0704.0628", + archivePrefix = "arXiv", + primaryClass = "gr-qc", + doi = "10.1103/PhysRevD.76.024008", + url = "https://doi.org/10.1103/PhysRevD.76.024008", + journal = "Phys. Rev. D", + volume = "76", + pages = "024008", + year = "2007" +} + @article{Kennedy2003, author = "Kennedy, Christopher A. and Carpenter, Mark H.", title = "Additive {Runge-Kutta} schemes for convection-diffusion-reaction @@ -2006,6 +2051,22 @@ @article{Muhlberger2014pja year = 2014 } +@article{Mundim2010hu, + author = "Mundim, Bruno C. and Kelly, Bernard J. and Zlochower, + Yosef and Nakano, Hiroyuki and Campanelli, Manuela", + editor = "Lehner, Luis and Pfeiffer, Harald P. and Poisson, E.", + title = "{Hybrid black-hole binary initial data}", + eprint = "1012.0886", + archivePrefix = "arXiv", + primaryClass = "gr-qc", + doi = "10.1088/0264-9381/28/13/134003", + url = "https://doi.org/10.1088/0264-9381/28/13/134003", + journal = "Class. Quant. Grav.", + volume = "28", + pages = "134003", + year = "2011" +} + @article{Nee2024bur, author = {Nee, Peter James and Lara, Guillermo and Pfeiffer, Harald P. and Vu, Nils L.}, @@ -2620,6 +2681,21 @@ @book{ThorneBlandford2017 publisher = "Princeton University Press", } +@article{Tichy2002ec, + author = "Tichy, Wolfgang and Bruegmann, Bernd and Campanelli, Manuela + and Diener, Peter", + title = "{Binary black hole initial data for numerical general relativity + based on postNewtonian data}", + eprint = "gr-qc/0207011", + archivePrefix = "arXiv", + reportNumber = "UTBRG-2002-06, AEI-2002-048", + doi = "10.1103/PhysRevD.67.064008", + journal = "Phys. Rev. D", + volume = "67", + pages = "064008", + year = "2003" +} + @article{Toro1994, author = "Toro, E. F. and Spruce, M. and Speares, W.", title = "{Restoration of the Contact Surface in the HLL–Riemann diff --git a/src/Elliptic/Executables/Xcts/SolveXcts.hpp b/src/Elliptic/Executables/Xcts/SolveXcts.hpp index 48161a5e0125..b9e49366da4f 100644 --- a/src/Elliptic/Executables/Xcts/SolveXcts.hpp +++ b/src/Elliptic/Executables/Xcts/SolveXcts.hpp @@ -43,6 +43,7 @@ #include "ParallelAlgorithms/LinearSolver/Multigrid/ElementsAllocator.hpp" #include "ParallelAlgorithms/LinearSolver/Multigrid/Tags.hpp" #include "PointwiseFunctions/AnalyticData/Xcts/Binary.hpp" +#include "PointwiseFunctions/AnalyticData/Xcts/BinaryWithGravitationalWaves.hpp" #include "PointwiseFunctions/AnalyticSolutions/Xcts/Factory.hpp" #include "PointwiseFunctions/Hydro/LowerSpatialFourVelocity.hpp" #include "PointwiseFunctions/Hydro/Tags.hpp" @@ -118,6 +119,9 @@ struct Metavariables { Xcts::Solutions::all_analytic_solutions, Xcts::AnalyticData::Binary, + Xcts::AnalyticData::BinaryWithGravitationalWaves< + elliptic::analytic_data::AnalyticSolution, + Xcts::Solutions::all_analytic_solutions>, elliptic::analytic_data::NumericData>; using factory_classes = tmpl::map< diff --git a/src/Elliptic/Systems/Xcts/BoundaryConditions/CMakeLists.txt b/src/Elliptic/Systems/Xcts/BoundaryConditions/CMakeLists.txt index 1f6a4ad21274..b97b6644b4e4 100644 --- a/src/Elliptic/Systems/Xcts/BoundaryConditions/CMakeLists.txt +++ b/src/Elliptic/Systems/Xcts/BoundaryConditions/CMakeLists.txt @@ -10,6 +10,7 @@ spectre_target_sources( PRIVATE ApparentHorizon.cpp Flatness.cpp + SuperposedBoostedBinary.cpp Robin.cpp ) @@ -20,6 +21,7 @@ spectre_target_headers( ApparentHorizon.hpp Factory.hpp Flatness.hpp + SuperposedBoostedBinary.hpp Robin.hpp ) diff --git a/src/Elliptic/Systems/Xcts/BoundaryConditions/Factory.hpp b/src/Elliptic/Systems/Xcts/BoundaryConditions/Factory.hpp index 71a85cabd3eb..88f4d6b5d0c7 100644 --- a/src/Elliptic/Systems/Xcts/BoundaryConditions/Factory.hpp +++ b/src/Elliptic/Systems/Xcts/BoundaryConditions/Factory.hpp @@ -7,13 +7,17 @@ #include "Elliptic/Systems/Xcts/BoundaryConditions/ApparentHorizon.hpp" #include "Elliptic/Systems/Xcts/BoundaryConditions/Flatness.hpp" #include "Elliptic/Systems/Xcts/BoundaryConditions/Robin.hpp" +#include "Elliptic/Systems/Xcts/BoundaryConditions/SuperposedBoostedBinary.hpp" +#include "PointwiseFunctions/AnalyticSolutions/Xcts/Factory.hpp" +#include "PointwiseFunctions/InitialDataUtilities/AnalyticSolution.hpp" #include "Utilities/TMPL.hpp" namespace Xcts::BoundaryConditions { template -using standard_boundary_conditions = - tmpl::list, - Flatness, - Robin, - ApparentHorizon>; +using standard_boundary_conditions = tmpl::list< + elliptic::BoundaryConditions::AnalyticSolution, + Flatness, Robin, + ApparentHorizon, + SuperposedBoostedBinary>; } diff --git a/src/Elliptic/Systems/Xcts/BoundaryConditions/SuperposedBoostedBinary.cpp b/src/Elliptic/Systems/Xcts/BoundaryConditions/SuperposedBoostedBinary.cpp new file mode 100644 index 000000000000..23f2e2c3dbc1 --- /dev/null +++ b/src/Elliptic/Systems/Xcts/BoundaryConditions/SuperposedBoostedBinary.cpp @@ -0,0 +1,300 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Elliptic/Systems/Xcts/BoundaryConditions/SuperposedBoostedBinary.hpp" + +#include + +#include +#include +#include +#include +#include + +#include "DataStructures/DataVector.hpp" +#include "DataStructures/Tensor/EagerMath/DeterminantAndInverse.hpp" +#include "DataStructures/Tensor/IndexType.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "DataStructures/Tensor/TypeAliases.hpp" +#include "Elliptic/BoundaryConditions/BoundaryConditionType.hpp" +#include "Elliptic/Systems/Xcts/Tags.hpp" +#include "PointwiseFunctions/AnalyticSolutions/Xcts/Factory.hpp" +#include "PointwiseFunctions/AnalyticSolutions/Xcts/Schwarzschild.hpp" +#include "PointwiseFunctions/GeneralRelativity/Lapse.hpp" +#include "PointwiseFunctions/GeneralRelativity/Shift.hpp" +#include "PointwiseFunctions/GeneralRelativity/SpacetimeMetric.hpp" +#include "PointwiseFunctions/GeneralRelativity/SpatialMetric.hpp" +#include "PointwiseFunctions/InitialDataUtilities/AnalyticSolution.hpp" +#include "PointwiseFunctions/SpecialRelativity/LorentzBoostMatrix.hpp" +#include "Utilities/CallWithDynamicType.hpp" +#include "Utilities/ConstantExpressions.hpp" +#include "Utilities/Gsl.hpp" +#include "Utilities/TaggedTuple.hpp" + +namespace Xcts::BoundaryConditions { + +namespace { + +template +void implement_apply_dirichlet( + const gsl::not_null*> conformal_factor_minus_one, + const gsl::not_null*> + lapse_times_conformal_factor_minus_one, + const gsl::not_null*> shift_excess, + const std::array>, 2>& + superposed_objects, + const std::array& xcoords, const std::array& masses, + const std::array& momentum_left, + const std::array& momentum_right, const double y_offset, + const double z_offset, const tnsr::I& x) { + using analytic_tags = + tmpl::list, + Xcts::Tags::LapseTimesConformalFactorMinusOne, + Xcts::Tags::ShiftExcess>; + + std::array, 2> x_isolated{{x, x}}; + const std::array, 2> coords_isolated{ + {{{xcoords[0], y_offset, z_offset}}, {{xcoords[1], y_offset, z_offset}}}}; + for (size_t i = 0; i < 2; ++i) { + for (size_t dim = 0; dim < 3; dim++) { + gsl::at(x_isolated, i).get(dim) -= + gsl::at(gsl::at(coords_isolated, i), dim); + } + } + std::array, 2> x_unboosted{{x, x}}; + sr::lorentz_boost(make_not_null(&(gsl::at(x_unboosted, 0))), + gsl::at(x_isolated, 0), 0., momentum_left / masses[0]); + sr::lorentz_boost(make_not_null(&(gsl::at(x_unboosted, 1))), + gsl::at(x_isolated, 1), 0., momentum_right / masses[1]); + + const auto solution_left = + call_with_dynamic_type, + IsolatedObjectClasses>( + (*(superposed_objects[0])).get(), + [&x_unboosted](const auto* const local_solution) { + return local_solution->variables(gsl::at(x_unboosted, 0), + analytic_tags{}); + }); + const auto solution_right = + call_with_dynamic_type, + IsolatedObjectClasses>( + (*(superposed_objects[1])).get(), + [&x_unboosted](const auto* const local_solution) { + return local_solution->variables(gsl::at(x_unboosted, 1), + analytic_tags{}); + }); + Scalar conformal_factor_minus_one_left = + get>(solution_left); + Scalar conformal_factor_minus_one_right = + get>(solution_right); + Scalar lapse_times_conformal_factor_minus_one_left = + get>( + solution_left); + Scalar lapse_times_conformal_factor_minus_one_right = + get>( + solution_right); + tnsr::I shift_excess_left = + get>( + solution_left); + tnsr::I shift_excess_right = + get>( + solution_right); + + // Boosted spacetime left + Scalar conformal_factor_left{get<0>(shift_excess_left).size()}; + get(conformal_factor_left) = 1. + get(conformal_factor_minus_one_left); + Scalar lapse_left{get<0>(shift_excess_left).size()}; + get(lapse_left) = (1. + get(lapse_times_conformal_factor_minus_one_left)) / + get(conformal_factor_left); + + tnsr::ii spatial_metric_left{get<0>(shift_excess_left).size()}; + std::fill(spatial_metric_left.begin(), spatial_metric_left.end(), 0.); + for (size_t i = 0; i < 3; ++i) { + spatial_metric_left.get(i, i) = square(square(get(conformal_factor_left))); + } + + tnsr::aa spacetime_metric_left{ + get<0>(shift_excess_left).size()}; + gr::spacetime_metric(make_not_null(&spacetime_metric_left), lapse_left, + shift_excess_left, spatial_metric_left); + + const tnsr::Ab lorentz_boost_matrix_double_left = + sr::lorentz_boost_matrix(-(momentum_left / masses[0])); + tnsr::Ab lorentz_boost_matrix_left{ + get<0>(shift_excess_left).size()}; + for (size_t i = 0; i < 4; ++i) { + for (size_t j = 0; j < 4; ++j) { + lorentz_boost_matrix_left.get(i, j) = + lorentz_boost_matrix_double_left.get(i, j); + } + } + + tnsr::aa spacetime_metric_boosted_left{ + get<0>(shift_excess_left).size()}; + for (size_t i = 0; i < 4; ++i) { + for (size_t j = 0; j <= i; ++j) { + spacetime_metric_boosted_left.get(i, j) = 0.; + for (size_t k = 0; k < 4; ++k) { + for (size_t l = 0; l < 4; ++l) { + spacetime_metric_boosted_left.get(i, j) += + lorentz_boost_matrix_left.get(k, i) * + lorentz_boost_matrix_left.get(l, j) * + spacetime_metric_left.get(k, l); + } + } + } + } + + // Boosted spacetime right + Scalar conformal_factor_right{get<0>(shift_excess_right).size()}; + get(conformal_factor_right) = 1. + get(conformal_factor_minus_one_right); + Scalar lapse_right{get<0>(shift_excess_right).size()}; + get(lapse_right) = (1. + get(lapse_times_conformal_factor_minus_one_right)) / + get(conformal_factor_right); + + tnsr::ii spatial_metric_right{ + get<0>(shift_excess_right).size()}; + std::fill(spatial_metric_right.begin(), spatial_metric_right.end(), 0.); + for (size_t i = 0; i < 3; ++i) { + spatial_metric_right.get(i, i) = + square(square(get(conformal_factor_right))); + } + + tnsr::aa spacetime_metric_right{ + get<0>(shift_excess_right).size()}; + gr::spacetime_metric(make_not_null(&spacetime_metric_right), lapse_right, + shift_excess_right, spatial_metric_right); + + const tnsr::Ab lorentz_boost_matrix_double_right = + sr::lorentz_boost_matrix(-(momentum_right / masses[1])); + tnsr::Ab lorentz_boost_matrix_right{ + get<0>(shift_excess_right).size()}; + for (size_t i = 0; i < 4; ++i) { + for (size_t j = 0; j < 4; ++j) { + lorentz_boost_matrix_right.get(i, j) = + lorentz_boost_matrix_double_right.get(i, j); + } + } + + tnsr::aa spacetime_metric_boosted_right{ + get<0>(shift_excess_right).size()}; + for (size_t i = 0; i < 4; ++i) { + for (size_t j = 0; j <= i; ++j) { + spacetime_metric_boosted_right.get(i, j) = 0.; + for (size_t k = 0; k < 4; ++k) { + for (size_t l = 0; l < 4; ++l) { + spacetime_metric_boosted_right.get(i, j) += + lorentz_boost_matrix_right.get(k, i) * + lorentz_boost_matrix_right.get(l, j) * + spacetime_metric_right.get(k, l); + } + } + } + } + + // Superposed lapse and shift (no longer conformally flat + // so conformal factor is 1. and "boosted spatial metric" should be on + // background) + get(*conformal_factor_minus_one) = 0.; + + const auto conformal_metric_left = + gr::spatial_metric(spacetime_metric_boosted_left); + const auto inv_conformal_metric_left = + determinant_and_inverse(conformal_metric_left).second; + const auto boosted_shift_left = + gr::shift(spacetime_metric_boosted_left, inv_conformal_metric_left); + const Scalar boosted_lapse_left = + gr::lapse(boosted_shift_left, spacetime_metric_boosted_left); + + const auto conformal_metric_right = + gr::spatial_metric(spacetime_metric_boosted_right); + const auto inv_conformal_metric_right = + determinant_and_inverse(conformal_metric_right).second; + const auto boosted_shift_right = + gr::shift(spacetime_metric_boosted_right, inv_conformal_metric_right); + const Scalar boosted_lapse_right = + gr::lapse(boosted_shift_right, spacetime_metric_boosted_right); + + get(*lapse_times_conformal_factor_minus_one) = + (get(boosted_lapse_left) * get(boosted_lapse_right) - 1.); + + for (size_t i = 0; i < 3; ++i) { + shift_excess->get(i) = + boosted_shift_left.get(i) + boosted_shift_right.get(i); + } +} + +} // namespace + +template +void SuperposedBoostedBinary::apply( + const gsl::not_null*> conformal_factor_minus_one, + const gsl::not_null*> + lapse_times_conformal_factor_minus_one, + const gsl::not_null*> shift_excess, + const gsl::not_null< + Scalar*> /*n_dot_conformal_factor_gradient*/, + const gsl::not_null*> + /*n_dot_lapse_times_conformal_factor_gradient*/, + const gsl::not_null*> + /*n_dot_longitudinal_shift_excess*/, + const tnsr::i& /*deriv_conformal_factor_correction1*/, + const tnsr::i& /*deriv_lapse_times_conformal_factor_correction1*/, + const tnsr::iJ& /*deriv_shift_excess_correction1*/, + const tnsr::I& x) const { + implement_apply_dirichlet( + conformal_factor_minus_one, lapse_times_conformal_factor_minus_one, + shift_excess, superposed_objects_, xcoords_, masses_, momentum_left_, + momentum_right_, y_offset_, z_offset_, x); +} + +template +void SuperposedBoostedBinary:: + apply_linearized( + const gsl::not_null*> conformal_factor_correction, + const gsl::not_null*> + lapse_times_conformal_factor_correction, + const gsl::not_null*> shift_excess_correction, + const gsl::not_null*> + /*n_dot_conformal_factor_gradient_correction*/, + const gsl::not_null*> + /*n_dot_lapse_times_conformal_factor_gradient_correction*/, + const gsl::not_null*> + /*n_dot_longitudinal_shift_excess_correction*/, + const tnsr::i& /*deriv_conformal_factor_correction*/, + const tnsr::i& + /*deriv_lapse_times_conformal_factor_correction*/, + const tnsr::iJ& /*deriv_shift_excess_correction*/) + const { + get(*conformal_factor_correction) = 0.; + get(*lapse_times_conformal_factor_correction) = 0.; + std::fill(shift_excess_correction->begin(), shift_excess_correction->end(), + 0.); +} + +template +bool operator==(const SuperposedBoostedBinary& /*lhs*/, + const SuperposedBoostedBinary& /*rhs*/) { + return true; +} + +template +bool operator!=(const SuperposedBoostedBinary& lhs, + const SuperposedBoostedBinary& rhs) { + return not(lhs == rhs); +} + +template class SuperposedBoostedBinary< + elliptic::analytic_data::AnalyticSolution, + Xcts::Solutions::all_analytic_solutions>; + +template class SuperposedBoostedBinary< + elliptic::analytic_data::AnalyticSolution, + tmpl::list>; + +} // namespace Xcts::BoundaryConditions diff --git a/src/Elliptic/Systems/Xcts/BoundaryConditions/SuperposedBoostedBinary.hpp b/src/Elliptic/Systems/Xcts/BoundaryConditions/SuperposedBoostedBinary.hpp new file mode 100644 index 000000000000..72f8216593fc --- /dev/null +++ b/src/Elliptic/Systems/Xcts/BoundaryConditions/SuperposedBoostedBinary.hpp @@ -0,0 +1,239 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include + +#include +#include +#include +#include +#include + +#include "DataStructures/Tensor/IndexType.hpp" +#include "DataStructures/Tensor/TypeAliases.hpp" +#include "Domain/BoundaryConditions/BoundaryCondition.hpp" +#include "Domain/Tags.hpp" +#include "Domain/Tags/FaceNormal.hpp" +#include "Elliptic/BoundaryConditions/BoundaryCondition.hpp" +#include "Elliptic/BoundaryConditions/BoundaryConditionType.hpp" +#include "Elliptic/Systems/Xcts/FluxesAndSources.hpp" +#include "Options/Context.hpp" +#include "Options/ParseError.hpp" +#include "Options/String.hpp" +#include "Utilities/Gsl.hpp" +#include "Utilities/Serialization/CharmPupable.hpp" + +/// \cond +class DataVector; +/// \endcond + +namespace Xcts::BoundaryConditions { + +/*! + * \brief Impose supperposed boosted binary system on the boundary. + * + * This takes two isolated objects and after applying a boost to each of them, + * superposes them. The superposed system is then imposed on the boundary, with + * Dirichlet boundary conditions. + * + */ + +template +class SuperposedBoostedBinary + : public elliptic::BoundaryConditions::BoundaryCondition<3> { + private: + using Base = elliptic::BoundaryConditions::BoundaryCondition<3>; + + public: + struct XCoords { + static constexpr Options::String help = + "The coordinates on the x-axis where the two objects are placed."; + using type = std::array; + }; + struct Masses { + static constexpr Options::String help = + "The masses of each object, first left and second right"; + using type = std::array; + }; + struct MomentumLeft { + static constexpr Options::String help = + "The momentum assigned to the left object."; + using type = std::array; + }; + struct MomentumRight { + static constexpr Options::String help = + "The momentum assigned to the right object."; + using type = std::array; + }; + struct CenterOfMassOffset { + static constexpr Options::String help = { + "Offset in the y and z axes applied to both objects in order to " + "control the center of mass."}; + using type = std::array; + }; + struct ObjectLeft { + static constexpr Options::String help = + "The object placed on the negative x-axis."; + using type = std::unique_ptr; + }; + struct ObjectRight { + static constexpr Options::String help = + "The object placed on the positive x-axis."; + using type = std::unique_ptr; + }; + + using options = tmpl::list; + static constexpr Options::String help = + "Impose supperposed boosted binary system on the boundary."; + + SuperposedBoostedBinary() = default; + SuperposedBoostedBinary(const SuperposedBoostedBinary&) = delete; + SuperposedBoostedBinary& operator=(const SuperposedBoostedBinary&) = delete; + SuperposedBoostedBinary(SuperposedBoostedBinary&&) = default; + SuperposedBoostedBinary& operator=(SuperposedBoostedBinary&&) = default; + ~SuperposedBoostedBinary() override = default; + + /// \cond + explicit SuperposedBoostedBinary(CkMigrateMessage* m) : Base(m) {} + using PUP::able::register_constructor; + WRAPPED_PUPable_decl_template(SuperposedBoostedBinary); + /// \endcond + + std::unique_ptr get_clone() + const override { + const std::array center_of_mass_offset = { + {y_offset_, z_offset_}}; + return std::make_unique( + xcoords_, masses_, momentum_left_, momentum_right_, + center_of_mass_offset, + superposed_objects_[0].has_value() + ? std::make_optional( + deserialize>( + serialize(superposed_objects_[0].value()).data())) + : std::nullopt, + superposed_objects_[1].has_value() + ? std::make_optional( + deserialize>( + serialize(superposed_objects_[1].value()).data())) + : std::nullopt); + } + + SuperposedBoostedBinary( + const std::array xcoords, const std::array masses, + const std::array momentum_left, + const std::array momentum_right, + const std::array center_of_mass_offset, + std::optional> object_left, + std::optional> object_right, + const Options::Context& context = {}) + : xcoords_(xcoords), + masses_(masses), + momentum_left_(momentum_left), + momentum_right_(momentum_right), + y_offset_(center_of_mass_offset[0]), + z_offset_(center_of_mass_offset[1]), + superposed_objects_({std::move(object_left), std::move(object_right)}) { + if (masses_[0] <= 0. || masses_[1] <= 0.) { + PARSE_ERROR(context, "The masses must be positive."); + } + if (xcoords_[0] >= xcoords_[1]) { + PARSE_ERROR(context, "Specify 'XCoords' ascending from left to right."); + } + } + + std::vector boundary_condition_types() + const override { + return {// Conformal factor + elliptic::BoundaryConditionType::Dirichlet, + // Lapse times conformal factor + elliptic::BoundaryConditionType::Dirichlet, + // Shift + elliptic::BoundaryConditionType::Dirichlet, + elliptic::BoundaryConditionType::Dirichlet, + elliptic::BoundaryConditionType::Dirichlet}; + } + + using argument_tags = + tmpl::flatten>>; + using volume_tags = tmpl::list<>; + + void apply( + gsl::not_null*> conformal_factor_minus_one, + gsl::not_null*> lapse_times_conformal_factor_minus_one, + gsl::not_null*> shift_excess, + gsl::not_null*> n_dot_conformal_factor_gradient, + gsl::not_null*> + n_dot_lapse_times_conformal_factor_gradient, + gsl::not_null*> n_dot_longitudinal_shift_excess, + const tnsr::i& deriv_conformal_factor_correction, + const tnsr::i& + deriv_lapse_times_conformal_factor_correction, + const tnsr::iJ& deriv_shift_excess_correction, + const tnsr::I& x) const; + + using argument_tags_linearized = tmpl::list<>; + using volume_tags_linearized = tmpl::list<>; + + void apply_linearized( + gsl::not_null*> conformal_factor_correction, + gsl::not_null*> + lapse_times_conformal_factor_correction, + gsl::not_null*> shift_excess_correction, + gsl::not_null*> + n_dot_conformal_factor_gradient_correction, + gsl::not_null*> + n_dot_lapse_times_conformal_factor_gradient_correction, + gsl::not_null*> + n_dot_longitudinal_shift_excess_correction, + const tnsr::i& deriv_conformal_factor_correction, + const tnsr::i& + deriv_lapse_times_conformal_factor_correction, + const tnsr::iJ& deriv_shift_excess_correction) const; + + // NOLINTNEXTLINE(google-runtime-references) + void pup(PUP::er& p) override { + Base::pup(p); + p | xcoords_; + p | masses_; + p | y_offset_; + p | z_offset_; + p | momentum_left_; + p | momentum_right_; + p | superposed_objects_; + } + + private: + std::array xcoords_{}; + std::array masses_{}; + std::array momentum_left_{}; + std::array momentum_right_{}; + double y_offset_{}; + double z_offset_{}; + std::array>, 2> + superposed_objects_{}; +}; + +template +bool operator==(const SuperposedBoostedBinary& lhs, + const SuperposedBoostedBinary& rhs); + +template +bool operator!=(const SuperposedBoostedBinary& lhs, + const SuperposedBoostedBinary& rhs); + +/// \cond +template +PUP::able::PUP_ID + SuperposedBoostedBinary::my_PUP_ID = // NOLINT + 0; // NOLINT +/// \endcond + +} // namespace Xcts::BoundaryConditions diff --git a/src/PointwiseFunctions/AnalyticData/Xcts/BinaryWithGravitationalWaves.cpp b/src/PointwiseFunctions/AnalyticData/Xcts/BinaryWithGravitationalWaves.cpp new file mode 100644 index 000000000000..0da205119327 --- /dev/null +++ b/src/PointwiseFunctions/AnalyticData/Xcts/BinaryWithGravitationalWaves.cpp @@ -0,0 +1,757 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "PointwiseFunctions/AnalyticData/Xcts/BinaryWithGravitationalWaves.hpp" + +#include + +#include +#include + +#include "DataStructures/DataBox/Prefixes.hpp" +#include "DataStructures/DataVector.hpp" +#include "DataStructures/Tensor/EagerMath/DeterminantAndInverse.hpp" +#include "DataStructures/Tensor/EagerMath/Magnitude.hpp" +#include "DataStructures/Tensor/EagerMath/RaiseOrLowerIndex.hpp" +#include "DataStructures/Tensor/EagerMath/Trace.hpp" +#include "DataStructures/Tensor/Expressions/TensorIndex.hpp" +#include "DataStructures/Tensor/IndexType.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "DataStructures/Tensor/TypeAliases.hpp" +#include "Elliptic/Systems/Xcts/Tags.hpp" +#include "NumericalAlgorithms/Spectral/Mesh.hpp" +#include "PointwiseFunctions/AnalyticData/Xcts/CommonVariables.hpp" +#include "PointwiseFunctions/GeneralRelativity/Christoffel.hpp" +#include "PointwiseFunctions/GeneralRelativity/Lapse.hpp" +#include "PointwiseFunctions/GeneralRelativity/Shift.hpp" +#include "PointwiseFunctions/GeneralRelativity/SpacetimeMetric.hpp" +#include "PointwiseFunctions/GeneralRelativity/SpatialMetric.hpp" +#include "PointwiseFunctions/GeneralRelativity/Tags.hpp" +#include "PointwiseFunctions/SpecialRelativity/LorentzBoostMatrix.hpp" +#include "Utilities/ConstantExpressions.hpp" +#include "Utilities/ContainerHelpers.hpp" +#include "Utilities/ErrorHandling/Assert.hpp" +#include "Utilities/ErrorHandling/Error.hpp" +#include "Utilities/Gsl.hpp" + +namespace Xcts::AnalyticData::detail { + +template +void BinaryWithGravitationalWavesVariables::operator()( + const gsl::not_null*> conformal_metric, + const gsl::not_null /*cache*/, + Xcts::Tags::ConformalMetric /*meta*/) const { + const auto& conformal_metric_aux = get_t_conformal_metric(present_time); + for (size_t i = 0; i < 3; ++i) { + for (size_t j = 0; j <= i; ++j) { + conformal_metric->get(i, j) = conformal_metric_aux.get(i, j); + } + } +} + +template +void BinaryWithGravitationalWavesVariables::operator()( + const gsl::not_null*> deriv_conformal_metric, + const gsl::not_null cache, + ::Tags::deriv, + tmpl::size_t, Frame::Inertial> /*meta*/) const { + ASSERT(mesh.has_value() and inv_jacobian.has_value(), + "Need a mesh and a Jacobian for numeric differentiation."); + if constexpr (std::is_same_v) { + const auto& conformal_metric = cache->get_var( + *this, Xcts::Tags::ConformalMetric{}); + partial_derivative(deriv_conformal_metric, conformal_metric, mesh->get(), + inv_jacobian->get()); + } else { + (void)deriv_conformal_metric; + (void)cache; + ERROR( + "Numeric differentiation only works with DataVectors because it needs " + "a grid."); + } +} + +template +void BinaryWithGravitationalWavesVariables::operator()( + const gsl::not_null*> trace_extrinsic_curvature, + const gsl::not_null /*cache*/, + gr::Tags::TraceExtrinsicCurvature /*meta*/) const { + ASSERT(mesh.has_value() and inv_jacobian.has_value(), + "Need a mesh and a Jacobian for numeric differentiation."); + if constexpr (std::is_same_v) { + get(*trace_extrinsic_curvature) = get(get_t_trace_extrinsic_curvature( + present_time, mesh->get(), inv_jacobian->get())); + } else { + (void)trace_extrinsic_curvature; + ERROR( + "Numeric differentiation only works with DataVectors because it needs " + "a grid."); + } +} + +template +void BinaryWithGravitationalWavesVariables::operator()( + const gsl::not_null*> dt_trace_extrinsic_curvature, + const gsl::not_null cache, + ::Tags::dt> /*meta*/) const { + const auto& trace_extrinsic_curvature = + cache->get_var(*this, gr::Tags::TraceExtrinsicCurvature{}); + DataType time_back(get(trace_extrinsic_curvature).size(), -time_displacement); + DataType time_back_two(get(trace_extrinsic_curvature).size(), + -2. * time_displacement); + DataType time_back_three(get(trace_extrinsic_curvature).size(), + -3. * time_displacement); + Scalar trace_extrinsic_curvature_back; + Scalar trace_extrinsic_curvature_back_two; + Scalar trace_extrinsic_curvature_back_three; + ASSERT(mesh.has_value() and inv_jacobian.has_value(), + "Need a mesh and a Jacobian for numeric differentiation."); + if constexpr (std::is_same_v) { + trace_extrinsic_curvature_back = get_t_trace_extrinsic_curvature( + time_back, mesh->get(), inv_jacobian->get()); + trace_extrinsic_curvature_back_two = get_t_trace_extrinsic_curvature( + time_back_two, mesh->get(), inv_jacobian->get()); + trace_extrinsic_curvature_back_three = get_t_trace_extrinsic_curvature( + time_back_three, mesh->get(), inv_jacobian->get()); + } else { + (void)dt_trace_extrinsic_curvature; + (void)cache; + ERROR( + "Numeric differentiation only works with DataVectors because it needs " + "a grid."); + } + // Third order + get(*dt_trace_extrinsic_curvature) = + (11. * get(trace_extrinsic_curvature) - + 18. * get(trace_extrinsic_curvature_back) + + 9. * get(trace_extrinsic_curvature_back_two) - + 2. * get(trace_extrinsic_curvature_back_three)) / + (6. * time_displacement); +} + +template +void BinaryWithGravitationalWavesVariables::operator()( + const gsl::not_null*> + longitudinal_shift_background_minus_dt_conformal_metric, + const gsl::not_null cache, + Xcts::Tags::LongitudinalShiftBackgroundMinusDtConformalMetric< + DataType, 3, Frame::Inertial> /*meta*/) const { + std::fill(longitudinal_shift_background_minus_dt_conformal_metric->begin(), + longitudinal_shift_background_minus_dt_conformal_metric->end(), 0.); + // LongitudinalShiftBackground is zero + // DtConformalMetric (finite difference 3rd order) + const auto& conformal_metric = cache->get_var( + *this, Xcts::Tags::ConformalMetric{}); + const auto& inv_conformal_metric = cache->get_var( + *this, + ::Xcts::Tags::InverseConformalMetric{}); + DataType time_back(get_size(x.get(0)), -time_displacement); + DataType time_back_two(get_size(x.get(0)), -2. * time_displacement); + DataType time_back_three(get_size(x.get(0)), -3. * time_displacement); + const auto conformal_metric_back = get_t_conformal_metric(time_back); + const auto conformal_metric_back_two = get_t_conformal_metric(time_back_two); + const auto conformal_metric_back_three = + get_t_conformal_metric(time_back_three); + + tnsr::ii dt_conformal_metric{get_size(x.get(0))}; + for (size_t i = 0; i < 3; ++i) { + for (size_t j = 0; j <= i; ++j) { + // Third order + dt_conformal_metric.get(i, j) = + (11. * conformal_metric.get(i, j) - + 18. * conformal_metric_back.get(i, j) + + 9. * conformal_metric_back_two.get(i, j) - + 2. * conformal_metric_back_three.get(i, j)) / + (6. * time_displacement); + } + } + + for (size_t i = 0; i < 3; ++i) { + for (size_t j = 0; j <= i; ++j) { + for (size_t k = 0; k < 3; ++k) { + for (size_t l = 0; l < 3; ++l) { + longitudinal_shift_background_minus_dt_conformal_metric->get(i, j) -= + inv_conformal_metric.get(i, k) * inv_conformal_metric.get(j, l) * + (dt_conformal_metric.get(k, l) - + (1. / 3.) * + get(trace(dt_conformal_metric, inv_conformal_metric)) * + conformal_metric.get(k, l)); + } + } + } + } +} + +template +void BinaryWithGravitationalWavesVariables::operator()( + const gsl::not_null*> conformal_factor_minus_one, + const gsl::not_null /*cache*/, + Xcts::Tags::ConformalFactorMinusOne /*meta*/) const { + get(*conformal_factor_minus_one) = 0.; +} + +template +void BinaryWithGravitationalWavesVariables::operator()( + const gsl::not_null*> + lapse_times_conformal_factor_minus_one, + const gsl::not_null cache, + Xcts::Tags::LapseTimesConformalFactorMinusOne /*meta*/) const { + const auto& conformal_factor_minus_one = + cache->get_var(*this, Xcts::Tags::ConformalFactorMinusOne{}); + const auto lapse = get_t_lapse(present_time); + get(*lapse_times_conformal_factor_minus_one) = + get(lapse) * (get(conformal_factor_minus_one) + 1.) - 1.; +} + +template +void BinaryWithGravitationalWavesVariables::operator()( + const gsl::not_null*> shift_excess, + const gsl::not_null /*cache*/, + Xcts::Tags::ShiftExcess /*meta*/) const { + std::fill(shift_excess->begin(), shift_excess->end(), 0.); + const auto shift_excess_aux = get_t_shift(present_time); + for (size_t i = 0; i < 3; ++i) { + shift_excess->get(i) += shift_excess_aux.get(i); + } +} + +// Private functions + +template +Scalar BinaryWithGravitationalWavesVariables:: + get_t_trace_extrinsic_curvature( + DataType t, Mesh<3> local_mesh, + InverseJacobian + local_inv_jacobian) const { + tnsr::ii extrinsic_curvature{t.size()}; + std::fill(extrinsic_curvature.begin(), extrinsic_curvature.end(), 0.); + + DataType time_back(get_size(x.get(0)), t[0] - time_displacement); + DataType time_back_two(get_size(x.get(0)), t[0] - 2. * time_displacement); + DataType time_back_three(get_size(x.get(0)), t[0] - 3. * time_displacement); + const auto conformal_metric = get_t_conformal_metric(t); + const auto inv_conformal_metric = + determinant_and_inverse(conformal_metric).second; + const auto conformal_metric_back = get_t_conformal_metric(time_back); + const auto conformal_metric_back_two = get_t_conformal_metric(time_back_two); + const auto conformal_metric_back_three = + get_t_conformal_metric(time_back_three); + const auto lapse = get_t_lapse(t); + + for (size_t i = 0; i < 3; ++i) { + for (size_t j = 0; j <= i; ++j) { + // Third order + extrinsic_curvature.get(i, j) += + (11. * conformal_metric.get(i, j) - + 18. * conformal_metric_back.get(i, j) + + 9. * conformal_metric_back_two.get(i, j) - + 2. * conformal_metric_back_three.get(i, j)) / + (6. * time_displacement * (-2.) * get(lapse)); + } + } + + const auto shift = get_t_shift(t); + const auto shift_down = raise_or_lower_index(shift, conformal_metric); + const auto deriv_shift = + partial_derivative(shift_down, local_mesh, local_inv_jacobian); + const auto deriv_conformal_metric = + partial_derivative(conformal_metric, local_mesh, local_inv_jacobian); + const auto christoffel_second_kind = + gr::christoffel_second_kind(deriv_conformal_metric, inv_conformal_metric); + const auto covariant_deriv_shift_contribution = tenex::evaluate( + deriv_shift(ti::i, ti::j) - + christoffel_second_kind(ti::K, ti::i, ti::j) * shift_down(ti::k) + + deriv_shift(ti::j, ti::i) - + christoffel_second_kind(ti::K, ti::j, ti::i) * shift_down(ti::k)); + + for (size_t i = 0; i < 3; ++i) { + for (size_t j = 0; j <= i; ++j) { + extrinsic_curvature.get(i, j) += + covariant_deriv_shift_contribution.get(i, j) / (2. * get(lapse)); + } + } + return trace(extrinsic_curvature, inv_conformal_metric); +} + +template +tnsr::ii +BinaryWithGravitationalWavesVariables::get_t_conformal_metric( + DataType t) const { + const auto radiative_term = get_t_radiative_term(t); + const auto superposed_spacetime_metric_t = + get_t_superposed_spacetime_metric(t); + auto conformal_metric = gr::spatial_metric(superposed_spacetime_metric_t); + for (size_t i = 0; i < 3; ++i) { + for (size_t j = 0; j <= i; ++j) { + conformal_metric.get(i, j) += radiative_term.get(i, j); + } + } + return conformal_metric; +} + +template +tnsr::ii +BinaryWithGravitationalWavesVariables::get_t_radiative_term( + DataType t) const { + const auto distance_left_t = get_t_distance_left(t); + const auto distance_right_t = get_t_distance_right(t); + const auto near_zone_term_t = get_t_near_zone_term(t); + const auto present_term_t = get_t_present_term(t); + const auto past_term_t = get_t_past_term(t); + const auto integral_term_t = get_t_integral_term(t); + double turn_off = .5; + if (attenuation_parameter == 0) { + turn_off = 1.; + } + tnsr::ii radiative_term_t{t.size()}; + for (size_t i = 0; i < 3; ++i) { + for (size_t j = 0; j <= i; ++j) { + radiative_term_t.get(i, j) = + (turn_off + .5 * tanh(attenuation_parameter * + (get(distance_left_t) - attenuation_radius))) * + (turn_off + .5 * tanh(attenuation_parameter * + (get(distance_right_t) - attenuation_radius))) * + (near_zone_term_t.get(i, j) + present_term_t.get(i, j) + + past_term_t.get(i, j) + integral_term_t.get(i, j)); + } + } + return radiative_term_t; +} + +template +tnsr::ii +BinaryWithGravitationalWavesVariables::get_t_near_zone_term( + DataType t) const { + // Computation will be added in the future + tnsr::ii near_zone_term_t{t.size()}; + std::fill(near_zone_term_t.begin(), near_zone_term_t.end(), 0.); + return near_zone_term_t; +} + +template +tnsr::ii +BinaryWithGravitationalWavesVariables::get_t_present_term( + DataType t) const { + // Computation will be added in the future + tnsr::ii present_term_t{t.size()}; + std::fill(present_term_t.begin(), present_term_t.end(), 0.); + return present_term_t; +} + +template +tnsr::ii +BinaryWithGravitationalWavesVariables::get_t_past_term( + DataType t) const { + // Computation will be added in the future + tnsr::ii past_term_t{t.size()}; + std::fill(past_term_t.begin(), past_term_t.end(), 0.); + return past_term_t; +} + +template +tnsr::ii +BinaryWithGravitationalWavesVariables::get_t_integral_term( + DataType t) const { + // Computation will be added in the future + tnsr::ii integral_term_t{t.size()}; + std::fill(integral_term_t.begin(), integral_term_t.end(), 0.); + return integral_term_t; +} + +template +Scalar BinaryWithGravitationalWavesVariables::get_t_lapse( + DataType t) const { + Scalar lapse_t{t.size()}; + std::fill(lapse_t.begin(), lapse_t.end(), 0.); + + // PN Lapse in the far zone + const auto distance_left_t = get_t_distance_left(t); + const auto distance_right_t = get_t_distance_right(t); + const auto separation_t = get_t_separation(t); + const auto momentum_left_t = get_t_momentum_left(t); + const auto momentum_right_t = get_t_momentum_right(t); + get(lapse_t) = + 1 - mass_left / get(distance_left_t) - + mass_right / get(distance_right_t) + + (mass_left * mass_right) / + (get(distance_left_t) * get(distance_right_t)) + + (mass_left * mass_right) / (get(distance_left_t) * get(separation_t)) + + (mass_left * mass_right) / (get(distance_right_t) * get(separation_t)) + + 0.5 * (square(mass_left / get(distance_left_t)) + + square(mass_right / get(distance_right_t)) - + 3. * (get(dot_product(momentum_left_t, momentum_left_t)) / + (get(distance_left_t) * mass_left) + + get(dot_product(momentum_right_t, momentum_right_t)) / + (get(distance_right_t) * mass_right))); + double turn_off = .5; + if (attenuation_parameter == 0) { + turn_off = 1.; + } + get(lapse_t) *= + (turn_off + .5 * tanh(attenuation_parameter * + (get(distance_left_t) - attenuation_radius))) * + (turn_off + .5 * tanh(attenuation_parameter * + (get(distance_right_t) - attenuation_radius))); + + // Boosted lapse in the inner zone + const auto spacetime_metric_left = get_t_boosted_spacetime_metric_left(t); + const auto conformal_metric_left = gr::spatial_metric(spacetime_metric_left); + const auto inv_conformal_metric_left = + determinant_and_inverse(conformal_metric_left).second; + const auto shift_left = + gr::shift(spacetime_metric_left, inv_conformal_metric_left); + const auto lapse_left = gr::lapse(shift_left, spacetime_metric_left); + const auto spacetime_metric_right = get_t_boosted_spacetime_metric_right(t); + const auto conformal_metric_right = + gr::spatial_metric(spacetime_metric_right); + const auto inv_conformal_metric_right = + determinant_and_inverse(conformal_metric_right).second; + const auto shift_right = + gr::shift(spacetime_metric_right, inv_conformal_metric_right); + const auto lapse_right = gr::lapse(shift_right, spacetime_metric_right); + get(lapse_t) += + (1. - + (turn_off + .5 * tanh(attenuation_parameter * + (get(distance_left_t) - attenuation_radius))) * + (turn_off + + .5 * tanh(attenuation_parameter * + (get(distance_right_t) - attenuation_radius)))) * + get(lapse_left) * get(lapse_right); + return lapse_t; +} + +template +tnsr::I +BinaryWithGravitationalWavesVariables::get_t_shift(DataType t) const { + tnsr::I shift_t{t.size()}; + std::fill(shift_t.begin(), shift_t.end(), 0.); + + // PN shift in the far zone + const auto distance_left_t = get_t_distance_left(present_time); + const auto distance_right_t = get_t_distance_right(present_time); + const auto separation_t = get_t_separation(present_time); + const auto momentum_left_t = get_t_momentum_left(present_time); + const auto momentum_right_t = get_t_momentum_right(present_time); + const auto normal_left_t = get_t_normal_left(present_time); + const auto normal_right_t = get_t_normal_right(present_time); + for (size_t i = 0; i < 3; ++i) { + shift_t.get(i) -= 4. * (momentum_left_t.get(i) / get(distance_left_t) + + momentum_right_t.get(i) / get(distance_right_t)); + for (size_t j = 0; j < 3; ++j) { + shift_t.get(i) += .5 * momentum_left_t.get(j) * + (-normal_left_t.get(i) * normal_left_t.get(j) / + get(distance_left_t)) + + .5 * momentum_right_t.get(j) * + (-normal_right_t.get(i) * normal_right_t.get(j) / + get(distance_right_t)); + } + shift_t.get(i) += 0.5 * momentum_left_t.get(i) / get(distance_left_t) + + 0.5 * momentum_right_t.get(i) / get(distance_right_t); + } + + double turn_off = .5; + if (attenuation_parameter == 0) { + turn_off = 1.; + } + for (size_t i = 0; i < 3; ++i) { + shift_t.get(i) *= + (turn_off + .5 * tanh(attenuation_parameter * + (get(distance_left_t) - attenuation_radius))) * + (turn_off + .5 * tanh(attenuation_parameter * + (get(distance_right_t) - attenuation_radius))); + } + + // Boosted shifts in the inner zone + const auto spacetime_metric_left = get_t_boosted_spacetime_metric_left(t); + const auto conformal_metric_left = gr::spatial_metric(spacetime_metric_left); + const auto inv_conformal_metric_left = + determinant_and_inverse(conformal_metric_left).second; + const auto shift_left = + gr::shift(spacetime_metric_left, inv_conformal_metric_left); + const auto lapse_left = gr::lapse(shift_left, spacetime_metric_left); + const auto spacetime_metric_right = get_t_boosted_spacetime_metric_right(t); + const auto conformal_metric_right = + gr::spatial_metric(spacetime_metric_right); + const auto inv_conformal_metric_right = + determinant_and_inverse(conformal_metric_right).second; + const auto shift_right = + gr::shift(spacetime_metric_right, inv_conformal_metric_right); + const auto lapse_right = gr::lapse(shift_right, spacetime_metric_right); + for (size_t i = 0; i < 3; ++i) { + shift_t.get(i) += + (1. - + (turn_off + .5 * tanh(attenuation_parameter * + (get(distance_left_t) - attenuation_radius))) * + (turn_off + + .5 * tanh(attenuation_parameter * + (get(distance_right_t) - attenuation_radius)))) * + (shift_left.get(i) + shift_right.get(i)); + } + return shift_t; +} + +template +Scalar +BinaryWithGravitationalWavesVariables::get_t_distance_left( + const DataType time) const { + // later it will depend on t by interpolation of past history + tnsr::I v = x; + const std::array position_left = {xcoords[0], offset[0], + offset[1]}; + for (size_t i = 0; i < time.size(); ++i) { + for (size_t j = 0; j < 3; ++j) { + v.get(j)[i] = x.get(j)[i] - gsl::at(position_left, j); + } + } + return magnitude(v); +} + +template +Scalar +BinaryWithGravitationalWavesVariables::get_t_distance_right( + const DataType time) const { + // later it will depend on t by interpolation of past history + tnsr::I v = x; + const std::array position_right = {xcoords[1], offset[0], + offset[1]}; + for (size_t i = 0; i < time.size(); ++i) { + for (size_t j = 0; j < 3; ++j) { + v.get(j)[i] = x.get(j)[i] - gsl::at(position_right, j); + } + } + return magnitude(v); +} + +template +Scalar +BinaryWithGravitationalWavesVariables::get_t_separation( + const DataType time) const { + // later it will depend on t by interpolation of past history + tnsr::I v = x; + const std::array separation = {xcoords[1] - xcoords[0], 0., 0.}; + for (size_t i = 0; i < time.size(); ++i) { + for (size_t j = 0; j < 3; ++j) { + v.get(j)[i] = gsl::at(separation, j); + } + } + return magnitude(v); +} + +template +tnsr::I +BinaryWithGravitationalWavesVariables::get_t_momentum_left( + const DataType time) const { + // later it will depend on t by interpolation of past history + tnsr::I v = x; + for (size_t i = 0; i < time.size(); ++i) { + for (size_t j = 0; j < 3; ++j) { + v.get(j)[i] = gsl::at(momentum_left, j); + } + } + return v; +} + +template +tnsr::I +BinaryWithGravitationalWavesVariables::get_t_momentum_right( + const DataType time) const { + // later it will depend on t by interpolation of past history + tnsr::I v = x; + for (size_t i = 0; i < time.size(); ++i) { + for (size_t j = 0; j < 3; ++j) { + v.get(j)[i] = gsl::at(momentum_right, j); + } + } + return v; +} + +template +tnsr::I +BinaryWithGravitationalWavesVariables::get_t_normal_left( + const DataType time) const { + tnsr::I v = x; + const Scalar distance_left = get_t_distance_left(time); + const std::array position_left = {xcoords[0], offset[0], + offset[1]}; + for (size_t i = 0; i < time.size(); ++i) { + for (size_t j = 0; j < 3; ++j) { + v.get(j)[i] = + (x.get(j)[i] - gsl::at(position_left, j)) / get(distance_left)[i]; + } + } + return v; +} + +template +tnsr::I +BinaryWithGravitationalWavesVariables::get_t_normal_right( + const DataType time) const { + tnsr::I v = x; + const Scalar distance_right = get_t_distance_right(time); + const std::array position_right = {xcoords[1], offset[0], + offset[1]}; + for (size_t i = 0; i < time.size(); ++i) { + for (size_t j = 0; j < 3; ++j) { + v.get(j)[i] = + (x.get(j)[i] - gsl::at(position_right, j)) / get(distance_right)[i]; + } + } + + return v; +} + +template +tnsr::I +BinaryWithGravitationalWavesVariables::get_t_normal_lr( + const DataType time) const { + tnsr::I v = x; + std::array separation = {xcoords[1] - xcoords[0], 0., 0.}; + Scalar separation_norm = get_t_separation(time); + for (size_t i = 0; i < time.size(); ++i) { + for (size_t j = 0; j < 3; ++j) { + v.get(j)[i] = gsl::at(separation, j) / get(separation_norm)[i]; + } + } + return v; +} + +template +tnsr::aa BinaryWithGravitationalWavesVariables< + DataType>::get_t_boosted_spacetime_metric_left(DataType t) const { + const auto& lapse_times_conformal_factor_minus_one = + get>(boost_vars[0]); + const auto& conformal_factor_minus_one = + get>(boost_vars[0]); + const auto& shift = + get>(boost_vars[0]); + + Scalar conformal_factor{t.size()}; + get(conformal_factor) = 1. + get(conformal_factor_minus_one); + Scalar lapse{t.size()}; + get(lapse) = (1. + get(lapse_times_conformal_factor_minus_one)) / + get(conformal_factor); + + tnsr::ii spatial_metric{t.size()}; + std::fill(spatial_metric.begin(), spatial_metric.end(), 0.); + for (size_t i = 0; i < 3; ++i) { + spatial_metric.get(i, i) = square(square(get(conformal_factor))); + } + + tnsr::aa spacetime_metric{t.size()}; + gr::spacetime_metric(make_not_null(&spacetime_metric), lapse, shift, + spatial_metric); + + const tnsr::I momentum_Datavector = get_t_momentum_left(t); + const std::array boost_velocity = { + -get<0>(momentum_Datavector)[0] / mass_left, + -get<1>(momentum_Datavector)[0] / mass_left, + -get<2>(momentum_Datavector)[0] / mass_left}; + const tnsr::Ab lorentz_boost_matrix_double = + sr::lorentz_boost_matrix(boost_velocity); + tnsr::Ab lorentz_boost_matrix{t.size()}; + for (size_t i = 0; i < 4; ++i) { + for (size_t j = 0; j < 4; ++j) { + lorentz_boost_matrix.get(i, j) = lorentz_boost_matrix_double.get(i, j); + } + } + + tnsr::aa spacetime_metric_boosted{t.size()}; + for (size_t i = 0; i < 4; ++i) { + for (size_t j = 0; j <= i; ++j) { + spacetime_metric_boosted.get(i, j) = 0.; + for (size_t k = 0; k < 4; ++k) { + for (size_t l = 0; l < 4; ++l) { + spacetime_metric_boosted.get(i, j) += lorentz_boost_matrix.get(k, i) * + lorentz_boost_matrix.get(l, j) * + spacetime_metric.get(k, l); + } + } + } + } + + return spacetime_metric_boosted; +} + +template +tnsr::aa BinaryWithGravitationalWavesVariables< + DataType>::get_t_boosted_spacetime_metric_right(DataType t) const { + const auto& lapse_times_conformal_factor_minus_one = + get>(boost_vars[1]); + const auto& conformal_factor_minus_one = + get>(boost_vars[1]); + const auto& shift = + get>(boost_vars[1]); + + Scalar conformal_factor{t.size()}; + get(conformal_factor) = 1. + get(conformal_factor_minus_one); + Scalar lapse{t.size()}; + get(lapse) = (1. + get(lapse_times_conformal_factor_minus_one)) / + get(conformal_factor); + + tnsr::ii spatial_metric{t.size()}; + std::fill(spatial_metric.begin(), spatial_metric.end(), 0.); + for (size_t i = 0; i < 3; ++i) { + spatial_metric.get(i, i) = square(square(get(conformal_factor))); + } + + tnsr::aa spacetime_metric{t.size()}; + gr::spacetime_metric(make_not_null(&spacetime_metric), lapse, shift, + spatial_metric); + + const tnsr::I momentum_Datavector = get_t_momentum_right(t); + const std::array boost_velocity = { + -get<0>(momentum_Datavector)[0] / mass_right, + -get<1>(momentum_Datavector)[0] / mass_right, + -get<2>(momentum_Datavector)[0] / mass_right}; + const tnsr::Ab lorentz_boost_matrix_double = + sr::lorentz_boost_matrix(boost_velocity); + tnsr::Ab lorentz_boost_matrix{t.size()}; + for (size_t i = 0; i < 4; ++i) { + for (size_t j = 0; j < 4; ++j) { + lorentz_boost_matrix.get(i, j) = lorentz_boost_matrix_double.get(i, j); + } + } + + tnsr::aa spacetime_metric_boosted{t.size()}; + for (size_t i = 0; i < 4; ++i) { + for (size_t j = 0; j <= i; ++j) { + spacetime_metric_boosted.get(i, j) = 0.; + for (size_t k = 0; k < 4; ++k) { + for (size_t l = 0; l < 4; ++l) { + spacetime_metric_boosted.get(i, j) += lorentz_boost_matrix.get(k, i) * + lorentz_boost_matrix.get(l, j) * + spacetime_metric.get(k, l); + } + } + } + } + + return spacetime_metric_boosted; +} + +template +tnsr::aa BinaryWithGravitationalWavesVariables< + DataType>::get_t_superposed_spacetime_metric(DataType t) const { + const auto spacetime_metric_left = get_t_boosted_spacetime_metric_left(t); + const auto spacetime_metric_right = get_t_boosted_spacetime_metric_right(t); + tnsr::aa superposed_spacetime_metric{t.size()}; + for (size_t i = 0; i < 4; ++i) { + for (size_t j = 0; j <= i; ++j) { + superposed_spacetime_metric.get(i, j) = + spacetime_metric_left.get(i, j) + spacetime_metric_right.get(i, j); + } + superposed_spacetime_metric.get(i, i) -= 1.; + } + + superposed_spacetime_metric.get(0, 0) += 2.; + + return superposed_spacetime_metric; +} + +template class BinaryWithGravitationalWavesVariables; + +} // namespace Xcts::AnalyticData::detail + +template class Xcts::AnalyticData::CommonVariables< + DataVector, typename Xcts::AnalyticData::detail:: + BinaryWithGravitationalWavesVariables::Cache>; diff --git a/src/PointwiseFunctions/AnalyticData/Xcts/BinaryWithGravitationalWaves.hpp b/src/PointwiseFunctions/AnalyticData/Xcts/BinaryWithGravitationalWaves.hpp new file mode 100644 index 000000000000..960bf2e9f1e5 --- /dev/null +++ b/src/PointwiseFunctions/AnalyticData/Xcts/BinaryWithGravitationalWaves.hpp @@ -0,0 +1,692 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include + +#include +#include +#include +#include + +#include "DataStructures/CachedTempBuffer.hpp" +#include "DataStructures/DataBox/Prefixes.hpp" +#include "DataStructures/Tensor/IndexType.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "DataStructures/Tensor/TypeAliases.hpp" +#include "Elliptic/Systems/Xcts/Tags.hpp" +#include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp" +#include "Options/Context.hpp" +#include "Options/ParseError.hpp" +#include "Options/String.hpp" +#include "PointwiseFunctions/AnalyticData/Xcts/CommonVariables.hpp" +#include "PointwiseFunctions/AnalyticSolutions/Xcts/Flatness.hpp" +#include "PointwiseFunctions/AnalyticSolutions/Xcts/Schwarzschild.hpp" +#include "PointwiseFunctions/GeneralRelativity/Tags.hpp" +#include "PointwiseFunctions/GeneralRelativity/Tags/Conformal.hpp" +#include "PointwiseFunctions/InitialDataUtilities/Background.hpp" +#include "PointwiseFunctions/InitialDataUtilities/InitialGuess.hpp" +#include "PointwiseFunctions/SpecialRelativity/LorentzBoostMatrix.hpp" +#include "Utilities/CallWithDynamicType.hpp" +#include "Utilities/Requires.hpp" +#include "Utilities/Serialization/CharmPupable.hpp" +#include "Utilities/TMPL.hpp" +#include "Utilities/TaggedTuple.hpp" + +/// \cond +namespace PUP { +class er; +} // namespace PUP +/// \endcond + +namespace Xcts::AnalyticData { + +namespace detail { + +template +using BinaryWithGravitationalWavesVariablesCache = + cached_temp_buffer_from_typelist, + tmpl::list< + ::Tags::deriv< + Xcts::Tags::ShiftBackground, + tmpl::size_t<3>, Frame::Inertial>, + gr::Tags::Conformal, 0>, + gr::Tags::Conformal, 0>, + gr::Tags::Conformal, 0>, + // For initial guesses + Xcts::Tags::ConformalFactorMinusOne, + Xcts::Tags::LapseTimesConformalFactorMinusOne, + Xcts::Tags::ShiftExcess>, + hydro_tags>>; + +template +struct BinaryWithGravitationalWavesVariables + : CommonVariables> { + static constexpr size_t Dim = 3; + using Cache = BinaryWithGravitationalWavesVariablesCache; + using Base = + CommonVariables>; + using Base::operator(); + + using superposed_tags = tmpl::append< + tmpl::list< + Xcts::Tags::ConformalMetric, + gr::Tags::Conformal, 0>, + gr::Tags::Conformal, 0>, + gr::Tags::Conformal, 0>>, + hydro_tags>; + + using boost_tags = tmpl::append< + tmpl::list, + Xcts::Tags::ConformalFactorMinusOne, + Xcts::Tags::LapseTimesConformalFactorMinusOne, + Xcts::Tags::ShiftExcess>>; + + BinaryWithGravitationalWavesVariables( + std::optional>> local_mesh, + std::optional>> + local_inv_jacobian, + const tnsr::I& local_x, + const std::array& local_masses, + const std::array& local_xcoords, + const std::array& local_momentum_left, + const std::array& local_momentum_right, + const double local_y_offset, const double local_z_offset, + const double local_attenuation_parameter, + const double local_attenuation_radius, + std::array, 2> local_x_isolated, + tuples::tagged_tuple_from_typelist local_flat_vars, + std::array, 2> + local_isolated_vars, + std::array, 2> + local_boost_vars) + : Base(std::move(local_mesh), std::move(local_inv_jacobian)), + mesh(std::move(local_mesh)), + inv_jacobian(std::move(local_inv_jacobian)), + x(local_x), + mass_left(local_masses[0]), + mass_right(local_masses[1]), + xcoords(local_xcoords), + momentum_left(local_momentum_left), + momentum_right(local_momentum_right), + offset({local_y_offset, local_z_offset}), + attenuation_parameter(local_attenuation_parameter), + attenuation_radius(local_attenuation_radius), + x_isolated(std::move(local_x_isolated)), + flat_vars(std::move(local_flat_vars)), + isolated_vars(std::move(local_isolated_vars)), + boost_vars(std::move(local_boost_vars)) {} + + std::optional>> mesh; + std::optional>> + inv_jacobian; + + tnsr::I x; + double mass_left; + double mass_right; + std::array xcoords; + std::array momentum_left; + std::array momentum_right; + std::array offset; + double attenuation_parameter; + double attenuation_radius; + std::array, 2> x_isolated; + tuples::tagged_tuple_from_typelist flat_vars; + std::array, 2> + isolated_vars; + std::array, 2> boost_vars; + + double time_displacement = 0.01; + DataType present_time = make_with_value(x, 1.); + + template > = nullptr> + void superposition(gsl::not_null superposed_var, + gsl::not_null /*cache*/, Tag /*meta*/) const { + for (size_t i = 0; i < superposed_var->size(); ++i) { + (*superposed_var)[i] = get(isolated_vars[0])[i] + + get(isolated_vars[1])[i] - + get(flat_vars)[i]; + } + } + + void operator()(gsl::not_null*> conformal_metric, + gsl::not_null cache, + Xcts::Tags::ConformalMetric + meta) const override; + void operator()( + gsl::not_null*> deriv_conformal_metric, + gsl::not_null cache, + ::Tags::deriv, + tmpl::size_t, Frame::Inertial> + meta) const override; + void operator()( + gsl::not_null*> trace_extrinsic_curvature, + gsl::not_null cache, + gr::Tags::TraceExtrinsicCurvature meta) const override; + void operator()(gsl::not_null*> dt_trace_extrinsic_curvature, + gsl::not_null cache, + ::Tags::dt> meta) + const override; + void operator()( + gsl::not_null*> shift_background, + gsl::not_null /*cache*/, + Xcts::Tags::ShiftBackground /*meta*/) + const override { + std::fill(shift_background->begin(), shift_background->end(), 0.); + } + void operator()( + gsl::not_null*> deriv_shift_background, + gsl::not_null /*cache*/, + ::Tags::deriv, + tmpl::size_t, Frame::Inertial> /*meta*/) const { + std::fill(deriv_shift_background->begin(), deriv_shift_background->end(), + 0.); + } + void operator()(gsl::not_null*> + longitudinal_shift_background_minus_dt_conformal_metric, + gsl::not_null cache, + Xcts::Tags::LongitudinalShiftBackgroundMinusDtConformalMetric< + DataType, Dim, Frame::Inertial> /*meta*/) const override; + void operator()( + const gsl::not_null*> conformal_energy_density, + const gsl::not_null cache, + gr::Tags::Conformal, 0> meta) const { + superposition(conformal_energy_density, cache, meta); + } + void operator()( + const gsl::not_null*> conformal_stress_trace, + const gsl::not_null cache, + gr::Tags::Conformal, 0> meta) const { + superposition(conformal_stress_trace, cache, meta); + } + void operator()( + const gsl::not_null*> conformal_momentum_density, + const gsl::not_null cache, + gr::Tags::Conformal, 0> meta) + const { + superposition(conformal_momentum_density, cache, meta); + } + void operator()(gsl::not_null*> conformal_factor_minus_one, + gsl::not_null cache, + Xcts::Tags::ConformalFactorMinusOne meta) const; + void operator()( + gsl::not_null*> lapse_times_conformal_factor_minus_one, + gsl::not_null cache, + Xcts::Tags::LapseTimesConformalFactorMinusOne meta) const; + void operator()( + gsl::not_null*> shift_excess, + gsl::not_null cache, + Xcts::Tags::ShiftExcess meta) const; + void operator()(const gsl::not_null*> rest_mass_density, + const gsl::not_null cache, + hydro::Tags::RestMassDensity meta) const { + superposition(rest_mass_density, cache, meta); + } + void operator()(const gsl::not_null*> specific_enthalpy, + const gsl::not_null cache, + hydro::Tags::SpecificEnthalpy meta) const { + superposition(specific_enthalpy, cache, meta); + } + void operator()(const gsl::not_null*> pressure, + const gsl::not_null cache, + hydro::Tags::Pressure meta) const { + superposition(pressure, cache, meta); + } + void operator()(const gsl::not_null*> spatial_velocity, + const gsl::not_null cache, + hydro::Tags::SpatialVelocity meta) const { + superposition(spatial_velocity, cache, meta); + } + void operator()(const gsl::not_null*> lorentz_factor, + const gsl::not_null cache, + hydro::Tags::LorentzFactor meta) const { + superposition(lorentz_factor, cache, meta); + } + void operator()(const gsl::not_null*> magnetic_field, + const gsl::not_null cache, + hydro::Tags::MagneticField meta) const { + superposition(magnetic_field, cache, meta); + } + + private: + // void interpolate_past_history(); + // DataType find_retarded_time_left(DataType t0) const; + // DataType find_retarded_time_right(DataType t0) const; + Scalar get_t_distance_left(DataType t) const; + Scalar get_t_distance_right(DataType t) const; + Scalar get_t_separation(DataType t) const; + tnsr::I get_t_momentum_left(DataType t) const; + tnsr::I get_t_momentum_right(DataType t) const; + tnsr::I get_t_normal_left(DataType t) const; + tnsr::I get_t_normal_right(DataType t) const; + tnsr::I get_t_normal_lr(DataType t) const; + // DataType integrate_term(DataType t, size_t i, size_t j, int left_right, + // double t0) const; + Scalar get_t_trace_extrinsic_curvature( + DataType t, Mesh<3> local_mesh, + InverseJacobian + local_inv_jacobian) const; + tnsr::ii get_t_conformal_metric(DataType t) const; + tnsr::ii get_t_radiative_term(DataType t) const; + tnsr::ii get_t_near_zone_term(DataType t) const; + tnsr::ii get_t_present_term(DataType t) const; + tnsr::ii get_t_past_term(DataType t) const; + tnsr::ii get_t_integral_term(DataType t) const; + Scalar get_t_lapse(DataType t) const; + tnsr::I get_t_shift(DataType t) const; + tnsr::aa get_t_boosted_spacetime_metric_left(DataType t) const; + tnsr::aa get_t_boosted_spacetime_metric_right(DataType t) const; + tnsr::aa get_t_superposed_spacetime_metric(DataType t) const; +}; +} // namespace detail + +/*! + * \brief Binary initial data with realistic wave background, + * constructed in Post-Newtonian approximations. + * + * The main goal of this implementation is to improve the extracted + * wave forms, for example, by minimizing junk radiation. + * The data is only valid for black holes without spin. Even so, there is some + * work done to describe such systems that could later be implemented. + * The objects are constructed from a superposition of two isolated objects that + * are boosted with respect to each other. The radiative data is constructed + * from Post-Newtonian expansions for the inspiral phase, in orders of + * \f$\epsilon = 1/c\f$, in \cite Jaranowski1997ky. In ADMTT gauge it is + * possible to get the 3-metric as \f$\gamma^{PN}_{ij} = \psi^{4}_{PN} + * \delta_{ij} + h^{TT}_{ij}\f$ where \f$h^{TT}_{ij}\f$ is the radiative part + * and the non-radiative Post-Newtonian conformal factor is given by + * + * \f{equation}{ + * \psi_{PN} = 1 + \sum_{a=1}^{2} \frac{E_a}{2 r_a} + O(\epsilon^6) + * \f} + * + * and + * + * \f{equation}{E_a = (\epsilon^2) m_a + (\epsilon^4) \Bigr(\frac{p_a^2}{2 m_a} + * - \frac{m_1 m_2}{2 r_{12}}\Bigr) \f} + * + * with \f$\vec{p}_a\f$ the linear momentum, \f$r_a\f$ the distance to each + * black hole center of mass from the point of calculation and \f$r_{12}\f$ + * separation between the two black holes and \f$m_a\f$ is the mass of each + * black hole. Near each black hole, the 3-metric can be approximated by the + * Schwarzschild 3-metric in isotropic coordinates. + * + * In \cite Mundim2010hu, the radiative term \f$h^{TT}_{ij}\f$ is decomposed + * into two parts, a near-zone that is only valid close to the black holes and a + * remainder that makes corrections far from the black holes, \f$h^{TT}_{ij} = + * h^{TT\ (NZ)}_{ij} + h^{TT\ (remainder)}_{ij} + O(\epsilon^5)\f$. The + * near-zone term is given by $h^{TT\ (NZ)}_{ij} = (\epsilon^4) h^{TT}_{(4)ij} + + * (\epsilon^5) h^{TT}_{(5)ij}$, with + * + * \f{align}{ + * h^{TT\ i j}_{(4)} &= \frac{1}{4} \sum_a \frac{1}{m_a r_a} \Bigr\{ + * [p_a^2-5(\hat{n}_a \cdot \vec{p}_a)^2] \delta^{i j}+2 p_a^i p_a^j + * +[3(\hat{n}_a \cdot \vec{p}_a)^2-5p_a^2] n_a^i n_a^j +12(\hat{n}_a \cdot + * \vec{p}_a) n_a^{(i} p_a^{j)} \Bigr\} \nonumber \\ + * &+\frac{1}{8} \sum_a \sum_{b \neq a} m_a m_b \Bigr\{-\frac{32}{s_{a + * b}}(\frac{1}{r_{a b}}+\frac{1}{s_{a b}}) n_{a b}^i n_{a + * b}^j+2(\frac{r_a+r_b}{r_{a b}^3}+\frac{12}{s_{a b}^2}) n_a^i + * n_b^j+32(\frac{2}{s_{a b}^2}-\frac{1}{r_{a b}^2}) n_a^{(i} n_{a b}^{j)} + * \label{eq:near_zone_term} \\ + * &+[\frac{5}{r_{a b} r_a}-\frac{1}{r_{a b}^3}(\frac{r_b^2}{r_a}+3 + * r_a)-\frac{8}{s_{a b}}(\frac{1}{r_a}+\frac{1}{s_{a b}})] n_a^i n_a^j+[5 + * \frac{r_a}{r_{a b}^3}(\frac{r_a}{r_b}-1)-\frac{17}{r_{a b} r_a}+\frac{4}{r_a + * r_b}+\frac{8}{s_{a b}}(\frac{1}{r_a}+\frac{4}{r_{a b}})] \delta^{i j}\Bigr\}, + * \nonumber \f} + * + * where \f$\hat{n}_a\f$ is the unit normal vector pointing to the black hole + * center of mass, \f$\hat{n}_{ab}\f$ is the unit normal vector pointing from + * black hole \f$a\f$ to black hole \f$b\f$ and \f$s_{ab} = r_a + r_b + + * r_{ab}\f$. The term \f$h^{TT}_{(5)ij}\f$ is a spatially constant field that + * just varies in time, for initial data we can choose an initial time such that + * \f$h^{TT}_{(5)ij} = 0\f$. + * + * Looking at \cite Kelly2007uc, the remainder term in itself is decomposed in + * general computations for specific vectors as + * + * \f{equation}{ + * h^{TT\ (remainder)}_{ij} = H^{TT\ 1}_{ij} \Bigr[ + * \frac{\vec{p_1}}{\sqrt{m_1}}\Bigr] + H^{TT\ 2}_{ij} \Bigr[ + * \frac{\vec{p_2}}{\sqrt{m_2}}\Bigr] - H^{TT\ 1}_{ij} \Bigr[ \sqrt{\frac{m_1 + * m_2}{2 r_{12}}} \hat{n_{12}} \Bigr] - H^{TT\ 2}_{ij} \Bigr[ \sqrt{\frac{m_1 + * m_2}{2 r_{12}}} \hat{n_{12}} \Bigr], \f} + * + * each of this is composed of three different computations: one computed at + * present time \f$t\f$, other at retarded time \f$t_{a}^{r}\f$ defined by \f$t + * - t_{a}^{r} - r_a(t_{a}^{r}) = 0\f$ and the last is an integral between the + * two times: + * + * \f{equation}{ + * H^{TT\ a}_{ij} [ \vec{u} ] = H^{TT\ a}_{ij} [ \vec{u} ; t] + H^{TT\ a}_{ij} [ + * \vec{u} ; t^{r}_a] + H^{TT\ a}_{ij} [ \vec{u} ; t_{a}^{r} \to t]. \f} + * + * Explicitly they are + * + * \f{equation}{ + * H^{TT\ a}_{ij} [ \vec{u} ; t] = -\frac{1}{4 r_a(t)} \Bigr\{ [u^2 - 5(\vec{u} + * \cdot \hat{n}_a)^2] \delta_{ij} + 2 u^iu^j + 3(\vec{u}\cdot\hat{n}_a)^2 - 5 + * u^2] n_a^i n_a^j + 12 (\vec{u} \cdot \hat{n}_a) u^{(i}n_a^{j)}\Bigr\}_t, + * \label{eq:present_term} \f} + * + * \f{equation}{ + * H^{TT\ a}_{ij} [ \vec{u} ; t^{r}_a] = \frac{1}{r_a(t^{r}_a)} \Bigr\{ [-2u^2 + * + 2(\vec{u} \cdot \hat{n}_a)^2] \delta^{ij} + 4u^iu^j + [2 u^2 + 2 (\vec{u} + * \cdot \hat{n}_a)^2 ] n_a^i n_a^j - 8(\vec{u}\cdot\hat{n}_a) u^{(i}n_a^{j)} + * \Bigr\}_{t^{r}_a}, + * \label{eq:retarded_term} \f} + * + * and + * + * \f{align}{ + * H^{TT\ a}_{ij} [ \vec{u} ; t^{r}_a \to t] &= \nonumber \\ + * &- \int^t_{t^{r}_a} d\tau \frac{(t-\tau)}{r_a(\tau)^3} \Bigr\{ [-5u^2 + + * 9(\vec{u} \cdot \hat{n}_a)^2] \delta^{ij} + 6u^iu^j - 12 + * (\vec{u}\cdot\hat{n}_a) u^{(i}n_a^{j)} + [9 u^2 - 15(\vec{u} \cdot + * \hat{n}_a)^2 ] n_a^i n_a^j\Bigr\} \label{eq:integral_term} \\ + * &- \int^t_{t^{r}_a} d\tau \frac{(t-\tau)^3}{r_a(\tau)^5} \Bigr\{ [u^2 - + * 5(\vec{u} \cdot \hat{n}_a)^2] \delta^{ij} + 2 u^iu^j - 20 + * (\vec{u}\cdot\hat{n}_a) u^{(i}n_a^{j)} + [-5 u^2 + 35(\vec{u} \cdot + * \hat{n}_a)^2 ] n_a^i n_a^j\Bigr\}. \nonumber \f} + * + * \warning The radiative terms, equations \f$\eqref{eq:near_zone_term}\f$, + * \f$\eqref{eq:present_term}\f$, \f$\eqref{eq:retarded_term}\f$ and + * \f$\eqref{eq:integral_term}\f$, are not implemented yet. Instead these terms + * are set to zero. + * + * With this the whole spatial metric is computed up to \f$2PN\f$ order and the + * radiative term agrees well with quadrupole predictions. In \cite Tichy2002ec, + * the extrinsic curvature is given up to \f$2.5PN\f$ order by + * + * \f{equation}{ + * K^{ij}_{PN} = - \psi^{-10}_{PN} \Bigr[ (\epsilon^3) \tilde{\pi}_{(3)}^{ij} + + * (\epsilon^5) \frac{1}{2} \dot{h}^{TT}_{(4)ij} + (\epsilon^5) (\phi_{(2)} + * \tilde{\pi}_{(3)}^{ij})^{TT} \Bigr] + O(\epsilon^6). \f} + * + * where + * + * \f{equation}{ + * \tilde{\pi}_{(3)}^{i j}=\frac{1}{16 \pi} \sum_a p_{a}^k\{-\delta_{i + * j}(\frac{1}{r_a})_{, k}+2[\delta_{i k}(\frac{1}{r_a})_{, j}+\delta_{j + * k}(\frac{1}{r_a})_{, i}]-\frac{1}{2} r_{a, i j k}\}. \f} + * + * + * To be able to calculate equations \f$\eqref{eq:retarded_term}\f$ and + * \f$\eqref{eq:integral_term}\f$ we need to look into the past history + * of the binary at least up to the time were the generated wave can reach the + * furthest point on the grid. To do so we must evolve the binary backward in + * time. Because we are only looking into the inspiral phase we can follow a + * simple Hamiltonian evolution computed in Post-Newtonian orders. The + * equations to be solved are + * + * \f{equation}{ + * \frac{d X^i}{d t}=\frac{\partial H}{\partial P_i} + * \f} + * + * and + * + * \f{equation}{ + * \frac{d P_i}{d t}=-\frac{\partial H}{\partial X^i}+F_i, + * \f} + * + * where $H$ is the Post-Newtonian Hamiltonian, $X^i$ is the separation + * vector between the two particles, $P_i$ is the momentum of one particle + * in the center of mass frame and $F_i$ is the radiation-reaction flux term. + * The Post-Newtonian Hamiltonian is given in \cite Buonanno2005xu. + */ +template +class BinaryWithGravitationalWaves + : public elliptic::analytic_data::Background, + public elliptic::analytic_data::InitialGuess { + public: + struct XCoords { + static constexpr Options::String help = + "The coordinates on the x-axis where the two objects are placed."; + using type = std::array; + }; + struct Masses { + static constexpr Options::String help = + "The mass of each object, first left and second right."; + using type = std::array; + }; + struct MomentumLeft { + static constexpr Options::String help = + "The momentum assigned to the left object."; + using type = std::array; + }; + struct MomentumRight { + static constexpr Options::String help = + "The momentum assigned to the right object."; + using type = std::array; + }; + struct CenterOfMassOffset { + static constexpr Options::String help = { + "Offset in the y and z axes applied to both objects in order to " + "control the center of mass."}; + using type = std::array; + }; + struct ObjectLeft { + static constexpr Options::String help = + "The object placed on the negative x-axis."; + using type = std::unique_ptr; + }; + struct ObjectRight { + static constexpr Options::String help = + "The object placed on the positive x-axis."; + using type = std::unique_ptr; + }; + struct AttenuationParameter { + static constexpr Options::String help = + "The parameter controlling the transition width of the attenuation " + "function."; + using type = double; + }; + struct AttenuationRadius { + static constexpr Options::String help = + "The parameter controlling the transition center of the attenuation " + "function."; + using type = double; + }; + + using options = tmpl::list; + static constexpr Options::String help = + "Binary black hole initial data with realistic wave background, " + "constructed in Post-Newtonian approximations. "; + + BinaryWithGravitationalWaves() = default; + BinaryWithGravitationalWaves(const BinaryWithGravitationalWaves&) = delete; + BinaryWithGravitationalWaves& operator=(const BinaryWithGravitationalWaves&) = + delete; + BinaryWithGravitationalWaves(BinaryWithGravitationalWaves&&) = default; + BinaryWithGravitationalWaves& operator=(BinaryWithGravitationalWaves&&) = + default; + ~BinaryWithGravitationalWaves() override = default; + + BinaryWithGravitationalWaves( + const std::array xcoords, const std::array masses, + const std::array momentum_left, + const std::array momentum_right, + const std::array center_of_mass_offset, + std::unique_ptr object_left, + std::unique_ptr object_right, + const double attenuation_parameter, const double attenuation_radius, + const Options::Context& context = {}) + : xcoords_(xcoords), + masses_(masses), + momentum_left_(momentum_left), + momentum_right_(momentum_right), + y_offset_(center_of_mass_offset[0]), + z_offset_(center_of_mass_offset[1]), + superposed_objects_({std::move(object_left), std::move(object_right)}), + attenuation_parameter_(attenuation_parameter), + attenuation_radius_(attenuation_radius) { + if (masses_[0] <= 0. || masses_[1] <= 0.) { + PARSE_ERROR(context, "The masses must be positive."); + } + if (xcoords_[0] >= xcoords_[1]) { + PARSE_ERROR(context, "Specify 'XCoords' ascending from left to right."); + } + if (attenuation_parameter_ < 0.) { + PARSE_ERROR(context, "'AttenuationParameter' must be positive."); + } + if (attenuation_radius_ <= 0.) { + PARSE_ERROR(context, "'AttenuationRadius' must be positive."); + } + } + + explicit BinaryWithGravitationalWaves(CkMigrateMessage* m) + : elliptic::analytic_data::Background(m), + elliptic::analytic_data::InitialGuess(m) {} + using PUP::able::register_constructor; + WRAPPED_PUPable_decl_template(BinaryWithGravitationalWaves); + + template + tuples::TaggedTuple variables( + const tnsr::I& x, + tmpl::list /*meta*/) const { + return variables_impl(x, std::nullopt, std::nullopt, + tmpl::list{}); + } + template + tuples::TaggedTuple variables( + const tnsr::I& x, const Mesh<3>& mesh, + const InverseJacobian& inv_jacobian, + tmpl::list /*meta*/) const { + return variables_impl(x, mesh, inv_jacobian, + tmpl::list{}); + } + + // NOLINTNEXTLINE + void pup(PUP::er& p) override { + elliptic::analytic_data::Background::pup(p); + elliptic::analytic_data::InitialGuess::pup(p); + p | xcoords_; + p | masses_; + p | y_offset_; + p | z_offset_; + p | momentum_left_; + p | momentum_right_; + p | superposed_objects_; + p | attenuation_parameter_; + p | attenuation_radius_; + } + + /// Coordinates of the objects, ascending left to right + const std::array& x_coords() const { return xcoords_; } + /// The momentum of the left object. + const std::array& momentum_left() const { return momentum_left_; } + /// The momentum of the right object. + const std::array& momentum_right() const { + return momentum_right_; + } + /// Offset in y and z coordinates of the objects + double y_offset() const { return y_offset_; } + double z_offset() const { return z_offset_; } + /// The two objects. First entry is the left object, second entry is the right + /// object. + const std::array, 2>& superposed_objects() + const { + return superposed_objects_; + } + double attenuation_parameter() const { return attenuation_parameter_; } + double attenuation_radius() const { return attenuation_radius_; } + + private: + std::array xcoords_{}; + std::array masses_{}; + std::array momentum_left_{}; + std::array momentum_right_{}; + double y_offset_{}; + double z_offset_{}; + std::array, 2> superposed_objects_{}; + Xcts::Solutions::Flatness flatness_{}; + double attenuation_parameter_{}; + double attenuation_radius_{}; + + template + tuples::TaggedTuple variables_impl( + const tnsr::I& x, + std::optional>> mesh, + std::optional>> + inv_jacobian, + tmpl::list /*meta*/) const { + std::array, 2> x_isolated{{x, x}}; + const std::array, 2> coords_isolated{ + {{{xcoords_[0], y_offset_, z_offset_}}, + {{xcoords_[1], y_offset_, z_offset_}}}}; + // Possible optimization: Only retrieve those superposed tags from the + // isolated solutions that are actually needed. This needs some dependency + // logic, because some of the non-superposed tags depend on superposed tags. + using VarsComputer = + detail::BinaryWithGravitationalWavesVariables; + using requested_superposed_tags = typename VarsComputer::superposed_tags; + std::array, 2> + isolated_vars; + for (size_t i = 0; i < 2; ++i) { + for (size_t dim = 0; dim < 3; dim++) { + gsl::at(x_isolated, i).get(dim) -= gsl::at(coords_isolated, i)[dim]; + } + gsl::at(isolated_vars, i) = get_isolated_vars( + *gsl::at(superposed_objects_, i), gsl::at(x_isolated, i)); + } + auto flat_vars = flatness_.variables(x, requested_superposed_tags{}); + using requested_boost_tags = typename VarsComputer::boost_tags; + std::array, 2> + boost_vars; + std::array, 2> x_unboosted{{x, x}}; + sr::lorentz_boost(make_not_null(&(gsl::at(x_unboosted, 0))), + gsl::at(x_isolated, 0), 0., momentum_left_ / masses_[0]); + sr::lorentz_boost(make_not_null(&(gsl::at(x_unboosted, 1))), + gsl::at(x_isolated, 1), 0., momentum_right_ / masses_[1]); + for (size_t i = 0; i < 2; ++i) { + gsl::at(boost_vars, i) = get_isolated_vars( + *gsl::at(superposed_objects_, i), gsl::at(x_unboosted, i)); + } + + typename VarsComputer::Cache cache{get_size(*x.begin())}; + const VarsComputer computer{std::move(mesh), + std::move(inv_jacobian), + x, + masses_, + xcoords_, + momentum_left_, + momentum_right_, + y_offset_, + z_offset_, + attenuation_parameter_, + attenuation_radius_, + std::move(x_isolated), + std::move(flat_vars), + std::move(isolated_vars), + std::move(boost_vars)}; + return {cache.get_var(computer, RequestedTags{})...}; + } + + template + tuples::tagged_tuple_from_typelist get_isolated_vars( + const IsolatedObjectBase& isolated_object, const Args&... args) const { + return call_with_dynamic_type, + IsolatedObjectClasses>( + &isolated_object, [&args...](const auto* const derived) { + return derived->variables(args..., TagsList{}); + }); + } +}; + +/// \cond +template +PUP::able::PUP_ID BinaryWithGravitationalWaves< + IsolatedObjectBase, IsolatedObjectClasses>::my_PUP_ID = 0; // NOLINT +/// \endcond + +} // namespace Xcts::AnalyticData diff --git a/src/PointwiseFunctions/AnalyticData/Xcts/CMakeLists.txt b/src/PointwiseFunctions/AnalyticData/Xcts/CMakeLists.txt index 7085faf0a67a..43b8934af91d 100644 --- a/src/PointwiseFunctions/AnalyticData/Xcts/CMakeLists.txt +++ b/src/PointwiseFunctions/AnalyticData/Xcts/CMakeLists.txt @@ -9,6 +9,7 @@ spectre_target_sources( ${LIBRARY} PRIVATE Binary.cpp + BinaryWithGravitationalWaves.cpp ) spectre_target_headers( @@ -17,6 +18,7 @@ spectre_target_headers( HEADERS AnalyticData.hpp Binary.hpp + BinaryWithGravitationalWaves.hpp CommonVariables.hpp CommonVariables.tpp ) @@ -31,6 +33,7 @@ target_link_libraries( LinearOperators Options Parallel + SpecialRelativity Spectral Utilities PRIVATE diff --git a/tests/Unit/Elliptic/Systems/Xcts/BoundaryConditions/CMakeLists.txt b/tests/Unit/Elliptic/Systems/Xcts/BoundaryConditions/CMakeLists.txt index 3a60b3aff691..687d104b869a 100644 --- a/tests/Unit/Elliptic/Systems/Xcts/BoundaryConditions/CMakeLists.txt +++ b/tests/Unit/Elliptic/Systems/Xcts/BoundaryConditions/CMakeLists.txt @@ -7,6 +7,7 @@ set(LIBRARY_SOURCES Test_ApparentHorizon.cpp Test_Flatness.cpp Test_Robin.cpp + Test_SuperposedBoostedBinary.cpp ) add_test_library(${LIBRARY} "${LIBRARY_SOURCES}") diff --git a/tests/Unit/Elliptic/Systems/Xcts/BoundaryConditions/SuperposedBoostedBinary.py b/tests/Unit/Elliptic/Systems/Xcts/BoundaryConditions/SuperposedBoostedBinary.py new file mode 100644 index 000000000000..7ecb5ebb0e85 --- /dev/null +++ b/tests/Unit/Elliptic/Systems/Xcts/BoundaryConditions/SuperposedBoostedBinary.py @@ -0,0 +1,161 @@ +# Distributed under the MIT License. +# See LICENSE.txt for details. + +import numpy as np + +x_coords = [-5.0, 6.0] +momentum_left = [0.01, 0.01, 0.01] +momentum_right = [-0.01, -0.01, -0.01] +y_offset = 0.02 +z_offset = 0.01 +masses = [1.1, 0.43] + + +def x_left(x): + x_left = np.array(x) + x_left[0] -= x_coords[0] + x_left[1] -= y_offset + x_left[2] -= z_offset + + velocity = np.array(momentum_left) / masses[0] + beta2 = np.dot(velocity, velocity) + gamma = 1.0 / np.sqrt(1.0 - beta2) + + boost_matrix = np.zeros((4, 4)) + boost_matrix[0, 0] = gamma + boost_matrix[0, 1:] = -gamma * velocity + boost_matrix[1:, 0] = -gamma * velocity + boost_matrix[1:, 1:] = ( + np.identity(3) + (gamma - 1.0) * np.outer(velocity, velocity) / beta2 + ) + + boosted_vector = np.dot(boost_matrix, np.append(0.0, x_left)) + + return boosted_vector[1:] + + +def x_right(x): + x_right = np.array(x) + x_right[0] -= x_coords[1] + x_right[1] -= y_offset + x_right[2] -= z_offset + + velocity = np.array(momentum_right) / masses[1] + beta2 = np.dot(velocity, velocity) + gamma = 1.0 / np.sqrt(1.0 - beta2) + + boost_matrix = np.zeros((4, 4)) + boost_matrix[0, 0] = gamma + boost_matrix[0, 1:] = -gamma * velocity + boost_matrix[1:, 0] = -gamma * velocity + boost_matrix[1:, 1:] = ( + np.identity(3) + (gamma - 1.0) * np.outer(velocity, velocity) / beta2 + ) + + boosted_vector = np.dot(boost_matrix, np.append(0.0, x_right)) + + return boosted_vector[1:] + + +def spacetime_left(x): + r = np.linalg.norm(x_left(x)) + conformalfactor = 1.0 + 0.5 * masses[0] / r + lapse = (1.0 - 0.5 * masses[0] / r) / conformalfactor + shift = np.zeros(3) + spatial_metric = np.identity(3) * pow(conformalfactor, 4) + + spacetime_metric = np.zeros((4, 4)) + spacetime_metric[0, 0] = -(lapse**2) + np.dot( + shift, np.dot(spatial_metric, shift) + ) + spacetime_metric[0, 1:] = np.dot(spatial_metric, shift) + spacetime_metric[1:, 0] = spacetime_metric[0, 1:] + spacetime_metric[1:, 1:] = spatial_metric + + return spacetime_metric + + +def spacetime_right(x): + r = np.linalg.norm(x_right(x)) + conformalfactor = 1.0 + 0.5 * masses[1] / r + lapse = (1.0 - 0.5 * masses[1] / r) / conformalfactor + shift = np.zeros(3) + spatial_metric = np.identity(3) * pow(conformalfactor, 4) + + spacetime_metric = np.zeros((4, 4)) + spacetime_metric[0, 0] = -(lapse**2) + np.dot( + shift, np.dot(spatial_metric, shift) + ) + spacetime_metric[0, 1:] = np.dot(spatial_metric, shift) + spacetime_metric[1:, 0] = spacetime_metric[0, 1:] + spacetime_metric[1:, 1:] = spatial_metric + + return spacetime_metric + + +def boost_spacetime_metric(spacetime_metric, velocity): + beta2 = np.dot(velocity, velocity) + gamma = 1.0 / np.sqrt(1.0 - beta2) + + boost_matrix = np.zeros((4, 4)) + boost_matrix[0, 0] = gamma + boost_matrix[0, 1:] = -gamma * velocity + boost_matrix[1:, 0] = -gamma * velocity + boost_matrix[1:, 1:] = ( + np.identity(3) + (gamma - 1.0) * np.outer(velocity, velocity) / beta2 + ) + + return np.dot(boost_matrix, np.dot(spacetime_metric, boost_matrix.T)) + + +def shift(spacetime_metric): + spatial_metric = spacetime_metric[1:, 1:] + inverse_spatial_metric = np.linalg.inv(spatial_metric) + g_jt = spacetime_metric[1:, 0] + + shift = np.dot(inverse_spatial_metric, g_jt) + + return shift + + +def lapse(spacetime_metric): + spatial_metric = spacetime_metric[1:, 1:] + inverse_spatial_metric = np.linalg.inv(spatial_metric) + g_jt = spacetime_metric[1:, 0] + shift = np.dot(inverse_spatial_metric, g_jt) + beta_g_it = np.dot(shift, g_jt) + g_tt = spacetime_metric[0, 0] + + lapse = np.sqrt(beta_g_it - g_tt) + + return lapse + + +def conformal_factor_minus_one(x): + return 0.0 + + +def lapse_times_conformal_factor_minus_one(x): + boosted_spacetime_metric_left = boost_spacetime_metric( + spacetime_left(x), np.array(momentum_left) / masses[0] + ) + boosted_spacetime_metric_right = boost_spacetime_metric( + spacetime_right(x), np.array(momentum_right) / masses[1] + ) + lapse_left = lapse(boosted_spacetime_metric_left) + lapse_right = lapse(boosted_spacetime_metric_right) + + return lapse_left * lapse_right - 1.0 + + +def shift_excess(x): + boosted_spacetime_metric_left = boost_spacetime_metric( + spacetime_left(x), np.array(momentum_left) / masses[0] + ) + boosted_spacetime_metric_right = boost_spacetime_metric( + spacetime_right(x), np.array(momentum_right) / masses[1] + ) + shift_left = shift(boosted_spacetime_metric_left) + shift_right = shift(boosted_spacetime_metric_right) + + return shift_left + shift_right diff --git a/tests/Unit/Elliptic/Systems/Xcts/BoundaryConditions/Test_SuperposedBoostedBinary.cpp b/tests/Unit/Elliptic/Systems/Xcts/BoundaryConditions/Test_SuperposedBoostedBinary.cpp new file mode 100644 index 000000000000..8295b875fce0 --- /dev/null +++ b/tests/Unit/Elliptic/Systems/Xcts/BoundaryConditions/Test_SuperposedBoostedBinary.cpp @@ -0,0 +1,181 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Framework/TestingFramework.hpp" + +#include +#include +#include +#include + +#include "DataStructures/DataBox/DataBox.hpp" +#include "DataStructures/DataVector.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "Domain/Structure/Direction.hpp" +#include "Domain/Tags.hpp" +#include "Domain/Tags/FaceNormal.hpp" +#include "Elliptic/BoundaryConditions/ApplyBoundaryCondition.hpp" +#include "Elliptic/BoundaryConditions/BoundaryCondition.hpp" +#include "Elliptic/BoundaryConditions/BoundaryConditionType.hpp" +#include "Elliptic/Systems/Xcts/BoundaryConditions/SuperposedBoostedBinary.hpp" +#include "Elliptic/Systems/Xcts/FluxesAndSources.hpp" +#include "Framework/CheckWithRandomValues.hpp" +#include "Framework/Pypp.hpp" +#include "Framework/SetupLocalPythonEnvironment.hpp" +#include "Framework/TestCreation.hpp" +#include "Framework/TestHelpers.hpp" +#include "NumericalAlgorithms/DiscontinuousGalerkin/NormalDotFlux.hpp" +#include "PointwiseFunctions/AnalyticSolutions/Xcts/Schwarzschild.hpp" +#include "PointwiseFunctions/InitialDataUtilities/AnalyticSolution.hpp" +#include "PointwiseFunctions/Xcts/LongitudinalOperator.hpp" +#include "Utilities/Gsl.hpp" +#include "Utilities/TMPL.hpp" + +namespace Xcts::BoundaryConditions { + +namespace { + +const std::string py_module{ + "Elliptic.Systems.Xcts.BoundaryConditions.SuperposedBoostedBinary"}; + +tnsr::i make_spherical_face_normal_flat_cartesian( + tnsr::I x, const std::array& center) { + for (size_t d = 0; d < 3; ++d) { + x.get(d) -= gsl::at(center, d); + } + Scalar euclidean_radius = magnitude(x); + tnsr::i face_normal{x.begin()->size()}; + get<0>(face_normal) = -get<0>(x) / get(euclidean_radius); + get<1>(face_normal) = -get<1>(x) / get(euclidean_radius); + get<2>(face_normal) = -get<2>(x) / get(euclidean_radius); + return face_normal; +} + +struct FactoryMetavars { + struct factory_creation + : tt::ConformsTo { + using factory_classes = + tmpl::map, + tmpl::list>>>, + tmpl::pair>>; + }; +}; + +template +void test_suite(const std::string& options_string) { + INFO("Test factory-creation"); + using IsolatedObjectBase = elliptic::analytic_data::AnalyticSolution; + using IsolatedObjectClasses = tmpl::list; + register_classes_with_charm(); + register_factory_classes_with_charm(); + const auto created = TestHelpers::test_creation< + std::unique_ptr>, + FactoryMetavars>(options_string); + const auto serialized = serialize_and_deserialize(created); + REQUIRE(dynamic_cast*>( + serialized.get()) != nullptr); + const auto& boundary_condition = dynamic_cast&>(*serialized); + { + INFO("Properties"); + CHECK(boundary_condition.boundary_condition_types() == + std::vector{ + 5, elliptic::BoundaryConditionType::Dirichlet}); + } + { + MAKE_GENERATOR(gen); + std::uniform_real_distribution<> dist(-1., 1.); + const size_t num_points = 3; + const auto direction = Direction<3>::upper_zeta(); + const std::array center{{0., 0., 0.}}; + const auto x = make_with_random_values>( + make_not_null(&gen), make_not_null(&dist), num_points); + const auto face_normal = + make_spherical_face_normal_flat_cartesian(x, center); + const auto box = db::create>, + domain::Tags::Faces<3, domain::Tags::FaceNormal<3, Frame::Inertial>>>>( + DirectionMap<3, tnsr::I>{{direction, x}}, + DirectionMap<3, tnsr::i>{{direction, face_normal}}); + auto conformal_factor_minus_one = + make_with_random_values>( + make_not_null(&gen), make_not_null(&dist), num_points); + auto lapse_times_conformal_factor_minus_one = + make_with_random_values>( + make_not_null(&gen), make_not_null(&dist), num_points); + auto shift_excess = make_with_random_values>( + make_not_null(&gen), make_not_null(&dist), num_points); + auto n_dot_conformal_factor_gradient = + make_with_random_values>( + make_not_null(&gen), make_not_null(&dist), num_points); + auto n_dot_lapse_times_conformal_factor_gradient = + make_with_random_values>( + make_not_null(&gen), make_not_null(&dist), num_points); + auto n_dot_longitudinal_shift_excess = + make_with_random_values>( + make_not_null(&gen), make_not_null(&dist), num_points); + const tnsr::i deriv_conformal_factor{ + num_points, std::numeric_limits::signaling_NaN()}; + const tnsr::i deriv_lapse_times_conformal_factor{ + num_points, std::numeric_limits::signaling_NaN()}; + const tnsr::iJ deriv_shift_excess{ + num_points, std::numeric_limits::signaling_NaN()}; + + elliptic::apply_boundary_condition< + Linearized, void, + tmpl::list>>( + boundary_condition, box, direction, + make_not_null(&conformal_factor_minus_one), + make_not_null(&lapse_times_conformal_factor_minus_one), + make_not_null(&shift_excess), + make_not_null(&n_dot_conformal_factor_gradient), + make_not_null(&n_dot_lapse_times_conformal_factor_gradient), + make_not_null(&n_dot_longitudinal_shift_excess), deriv_conformal_factor, + deriv_lapse_times_conformal_factor, deriv_shift_excess); + + const auto expected_conformal_factor_minus_one = + pypp::call>(py_module, "conformal_factor_minus_one", + x); + const auto expected_lapse_times_conformal_factor_minus_one = + pypp::call>( + py_module, "lapse_times_conformal_factor_minus_one", x); + const auto expected_shift_excess = + pypp::call>(py_module, "shift_excess", x); + + CHECK_ITERABLE_APPROX(get(conformal_factor_minus_one), + get(expected_conformal_factor_minus_one)); + CHECK_ITERABLE_APPROX(get(lapse_times_conformal_factor_minus_one), + get(expected_lapse_times_conformal_factor_minus_one)); + CHECK_ITERABLE_APPROX(get<0>(shift_excess), get<0>(expected_shift_excess)); + CHECK_ITERABLE_APPROX(get<1>(shift_excess), get<1>(expected_shift_excess)); + CHECK_ITERABLE_APPROX(get<2>(shift_excess), get<2>(expected_shift_excess)); + } +} +} // namespace + +SPECTRE_TEST_CASE("Unit.Xcts.BoundaryConditions.SuperposedBoostedBinary", + "[Unit][Elliptic]") { + const pypp::SetupLocalPythonEnvironment local_python_env(""); + test_suite( + "SuperposedBoostedBinary:\n" + " XCoords: [-5., 6.]\n" + " Masses: [1.1, 0.43]\n" + " MomentumLeft: [0.01, 0.01, 0.01]\n" + " MomentumRight: [-0.01, -0.01, -0.01]\n" + " CenterOfMassOffset: [0.02, 0.01]\n" + " ObjectLeft:\n" + " Schwarzschild:\n" + " Mass: 1.1\n" + " Coordinates: Isotropic\n" + " ObjectRight:\n" + " Schwarzschild:\n" + " Mass: 0.43\n" + " Coordinates: Isotropic"); +} + +} // namespace Xcts::BoundaryConditions diff --git a/tests/Unit/PointwiseFunctions/AnalyticData/Xcts/BinaryWithGravitationalWaves.py b/tests/Unit/PointwiseFunctions/AnalyticData/Xcts/BinaryWithGravitationalWaves.py new file mode 100644 index 000000000000..435fb5106692 --- /dev/null +++ b/tests/Unit/PointwiseFunctions/AnalyticData/Xcts/BinaryWithGravitationalWaves.py @@ -0,0 +1,201 @@ +# Distributed under the MIT License. +# See LICENSE.txt for details. + +import numpy as np + +x_coords = [-5.0, 6.0] +momentum_left = [-0.01, -0.01, -0.01] +momentum_right = [0.01, 0.01, 0.01] +y_offset = 0.02 +z_offset = 0.01 +masses = [1.1, 0.43] +attenuation_parameter = 1.1 +attenuation_radius = 3.6 + + +def x_left(x): + x_left = np.array(x) + x_left[0] -= x_coords[0] + x_left[1] -= y_offset + x_left[2] -= z_offset + + velocity = np.array(momentum_left) / masses[0] + beta2 = np.dot(velocity, velocity) + gamma = 1.0 / np.sqrt(1.0 - beta2) + + boost_matrix = np.zeros((4, 4)) + boost_matrix[0, 0] = gamma + boost_matrix[0, 1:] = -gamma * velocity + boost_matrix[1:, 0] = -gamma * velocity + boost_matrix[1:, 1:] = ( + np.identity(3) + (gamma - 1.0) * np.outer(velocity, velocity) / beta2 + ) + + boosted_vector = np.dot(boost_matrix, np.append(0.0, x_left)) + + return boosted_vector[1:] + + +def x_right(x): + x_right = np.array(x) + x_right[0] -= x_coords[1] + x_right[1] -= y_offset + x_right[2] -= z_offset + + velocity = np.array(momentum_right) / masses[1] + beta2 = np.dot(velocity, velocity) + gamma = 1.0 / np.sqrt(1.0 - beta2) + + boost_matrix = np.zeros((4, 4)) + boost_matrix[0, 0] = gamma + boost_matrix[0, 1:] = -gamma * velocity + boost_matrix[1:, 0] = -gamma * velocity + boost_matrix[1:, 1:] = ( + np.identity(3) + (gamma - 1.0) * np.outer(velocity, velocity) / beta2 + ) + + boosted_vector = np.dot(boost_matrix, np.append(0.0, x_right)) + + return boosted_vector[1:] + + +def spacetime_left(x): + r = np.linalg.norm(x_left(x)) + conformalfactor = 1.0 + 0.5 * masses[0] / r + lapse = (1.0 - 0.5 * masses[0] / r) / conformalfactor + shift = np.zeros(3) + spatial_metric = np.identity(3) * pow(conformalfactor, 4) + + spacetime_metric = np.zeros((4, 4)) + spacetime_metric[0, 0] = -(lapse**2) + np.dot( + shift, np.dot(spatial_metric, shift) + ) + spacetime_metric[0, 1:] = np.dot(spatial_metric, shift) + spacetime_metric[1:, 0] = spacetime_metric[0, 1:] + spacetime_metric[1:, 1:] = spatial_metric + + return spacetime_metric + + +def spacetime_right(x): + r = np.linalg.norm(x_right(x)) + conformalfactor = 1.0 + 0.5 * masses[1] / r + lapse = (1.0 - 0.5 * masses[1] / r) / conformalfactor + shift = np.zeros(3) + spatial_metric = np.identity(3) * pow(conformalfactor, 4) + + spacetime_metric = np.zeros((4, 4)) + spacetime_metric[0, 0] = -(lapse**2) + np.dot( + shift, np.dot(spatial_metric, shift) + ) + spacetime_metric[0, 1:] = np.dot(spatial_metric, shift) + spacetime_metric[1:, 0] = spacetime_metric[0, 1:] + spacetime_metric[1:, 1:] = spatial_metric + + return spacetime_metric + + +def boost_spacetime_metric(spacetime_metric, velocity): + beta2 = np.dot(velocity, velocity) + gamma = 1.0 / np.sqrt(1.0 - beta2) + + boost_matrix = np.zeros((4, 4)) + boost_matrix[0, 0] = gamma + boost_matrix[0, 1:] = -gamma * velocity + boost_matrix[1:, 0] = -gamma * velocity + boost_matrix[1:, 1:] = ( + np.identity(3) + (gamma - 1.0) * np.outer(velocity, velocity) / beta2 + ) + + return np.dot(boost_matrix, np.dot(spacetime_metric, boost_matrix.T)) + + +def shift(spacetime_metric): + spatial_metric = spacetime_metric[1:, 1:] + inverse_spatial_metric = np.linalg.inv(spatial_metric) + g_jt = spacetime_metric[1:, 0] + + shift = np.dot(inverse_spatial_metric, g_jt) + + return shift + + +def lapse(spacetime_metric): + spatial_metric = spacetime_metric[1:, 1:] + inverse_spatial_metric = np.linalg.inv(spatial_metric) + g_jt = spacetime_metric[1:, 0] + shift = np.dot(inverse_spatial_metric, g_jt) + beta_g_it = np.dot(shift, g_jt) + g_tt = spacetime_metric[0, 0] + + lapse = np.sqrt(beta_g_it - g_tt) + + return lapse + + +def superposed_spacetime_metric(x): + spacetime_metric_left = spacetime_left(x) + spacetime_metric_right = spacetime_right(x) + boosted_spacetime_metric_left = boost_spacetime_metric( + spacetime_metric_left, np.array(momentum_left) / masses[0] + ) + boosted_spacetime_metric_right = boost_spacetime_metric( + spacetime_metric_right, np.array(momentum_right) / masses[1] + ) + + spatial_metric_left = boosted_spacetime_metric_left[1:, 1:] + spatial_metric_right = boosted_spacetime_metric_right[1:, 1:] + shift_left = shift(boosted_spacetime_metric_left) + shift_right = shift(boosted_spacetime_metric_right) + lapse_left = lapse(boosted_spacetime_metric_left) + lapse_right = lapse(boosted_spacetime_metric_right) + + superposed_spatial_metric = ( + spatial_metric_left + spatial_metric_right - np.identity(3) + ) + superposed_shift = shift_left + shift_right + superposed_lapse = lapse_left * lapse_right + + superposed_spacetime_metric = np.zeros((4, 4)) + superposed_spacetime_metric[0, 0] = -(superposed_lapse**2) + np.dot( + superposed_shift, np.dot(superposed_spatial_metric, superposed_shift) + ) + superposed_spacetime_metric[0, 1:] = np.dot( + superposed_spatial_metric, superposed_shift + ) + superposed_spacetime_metric[1:, 0] = superposed_spacetime_metric[0, 1:] + superposed_spacetime_metric[1:, 1:] = superposed_spatial_metric + + return superposed_spacetime_metric + + +def conformal_metric_bbh_isotropic(x): + return superposed_spacetime_metric(x)[1:, 1:] + + +def inv_conformal_metric_bbh_isotropic(x): + return np.linalg.inv(conformal_metric_bbh_isotropic(x)) + + +def shift_background(x): + return np.zeros(3) + + +def longitudinal_shift_background_bbh_isotropic(x): + return np.zeros((3, 3)) + + +def conformal_factor_minus_one_bbh_isotropic(x): + return 0.0 + + +def energy_density_bbh_isotropic(x): + return 0.0 + + +def stress_trace_bbh_isotropic(x): + return 0.0 + + +def momentum_density_bbh_isotropic(x): + return np.zeros(3) diff --git a/tests/Unit/PointwiseFunctions/AnalyticData/Xcts/CMakeLists.txt b/tests/Unit/PointwiseFunctions/AnalyticData/Xcts/CMakeLists.txt index ba305bc1313a..f0c86d72a506 100644 --- a/tests/Unit/PointwiseFunctions/AnalyticData/Xcts/CMakeLists.txt +++ b/tests/Unit/PointwiseFunctions/AnalyticData/Xcts/CMakeLists.txt @@ -5,6 +5,7 @@ set(LIBRARY "Test_XctsAnalyticData") set(LIBRARY_SOURCES Test_Binary.cpp + Test_BinaryWithGravitationalWaves.cpp ) add_test_library(${LIBRARY} "${LIBRARY_SOURCES}") diff --git a/tests/Unit/PointwiseFunctions/AnalyticData/Xcts/Test_BinaryWithGravitationalWaves.cpp b/tests/Unit/PointwiseFunctions/AnalyticData/Xcts/Test_BinaryWithGravitationalWaves.cpp new file mode 100644 index 000000000000..8f5190950d30 --- /dev/null +++ b/tests/Unit/PointwiseFunctions/AnalyticData/Xcts/Test_BinaryWithGravitationalWaves.cpp @@ -0,0 +1,152 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Framework/TestingFramework.hpp" + +#include +#include +#include + +#include "DataStructures/DataBox/PrefixHelpers.hpp" +#include "DataStructures/DataBox/Prefixes.hpp" +#include "DataStructures/DataVector.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "Elliptic/Systems/Xcts/Tags.hpp" +#include "Framework/CheckWithRandomValues.hpp" +#include "Framework/SetupLocalPythonEnvironment.hpp" +#include "Framework/TestCreation.hpp" +#include "Framework/TestHelpers.hpp" +#include "PointwiseFunctions/AnalyticData/Xcts/BinaryWithGravitationalWaves.hpp" +#include "PointwiseFunctions/AnalyticSolutions/Xcts/Schwarzschild.hpp" +#include "PointwiseFunctions/InitialDataUtilities/AnalyticSolution.hpp" +#include "Utilities/Serialization/RegisterDerivedClassesWithCharm.hpp" +#include "Utilities/TMPL.hpp" +#include "Utilities/TaggedTuple.hpp" + +namespace Xcts::AnalyticData { +namespace { + +using test_tags = tmpl::list< + Tags::ConformalMetric, + Tags::InverseConformalMetric, + Tags::ShiftBackground, + Tags::LongitudinalShiftBackgroundMinusDtConformalMetric, + Tags::ConformalFactorMinusOne, + gr::Tags::Conformal, 0>, + gr::Tags::Conformal, 0>, + gr::Tags::Conformal, 0>>; + +template +struct BinaryWithGravitationalWavesProxy { + tuples::tagged_tuple_from_typelist test_variables( + const tnsr::I& x) const { + return BinaryWithGravitationalWaves_ptr->variables(x, test_tags{}); + } + + BinaryWithGravitationalWaves* + BinaryWithGravitationalWaves_ptr; +}; + +struct Metavariables { + struct factory_creation + : tt::ConformsTo { + using factory_classes = + tmpl::map>>>, + tmpl::pair>>; + }; +}; + +void test_data(const std::array& x_coords, + const std::array& masses, + const std::array& momentum_left, + const std::array& momentum_right, + const std::array& center_of_mass_offset, + const double& attenuation_parameter, + const double& attenuation_radius, + const std::string& py_functions_suffix, + const std::string& options_string) { + using IsolatedObjectBase = elliptic::analytic_data::AnalyticSolution; + using IsolatedObjectClasses = tmpl::list; + register_classes_with_charm(); + const auto created = TestHelpers::test_creation< + std::unique_ptr, Metavariables>( + options_string); + REQUIRE(dynamic_cast*>(created.get()) != + nullptr); + const auto& derived = dynamic_cast&>(*created); + auto BinaryWithGravitationalWaves = serialize_and_deserialize(derived); + { + INFO("Properties"); + CHECK(BinaryWithGravitationalWaves.x_coords() == x_coords); + CHECK(BinaryWithGravitationalWaves.y_offset() == center_of_mass_offset[0]); + CHECK(BinaryWithGravitationalWaves.z_offset() == center_of_mass_offset[1]); + CHECK(BinaryWithGravitationalWaves.momentum_left() == momentum_left); + CHECK(BinaryWithGravitationalWaves.momentum_right() == momentum_right); + const auto& superposed_objects = + BinaryWithGravitationalWaves.superposed_objects(); + CHECK(dynamic_cast( + *superposed_objects[0]) + .mass() == masses[0]); + CHECK(dynamic_cast( + *superposed_objects[1]) + .mass() == masses[1]); + CHECK(BinaryWithGravitationalWaves.attenuation_parameter() == + attenuation_parameter); + CHECK(BinaryWithGravitationalWaves.attenuation_radius() == + attenuation_radius); + } + { + const BinaryWithGravitationalWavesProxy + proxy{&BinaryWithGravitationalWaves}; + pypp::check_with_random_values<1>( + &BinaryWithGravitationalWavesProxy< + IsolatedObjectBase, IsolatedObjectClasses>::test_variables, + proxy, "BinaryWithGravitationalWaves", + {"conformal_metric_" + py_functions_suffix, + "inv_conformal_metric_" + py_functions_suffix, "shift_background", + "longitudinal_shift_background_" + py_functions_suffix, + "conformal_factor_minus_one_" + py_functions_suffix, + "energy_density_" + py_functions_suffix, + "stress_trace_" + py_functions_suffix, + "momentum_density_" + py_functions_suffix}, + {{{x_coords[0] * 2, x_coords[1] * 2}}}, std::make_tuple(), + DataVector(5)); + } +} + +} // namespace + +SPECTRE_TEST_CASE( + "Unit.PointwiseFunctions.AnalyticData.Xcts.BinaryWithGravitationalWaves", + "[PointwiseFunctions][Unit]") { + const pypp::SetupLocalPythonEnvironment local_python_env{ + "PointwiseFunctions/AnalyticData/Xcts"}; + test_data({{-5., 6.}}, {{1.1, 0.43}}, {{-0.01, -0.01, -0.01}}, + {{0.01, 0.01, 0.01}}, {{0.02, 0.01}}, 1.1, 3.6, "bbh_isotropic", + "BinaryWithGravitationalWaves:\n" + " XCoords: [-5., 6.]\n" + " Masses: [1.1, 0.43]\n" + " MomentumLeft: [-0.01, -0.01, -0.01]\n" + " MomentumRight: [0.01, 0.01, 0.01]\n" + " CenterOfMassOffset: [0.02, 0.01]\n" + " ObjectLeft:\n" + " Schwarzschild:\n" + " Mass: 1.1\n" + " Coordinates: Isotropic\n" + " ObjectRight:\n" + " Schwarzschild:\n" + " Mass: 0.43\n" + " Coordinates: Isotropic\n" + " AttenuationParameter: 1.1\n" + " AttenuationRadius: 3.6"); +} + +} // namespace Xcts::AnalyticData From 1b1e6eb6e475c8e1d1d0a1727748857a6191705d Mon Sep 17 00:00:00 2001 From: Megum Rebelo Date: Thu, 6 Feb 2025 17:08:47 +0000 Subject: [PATCH 2/3] Fixup: comments from PR by Hannes --- docs/References.bib | 25 +- .../SuperposedBoostedBinary.hpp | 15 +- .../Xcts/BinaryWithGravitationalWaves.cpp | 92 +++---- .../Xcts/BinaryWithGravitationalWaves.hpp | 13 +- .../Xcts/BinaryWithGravitationalWaves.yaml | 246 ++++++++++++++++++ .../Test_BinaryWithGravitationalWaves.cpp | 4 +- 6 files changed, 326 insertions(+), 69 deletions(-) create mode 100644 tests/InputFiles/Xcts/BinaryWithGravitationalWaves.yaml diff --git a/docs/References.bib b/docs/References.bib index 9a68a4766643..f52e36145050 100644 --- a/docs/References.bib +++ b/docs/References.bib @@ -1360,9 +1360,9 @@ @Article{Hou2007 } @article{Jaranowski1997ky, - author = "Jaranowski, Piotr and Schaefer, Gerhard", + author = "{Jaranowski, Piotr and Sch\"afer, Gerhard}", title = "{Third postNewtonian higher order ADM Hamilton dynamics for - two-body point mass systems}", + two-body point-mass systems}", eprint = "gr-qc/9712075", archivePrefix = "arXiv", doi = "10.1103/PhysRevD.57.7274", @@ -2554,6 +2554,20 @@ @article{Stamm2010 volume = "79", } +@article{Steinhoff2008zr, + author = "Steinhoff, Jan and Schaefer, Gerhard and Hergt, Steven", + title = "{ADM canonical formalism for gravitating spinning objects}", + eprint = "0805.3136", + archivePrefix = "arXiv", + primaryClass = "gr-qc", + doi = "10.1103/PhysRevD.77.104018", + url = "https://doi.org/10.1103/PhysRevD.77.104018", + journal = "Phys. Rev. D", + volume = "77", + pages = "104018", + year = "2008" +} + @article{Stiller2016a, author = "Stiller, Jörg", title = "Robust multigrid for high-order discontinuous {Galerkin} @@ -2682,14 +2696,15 @@ @book{ThorneBlandford2017 } @article{Tichy2002ec, - author = "Tichy, Wolfgang and Bruegmann, Bernd and Campanelli, Manuela - and Diener, Peter", + author = "{Tichy, Wolfgang and Br\"ugmann, Bernd and Campanelli, Manuela + and Diener, Peter}", title = "{Binary black hole initial data for numerical general relativity - based on postNewtonian data}", + based on post-Newtonian data}", eprint = "gr-qc/0207011", archivePrefix = "arXiv", reportNumber = "UTBRG-2002-06, AEI-2002-048", doi = "10.1103/PhysRevD.67.064008", + url = {https://link.aps.org/doi/10.1103/PhysRevD.67.064008}, journal = "Phys. Rev. D", volume = "67", pages = "064008", diff --git a/src/Elliptic/Systems/Xcts/BoundaryConditions/SuperposedBoostedBinary.hpp b/src/Elliptic/Systems/Xcts/BoundaryConditions/SuperposedBoostedBinary.hpp index 72f8216593fc..635683f6fd3d 100644 --- a/src/Elliptic/Systems/Xcts/BoundaryConditions/SuperposedBoostedBinary.hpp +++ b/src/Elliptic/Systems/Xcts/BoundaryConditions/SuperposedBoostedBinary.hpp @@ -32,10 +32,21 @@ class DataVector; namespace Xcts::BoundaryConditions { /*! - * \brief Impose supperposed boosted binary system on the boundary. + * \brief Impose superposed boosted binary system on the boundary. * * This takes two isolated objects and after applying a boost to each of them, - * superposes them. The superposed system is then imposed on the boundary, with + * the isolated lapses and shifts are extracted and superposed, following + * + * \f{align} + * \alpha = \alpha^1 \alpha^2 + * \text{ and } + * \beta^i = \beta^i_1 + \beta^i_2. + * \f} + * + * The conformal factor is set to 1 and the background spatial + * metric should contain the boosted isolated objects. + * + * The superposed system is then imposed on the boundary with * Dirichlet boundary conditions. * */ diff --git a/src/PointwiseFunctions/AnalyticData/Xcts/BinaryWithGravitationalWaves.cpp b/src/PointwiseFunctions/AnalyticData/Xcts/BinaryWithGravitationalWaves.cpp index 0da205119327..90ca469cd9bf 100644 --- a/src/PointwiseFunctions/AnalyticData/Xcts/BinaryWithGravitationalWaves.cpp +++ b/src/PointwiseFunctions/AnalyticData/Xcts/BinaryWithGravitationalWaves.cpp @@ -120,7 +120,7 @@ void BinaryWithGravitationalWavesVariables::operator()( "Numeric differentiation only works with DataVectors because it needs " "a grid."); } - // Third order + // Third order backwards finite difference get(*dt_trace_extrinsic_curvature) = (11. * get(trace_extrinsic_curvature) - 18. * get(trace_extrinsic_curvature_back) + @@ -156,7 +156,7 @@ void BinaryWithGravitationalWavesVariables::operator()( tnsr::ii dt_conformal_metric{get_size(x.get(0))}; for (size_t i = 0; i < 3; ++i) { for (size_t j = 0; j <= i; ++j) { - // Third order + // Third order backwards finite difference dt_conformal_metric.get(i, j) = (11. * conformal_metric.get(i, j) - 18. * conformal_metric_back.get(i, j) + @@ -166,18 +166,19 @@ void BinaryWithGravitationalWavesVariables::operator()( } } + const auto auxiliar_one = tenex::evaluate( + inv_conformal_metric(ti::I, ti::K) * inv_conformal_metric(ti::J, ti::L) * + dt_conformal_metric(ti::k, ti::l)); + const auto auxiliar_two = tenex::evaluate( + inv_conformal_metric(ti::I, ti::K) * inv_conformal_metric(ti::J, ti::L) * + conformal_metric(ti::k, ti::l)); + for (size_t i = 0; i < 3; ++i) { for (size_t j = 0; j <= i; ++j) { - for (size_t k = 0; k < 3; ++k) { - for (size_t l = 0; l < 3; ++l) { - longitudinal_shift_background_minus_dt_conformal_metric->get(i, j) -= - inv_conformal_metric.get(i, k) * inv_conformal_metric.get(j, l) * - (dt_conformal_metric.get(k, l) - - (1. / 3.) * - get(trace(dt_conformal_metric, inv_conformal_metric)) * - conformal_metric.get(k, l)); - } - } + longitudinal_shift_background_minus_dt_conformal_metric->get(i, j) -= + (auxiliar_one.get(i, j) - + (1. / 3.) * get(trace(dt_conformal_metric, inv_conformal_metric)) * + auxiliar_two.get(i, j)); } } } @@ -240,7 +241,7 @@ Scalar BinaryWithGravitationalWavesVariables:: for (size_t i = 0; i < 3; ++i) { for (size_t j = 0; j <= i; ++j) { - // Third order + // Third order backwards finite difference (over -2 alpha) extrinsic_curvature.get(i, j) += (11. * conformal_metric.get(i, j) - 18. * conformal_metric_back.get(i, j) + @@ -289,6 +290,22 @@ BinaryWithGravitationalWavesVariables::get_t_conformal_metric( return conformal_metric; } +template +DataType +BinaryWithGravitationalWavesVariables::get_t_attenuation_function( + DataType t) const { + const auto distance_left_t = get_t_distance_left(t); + const auto distance_right_t = get_t_distance_right(t); + double turn_off = .5; + if (attenuation_parameter == 0) { + turn_off = 1.; + } + return (turn_off + .5 * tanh(attenuation_parameter * + (get(distance_left_t) - attenuation_radius))) * + (turn_off + .5 * tanh(attenuation_parameter * + (get(distance_right_t) - attenuation_radius))); +} + template tnsr::ii BinaryWithGravitationalWavesVariables::get_t_radiative_term( @@ -299,18 +316,12 @@ BinaryWithGravitationalWavesVariables::get_t_radiative_term( const auto present_term_t = get_t_present_term(t); const auto past_term_t = get_t_past_term(t); const auto integral_term_t = get_t_integral_term(t); - double turn_off = .5; - if (attenuation_parameter == 0) { - turn_off = 1.; - } + const auto attenuation_t = get_t_attenuation_function(t); tnsr::ii radiative_term_t{t.size()}; for (size_t i = 0; i < 3; ++i) { for (size_t j = 0; j <= i; ++j) { radiative_term_t.get(i, j) = - (turn_off + .5 * tanh(attenuation_parameter * - (get(distance_left_t) - attenuation_radius))) * - (turn_off + .5 * tanh(attenuation_parameter * - (get(distance_right_t) - attenuation_radius))) * + attenuation_t * (near_zone_term_t.get(i, j) + present_term_t.get(i, j) + past_term_t.get(i, j) + integral_term_t.get(i, j)); } @@ -370,6 +381,7 @@ Scalar BinaryWithGravitationalWavesVariables::get_t_lapse( const auto separation_t = get_t_separation(t); const auto momentum_left_t = get_t_momentum_left(t); const auto momentum_right_t = get_t_momentum_right(t); + const auto attenuation_t = get_t_attenuation_function(t); get(lapse_t) = 1 - mass_left / get(distance_left_t) - mass_right / get(distance_right_t) + @@ -383,15 +395,7 @@ Scalar BinaryWithGravitationalWavesVariables::get_t_lapse( (get(distance_left_t) * mass_left) + get(dot_product(momentum_right_t, momentum_right_t)) / (get(distance_right_t) * mass_right))); - double turn_off = .5; - if (attenuation_parameter == 0) { - turn_off = 1.; - } - get(lapse_t) *= - (turn_off + .5 * tanh(attenuation_parameter * - (get(distance_left_t) - attenuation_radius))) * - (turn_off + .5 * tanh(attenuation_parameter * - (get(distance_right_t) - attenuation_radius))); + get(lapse_t) *= attenuation_t; // Boosted lapse in the inner zone const auto spacetime_metric_left = get_t_boosted_spacetime_metric_left(t); @@ -409,14 +413,7 @@ Scalar BinaryWithGravitationalWavesVariables::get_t_lapse( const auto shift_right = gr::shift(spacetime_metric_right, inv_conformal_metric_right); const auto lapse_right = gr::lapse(shift_right, spacetime_metric_right); - get(lapse_t) += - (1. - - (turn_off + .5 * tanh(attenuation_parameter * - (get(distance_left_t) - attenuation_radius))) * - (turn_off + - .5 * tanh(attenuation_parameter * - (get(distance_right_t) - attenuation_radius)))) * - get(lapse_left) * get(lapse_right); + get(lapse_t) += (1. - attenuation_t) * get(lapse_left) * get(lapse_right); return lapse_t; } @@ -434,6 +431,7 @@ BinaryWithGravitationalWavesVariables::get_t_shift(DataType t) const { const auto momentum_right_t = get_t_momentum_right(present_time); const auto normal_left_t = get_t_normal_left(present_time); const auto normal_right_t = get_t_normal_right(present_time); + const auto attenuation_t = get_t_attenuation_function(present_time); for (size_t i = 0; i < 3; ++i) { shift_t.get(i) -= 4. * (momentum_left_t.get(i) / get(distance_left_t) + momentum_right_t.get(i) / get(distance_right_t)); @@ -449,16 +447,8 @@ BinaryWithGravitationalWavesVariables::get_t_shift(DataType t) const { 0.5 * momentum_right_t.get(i) / get(distance_right_t); } - double turn_off = .5; - if (attenuation_parameter == 0) { - turn_off = 1.; - } for (size_t i = 0; i < 3; ++i) { - shift_t.get(i) *= - (turn_off + .5 * tanh(attenuation_parameter * - (get(distance_left_t) - attenuation_radius))) * - (turn_off + .5 * tanh(attenuation_parameter * - (get(distance_right_t) - attenuation_radius))); + shift_t.get(i) *= attenuation_t; } // Boosted shifts in the inner zone @@ -479,13 +469,7 @@ BinaryWithGravitationalWavesVariables::get_t_shift(DataType t) const { const auto lapse_right = gr::lapse(shift_right, spacetime_metric_right); for (size_t i = 0; i < 3; ++i) { shift_t.get(i) += - (1. - - (turn_off + .5 * tanh(attenuation_parameter * - (get(distance_left_t) - attenuation_radius))) * - (turn_off + - .5 * tanh(attenuation_parameter * - (get(distance_right_t) - attenuation_radius)))) * - (shift_left.get(i) + shift_right.get(i)); + (1. - attenuation_t) * (shift_left.get(i) + shift_right.get(i)); } return shift_t; } diff --git a/src/PointwiseFunctions/AnalyticData/Xcts/BinaryWithGravitationalWaves.hpp b/src/PointwiseFunctions/AnalyticData/Xcts/BinaryWithGravitationalWaves.hpp index 960bf2e9f1e5..4aa6959ed99b 100644 --- a/src/PointwiseFunctions/AnalyticData/Xcts/BinaryWithGravitationalWaves.hpp +++ b/src/PointwiseFunctions/AnalyticData/Xcts/BinaryWithGravitationalWaves.hpp @@ -273,6 +273,7 @@ struct BinaryWithGravitationalWavesVariables DataType t, Mesh<3> local_mesh, InverseJacobian local_inv_jacobian) const; + DataType get_t_attenuation_function(DataType t) const; tnsr::ii get_t_conformal_metric(DataType t) const; tnsr::ii get_t_radiative_term(DataType t) const; tnsr::ii get_t_near_zone_term(DataType t) const; @@ -294,12 +295,12 @@ struct BinaryWithGravitationalWavesVariables * The main goal of this implementation is to improve the extracted * wave forms, for example, by minimizing junk radiation. * The data is only valid for black holes without spin. Even so, there is some - * work done to describe such systems that could later be implemented. - * The objects are constructed from a superposition of two isolated objects that - * are boosted with respect to each other. The radiative data is constructed - * from Post-Newtonian expansions for the inspiral phase, in orders of - * \f$\epsilon = 1/c\f$, in \cite Jaranowski1997ky. In ADMTT gauge it is - * possible to get the 3-metric as \f$\gamma^{PN}_{ij} = \psi^{4}_{PN} + * work done to describe such systems that could later be implemented, \cite + * Steinhoff2008zr. The objects are constructed from a superposition of two + * isolated objects that are boosted with respect to each other. The radiative + * data is constructed from Post-Newtonian expansions for the inspiral phase, in + * orders of \f$\epsilon = 1/c\f$, in \cite Jaranowski1997ky. In ADMTT gauge it + * is possible to get the 3-metric as \f$\gamma^{PN}_{ij} = \psi^{4}_{PN} * \delta_{ij} + h^{TT}_{ij}\f$ where \f$h^{TT}_{ij}\f$ is the radiative part * and the non-radiative Post-Newtonian conformal factor is given by * diff --git a/tests/InputFiles/Xcts/BinaryWithGravitationalWaves.yaml b/tests/InputFiles/Xcts/BinaryWithGravitationalWaves.yaml new file mode 100644 index 000000000000..c0b1a667392f --- /dev/null +++ b/tests/InputFiles/Xcts/BinaryWithGravitationalWaves.yaml @@ -0,0 +1,246 @@ +# Distributed under the MIT License. +# See LICENSE.txt for details. + +Executable: SolveXcts +Testing: + Check: parse + Timeout: 30 + Priority: High +ExpectedOutput: + - BwGWReductions.h5 + - BwGWVolume0.h5 + +--- + +Parallelization: + ElementDistribution: NumGridPoints + +ResourceInfo: + AvoidGlobalProc0: false + Singletons: Auto + +Background: &background + BinaryWithGravitationalWaves: + Masses: [&mass_left .5, &mass_right .5] + XCoords: [&x_left -6.0, &x_right 6.0] + MomentumLeft: [0., &momentum_left -0.08, 0.] + MomentumRight: [0., &momentum_right 0.08, 0.] + CenterOfMassOffset: ¢er_offset [&y_center 0.0, &z_center 0.0] + ObjectLeft: &object_left + Schwarzschild: + Mass: *mass_left + Coordinates: MaximalIsotropic + ObjectRight: &object_right + Schwarzschild: + Mass: *mass_right + Coordinates: MaximalIsotropic + AttenuationParameter: 1. + AttenuationRadius: 100. + +InitialGuess: *background + +RandomizeInitialGuess: None + +DomainCreator: + BinaryCompactObject: + ObjectA: + InnerRadius: .3 + OuterRadius: 4. + XCoord: *x_right + Interior: + ExciseWithBoundaryCondition: + SuperposedBoostedBinary: + Masses: [*mass_left, *mass_right] + XCoords: [*x_left, *x_right] + MomentumLeft: [0., *momentum_left, 0.] + MomentumRight: [0., *momentum_right, 0.] + CenterOfMassOffset: *center_offset + ObjectLeft: *object_left + ObjectRight: *object_right + UseLogarithmicMap: True + ObjectB: + InnerRadius: .3 + OuterRadius: 4. + XCoord: *x_left + Interior: + ExciseWithBoundaryCondition: + SuperposedBoostedBinary: + Masses: [*mass_left, *mass_right] + XCoords: [*x_left, *x_right] + MomentumLeft: [0., *momentum_left, 0.] + MomentumRight: [0., *momentum_right, 0.] + CenterOfMassOffset: *center_offset + ObjectLeft: *object_left + ObjectRight: *object_right + UseLogarithmicMap: True + CenterOfMassOffset: *center_offset + CubeScale: 1.0 + Envelope: + Radius: &outer_shell_inner_radius 60. + RadialDistribution: Logarithmic + OuterShell: + Radius: &outer_radius 70. + RadialDistribution: &outer_shell_distribution Linear + OpeningAngle: 120.0 + BoundaryCondition: Robin + UseEquiangularMap: True + InitialRefinement: + ObjectAShell: [0, 0, 0] + ObjectBShell: [0, 0, 0] + ObjectACube: [0, 0, 0] + ObjectBCube: [0, 0, 0] + Envelope: [0, 0, 1] + OuterShell: [0, 0, 2] + # This p-refinement represents a crude manual optimization of the domain. We + # will need AMR to optimize the domain further. + InitialGridPoints: + ObjectAShell: [6, 6, 10] + ObjectBShell: [6, 6, 10] + ObjectACube: [6, 6, 7] + ObjectBCube: [6, 6, 7] + Envelope: [6, 6, 6] + OuterShell: [6, 6, 6] + TimeDependentMaps: None + +Amr: + Verbosity: Verbose + Criteria: [] + Policies: + EnforceTwoToOneBalanceInNormalDirection: true + Isotropy: Anisotropic + Limits: + NumGridPoints: Auto + RefinementLevel: Auto + ErrorBeyondLimits: False + Iterations: 1 + +PhaseChangeAndTriggers: [] + +Discretization: + DiscontinuousGalerkin: + PenaltyParameter: 1.5 + Massive: True + Quadrature: GaussLobatto + Formulation: WeakInertial + +Observers: + VolumeFileName: "BwGWVolume" + ReductionFileName: "BwGWReductions" + +NonlinearSolver: + NewtonRaphson: + ConvergenceCriteria: + MaxIterations: 20 + RelativeResidual: 1.e-10 + AbsoluteResidual: 1.e-11 + SufficientDecrease: 1.e-4 + MaxGlobalizationSteps: 40 + DampingFactor: 1. + Verbosity: Verbose + +LinearSolver: + Gmres: + ConvergenceCriteria: + MaxIterations: 100 + RelativeResidual: 1.e-4 + AbsoluteResidual: 1.e-14 + Verbosity: Quiet + + Multigrid: + Iterations: 1 + MaxLevels: Auto + PreSmoothing: True + PostSmoothingAtBottom: True + Verbosity: Silent + OutputVolumeData: False + + SchwarzSmoother: + MaxOverlap: 2 + Iterations: 3 + Verbosity: Silent + SubdomainSolver: + Gmres: + ConvergenceCriteria: + MaxIterations: 3 + RelativeResidual: 1.e-4 + AbsoluteResidual: 1.e-10 + Verbosity: Silent + Restart: None + Preconditioner: + MinusLaplacian: + Solver: + ExplicitInverse: + WriteMatrixToFile: None + BoundaryConditions: Auto + SkipResets: True + ObservePerCoreReductions: False + +RadiallyCompressedCoordinates: + InnerRadius: *outer_shell_inner_radius + OuterRadius: *outer_radius + Compression: *outer_shell_distribution + +EventsAndTriggersAtIterations: + - Trigger: + EveryNIterations: + N: 1000000000 + Offset: 0 + Events: + - ObserveFields: + SubfileName: BackgroundData + VariablesToObserve: + - Conformal(SpatialMetric) + - TraceExtrinsicCurvature + - Lapse + - Shift + - ShiftBackground + - LongitudinalShiftBackgroundMinusDtConformalMetric + - dt(TraceExtrinsicCurvature) + - HamiltonianConstraint + - MomentumConstraint + BlocksToObserve: All + InterpolateToMesh: None + CoordinatesFloatingPointType: Double + FloatingPointTypes: [Double] + - Trigger: HasConverged + Events: + - ObserveNorms: + SubfileName: Norms + TensorsToObserve: + - Name: ConformalFactor + NormType: Max + Components: Individual + - Name: Lapse + NormType: Min + Components: Individual + - Name: Magnitude(ShiftExcess) + NormType: Max + Components: Individual + - Name: HamiltonianConstraint + NormType: L2Norm + Components: Individual + - Name: MomentumConstraint + NormType: L2Norm + Components: Individual + - ObserveFields: + SubfileName: VolumeData + VariablesToObserve: + - ConformalFactor + - Lapse + - Shift + - ShiftExcess + - SpatialMetric + - InverseSpatialMetric + - SpatialChristoffelSecondKind + - ExtrinsicCurvature + - SpatialRicci + - RadiallyCompressedCoordinates + - HamiltonianConstraint + - MomentumConstraint + BlocksToObserve: All + InterpolateToMesh: None + CoordinatesFloatingPointType: Double + FloatingPointTypes: [Double] + - ObserveAdmIntegrals + + diff --git a/tests/Unit/PointwiseFunctions/AnalyticData/Xcts/Test_BinaryWithGravitationalWaves.cpp b/tests/Unit/PointwiseFunctions/AnalyticData/Xcts/Test_BinaryWithGravitationalWaves.cpp index 8f5190950d30..5545a4566e14 100644 --- a/tests/Unit/PointwiseFunctions/AnalyticData/Xcts/Test_BinaryWithGravitationalWaves.cpp +++ b/tests/Unit/PointwiseFunctions/AnalyticData/Xcts/Test_BinaryWithGravitationalWaves.cpp @@ -41,11 +41,11 @@ template struct BinaryWithGravitationalWavesProxy { tuples::tagged_tuple_from_typelist test_variables( const tnsr::I& x) const { - return BinaryWithGravitationalWaves_ptr->variables(x, test_tags{}); + return binary_with_gravitational_waves->variables(x, test_tags{}); } BinaryWithGravitationalWaves* - BinaryWithGravitationalWaves_ptr; + binary_with_gravitational_waves; }; struct Metavariables { From 0b2ba94fd01caa607eb2b66fd71befc31434d88b Mon Sep 17 00:00:00 2001 From: Megum Rebelo Date: Sun, 9 Feb 2025 21:58:59 +0000 Subject: [PATCH 3/3] Fixup: last two comments in PR --- docs/References.bib | 2 +- .../Xcts/BinaryWithGravitationalWaves.cpp | 5 +---- .../Xcts/BinaryWithGravitationalWaves.hpp | 17 +++++++++-------- 3 files changed, 11 insertions(+), 13 deletions(-) diff --git a/docs/References.bib b/docs/References.bib index f52e36145050..f35ed36e22b9 100644 --- a/docs/References.bib +++ b/docs/References.bib @@ -2555,7 +2555,7 @@ @article{Stamm2010 } @article{Steinhoff2008zr, - author = "Steinhoff, Jan and Schaefer, Gerhard and Hergt, Steven", + author = "{Steinhoff, Jan and Sch\"afer, Gerhard and Hergt, Steven}", title = "{ADM canonical formalism for gravitating spinning objects}", eprint = "0805.3136", archivePrefix = "arXiv", diff --git a/src/PointwiseFunctions/AnalyticData/Xcts/BinaryWithGravitationalWaves.cpp b/src/PointwiseFunctions/AnalyticData/Xcts/BinaryWithGravitationalWaves.cpp index 90ca469cd9bf..47575e58822e 100644 --- a/src/PointwiseFunctions/AnalyticData/Xcts/BinaryWithGravitationalWaves.cpp +++ b/src/PointwiseFunctions/AnalyticData/Xcts/BinaryWithGravitationalWaves.cpp @@ -169,16 +169,13 @@ void BinaryWithGravitationalWavesVariables::operator()( const auto auxiliar_one = tenex::evaluate( inv_conformal_metric(ti::I, ti::K) * inv_conformal_metric(ti::J, ti::L) * dt_conformal_metric(ti::k, ti::l)); - const auto auxiliar_two = tenex::evaluate( - inv_conformal_metric(ti::I, ti::K) * inv_conformal_metric(ti::J, ti::L) * - conformal_metric(ti::k, ti::l)); for (size_t i = 0; i < 3; ++i) { for (size_t j = 0; j <= i; ++j) { longitudinal_shift_background_minus_dt_conformal_metric->get(i, j) -= (auxiliar_one.get(i, j) - (1. / 3.) * get(trace(dt_conformal_metric, inv_conformal_metric)) * - auxiliar_two.get(i, j)); + inv_conformal_metric.get(i, j)); } } } diff --git a/src/PointwiseFunctions/AnalyticData/Xcts/BinaryWithGravitationalWaves.hpp b/src/PointwiseFunctions/AnalyticData/Xcts/BinaryWithGravitationalWaves.hpp index 4aa6959ed99b..dddccb189542 100644 --- a/src/PointwiseFunctions/AnalyticData/Xcts/BinaryWithGravitationalWaves.hpp +++ b/src/PointwiseFunctions/AnalyticData/Xcts/BinaryWithGravitationalWaves.hpp @@ -295,14 +295,15 @@ struct BinaryWithGravitationalWavesVariables * The main goal of this implementation is to improve the extracted * wave forms, for example, by minimizing junk radiation. * The data is only valid for black holes without spin. Even so, there is some - * work done to describe such systems that could later be implemented, \cite - * Steinhoff2008zr. The objects are constructed from a superposition of two - * isolated objects that are boosted with respect to each other. The radiative - * data is constructed from Post-Newtonian expansions for the inspiral phase, in - * orders of \f$\epsilon = 1/c\f$, in \cite Jaranowski1997ky. In ADMTT gauge it - * is possible to get the 3-metric as \f$\gamma^{PN}_{ij} = \psi^{4}_{PN} - * \delta_{ij} + h^{TT}_{ij}\f$ where \f$h^{TT}_{ij}\f$ is the radiative part - * and the non-radiative Post-Newtonian conformal factor is given by + * work done to describe such systems that could later be implemented + * \cite Steinhoff2008zr. The objects are constructed from a superposition of + * two isolated objects that are boosted with respect to each other. The + * radiative data is constructed from Post-Newtonian expansions for the + * inspiral phase, in orders of \f$\epsilon = 1/c\f$, in \cite Jaranowski1997ky. + * In ADMTT gauge it is possible to get the 3-metric as + * \f$\gamma^{PN}_{ij} = \psi^{4}_{PN} \delta_{ij} + h^{TT}_{ij}\f$ where + * \f$h^{TT}_{ij}\f$ is the radiative part and the non-radiative Post-Newtonian + * conformal factor is given by * * \f{equation}{ * \psi_{PN} = 1 + \sum_{a=1}^{2} \frac{E_a}{2 r_a} + O(\epsilon^6)