From 4921ffbbce011b6efb820bf4b16dc137d9c899ba Mon Sep 17 00:00:00 2001 From: Nils Vu Date: Wed, 8 Feb 2023 22:05:02 -0800 Subject: [PATCH 1/5] Add equatorial_radius accessor to RotatingStar --- .../AnalyticSolutions/RelativisticEuler/RotatingStar.hpp | 4 ++++ .../AnalyticSolutions/RelativisticEuler/Test_RotatingStar.cpp | 2 ++ 2 files changed, 6 insertions(+) diff --git a/src/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/RotatingStar.hpp b/src/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/RotatingStar.hpp index 016d57124649..0b119d6ed6f7 100644 --- a/src/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/RotatingStar.hpp +++ b/src/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/RotatingStar.hpp @@ -413,6 +413,10 @@ class RotatingStar : public virtual evolution::initial_data::InitialData, return equation_of_state_; } + double equatorial_radius() const { + return cst_solution_.equatorial_radius(); + } + protected: template using DerivLapse = ::Tags::deriv, tmpl::size_t<3>, diff --git a/tests/Unit/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/Test_RotatingStar.cpp b/tests/Unit/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/Test_RotatingStar.cpp index 2636d971293c..aeaac704d5ab 100644 --- a/tests/Unit/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/Test_RotatingStar.cpp +++ b/tests/Unit/PointwiseFunctions/AnalyticSolutions/RelativisticEuler/Test_RotatingStar.cpp @@ -59,6 +59,8 @@ SPECTRE_TEST_CASE( 100) .get_clone()))); + CHECK(solution.equatorial_radius() == approx(7.155891353887)); + // Near the center the finite difference derivatives cause "large" errors // (1e-8) and so the check fails. verify_solution(solution, {{1., 0., 0.}}, error_tolerance); From 2654cfb0262ca4b09046af18604b0ffaabb83dc1 Mon Sep 17 00:00:00 2001 From: Nils Vu Date: Tue, 7 Feb 2023 22:13:40 -0800 Subject: [PATCH 2/5] Remove an unused template parameter --- src/PointwiseFunctions/AnalyticSolutions/Xcts/TovStar.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/PointwiseFunctions/AnalyticSolutions/Xcts/TovStar.hpp b/src/PointwiseFunctions/AnalyticSolutions/Xcts/TovStar.hpp index e64a49006291..b3bafe637dd6 100644 --- a/src/PointwiseFunctions/AnalyticSolutions/Xcts/TovStar.hpp +++ b/src/PointwiseFunctions/AnalyticSolutions/Xcts/TovStar.hpp @@ -188,7 +188,6 @@ class TovStar : public elliptic::analytic_data::AnalyticSolution { TovStar& operator=(TovStar&&) = default; ~TovStar() = default; - template TovStar(double central_rest_mass_density, std::unique_ptr> equation_of_state, From 98c0d92b4daadd56008c7cdf9d99d866b87dedf4 Mon Sep 17 00:00:00 2001 From: Nils Vu Date: Tue, 7 Feb 2023 22:18:21 -0800 Subject: [PATCH 3/5] Switch Xcts::WrappedGr from inheritance to composition Avoids conflicts with overriding functions from different superclasses. In particular, evolution::InitialData and elliptic::AnalyticSolution define `get_clone` functions with different signatures. --- .../AnalyticSolutions/Xcts/WrappedGr.cpp | 6 ++- .../AnalyticSolutions/Xcts/WrappedGr.hpp | 44 ++++++++++++++----- .../Xcts/Test_HarmonicSchwarzschild.cpp | 4 +- .../AnalyticSolutions/Xcts/Test_Kerr.cpp | 6 +-- .../Xcts/Test_SphericalKerr.cpp | 6 +-- 5 files changed, 45 insertions(+), 21 deletions(-) diff --git a/src/PointwiseFunctions/AnalyticSolutions/Xcts/WrappedGr.cpp b/src/PointwiseFunctions/AnalyticSolutions/Xcts/WrappedGr.cpp index c93b1a4808e2..6fb98b69dff9 100644 --- a/src/PointwiseFunctions/AnalyticSolutions/Xcts/WrappedGr.cpp +++ b/src/PointwiseFunctions/AnalyticSolutions/Xcts/WrappedGr.cpp @@ -299,8 +299,10 @@ template class WrappedGrVariables; } // namespace detail -template -PUP::able::PUP_ID WrappedGr::my_PUP_ID = 0; // NOLINT +template +PUP::able::PUP_ID + WrappedGr>::my_PUP_ID = + 0; // NOLINT } // namespace Xcts::Solutions diff --git a/src/PointwiseFunctions/AnalyticSolutions/Xcts/WrappedGr.hpp b/src/PointwiseFunctions/AnalyticSolutions/Xcts/WrappedGr.hpp index fc5f557b355e..32c09f55cf60 100644 --- a/src/PointwiseFunctions/AnalyticSolutions/Xcts/WrappedGr.hpp +++ b/src/PointwiseFunctions/AnalyticSolutions/Xcts/WrappedGr.hpp @@ -211,9 +211,13 @@ struct WrappedGrVariables * initial guess than flatness, such as a superposition of Kerr solutions for * black-hole binary initial data. */ -template -class WrappedGr : public elliptic::analytic_data::AnalyticSolution, - public GrSolution { +template +class WrappedGr; + +template +class WrappedGr> + : public elliptic::analytic_data::AnalyticSolution { public: static constexpr size_t Dim = 3; @@ -221,21 +225,36 @@ class WrappedGr : public elliptic::analytic_data::AnalyticSolution, static constexpr Options::String help = GrSolution::help; static std::string name() { return pretty_type::name(); } - using GrSolution::GrSolution; + WrappedGr() = default; + WrappedGr(const WrappedGr&) = default; + WrappedGr& operator=(const WrappedGr&) = default; + WrappedGr(WrappedGr&&) = default; + WrappedGr& operator=(WrappedGr&&) = default; + ~WrappedGr() = default; + + WrappedGr(typename GrSolutionOptions::type... gr_solution_options) + : gr_solution_(std::move(gr_solution_options)...) {} + + const GrSolution& gr_solution() const { return gr_solution_; } + + /// \cond + explicit WrappedGr(CkMigrateMessage* m) + : elliptic::analytic_data::AnalyticSolution(m) {} using PUP::able::register_constructor; - WRAPPED_PUPable_decl_template(WrappedGr); + WRAPPED_PUPable_decl_template(WrappedGr); std::unique_ptr get_clone() const override { - return std::make_unique>(*this); + return std::make_unique(*this); } + /// \endcond template tuples::TaggedTuple variables( const tnsr::I& x, tmpl::list /*meta*/) const { const auto gr_solution = - GrSolution::variables(x, std::numeric_limits::signaling_NaN(), - detail::gr_solution_vars{}); + gr_solution_.variables(x, std::numeric_limits::signaling_NaN(), + detail::gr_solution_vars{}); using VarsComputer = detail::WrappedGrVariables; const size_t num_points = get_size(*x.begin()); typename VarsComputer::Cache cache{num_points}; @@ -250,8 +269,8 @@ class WrappedGr : public elliptic::analytic_data::AnalyticSolution, Frame::Inertial>& inv_jacobian, tmpl::list /*meta*/) const { const auto gr_solution = - GrSolution::variables(x, std::numeric_limits::signaling_NaN(), - detail::gr_solution_vars{}); + gr_solution_.variables(x, std::numeric_limits::signaling_NaN(), + detail::gr_solution_vars{}); using VarsComputer = detail::WrappedGrVariables; const size_t num_points = get_size(*x.begin()); typename VarsComputer::Cache cache{num_points}; @@ -260,9 +279,12 @@ class WrappedGr : public elliptic::analytic_data::AnalyticSolution, } void pup(PUP::er& p) override { - GrSolution::pup(p); elliptic::analytic_data::AnalyticSolution::pup(p); + p | gr_solution_; } + + private: + GrSolution gr_solution_; }; } // namespace Xcts::Solutions diff --git a/tests/Unit/PointwiseFunctions/AnalyticSolutions/Xcts/Test_HarmonicSchwarzschild.cpp b/tests/Unit/PointwiseFunctions/AnalyticSolutions/Xcts/Test_HarmonicSchwarzschild.cpp index 2164ad3b1f87..48d1b8e084ae 100644 --- a/tests/Unit/PointwiseFunctions/AnalyticSolutions/Xcts/Test_HarmonicSchwarzschild.cpp +++ b/tests/Unit/PointwiseFunctions/AnalyticSolutions/Xcts/Test_HarmonicSchwarzschild.cpp @@ -31,8 +31,8 @@ void test_solution(const double mass, const std::array& center, const auto& solution = dynamic_cast(*created); { INFO("Properties"); - CHECK(solution.mass() == mass); - CHECK(solution.center() == center); + CHECK(solution.gr_solution().mass() == mass); + CHECK(solution.gr_solution().center() == center); } { INFO("Semantics"); diff --git a/tests/Unit/PointwiseFunctions/AnalyticSolutions/Xcts/Test_Kerr.cpp b/tests/Unit/PointwiseFunctions/AnalyticSolutions/Xcts/Test_Kerr.cpp index 12a8eb7c9684..1e70c4d4e416 100644 --- a/tests/Unit/PointwiseFunctions/AnalyticSolutions/Xcts/Test_Kerr.cpp +++ b/tests/Unit/PointwiseFunctions/AnalyticSolutions/Xcts/Test_Kerr.cpp @@ -33,9 +33,9 @@ void test_solution(const double mass, const std::array spin, const auto& solution = dynamic_cast(*created); { INFO("Properties"); - CHECK(solution.mass() == mass); - CHECK(solution.dimensionless_spin() == spin); - CHECK(solution.center() == center); + CHECK(solution.gr_solution().mass() == mass); + CHECK(solution.gr_solution().dimensionless_spin() == spin); + CHECK(solution.gr_solution().center() == center); } { INFO("Semantics"); diff --git a/tests/Unit/PointwiseFunctions/AnalyticSolutions/Xcts/Test_SphericalKerr.cpp b/tests/Unit/PointwiseFunctions/AnalyticSolutions/Xcts/Test_SphericalKerr.cpp index e7dc14692639..eb0b3d7e3ebc 100644 --- a/tests/Unit/PointwiseFunctions/AnalyticSolutions/Xcts/Test_SphericalKerr.cpp +++ b/tests/Unit/PointwiseFunctions/AnalyticSolutions/Xcts/Test_SphericalKerr.cpp @@ -33,9 +33,9 @@ void test_solution(const double mass, const std::array spin, const auto& solution = dynamic_cast(*created); { INFO("Properties"); - CHECK(solution.mass() == mass); - CHECK(solution.dimensionless_spin() == spin); - CHECK(solution.center() == center); + CHECK(solution.gr_solution().mass() == mass); + CHECK(solution.gr_solution().dimensionless_spin() == spin); + CHECK(solution.gr_solution().center() == center); } { INFO("Semantics"); From 03a940ede7798b4ccb4cbe1fffa65f0a82961e14 Mon Sep 17 00:00:00 2001 From: Nils Vu Date: Wed, 8 Feb 2023 11:38:57 -0800 Subject: [PATCH 4/5] Support AnalyticData in Xcts::WrappedGr --- .../AnalyticSolutions/Xcts/WrappedGr.hpp | 27 ++++++++++++++----- 1 file changed, 21 insertions(+), 6 deletions(-) diff --git a/src/PointwiseFunctions/AnalyticSolutions/Xcts/WrappedGr.hpp b/src/PointwiseFunctions/AnalyticSolutions/Xcts/WrappedGr.hpp index 32c09f55cf60..ef582209c5ba 100644 --- a/src/PointwiseFunctions/AnalyticSolutions/Xcts/WrappedGr.hpp +++ b/src/PointwiseFunctions/AnalyticSolutions/Xcts/WrappedGr.hpp @@ -16,6 +16,7 @@ #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp" #include "Options/Options.hpp" #include "Parallel/CharmPupable.hpp" +#include "PointwiseFunctions/AnalyticSolutions/AnalyticSolution.hpp" #include "PointwiseFunctions/AnalyticSolutions/Xcts/CommonVariables.hpp" #include "PointwiseFunctions/GeneralRelativity/Tags.hpp" #include "PointwiseFunctions/GeneralRelativity/Tags/Conformal.hpp" @@ -252,9 +253,16 @@ class WrappedGr> tuples::TaggedTuple variables( const tnsr::I& x, tmpl::list /*meta*/) const { - const auto gr_solution = - gr_solution_.variables(x, std::numeric_limits::signaling_NaN(), - detail::gr_solution_vars{}); + tuples::tagged_tuple_from_typelist> + gr_solution; + if constexpr (is_analytic_solution_v) { + gr_solution = gr_solution_.variables( + x, std::numeric_limits::signaling_NaN(), + detail::gr_solution_vars{}); + } else { + gr_solution = + gr_solution_.variables(x, detail::gr_solution_vars{}); + } using VarsComputer = detail::WrappedGrVariables; const size_t num_points = get_size(*x.begin()); typename VarsComputer::Cache cache{num_points}; @@ -268,9 +276,16 @@ class WrappedGr> const InverseJacobian& inv_jacobian, tmpl::list /*meta*/) const { - const auto gr_solution = - gr_solution_.variables(x, std::numeric_limits::signaling_NaN(), - detail::gr_solution_vars{}); + tuples::tagged_tuple_from_typelist> + gr_solution; + if constexpr (is_analytic_solution_v) { + gr_solution = gr_solution_.variables( + x, std::numeric_limits::signaling_NaN(), + detail::gr_solution_vars{}); + } else { + gr_solution = + gr_solution_.variables(x, detail::gr_solution_vars{}); + } using VarsComputer = detail::WrappedGrVariables; const size_t num_points = get_size(*x.begin()); typename VarsComputer::Cache cache{num_points}; From 08a28759446b5e13ce42f92b5b9767d832f9f0db Mon Sep 17 00:00:00 2001 From: Nils Vu Date: Wed, 8 Feb 2023 22:16:52 -0800 Subject: [PATCH 5/5] Extend Xcts::WrappedGr to GRMHD Allows to solve the Einstein constraints for any fixed GRMHD profile, such as a magnetized TOV star or a CCSN profile. Not that this does _not_ solve for the hydro sector, only for the gravity sector given the matter content. --- docs/References.bib | 17 ++ .../AnalyticSolutions/Xcts/CMakeLists.txt | 1 + .../AnalyticSolutions/Xcts/Factory.hpp | 14 +- .../AnalyticSolutions/Xcts/WrappedGr.cpp | 275 ++++++++++++++---- .../AnalyticSolutions/Xcts/WrappedGr.hpp | 111 ++++++- .../AnalyticSolutions/Xcts/CMakeLists.txt | 2 + .../Xcts/Test_RotatingStar.cpp | 58 ++++ 7 files changed, 403 insertions(+), 75 deletions(-) create mode 100644 tests/Unit/PointwiseFunctions/AnalyticSolutions/Xcts/Test_RotatingStar.cpp diff --git a/docs/References.bib b/docs/References.bib index 435806b3ee3b..ee5c5905c475 100644 --- a/docs/References.bib +++ b/docs/References.bib @@ -2062,6 +2062,23 @@ @article{Szilagyi2014fna year = "2014" } +@article{Tacik2016zal, + author = {Tacik, Nick and Foucart, Francois and Pfeiffer, Harald P. + and Muhlberger, Curran and Kidder, Lawrence E. and Scheel, + Mark A. and Szil\'agyi, B.}, + title = {Initial data for black hole-neutron star binaries, with + rotating stars}, + eprint = {1607.07962}, + archiveprefix = {arXiv}, + primaryclass = {gr-qc}, + doi = {10.1088/0264-9381/33/22/225012}, + journal = {Class. Quant. Grav.}, + volume = {33}, + number = {22}, + pages = {225012}, + year = {2016} +} + @article{Teukolsky2015ega, author = "Teukolsky, Saul A.", title = "Formulation of discontinuous {Galerkin} methods for diff --git a/src/PointwiseFunctions/AnalyticSolutions/Xcts/CMakeLists.txt b/src/PointwiseFunctions/AnalyticSolutions/Xcts/CMakeLists.txt index ea1fe5d73914..38dadbd4925f 100644 --- a/src/PointwiseFunctions/AnalyticSolutions/Xcts/CMakeLists.txt +++ b/src/PointwiseFunctions/AnalyticSolutions/Xcts/CMakeLists.txt @@ -37,6 +37,7 @@ target_link_libraries( ErrorHandling GeneralRelativity GeneralRelativitySolutions + GrMhdAnalyticData InitialDataUtilities Options Parallel diff --git a/src/PointwiseFunctions/AnalyticSolutions/Xcts/Factory.hpp b/src/PointwiseFunctions/AnalyticSolutions/Xcts/Factory.hpp index 3889744c399f..5cf1276fa219 100644 --- a/src/PointwiseFunctions/AnalyticSolutions/Xcts/Factory.hpp +++ b/src/PointwiseFunctions/AnalyticSolutions/Xcts/Factory.hpp @@ -3,9 +3,12 @@ #pragma once +#include "PointwiseFunctions/AnalyticData/GrMhd/CcsnCollapse.hpp" +#include "PointwiseFunctions/AnalyticData/GrMhd/MagnetizedTovStar.hpp" #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/HarmonicSchwarzschild.hpp" #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrSchild.hpp" #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/SphericalKerrSchild.hpp" +#include "PointwiseFunctions/AnalyticSolutions/RelativisticEuler/RotatingStar.hpp" #include "PointwiseFunctions/AnalyticSolutions/Xcts/Flatness.hpp" #include "PointwiseFunctions/AnalyticSolutions/Xcts/Schwarzschild.hpp" #include "PointwiseFunctions/AnalyticSolutions/Xcts/TovStar.hpp" @@ -18,6 +21,15 @@ namespace Solutions { using all_analytic_solutions = tmpl::list, WrappedGr, Schwarzschild, - WrappedGr, TovStar>; + WrappedGr, TovStar, + WrappedGrMhd, + // The following are only approximate solutions to the XCTS + // equations. We list them here (as opposed to AnalyticData) to + // avoid a cyclic dependency with the `WrappedGr` class, and also + // for convenience: by treating them as analytic solutions we can + // easily measure the difference between the approximate + // (analytic) solution and the numerical solution. + WrappedGrMhd, + WrappedGrMhd>; } // namespace Solutions } // namespace Xcts diff --git a/src/PointwiseFunctions/AnalyticSolutions/Xcts/WrappedGr.cpp b/src/PointwiseFunctions/AnalyticSolutions/Xcts/WrappedGr.cpp index 6fb98b69dff9..7ad135a363a8 100644 --- a/src/PointwiseFunctions/AnalyticSolutions/Xcts/WrappedGr.cpp +++ b/src/PointwiseFunctions/AnalyticSolutions/Xcts/WrappedGr.cpp @@ -12,14 +12,19 @@ #include "Elliptic/Systems/Xcts/Tags.hpp" #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp" #include "Options/Options.hpp" +#include "PointwiseFunctions/AnalyticData/GrMhd/CcsnCollapse.hpp" +#include "PointwiseFunctions/AnalyticData/GrMhd/MagnetizedTovStar.hpp" #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/HarmonicSchwarzschild.hpp" #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/KerrSchild.hpp" #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/SphericalKerrSchild.hpp" +#include "PointwiseFunctions/AnalyticSolutions/RelativisticEuler/RotatingStar.hpp" #include "PointwiseFunctions/AnalyticSolutions/Xcts/CommonVariables.tpp" #include "PointwiseFunctions/Elasticity/Strain.hpp" #include "PointwiseFunctions/GeneralRelativity/IndexManipulation.hpp" #include "PointwiseFunctions/GeneralRelativity/Tags.hpp" #include "PointwiseFunctions/GeneralRelativity/Tags/Conformal.hpp" +#include "PointwiseFunctions/Hydro/ComovingMagneticField.hpp" +#include "PointwiseFunctions/Hydro/StressEnergy.hpp" #include "PointwiseFunctions/Xcts/LongitudinalOperator.hpp" #include "Utilities/ConstantExpressions.hpp" #include "Utilities/GenerateInstantiations.hpp" @@ -28,8 +33,8 @@ namespace Xcts::Solutions { namespace detail { -template -void WrappedGrVariables::operator()( +template +void WrappedGrVariables::operator()( const gsl::not_null*> spatial_metric, const gsl::not_null /*cache*/, gr::Tags::SpatialMetric /*meta*/) const { @@ -37,8 +42,8 @@ void WrappedGrVariables::operator()( get>(gr_solution); } -template -void WrappedGrVariables::operator()( +template +void WrappedGrVariables::operator()( const gsl::not_null*> inv_spatial_metric, const gsl::not_null /*cache*/, gr::Tags::InverseSpatialMetric /*meta*/) @@ -48,8 +53,8 @@ void WrappedGrVariables::operator()( gr_solution); } -template -void WrappedGrVariables::operator()( +template +void WrappedGrVariables::operator()( const gsl::not_null*> deriv_spatial_metric, const gsl::not_null /*cache*/, ::Tags::deriv, @@ -59,8 +64,8 @@ void WrappedGrVariables::operator()( tmpl::size_t, Frame::Inertial>>(gr_solution); } -template -void WrappedGrVariables::operator()( +template +void WrappedGrVariables::operator()( const gsl::not_null*> conformal_metric, const gsl::not_null cache, Xcts::Tags::ConformalMetric /*meta*/) @@ -76,8 +81,8 @@ void WrappedGrVariables::operator()( } } -template -void WrappedGrVariables::operator()( +template +void WrappedGrVariables::operator()( const gsl::not_null*> inv_conformal_metric, const gsl::not_null cache, Xcts::Tags::InverseConformalMetric /*meta*/) @@ -94,8 +99,8 @@ void WrappedGrVariables::operator()( } } -template -void WrappedGrVariables::operator()( +template +void WrappedGrVariables::operator()( const gsl::not_null*> deriv_conformal_metric, const gsl::not_null cache, ::Tags::deriv, @@ -122,8 +127,8 @@ void WrappedGrVariables::operator()( } } -template -void WrappedGrVariables::operator()( +template +void WrappedGrVariables::operator()( const gsl::not_null*> trace_extrinsic_curvature, const gsl::not_null /*cache*/, gr::Tags::TraceExtrinsicCurvature /*meta*/) const { @@ -136,24 +141,24 @@ void WrappedGrVariables::operator()( trace(trace_extrinsic_curvature, extrinsic_curvature, inv_spatial_metric); } -template -void WrappedGrVariables::operator()( +template +void WrappedGrVariables::operator()( const gsl::not_null*> dt_trace_extrinsic_curvature, const gsl::not_null /*cache*/, ::Tags::dt> /*meta*/) const { get(*dt_trace_extrinsic_curvature) = 0.; } -template -void WrappedGrVariables::operator()( +template +void WrappedGrVariables::operator()( const gsl::not_null*> conformal_factor, const gsl::not_null /*cache*/, Xcts::Tags::ConformalFactor /*meta*/) const { get(*conformal_factor) = 1.; } -template -void WrappedGrVariables::operator()( +template +void WrappedGrVariables::operator()( const gsl::not_null*> conformal_factor_gradient, const gsl::not_null /*cache*/, ::Tags::deriv, tmpl::size_t, @@ -162,16 +167,16 @@ void WrappedGrVariables::operator()( conformal_factor_gradient->end(), 0.); } -template -void WrappedGrVariables::operator()( +template +void WrappedGrVariables::operator()( const gsl::not_null*> lapse, const gsl::not_null /*cache*/, gr::Tags::Lapse /*meta*/) const { *lapse = get>(gr_solution); } -template -void WrappedGrVariables::operator()( +template +void WrappedGrVariables::operator()( const gsl::not_null*> lapse_times_conformal_factor, const gsl::not_null cache, Xcts::Tags::LapseTimesConformalFactor /*meta*/) const { @@ -182,8 +187,8 @@ void WrappedGrVariables::operator()( get(*lapse_times_conformal_factor) *= get(conformal_factor); } -template -void WrappedGrVariables::operator()( +template +void WrappedGrVariables::operator()( const gsl::not_null*> lapse_times_conformal_factor_gradient, const gsl::not_null cache, @@ -205,8 +210,8 @@ void WrappedGrVariables::operator()( } } -template -void WrappedGrVariables::operator()( +template +void WrappedGrVariables::operator()( const gsl::not_null*> shift_background, const gsl::not_null /*cache*/, Xcts::Tags::ShiftBackground /*meta*/) @@ -214,8 +219,8 @@ void WrappedGrVariables::operator()( std::fill(shift_background->begin(), shift_background->end(), 0.); } -template -void WrappedGrVariables::operator()( +template +void WrappedGrVariables::operator()( const gsl::not_null*> longitudinal_shift_background_minus_dt_conformal_metric, const gsl::not_null /*cache*/, @@ -225,8 +230,8 @@ void WrappedGrVariables::operator()( longitudinal_shift_background_minus_dt_conformal_metric->end(), 0.); } -template -void WrappedGrVariables::operator()( +template +void WrappedGrVariables::operator()( const gsl::not_null*> shift_excess, const gsl::not_null /*cache*/, Xcts::Tags::ShiftExcess /*meta*/) const { @@ -234,8 +239,8 @@ void WrappedGrVariables::operator()( get>(gr_solution); } -template -void WrappedGrVariables::operator()( +template +void WrappedGrVariables::operator()( const gsl::not_null*> shift_strain, const gsl::not_null cache, Xcts::Tags::ShiftStrain /*meta*/) const { @@ -258,8 +263,8 @@ void WrappedGrVariables::operator()( shift); } -template -void WrappedGrVariables::operator()( +template +void WrappedGrVariables::operator()( const gsl::not_null*> extrinsic_curvature, const gsl::not_null /*cache*/, gr::Tags::ExtrinsicCurvature<3, Frame::Inertial, DataType> /*meta*/) const { @@ -268,40 +273,175 @@ void WrappedGrVariables::operator()( gr_solution); } -template -void WrappedGrVariables::operator()( - const gsl::not_null*> energy_density, +#if defined(__GNUC__) && !defined(__clang__) +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wsuggest-attribute=noreturn" +#endif // defined(__GNUC__) && !defined(__clang__) + +template +void WrappedGrVariables::operator()( + [[maybe_unused]] const gsl::not_null*> + magnetic_field_dot_spatial_velocity, const gsl::not_null /*cache*/, + hydro::Tags::MagneticFieldDotSpatialVelocity /*meta*/) const { + if constexpr (HasMhd) { + const auto& spatial_velocity = + get>( + hydro_solution); + const auto& magnetic_field = + get>( + hydro_solution); + const auto& spatial_metric = + get>( + gr_solution); + tenex::evaluate(magnetic_field_dot_spatial_velocity, + magnetic_field(ti::I) * spatial_velocity(ti::J) * + spatial_metric(ti::i, ti::j)); + } else { + ERROR( + "The 'MagneticFieldDotSpatialVelocity' should not be needed in vacuum " + "GR."); + } +} + +template +void WrappedGrVariables::operator()( + [[maybe_unused]] const gsl::not_null*> + comoving_magnetic_field_squared, + [[maybe_unused]] const gsl::not_null cache, + hydro::Tags::ComovingMagneticFieldSquared /*meta*/) const { + if constexpr (HasMhd) { + const auto& lorentz_factor = + get>(hydro_solution); + const auto& magnetic_field = + get>( + hydro_solution); + const auto& spatial_metric = + get>( + gr_solution); + const auto magnetic_field_squared = + tenex::evaluate(magnetic_field(ti::I) * magnetic_field(ti::J) * + spatial_metric(ti::i, ti::j)); + const auto& magnetic_field_dot_spatial_velocity = cache->get_var( + *this, hydro::Tags::MagneticFieldDotSpatialVelocity{}); + hydro::comoving_magnetic_field_squared( + comoving_magnetic_field_squared, magnetic_field_squared, + magnetic_field_dot_spatial_velocity, lorentz_factor); + } else { + ERROR( + "The 'ComovingMagneticFieldSquared' should not be needed in vacuum " + "GR."); + } +} + +#if defined(__GNUC__) && !defined(__clang__) +#pragma GCC diagnostic pop +#endif // defined(__GNUC__) && !defined(__clang__) + +template +void WrappedGrVariables::operator()( + const gsl::not_null*> energy_density, + [[maybe_unused]] const gsl::not_null cache, gr::Tags::Conformal, 0> /*meta*/) const { - std::fill(energy_density->begin(), energy_density->end(), 0.); + if constexpr (HasMhd) { + const auto& rest_mass_density = + get>(hydro_solution); + const auto& specific_enthalpy = + get>(hydro_solution); + const auto& pressure = get>(hydro_solution); + const auto& lorentz_factor = + get>(hydro_solution); + const auto& magnetic_field_dot_spatial_velocity = cache->get_var( + *this, hydro::Tags::MagneticFieldDotSpatialVelocity{}); + const auto& comoving_magnetic_field_squared = cache->get_var( + *this, hydro::Tags::ComovingMagneticFieldSquared{}); + hydro::energy_density(energy_density, rest_mass_density, specific_enthalpy, + pressure, lorentz_factor, + magnetic_field_dot_spatial_velocity, + comoving_magnetic_field_squared); + } else { + std::fill(energy_density->begin(), energy_density->end(), 0.); + } } -template -void WrappedGrVariables::operator()( +template +void WrappedGrVariables::operator()( const gsl::not_null*> stress_trace, - const gsl::not_null /*cache*/, + [[maybe_unused]] const gsl::not_null cache, gr::Tags::Conformal, 0> /*meta*/) const { - std::fill(stress_trace->begin(), stress_trace->end(), 0.); + if constexpr (HasMhd) { + const auto& rest_mass_density = + get>(hydro_solution); + const auto& specific_enthalpy = + get>(hydro_solution); + const auto& pressure = get>(hydro_solution); + const auto& spatial_velocity = + get>( + hydro_solution); + const auto& spatial_metric = + get>( + gr_solution); + const auto spatial_velocity_squared = + tenex::evaluate(spatial_velocity(ti::I) * spatial_velocity(ti::J) * + spatial_metric(ti::i, ti::j)); + const auto& lorentz_factor = + get>(hydro_solution); + const auto& magnetic_field_dot_spatial_velocity = cache->get_var( + *this, hydro::Tags::MagneticFieldDotSpatialVelocity{}); + const auto& comoving_magnetic_field_squared = cache->get_var( + *this, hydro::Tags::ComovingMagneticFieldSquared{}); + hydro::stress_trace(stress_trace, rest_mass_density, specific_enthalpy, + pressure, spatial_velocity_squared, lorentz_factor, + magnetic_field_dot_spatial_velocity, + comoving_magnetic_field_squared); + } else { + std::fill(stress_trace->begin(), stress_trace->end(), 0.); + } } -template -void WrappedGrVariables::operator()( +template +void WrappedGrVariables::operator()( const gsl::not_null*> momentum_density, - const gsl::not_null /*cache*/, + [[maybe_unused]] const gsl::not_null cache, gr::Tags::Conformal< gr::Tags::MomentumDensity, 0> /*meta*/) const { - std::fill(momentum_density->begin(), momentum_density->end(), 0.); + if constexpr (HasMhd) { + const auto& rest_mass_density = + get>(hydro_solution); + const auto& specific_enthalpy = + get>(hydro_solution); + const auto& spatial_velocity = + get>( + hydro_solution); + const auto& lorentz_factor = + get>(hydro_solution); + const auto& magnetic_field = + get>( + hydro_solution); + const auto& magnetic_field_dot_spatial_velocity = cache->get_var( + *this, hydro::Tags::MagneticFieldDotSpatialVelocity{}); + const auto& comoving_magnetic_field_squared = cache->get_var( + *this, hydro::Tags::ComovingMagneticFieldSquared{}); + hydro::momentum_density(momentum_density, rest_mass_density, + specific_enthalpy, spatial_velocity, lorentz_factor, + magnetic_field, magnetic_field_dot_spatial_velocity, + comoving_magnetic_field_squared); + } else { + std::fill(momentum_density->begin(), momentum_density->end(), 0.); + } } -template class WrappedGrVariables; -template class WrappedGrVariables; +template class WrappedGrVariables; +template class WrappedGrVariables; +template class WrappedGrVariables; +template class WrappedGrVariables; } // namespace detail -template +template PUP::able::PUP_ID - WrappedGr>::my_PUP_ID = + WrappedGr>::my_PUP_ID = 0; // NOLINT } // namespace Xcts::Solutions @@ -309,25 +449,32 @@ PUP::able::PUP_ID // Instantiate implementations for common variables template class Xcts::Solutions::CommonVariables< double, - typename Xcts::Solutions::detail::WrappedGrVariables::Cache>; + typename Xcts::Solutions::detail::WrappedGrVariables::Cache>; template class Xcts::Solutions::CommonVariables< - DataVector, - typename Xcts::Solutions::detail::WrappedGrVariables::Cache>; + DataVector, typename Xcts::Solutions::detail::WrappedGrVariables< + DataVector, false>::Cache>; template class Xcts::AnalyticData::CommonVariables< double, - typename Xcts::Solutions::detail::WrappedGrVariables::Cache>; + typename Xcts::Solutions::detail::WrappedGrVariables::Cache>; template class Xcts::AnalyticData::CommonVariables< - DataVector, - typename Xcts::Solutions::detail::WrappedGrVariables::Cache>; + DataVector, typename Xcts::Solutions::detail::WrappedGrVariables< + DataVector, false>::Cache>; #define STYPE(data) BOOST_PP_TUPLE_ELEM(0, data) -#define INSTANTIATE(_, data) \ - template class Xcts::Solutions::WrappedGr; +#define INSTANTIATE_GR(_, data) \ + template class Xcts::Solutions::WrappedGr; +#define INSTANTIATE_GRMHD(_, data) \ + template class Xcts::Solutions::WrappedGr; -GENERATE_INSTANTIATIONS(INSTANTIATE, (gr::Solutions::KerrSchild, - gr::Solutions::SphericalKerrSchild, - gr::Solutions::HarmonicSchwarzschild)) +GENERATE_INSTANTIATIONS(INSTANTIATE_GR, (gr::Solutions::KerrSchild, + gr::Solutions::SphericalKerrSchild, + gr::Solutions::HarmonicSchwarzschild)) +GENERATE_INSTANTIATIONS(INSTANTIATE_GRMHD, + (grmhd::AnalyticData::CcsnCollapse, + grmhd::AnalyticData::MagnetizedTovStar, + RelativisticEuler::Solutions::RotatingStar)) #undef STYPE -#undef INSTANTIATE +#undef INSTANTIATE_GR +#undef INSTANTIATE_GRMHD diff --git a/src/PointwiseFunctions/AnalyticSolutions/Xcts/WrappedGr.hpp b/src/PointwiseFunctions/AnalyticSolutions/Xcts/WrappedGr.hpp index ef582209c5ba..1f57bf875f6c 100644 --- a/src/PointwiseFunctions/AnalyticSolutions/Xcts/WrappedGr.hpp +++ b/src/PointwiseFunctions/AnalyticSolutions/Xcts/WrappedGr.hpp @@ -20,6 +20,7 @@ #include "PointwiseFunctions/AnalyticSolutions/Xcts/CommonVariables.hpp" #include "PointwiseFunctions/GeneralRelativity/Tags.hpp" #include "PointwiseFunctions/GeneralRelativity/Tags/Conformal.hpp" +#include "PointwiseFunctions/Hydro/Tags.hpp" #include "PointwiseFunctions/InitialDataUtilities/AnalyticSolution.hpp" #include "Utilities/Gsl.hpp" #include "Utilities/PrettyType.hpp" @@ -43,16 +44,27 @@ using gr_solution_vars = tmpl::list< tmpl::size_t, Frame::Inertial>, gr::Tags::ExtrinsicCurvature>; +template +using hydro_solution_vars = + tmpl::list, + hydro::Tags::SpecificEnthalpy, + hydro::Tags::Pressure, + hydro::Tags::SpatialVelocity, + hydro::Tags::LorentzFactor, + hydro::Tags::MagneticField>; + template using WrappedGrVariablesCache = cached_temp_buffer_from_typelist, + hydro::Tags::MagneticFieldDotSpatialVelocity, + hydro::Tags::ComovingMagneticFieldSquared, gr::Tags::Conformal, 0>, gr::Tags::Conformal, 0>, gr::Tags::Conformal< gr::Tags::MomentumDensity<3, Frame::Inertial, DataType>, 0>>>; -template +template struct WrappedGrVariables : CommonVariables> { static constexpr size_t Dim = 3; @@ -67,14 +79,19 @@ struct WrappedGrVariables local_inv_jacobian, const tnsr::I& local_x, const tuples::tagged_tuple_from_typelist>& - local_gr_solution) + local_gr_solution, + const tuples::tagged_tuple_from_typelist< + hydro_solution_vars>& local_hydro_solution) : Base(std::move(local_mesh), std::move(local_inv_jacobian)), x(local_x), - gr_solution(local_gr_solution) {} + gr_solution(local_gr_solution), + hydro_solution(local_hydro_solution) {} const tnsr::I& x; const tuples::tagged_tuple_from_typelist>& gr_solution; + const tuples::tagged_tuple_from_typelist>& + hydro_solution; void operator()( gsl::not_null*> conformal_metric, @@ -162,6 +179,14 @@ struct WrappedGrVariables gsl::not_null cache, gr::Tags::ExtrinsicCurvature<3, Frame::Inertial, DataType> /*meta*/) const override; + void operator()( + gsl::not_null*> magnetic_field_dot_spatial_velocity, + gsl::not_null cache, + hydro::Tags::MagneticFieldDotSpatialVelocity /*meta*/) const; + void operator()( + gsl::not_null*> comoving_magnetic_field_squared, + gsl::not_null cache, + hydro::Tags::ComovingMagneticFieldSquared /*meta*/) const; void operator()( gsl::not_null*> energy_density, gsl::not_null cache, @@ -211,13 +236,38 @@ struct WrappedGrVariables * production solves this is not an issue because we choose a much better * initial guess than flatness, such as a superposition of Kerr solutions for * black-hole binary initial data. + * + * \warning + * The computation of the XCTS matter source terms (energy density $\rho$, + * momentum density $S^i$, stress trace $S$) uses GR quantities (lapse $\alpha$, + * shift $\beta^i$, spatial metric $\gamma_{ij}$), which means these GR + * quantities are not treated dynamically in the source terms when solving the + * XCTS equations. If the GR quantities satisfy the Einstein constraints (as is + * the case if the `GrSolution` is actually a solution to the Einstein + * equations), then the XCTS solve will reproduce the GR quantities given the + * fixed sources computed here. However, if the GR quantities don't satisfy the + * Einstein constraints (e.g. because a magnetic field was added to the solution + * but ignored in the gravity sector, or because it is a hydrodynamic solution + * on a fixed background metric) then the XCTS solution will depend on our + * treatment of the source terms: fixing the source terms (the simple approach + * taken here) means we're making a choice of $W$ and $u^i$. This is what + * initial data codes usually do when they iterate back and forth between a + * hydro solve and an XCTS solve (e.g. see \cite Tacik2016zal). Alternatively, + * we could fix $v^i$ and compute $W$ and $u^i$ from $v^i$ and the dynamic + * metric variables at every step in the XCTS solver algorithm. This requires + * adding the source terms and their linearization to the XCTS equations, and + * could be interesting to explore. + * + * \tparam GrSolution Any solution to the Einstein constraint equations + * \tparam HasMhd Enable to compute matter source terms. Disable to set matter + * source terms to zero. */ -template class WrappedGr; -template -class WrappedGr> +template +class WrappedGr> : public elliptic::analytic_data::AnalyticSolution { public: static constexpr size_t Dim = 3; @@ -263,10 +313,24 @@ class WrappedGr> gr_solution = gr_solution_.variables(x, detail::gr_solution_vars{}); } - using VarsComputer = detail::WrappedGrVariables; + tuples::tagged_tuple_from_typelist< + detail::hydro_solution_vars> + hydro_solution; + if constexpr (HasMhd) { + if constexpr (is_analytic_solution_v) { + hydro_solution = gr_solution_.variables( + x, std::numeric_limits::signaling_NaN(), + detail::hydro_solution_vars{}); + } else { + hydro_solution = gr_solution_.variables( + x, detail::hydro_solution_vars{}); + } + } + using VarsComputer = detail::WrappedGrVariables; const size_t num_points = get_size(*x.begin()); typename VarsComputer::Cache cache{num_points}; - const VarsComputer computer{std::nullopt, std::nullopt, x, gr_solution}; + const VarsComputer computer{std::nullopt, std::nullopt, x, gr_solution, + hydro_solution}; return {cache.get_var(computer, RequestedTags{})...}; } @@ -286,10 +350,23 @@ class WrappedGr> gr_solution = gr_solution_.variables(x, detail::gr_solution_vars{}); } - using VarsComputer = detail::WrappedGrVariables; + tuples::tagged_tuple_from_typelist< + detail::hydro_solution_vars> + hydro_solution; + if constexpr (HasMhd) { + if constexpr (is_analytic_solution_v) { + hydro_solution = gr_solution_.variables( + x, std::numeric_limits::signaling_NaN(), + detail::hydro_solution_vars{}); + } else { + hydro_solution = gr_solution_.variables( + x, detail::hydro_solution_vars{}); + } + } + using VarsComputer = detail::WrappedGrVariables; const size_t num_points = get_size(*x.begin()); typename VarsComputer::Cache cache{num_points}; - VarsComputer computer{mesh, inv_jacobian, x, gr_solution}; + VarsComputer computer{mesh, inv_jacobian, x, gr_solution, hydro_solution}; return {cache.get_var(computer, RequestedTags{})...}; } @@ -299,7 +376,21 @@ class WrappedGr> } private: + friend bool operator==(const WrappedGr& lhs, + const WrappedGr& rhs) { + return lhs.gr_solution_ == rhs.gr_solution_; + } + GrSolution gr_solution_; }; +template +inline bool operator!=(const WrappedGr& lhs, + const WrappedGr& rhs) { + return not(lhs == rhs); +} + +template +using WrappedGrMhd = WrappedGr; + } // namespace Xcts::Solutions diff --git a/tests/Unit/PointwiseFunctions/AnalyticSolutions/Xcts/CMakeLists.txt b/tests/Unit/PointwiseFunctions/AnalyticSolutions/Xcts/CMakeLists.txt index 984545c2ce40..b8220048319d 100644 --- a/tests/Unit/PointwiseFunctions/AnalyticSolutions/Xcts/CMakeLists.txt +++ b/tests/Unit/PointwiseFunctions/AnalyticSolutions/Xcts/CMakeLists.txt @@ -8,6 +8,7 @@ set(LIBRARY_SOURCES Test_Flatness.cpp Test_HarmonicSchwarzschild.cpp Test_Kerr.cpp + Test_RotatingStar.cpp Test_Schwarzschild.cpp Test_SphericalKerr.cpp Test_TovStar.cpp @@ -26,6 +27,7 @@ target_link_libraries( CoordinateMaps DataStructures Domain + RelativisticEulerSolutions Spectral Utilities Xcts diff --git a/tests/Unit/PointwiseFunctions/AnalyticSolutions/Xcts/Test_RotatingStar.cpp b/tests/Unit/PointwiseFunctions/AnalyticSolutions/Xcts/Test_RotatingStar.cpp new file mode 100644 index 000000000000..ac4fc69d0246 --- /dev/null +++ b/tests/Unit/PointwiseFunctions/AnalyticSolutions/Xcts/Test_RotatingStar.cpp @@ -0,0 +1,58 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Framework/TestingFramework.hpp" + +#include + +#include "Elliptic/Systems/Xcts/Geometry.hpp" +#include "Framework/TestCreation.hpp" +#include "Framework/TestHelpers.hpp" +#include "Helpers/PointwiseFunctions/AnalyticSolutions/Xcts/VerifySolution.hpp" +#include "Informer/InfoFromBuild.hpp" +#include "PointwiseFunctions/AnalyticSolutions/RelativisticEuler/RotatingStar.hpp" +#include "PointwiseFunctions/AnalyticSolutions/Xcts/WrappedGr.hpp" + +namespace Xcts::Solutions { + +// [[Timeout, 20]] +SPECTRE_TEST_CASE("Unit.PointwiseFunctions.AnalyticSolutions.Xcts.RotatingStar", + "[PointwiseFunctions][Unit]") { + using RotatingStar = WrappedGrMhd; + const std::string rotns_filename = unit_test_src_path() + + "/PointwiseFunctions/AnalyticSolutions/" + "RelativisticEuler/RotatingStarId.dat"; + const auto created = TestHelpers::test_factory_creation< + elliptic::analytic_data::AnalyticSolution, RotatingStar>( + "RotatingStar:\n" + " RotNsFilename: " + + rotns_filename + + "\n" + " PolytropicConstant: 100.\n"); + REQUIRE(dynamic_cast(created.get()) != nullptr); + const auto& solution = dynamic_cast(*created); + { + INFO("Semantics"); + test_serialization(solution); + test_copy_semantics(solution); + auto move_solution = solution; + test_move_semantics(std::move(move_solution), solution); + } + { + INFO("Verify the solution solves the XCTS equations"); + const auto& equatorial_radius = solution.gr_solution().equatorial_radius(); + CAPTURE(equatorial_radius); + for (const double fraction_of_star_radius : {0.2, 1., 1.5}) { + CAPTURE(fraction_of_star_radius); + const double inner_radius = fraction_of_star_radius * equatorial_radius; + const double outer_radius = + (fraction_of_star_radius + 0.01) * equatorial_radius; + CAPTURE(inner_radius); + CAPTURE(outer_radius); + TestHelpers::Xcts::Solutions::verify_solution( + solution, {{0., 0., 0.}}, inner_radius, outer_radius, 2.e-3); + } + } +} + +} // namespace Xcts::Solutions