diff --git a/docs/References.bib b/docs/References.bib index aa24a3d0c32d..660f9140838b 100644 --- a/docs/References.bib +++ b/docs/References.bib @@ -182,6 +182,20 @@ @article{Baumgarte2006ug SLACcitation = "%%CITATION = GR-QC/0610120;%%" } +@article{Baumgarte2007ht, + author = "Baumgarte, Thomas W. and Naculich, Stephen G.", + title = "{Analytical representation of a black hole puncture solution}", + eprint = "gr-qc/0701037", + archivePrefix = "arXiv", + reportNumber = "BOW-PH-138", + doi = "10.1103/PhysRevD.75.067502", + url = {https://doi.org/10.1103/PhysRevD.75.067502}, + journal = "Phys. Rev. D", + volume = "75", + pages = "067502", + year = "2007" +} + @book{BaumgarteShapiro, author = {Baumgarte, Thomas W. and Shapiro, Stuart L.}, title = {Numerical Relativity: Solving {Einstein's} Equations on the @@ -885,6 +899,19 @@ @article{Dumbser2017okk year = "2018" } +@article{Estabrook1973ue, + author = "{Estabrook, Frank and Wahlquist, Hugo and Christensen, Steven and + DeWitt, Bryce and Smarr, Larry and Tsiang, Elaine}", + title = "{Maximally slicing a black hole}", + doi = "10.1103/PhysRevD.7.2814", + url = "https://doi.org/10.1103/PhysRevD.7.2814", + journal = "Phys. Rev. D", + volume = "7", + pages = "2814--2817", + year = "1973" +} + + @article{Etienne2008re, author = {Etienne, Zachariah B. and Liu, Yuk Tung and Shapiro, Stuart L. and Baumgarte, Thomas W.}, diff --git a/src/PointwiseFunctions/AnalyticSolutions/Xcts/Schwarzschild.cpp b/src/PointwiseFunctions/AnalyticSolutions/Xcts/Schwarzschild.cpp index db708f884df0..6ffaba27fd01 100644 --- a/src/PointwiseFunctions/AnalyticSolutions/Xcts/Schwarzschild.cpp +++ b/src/PointwiseFunctions/AnalyticSolutions/Xcts/Schwarzschild.cpp @@ -42,6 +42,8 @@ std::ostream& operator<<(std::ostream& os, return os << "PainleveGullstrand"; case SchwarzschildCoordinates::KerrSchildIsotropic: return os << "KerrSchildIsotropic"; + case SchwarzschildCoordinates::MaximalIsotropic: + return os << "MaximalIsotropic"; // LCOV_EXCL_START default: ERROR("Unknown SchwarzschildCoordinates"); @@ -62,13 +64,15 @@ Options::create_from_yaml::create< return Xcts::Solutions::SchwarzschildCoordinates::PainleveGullstrand; } else if ("KerrSchildIsotropic" == type_read) { return Xcts::Solutions::SchwarzschildCoordinates::KerrSchildIsotropic; + } else if ("MaximalIsotropic" == type_read) { + return Xcts::Solutions::SchwarzschildCoordinates::MaximalIsotropic; } PARSE_ERROR(options.context(), "Failed to convert \"" << type_read << "\" to Xcts::Solutions::SchwarzschildCoordinates. Must be " "one of 'Isotropic', 'PainleveGullstrand', " - "'KerrSchildIsotropic'."); + "'KerrSchildIsotropic', 'MaximalIsotropic'."); } namespace Xcts::Solutions { @@ -118,6 +122,60 @@ DataVector kerr_schild_areal_radius_from_isotropic( }, isotropic_radius, isotropic_radius + mass, 1.0e-12, 1.0e-15); } + +template +DataType maximal_isotropic_radius_from_areal(const DataType& areal_radius, + const double mass) { + return .25 * + (2. * areal_radius + mass + + sqrt(4. * square(areal_radius) + 4. * areal_radius * mass + + 3. * square(mass))) * + pow((4. + 3. * sqrt(2.)) * (2. * areal_radius - 3. * mass) / + (8. * areal_radius + 6. * mass + + 3. * sqrt(8. * square(areal_radius) + + 8. * areal_radius * mass + 6. * square(mass))), + 1. / sqrt(2.)); +} + +template +DataType maximal_isotropic_radius_from_areal_deriv( + const DataType& areal_radius, const DataType& isotropic_radius, + const double mass) { + // r_isotropic = (A/4) * B^(1/sqrt(2)) + const DataType S = sqrt(8. * square(areal_radius) + 8. * areal_radius * mass + + 6. * square(mass)); + const double C = 4. + 3. * sqrt(2.); + const DataType D = 8. * areal_radius + 6. * mass + 3. * S; + const DataType E = 2. * areal_radius - 3. * mass; + const DataType F = 2. * areal_radius + mass; + const DataType A = F + S / sqrt(2.); + const DataType B = C * E / D; + const DataType dAdR = 2. + (4. / sqrt(2.)) * F / S; + const DataType dBdR = C * (2. * D - E * (8. + 12. * F / S)) / square(D); + + return isotropic_radius * (dAdR / A + dBdR / (B * sqrt(2.))); +} + +template +DataType areal_radius_from_maximal_isotropic(const DataType& isotropic_radius, + const double mass) { + const auto residual = [&isotropic_radius, &mass](const auto areal_radius, + const size_t i = 0) { + if constexpr (simd::is_batch>::value) { + return maximal_isotropic_radius_from_areal(areal_radius, mass) - + simd::load_unaligned(&(get_element(isotropic_radius, i))); + } else { + return maximal_isotropic_radius_from_areal(areal_radius, mass) - + get_element(isotropic_radius, i); + } + }; + const auto lower_bound = + make_with_value(isotropic_radius, 1.5 * mass); + const auto upper_bound = make_with_value(isotropic_radius, 1e9); + return RootFinder::toms748(residual, lower_bound, upper_bound, 1.0e-12, + 1.0e-15); +} + } // namespace SchwarzschildImpl::SchwarzschildImpl( @@ -138,6 +196,8 @@ double SchwarzschildImpl::radius_at_horizon() const { return 2. * mass_; case SchwarzschildCoordinates::KerrSchildIsotropic: return kerr_schild_isotropic_radius_from_areal(2. * mass_, mass_); + case SchwarzschildCoordinates::MaximalIsotropic: + return maximal_isotropic_radius_from_areal(2. * mass_, mass_); // LCOV_EXCL_START default: ERROR("Missing case for SchwarzschildCoordinates"); @@ -172,13 +232,22 @@ void SchwarzschildVariables::operator()( const gsl::not_null*> areal_radius, const gsl::not_null cache, detail::Tags::ArealRadius /*meta*/) const { - ASSERT(coordinate_system == SchwarzschildCoordinates::KerrSchildIsotropic, - "The areal radius is only needed for 'KerrSchildIsotropic' " - "coordinates."); const auto& isotropic_radius = cache->get_var(*this, detail::Tags::Radius{}); - get(*areal_radius) = - kerr_schild_areal_radius_from_isotropic(get(isotropic_radius), mass); + switch (coordinate_system) { + case SchwarzschildCoordinates::KerrSchildIsotropic: + get(*areal_radius) = + kerr_schild_areal_radius_from_isotropic(get(isotropic_radius), mass); + break; + case SchwarzschildCoordinates::MaximalIsotropic: + get(*areal_radius) = + areal_radius_from_maximal_isotropic(get(isotropic_radius), mass); + break; + default: + ERROR( + "The areal radius is only needed for 'KerrSchildIsotropic' " + "or 'MaximalIsotropic' coordinates."); + } } template @@ -241,6 +310,10 @@ void SchwarzschildVariables::operator()( get(*trace_extrinsic_curvature) = 2. * mass * cube(lapse) / square(r) * (1. + 3. * mass / r); break; + } + case SchwarzschildCoordinates::MaximalIsotropic: { + get(*trace_extrinsic_curvature) = 0.; + break; } // LCOV_EXCL_START default: @@ -303,6 +376,11 @@ void SchwarzschildVariables::operator()( get<2>(*trace_extrinsic_curvature_gradient) = isotropic_prefactor * get<2>(x); break; + } + case SchwarzschildCoordinates::MaximalIsotropic: { + std::fill(trace_extrinsic_curvature_gradient->begin(), + trace_extrinsic_curvature_gradient->end(), 0.); + break; } // LCOV_EXCL_START default: @@ -332,6 +410,12 @@ void SchwarzschildVariables::operator()( get(cache->get_var(*this, Xcts::Tags::ConformalFactor{})); get(*conformal_factor_minus_one) = conformal_factor - 1.; break; + } + case SchwarzschildCoordinates::MaximalIsotropic: { + const auto& conformal_factor = + get(cache->get_var(*this, Xcts::Tags::ConformalFactor{})); + get(*conformal_factor_minus_one) = conformal_factor - 1.; + break; } // LCOV_EXCL_START default: @@ -362,6 +446,14 @@ void SchwarzschildVariables::operator()( get(cache->get_var(*this, gr::Tags::Lapse{})); get(*conformal_factor) = 2. * exp(1. / lapse - 1.) / (1. + 1. / lapse); break; + } + case SchwarzschildCoordinates::MaximalIsotropic: { + const auto r_areal = + get(cache->get_var(*this, detail::Tags::ArealRadius{})); + const auto& r_iso = + get(cache->get_var(*this, detail::Tags::Radius{})); + get(*conformal_factor) = sqrt(r_areal / r_iso); + break; } // LCOV_EXCL_START default: @@ -406,6 +498,22 @@ void SchwarzschildVariables::operator()( get<1>(*conformal_factor_gradient) = isotropic_prefactor * get<1>(x); get<2>(*conformal_factor_gradient) = isotropic_prefactor * get<2>(x); break; + } + case SchwarzschildCoordinates::MaximalIsotropic: { + const auto& r_iso = + get(cache->get_var(*this, detail::Tags::Radius{})); + const auto& r_areal = + get(cache->get_var(*this, detail::Tags::ArealRadius{})); + + const DataType drdR = + maximal_isotropic_radius_from_areal_deriv(r_areal, r_iso, mass); + const DataType dPsidR = .5 * sqrt(r_iso / r_areal) * + (1. / r_iso - (r_areal / square(r_iso)) * drdR); + const DataType prefactor = dPsidR / (drdR * r_iso); + get<0>(*conformal_factor_gradient) = prefactor * get<0>(x); + get<1>(*conformal_factor_gradient) = prefactor * get<1>(x); + get<2>(*conformal_factor_gradient) = prefactor * get<2>(x); + break; } // LCOV_EXCL_START default: @@ -437,6 +545,13 @@ void SchwarzschildVariables::operator()( get(*lapse_times_conformal_factor_minus_one) = lapse_times_conformal_factor - 1.; break; + } + case SchwarzschildCoordinates::MaximalIsotropic: { + const auto& lapse_times_conformal_factor = get(cache->get_var( + *this, Xcts::Tags::LapseTimesConformalFactor{})); + get(*lapse_times_conformal_factor_minus_one) = + lapse_times_conformal_factor - 1.; + break; } // LCOV_EXCL_START default: @@ -469,6 +584,14 @@ void SchwarzschildVariables::operator()( get(cache->get_var(*this, Xcts::Tags::ConformalFactor{})); get(*lapse_times_conformal_factor) = lapse * conformal_factor; break; + } + case SchwarzschildCoordinates::MaximalIsotropic: { + const auto& lapse = + get(cache->get_var(*this, gr::Tags::Lapse{})); + const auto& conformal_factor = + get(cache->get_var(*this, Xcts::Tags::ConformalFactor{})); + get(*lapse_times_conformal_factor) = lapse * conformal_factor; + break; } // LCOV_EXCL_START default: @@ -500,6 +623,14 @@ void SchwarzschildVariables::operator()( get(cache->get_var(*this, detail::Tags::ArealRadius{})); get(*lapse) = 1. / sqrt(1. + 2. * mass / r); break; + } + case SchwarzschildCoordinates::MaximalIsotropic: { + const auto& r_areal = + get(cache->get_var(*this, detail::Tags::ArealRadius{})); + get(*lapse) = + sqrt(1. - 2. * mass / r_areal + + 27. * square(square(mass)) / (16. * square(square(r_areal)))); + break; } // LCOV_EXCL_START default: @@ -544,6 +675,26 @@ void SchwarzschildVariables::operator()( get<1>(*deriv_lapse) = isotropic_prefactor * get<1>(x); get<2>(*deriv_lapse) = isotropic_prefactor * get<2>(x); break; + } + case SchwarzschildCoordinates::MaximalIsotropic: { + const auto& r_areal = + get(cache->get_var(*this, detail::Tags::ArealRadius{})); + const auto& r_iso = + get(cache->get_var(*this, detail::Tags::Radius{})); + const auto& lapse = + get(cache->get_var(*this, gr::Tags::Lapse{})); + + const DataType drdR = + maximal_isotropic_radius_from_areal_deriv(r_areal, r_iso, mass); + const DataType dLapsedR = + (1. / (2. * lapse)) * + (2. * mass / square(r_areal) - + 27. * square(square(mass)) / (4. * pow(r_areal, 5))); + const DataType prefactor = dLapsedR / (drdR * r_iso); + get<0>(*deriv_lapse) = prefactor * get<0>(x); + get<1>(*deriv_lapse) = prefactor * get<1>(x); + get<2>(*deriv_lapse) = prefactor * get<2>(x); + break; } // LCOV_EXCL_START default: @@ -599,6 +750,29 @@ void SchwarzschildVariables::operator()( get<2>(*lapse_times_conformal_factor_gradient) += isotropic_prefactor * get<2>(x); break; + } + case SchwarzschildCoordinates::MaximalIsotropic: { + const auto& conformal_factor = + get(cache->get_var(*this, Xcts::Tags::ConformalFactor{})); + const auto& conformal_factor_gradient = cache->get_var( + *this, ::Tags::deriv, + tmpl::size_t<3>, Frame::Inertial>{}); + const auto& lapse = + get(cache->get_var(*this, gr::Tags::Lapse{})); + const auto& deriv_lapse = cache->get_var( + *this, ::Tags::deriv, tmpl::size_t<3>, + Frame::Inertial>{}); + *lapse_times_conformal_factor_gradient = conformal_factor_gradient; + get<0>(*lapse_times_conformal_factor_gradient) *= lapse; + get<1>(*lapse_times_conformal_factor_gradient) *= lapse; + get<2>(*lapse_times_conformal_factor_gradient) *= lapse; + get<0>(*lapse_times_conformal_factor_gradient) += + conformal_factor * get<0>(deriv_lapse); + get<1>(*lapse_times_conformal_factor_gradient) += + conformal_factor * get<1>(deriv_lapse); + get<2>(*lapse_times_conformal_factor_gradient) += + conformal_factor * get<2>(deriv_lapse); + break; } // LCOV_EXCL_START default: @@ -660,6 +834,17 @@ void SchwarzschildVariables::operator()( get<1>(*shift_excess) *= isotropic_prefactor; get<2>(*shift_excess) *= isotropic_prefactor; break; + } + case SchwarzschildCoordinates::MaximalIsotropic: { + const auto& r_areal = + get(cache->get_var(*this, detail::Tags::ArealRadius{})); + get<0>(*shift_excess) = + .75 * sqrt(3) * square(mass) * get<0>(x) / pow(r_areal, 3); + get<1>(*shift_excess) = + .75 * sqrt(3) * square(mass) * get<1>(x) / pow(r_areal, 3); + get<2>(*shift_excess) = + .75 * sqrt(3) * square(mass) * get<2>(x) / pow(r_areal, 3); + break; } // LCOV_EXCL_START default: @@ -713,6 +898,30 @@ void SchwarzschildVariables::operator()( deriv_shift_excess->get(i, i) += diagonal_prefactor; } break; + } + case SchwarzschildCoordinates::MaximalIsotropic: { + const auto& r_areal = + get(cache->get_var(*this, detail::Tags::ArealRadius{})); + const auto& r_iso = + get(cache->get_var(*this, detail::Tags::Radius{})); + const DataType BetaRadial = + .75 * sqrt(3) * square(mass) * r_iso / pow(r_areal, 3); + + const DataType drdR = + maximal_isotropic_radius_from_areal_deriv(r_areal, r_iso, mass); + const DataType dBetaRadialdR = + .75 * sqrt(3) * square(mass) * + (drdR / pow(r_areal, 3) - 3. * r_iso / pow(r_areal, 4)); + for (size_t i = 0; i < 3; ++i) { + for (size_t j = 0; j < 3; ++j) { + deriv_shift_excess->get(i, j) = + (dBetaRadialdR / (drdR * square(r_iso)) - + BetaRadial / pow(r_iso, 3)) * + x.get(i) * x.get(j); + } + deriv_shift_excess->get(i, i) += BetaRadial / r_iso; + } + break; } // LCOV_EXCL_START default: @@ -761,6 +970,22 @@ void SchwarzschildVariables::operator()( longitudinal_shift_minus_dt_conformal_metric, cache->get_var(*this, gr::Tags::TraceExtrinsicCurvature{})); break; + } + case SchwarzschildCoordinates::MaximalIsotropic: { + // Background shift and \bar{u} are both zero + const auto& longitudinal_shift_minus_dt_conformal_metric = cache->get_var( + *this, + Xcts::Tags::LongitudinalShiftExcess{}); + Xcts::extrinsic_curvature( + extrinsic_curvature, + cache->get_var(*this, Xcts::Tags::ConformalFactor{}), + cache->get_var(*this, gr::Tags::Lapse{}), + cache->get_var( + *this, + Xcts::Tags::ConformalMetric{}), + longitudinal_shift_minus_dt_conformal_metric, + cache->get_var(*this, gr::Tags::TraceExtrinsicCurvature{})); + break; } // LCOV_EXCL_START default: diff --git a/src/PointwiseFunctions/AnalyticSolutions/Xcts/Schwarzschild.hpp b/src/PointwiseFunctions/AnalyticSolutions/Xcts/Schwarzschild.hpp index 2ffd965f712b..9f01c1d3b7f8 100644 --- a/src/PointwiseFunctions/AnalyticSolutions/Xcts/Schwarzschild.hpp +++ b/src/PointwiseFunctions/AnalyticSolutions/Xcts/Schwarzschild.hpp @@ -160,6 +160,54 @@ enum class SchwarzschildCoordinates { * \end{equation} */ KerrSchildIsotropic, + /*! + * \brief Maximal Isotropic (Horizon Penetrating) Schwarzschild coordinates + * + * Schwarzschild coordinates with a radial transformation such that the radius + * is isotropic and the coordinates are horizon penetrating. + * + * These arise from first choosing a family of time-independent, maximal + * slicings of the Schwarzschild spacetime and a slicing condition that give a + * unique solution with a limiting surface at \f$R=3M/2\f$ and horizon at + * \f$R=2M\f$ \cite Estabrook1973ue. The latter is then changed by the radial + * transformation \cite Baumgarte2007ht + * + * \f{equation} + * r = \frac{\left[2R + M + \sqrt{4R^2 + 4MR + 3M^2}\right]}{4} + * \times \left[ + * \frac{(4 + 3\sqrt{2})(2R - 3M)}{8R + 6M + 3\sqrt{8R^2 + 8MR + 6M^2}} + * \right]^{1/\sqrt{2}} + * \f} + * + * where \f$R\f$ is the canonical + * Schwarzschild radius, also referred to as "areal" radius because it is + * defined such that spheres with constant \f$R\f$ have the area \f$4\pi + * R^2\f$, and \f$r\f$ is the "isotropic" radius. In the isotropic + * radius the Schwarzschild spatial metric is conformally flat: + * + * \f{equation} + * \gamma_{ij}=\psi^4\eta_{ij} \quad \text{with conformal factor} \quad + * \psi=\sqrt{\frac{R}{r}}. + * \f} + * + * The lapse in the conformal radius is + * + * \f{equation} + * \alpha = \left(1 - \frac{2M}{R} + \frac{27M^4}{16R^4} \right)^{1/2} + * \f} + * + * and the shift is given by + * + * \f{equation} + * \beta^r = \frac{3 \sqrt{3} M^2}{4} \frac{r}{R^3} + * \f} + * + * The solution remains maximally sliced, i.e. \f$K=0\f$. And the horizon in + * these coordinates is at \f$r\approx 0.7793271080557972 M\f$ due to the + * radial transformation from \f$R=2M\f$. + * + */ + MaximalIsotropic, }; std::ostream& operator<<(std::ostream& os, SchwarzschildCoordinates coords); diff --git a/tests/Unit/PointwiseFunctions/AnalyticSolutions/Xcts/SchwarzschildMaximalIsotropic.py b/tests/Unit/PointwiseFunctions/AnalyticSolutions/Xcts/SchwarzschildMaximalIsotropic.py new file mode 100644 index 000000000000..51ba34523872 --- /dev/null +++ b/tests/Unit/PointwiseFunctions/AnalyticSolutions/Xcts/SchwarzschildMaximalIsotropic.py @@ -0,0 +1,250 @@ +# Distributed under the MIT License. +# See LICENSE.txt for details. + +import numpy as np +from numpy import exp, sqrt +from scipy.optimize import newton + +# Schwarzschild Maximal Isotropic coordinates + + +def isotropic_radius_from_areal(r_areal, mass): + # First part of the equation + part1 = ( + 2 * r_areal + + mass + + sqrt(4 * r_areal**2 + 4 * mass * r_areal + 3 * mass**2) + ) / 4.0 + + # Second part of the equation + part2_numerator = (4 + 3 * sqrt(2)) * (2 * r_areal - 3 * mass) + part2_denominator = ( + 8 * r_areal + + 6 * mass + + 3 * sqrt(8 * r_areal**2 + 8 * mass * r_areal + 6 * mass**2) + ) + part2 = (part2_numerator / part2_denominator) ** (1 / sqrt(2)) + + return part1 * part2 + + +def isotropic_radius_from_areal_deriv(r_areal, mass): + h = 1e-6 + term1 = -isotropic_radius_from_areal(r_areal + 2 * h, mass) + term2 = 8 * isotropic_radius_from_areal(r_areal + h, mass) + term3 = -8 * isotropic_radius_from_areal(r_areal - h, mass) + term4 = isotropic_radius_from_areal(r_areal - 2 * h, mass) + + return (term1 + term2 + term3 + term4) / (12 * h) + + +def areal_radius_from_isotropic(r_isotropic, mass): + def f(r_areal): + return isotropic_radius_from_areal(r_areal, mass) - r_isotropic + + def fprime(r_areal): + return isotropic_radius_from_areal_deriv(r_areal, mass) + + return newton(func=f, fprime=fprime, x0=r_isotropic, tol=1.0e-12) + + +def conformal_metric(x, mass): + return np.identity(3) + + +def inv_conformal_metric(x, mass): + return np.identity(3) + + +def deriv_conformal_metric(x, mass): + return np.zeros((3, 3, 3)) + + +def extrinsic_curvature_trace(x, mass): + return 0.0 + + +def extrinsic_curvature_trace_gradient(x, mass): + return np.zeros(3) + + +def conformal_factor(x, mass): + r_isotropic = np.linalg.norm(x) + r_areal = areal_radius_from_isotropic(r_isotropic, mass) + return sqrt(r_areal / r_isotropic) + + +def conformal_factor_gradient(x, mass): + h = 1e-6 + + x_forward_2h = np.array([x[0] + 2 * h, x[1], x[2]]) + x_forward_h = np.array([x[0] + h, x[1], x[2]]) + x_backward_h = np.array([x[0] - h, x[1], x[2]]) + x_backward_2h = np.array([x[0] - 2 * h, x[1], x[2]]) + dx = ( + 1 / 12 * conformal_factor(x_backward_2h, mass) + - 2 / 3 * conformal_factor(x_backward_h, mass) + + 2 / 3 * conformal_factor(x_forward_h, mass) + - 1 / 12 * conformal_factor(x_forward_2h, mass) + ) / h + + y_forward_2h = np.array([x[0], x[1] + 2 * h, x[2]]) + y_forward_h = np.array([x[0], x[1] + h, x[2]]) + y_backward_h = np.array([x[0], x[1] - h, x[2]]) + y_backward_2h = np.array([x[0], x[1] - 2 * h, x[2]]) + dy = ( + 1 / 12 * conformal_factor(y_backward_2h, mass) + - 2 / 3 * conformal_factor(y_backward_h, mass) + + 2 / 3 * conformal_factor(y_forward_h, mass) + - 1 / 12 * conformal_factor(y_forward_2h, mass) + ) / h + + z_forward_2h = np.array([x[0], x[1], x[2] + 2 * h]) + z_forward_h = np.array([x[0], x[1], x[2] + h]) + z_backward_h = np.array([x[0], x[1], x[2] - h]) + z_backward_2h = np.array([x[0], x[1], x[2] - 2 * h]) + dz = ( + 1 / 12 * conformal_factor(z_backward_2h, mass) + - 2 / 3 * conformal_factor(z_backward_h, mass) + + 2 / 3 * conformal_factor(z_forward_h, mass) + - 1 / 12 * conformal_factor(z_forward_2h, mass) + ) / h + + return np.array([dx, dy, dz]) + + +def lapse(x, mass): + r_isotropic = np.linalg.norm(x) + r_areal = areal_radius_from_isotropic(r_isotropic, mass) + return sqrt(1 - 2 * mass / r_areal + (27 / 16) * mass**4 / r_areal**4) + + +def lapse_times_conformal_factor(x, mass): + lapse_val = lapse(x, mass) + conformal_factor_val = conformal_factor(x, mass) + return lapse_val * conformal_factor_val + + +def lapse_times_conformal_factor_gradient(x, mass): + h = 1e-6 + + x_forward_2h = np.array([x[0] + 2 * h, x[1], x[2]]) + x_forward_h = np.array([x[0] + h, x[1], x[2]]) + x_backward_h = np.array([x[0] - h, x[1], x[2]]) + x_backward_2h = np.array([x[0] - 2 * h, x[1], x[2]]) + dx = ( + 1 / 12 * lapse_times_conformal_factor(x_backward_2h, mass) + - 2 / 3 * lapse_times_conformal_factor(x_backward_h, mass) + + 2 / 3 * lapse_times_conformal_factor(x_forward_h, mass) + - 1 / 12 * lapse_times_conformal_factor(x_forward_2h, mass) + ) / h + + y_forward_2h = np.array([x[0], x[1] + 2 * h, x[2]]) + y_forward_h = np.array([x[0], x[1] + h, x[2]]) + y_backward_h = np.array([x[0], x[1] - h, x[2]]) + y_backward_2h = np.array([x[0], x[1] - 2 * h, x[2]]) + dy = ( + 1 / 12 * lapse_times_conformal_factor(y_backward_2h, mass) + - 2 / 3 * lapse_times_conformal_factor(y_backward_h, mass) + + 2 / 3 * lapse_times_conformal_factor(y_forward_h, mass) + - 1 / 12 * lapse_times_conformal_factor(y_forward_2h, mass) + ) / h + + z_forward_2h = np.array([x[0], x[1], x[2] + 2 * h]) + z_forward_h = np.array([x[0], x[1], x[2] + h]) + z_backward_h = np.array([x[0], x[1], x[2] - h]) + z_backward_2h = np.array([x[0], x[1], x[2] - 2 * h]) + dz = ( + 1 / 12 * lapse_times_conformal_factor(z_backward_2h, mass) + - 2 / 3 * lapse_times_conformal_factor(z_backward_h, mass) + + 2 / 3 * lapse_times_conformal_factor(z_forward_h, mass) + - 1 / 12 * lapse_times_conformal_factor(z_forward_2h, mass) + ) / h + + return np.array([dx, dy, dz]) + + +def shift_background(x, mass): + return np.zeros(3) + + +def longitudinal_shift_background_minus_dt_conformal_metric(x, mass): + return np.zeros((3, 3)) + + +def shift(x, mass): + r_isotropic = np.linalg.norm(x) + r_areal = areal_radius_from_isotropic(r_isotropic, mass) + return ( + (0.75 * sqrt(3.0) * mass**2 * (r_isotropic / r_areal**3)) + * x + / r_isotropic + ) + + +def shift_strain(x, mass): + h = 1e-6 + n = len(x) + strain_tensor = np.zeros((n, n)) + + for i in range(n): + for j in range(n): + x_forward_2h = np.copy(x) + x_forward_h = np.copy(x) + x_backward_h = np.copy(x) + x_backward_2h = np.copy(x) + + x_forward_2h[j] += 2 * h + x_forward_h[j] += h + x_backward_h[j] -= h + x_backward_2h[j] -= 2 * h + + beta_forward_2h = shift(x_forward_2h, mass)[i] + beta_forward_h = shift(x_forward_h, mass)[i] + beta_backward_h = shift(x_backward_h, mass)[i] + beta_backward_2h = shift(x_backward_2h, mass)[i] + + strain_tensor[i, j] = ( + 1 / 12 * beta_backward_2h + - 2 / 3 * beta_backward_h + + 2 / 3 * beta_forward_h + - 1 / 12 * beta_forward_2h + ) / h + + return strain_tensor + + +def longitudinal_shift(x, mass): + B = shift_strain(x, mass) + return 2 * (B - 1.0 / 3.0 * np.identity(3) * np.trace(B)) + + +def shift_dot_extrinsic_curvature_trace_gradient(x, mass): + beta = shift(x, mass) + dK = extrinsic_curvature_trace_gradient(x, mass) + return np.dot(beta, dK) + + +def longitudinal_shift_minus_dt_conformal_metric_square(x, mass): + LB = longitudinal_shift(x, mass) + return np.einsum("ij,ij", LB, LB) + + +def longitudinal_shift_minus_dt_conformal_metric_over_lapse_square(x, mass): + LB = longitudinal_shift(x, mass) + return np.einsum("ij,ij", LB, LB) / lapse(x, mass) ** 2 + + +# Matter sources + + +def energy_density(x, mass): + return 0.0 + + +def stress_trace(x, mass): + return 0.0 + + +def momentum_density(x, mass): + return np.zeros(3) diff --git a/tests/Unit/PointwiseFunctions/AnalyticSolutions/Xcts/Test_Schwarzschild.cpp b/tests/Unit/PointwiseFunctions/AnalyticSolutions/Xcts/Test_Schwarzschild.cpp index dc74e272bf25..9d6b74c0b650 100644 --- a/tests/Unit/PointwiseFunctions/AnalyticSolutions/Xcts/Test_Schwarzschild.cpp +++ b/tests/Unit/PointwiseFunctions/AnalyticSolutions/Xcts/Test_Schwarzschild.cpp @@ -132,7 +132,8 @@ void test_solution(const double mass, "shift_dot_extrinsic_curvature_trace_gradient", "longitudinal_shift_minus_dt_conformal_metric_square", "longitudinal_shift_minus_dt_conformal_metric_over_lapse_square"}, - {{{inner_radius, outer_radius}}}, std::make_tuple(mass), DataVector(5)); + {{{inner_radius, outer_radius}}}, std::make_tuple(mass), DataVector(3), + 1e-8); } { INFO("Verify the solution solves the XCTS equations"); @@ -185,6 +186,18 @@ SPECTRE_TEST_CASE( "Schwarzschild:\n" " Mass: 0.8\n" " Coordinates: KerrSchildIsotropic"); + test_solution(1., SchwarzschildCoordinates::MaximalIsotropic, + 0.7793271080557972, 1.601972683625847, + "SchwarzschildMaximalIsotropic", + "Schwarzschild:\n" + " Mass: 1.\n" + " Coordinates: MaximalIsotropic"); + test_solution(0.8, SchwarzschildCoordinates::MaximalIsotropic, + 0.7793271080557972 * 0.8, 1.601972683625847, + "SchwarzschildMaximalIsotropic", + "Schwarzschild:\n" + " Mass: 0.8\n" + " Coordinates: MaximalIsotropic"); } } // namespace Xcts::Solutions