SpECTRE Documentation Coverage Report
Current view: top level - Elliptic/Systems/Xcts/BoundaryConditions - ApparentHorizon.hpp Hit Total Coverage
Commit: 3c2e9d3ed337bca2146eee9de07432e292a38c3a Lines: 1 44 2.3 %
Date: 2024-06-11 22:56:19
Legend: Lines: hit not hit

          Line data    Source code
       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <cstddef>
       7             : #include <memory>
       8             : #include <optional>
       9             : #include <pup.h>
      10             : #include <string>
      11             : #include <vector>
      12             : 
      13             : #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
      14             : #include "DataStructures/Tensor/TypeAliases.hpp"
      15             : #include "Domain/Tags.hpp"
      16             : #include "Domain/Tags/FaceNormal.hpp"
      17             : #include "Elliptic/BoundaryConditions/BoundaryCondition.hpp"
      18             : #include "Elliptic/BoundaryConditions/BoundaryConditionType.hpp"
      19             : #include "Elliptic/Systems/Xcts/Geometry.hpp"
      20             : #include "Elliptic/Systems/Xcts/Tags.hpp"
      21             : #include "Options/Auto.hpp"
      22             : #include "Options/Context.hpp"
      23             : #include "Options/String.hpp"
      24             : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
      25             : #include "PointwiseFunctions/InitialDataUtilities/AnalyticSolution.hpp"
      26             : #include "Utilities/Gsl.hpp"
      27             : #include "Utilities/MakeArray.hpp"
      28             : #include "Utilities/Serialization/CharmPupable.hpp"
      29             : #include "Utilities/TMPL.hpp"
      30             : 
      31             : /// \cond
      32             : class DataVector;
      33             : /// \endcond
      34             : 
      35           0 : namespace Xcts::BoundaryConditions {
      36             : 
      37             : /*!
      38             :  * \brief Impose the surface is a quasi-equilibrium apparent horizon.
      39             :  *
      40             :  * These boundary conditions on the conformal factor \f$\psi\f$, the lapse
      41             :  * \f$\alpha\f$ and the shift \f$\beta^i\f$ impose the surface is an apparent
      42             :  * horizon, i.e. that the expansion on the surface vanishes: \f$\Theta=0\f$.
      43             :  * Specifically, we impose:
      44             :  *
      45             :  * \f{align}
      46             :  * \label{eq:ah_psi}
      47             :  * \bar{s}^k\bar{D}_k\psi &= -\frac{\psi^3}{8\alpha}\bar{s}_i\bar{s}_j\left(
      48             :  * (\bar{L}\beta)^{ij} - \bar{u}^{ij}\right)
      49             :  * - \frac{\psi}{4}\bar{m}^{ij}\bar{\nabla}_i\bar{s}_j + \frac{1}{6}K\psi^3
      50             :  * \\
      51             :  * \label{eq:ah_beta}
      52             :  * \beta_\mathrm{excess}^i &= \frac{\alpha}{\psi^2}\bar{s}^i
      53             :  * + \epsilon_{ijk}\Omega^j x^k - \beta_\mathrm{background}^i
      54             :  * \f}
      55             :  *
      56             :  * following section 7.2 of \cite Pfeiffer2005zm, section 12.3.2 of
      57             :  * \cite BaumgarteShapiro or section II.B.1 of \cite Varma2018sqd. In these
      58             :  * equations \f$\bar{s}_i\f$ is the conformal surface normal to the apparent
      59             :  * horizon, \f$\bar{m}^{ij}=\bar{\gamma}^{ij}-\bar{s}^i\bar{s}^j\f$ is the
      60             :  * induced conformal surface metric (denoted \f$\tilde{h}^{ij}\f$ in
      61             :  * \cite Pfeiffer2005zm and \cite Varma2018sqd) and \f$\bar{D}\f$ is the
      62             :  * covariant derivative w.r.t. to the conformal metric \f$\bar{\gamma}_{ij}\f$.
      63             :  * Note that e.g. in \cite Varma2018sqd Eq. (16) appears the surface-normal
      64             :  * \f$s^i\f$, not the _conformal_ surface normal \f$\bar{s}^i = \psi^2 s^i\f$.
      65             :  * To incur a spin on the apparent horizon we can freely choose the rotational
      66             :  * parameters \f$\boldsymbol{\Omega}\f$. Note that for a Kerr solution with
      67             :  * dimensionless spin \f$\boldsymbol{\chi}\f$ the rotational parameters at the
      68             :  * outer horizon are \f$\boldsymbol{\Omega} =
      69             :  * -\frac{\boldsymbol{\chi}}{2r_+}\f$, where \f$r_+ / M = 1 + \sqrt{1 -
      70             :  * \chi^2}\f$ (see e.g. Eq. (8) in \cite Ossokine2015yla). \f$\epsilon_{ijk}\f$
      71             :  * is the flat-space Levi-Civita symbol.
      72             :  *
      73             :  * Note that the quasi-equilibrium conditions don't restrict the boundary
      74             :  * condition for the lapse. The choice for the lapse boundary condition is made
      75             :  * by the `Lapse` input-file options. Currently implemented choices are:
      76             :  * - A zero von-Neumann boundary condition:
      77             :  *   \f$\bar{s}^k\bar{D}_k(\alpha\psi) = 0\f$
      78             :  * - A Dirichlet boundary condition imposing the lapse of an analytic solution.
      79             :  *
      80             :  * \par Negative-expansion boundary conditions:
      81             :  * This class also supports negative-expansion boundary conditions following
      82             :  * section II.B.2 of \cite Varma2018sqd by taking a Kerr solution as additional
      83             :  * parameter and computing its expansion at the excision surface. Choosing an
      84             :  * excision surface _within_ the apparent horizon of the Kerr solution will
      85             :  * result in a negative expansion that is added to the boundary condition for
      86             :  * the conformal factor. Therefore, the excision surface will lie within an
      87             :  * apparent horizon. Specifically, we add the quantity
      88             :  * \f$\frac{\psi^3}{4}\Theta_\mathrm{Kerr}\f$ to (\f$\ref{eq:ah_psi}\f$), where
      89             :  * \f$\Theta_\mathrm{Kerr}\f$ is the expansion of the specified Kerr solution at
      90             :  * the excision surface, and we add the quantity \f$\epsilon = s_i
      91             :  * \beta_\mathrm{Kerr}^i - \alpha_\mathrm{Kerr}\f$ to the orthogonal part
      92             :  * \f$s_i\beta_\mathrm{excess}^i\f$ of (\f$\ref{eq:ah_beta}\f$).
      93             :  *
      94             :  * \note When negative-expansion boundary conditions are selected, the Kerr
      95             :  * solution gets evaluated at every boundary-condition application. This may
      96             :  * turn out to incur a significant computational cost, in which case a future
      97             :  * optimization might be to pre-compute the negative-expansion quantities at
      98             :  * initialization time and store them in the DataBox.
      99             :  */
     100             : template <Xcts::Geometry ConformalGeometry>
     101           1 : class ApparentHorizon
     102             :     : public elliptic::BoundaryConditions::BoundaryCondition<3> {
     103             :  private:
     104           0 :   using Base = elliptic::BoundaryConditions::BoundaryCondition<3>;
     105             : 
     106             :  public:
     107           0 :   static constexpr Options::String help =
     108             :       "Impose the boundary is a quasi-equilibrium apparent horizon.";
     109             : 
     110           0 :   struct Center {
     111           0 :     using type = std::array<double, 3>;
     112           0 :     static constexpr Options::String help =
     113             :         "The center of the excision surface representing the apparent-horizon "
     114             :         "surface";
     115             :   };
     116           0 :   struct Rotation {
     117           0 :     using type = std::array<double, 3>;
     118           0 :     static constexpr Options::String help =
     119             :         "The rotational parameters 'Omega' on the surface, which parametrize "
     120             :         "the spin of the black hole. The rotational parameters enter the "
     121             :         "Dirichlet boundary conditions for the shift in a term "
     122             :         "'Omega x (r - Center)', where 'r' are the coordinates on the surface.";
     123             :   };
     124           0 :   struct Lapse {
     125           0 :     using type = Options::Auto<
     126             :         std::unique_ptr<elliptic::analytic_data::AnalyticSolution>>;
     127           0 :     static constexpr Options::String help =
     128             :         "Specify an analytic solution to impose a Dirichlet condition on the "
     129             :         "lapse. The analytic solution will be evaluated at coordinates "
     130             :         "centered at the apparent horizon. "
     131             :         "Alternatively, set this option to 'None' "
     132             :         "to impose a zero von-Neumann boundary condition on the lapse. Note "
     133             :         "that the latter will not result in the standard Kerr-Schild slicing "
     134             :         "for a single black hole.";
     135             :   };
     136           0 :   struct NegativeExpansion {
     137           0 :     using type = Options::Auto<
     138             :         std::unique_ptr<elliptic::analytic_data::AnalyticSolution>,
     139             :         Options::AutoLabel::None>;
     140           0 :     static constexpr Options::String help =
     141             :         "Specify an analytic solution to impose its expansion at the excision "
     142             :         "surface. The analytic solution will be evaluated at coordinates "
     143             :         "centered at the apparent horizon. "
     144             :         "If the excision surface lies within the solution's "
     145             :         "apparent horizon, the imposed expansion will be negative and thus the "
     146             :         "excision surface will lie within an apparent horizon. Alternatively, "
     147             :         "set this option to 'None' to impose the expansion is zero at the "
     148             :         "excision surface, meaning the excision surface _is_ an apparent "
     149             :         "horizon.";
     150             :   };
     151             : 
     152           0 :   using options = tmpl::list<Center, Rotation, Lapse, NegativeExpansion>;
     153             : 
     154           0 :   ApparentHorizon() = default;
     155           0 :   ApparentHorizon(const ApparentHorizon&) = delete;
     156           0 :   ApparentHorizon& operator=(const ApparentHorizon&) = delete;
     157           0 :   ApparentHorizon(ApparentHorizon&&) = default;
     158           0 :   ApparentHorizon& operator=(ApparentHorizon&&) = default;
     159           0 :   ~ApparentHorizon() = default;
     160             : 
     161             :   /// \cond
     162             :   explicit ApparentHorizon(CkMigrateMessage* m) : Base(m) {}
     163             :   using PUP::able::register_constructor;
     164             :   WRAPPED_PUPable_decl_template(ApparentHorizon);
     165             :   /// \endcond
     166             : 
     167           0 :   std::unique_ptr<domain::BoundaryConditions::BoundaryCondition> get_clone()
     168             :       const override {
     169             :     return std::make_unique<ApparentHorizon>(
     170             :         center_, rotation_,
     171             :         solution_for_lapse_.has_value()
     172             :             ? std::make_optional(solution_for_lapse_.value()->get_clone())
     173             :             : std::nullopt,
     174             :         solution_for_negative_expansion_.has_value()
     175             :             ? std::make_optional(
     176             :                   solution_for_negative_expansion_.value()->get_clone())
     177             :             : std::nullopt);
     178             :   }
     179             : 
     180           0 :   ApparentHorizon(
     181             :       std::array<double, 3> center, std::array<double, 3> rotation,
     182             :       std::optional<std::unique_ptr<elliptic::analytic_data::AnalyticSolution>>
     183             :           solution_for_lapse,
     184             :       std::optional<std::unique_ptr<elliptic::analytic_data::AnalyticSolution>>
     185             :           solution_for_negative_expansion,
     186             :       const Options::Context& context = {});
     187             : 
     188           0 :   const std::array<double, 3>& center() const { return center_; }
     189           0 :   const std::array<double, 3>& rotation() const { return rotation_; }
     190             :   const std::optional<
     191             :       std::unique_ptr<elliptic::analytic_data::AnalyticSolution>>&
     192           0 :   solution_for_lapse() const {
     193             :     return solution_for_lapse_;
     194             :   }
     195             :   const std::optional<
     196             :       std::unique_ptr<elliptic::analytic_data::AnalyticSolution>>&
     197           0 :   solution_for_negative_expansion() const {
     198             :     return solution_for_negative_expansion_;
     199             :   }
     200             : 
     201           0 :   std::vector<elliptic::BoundaryConditionType> boundary_condition_types()
     202             :       const override {
     203             :     return {// Conformal factor
     204             :             elliptic::BoundaryConditionType::Neumann,
     205             :             // Lapse times conformal factor
     206             :             this->solution_for_lapse_.has_value()
     207             :                 ? elliptic::BoundaryConditionType::Dirichlet
     208             :                 : elliptic::BoundaryConditionType::Neumann,
     209             :             // Shift
     210             :             elliptic::BoundaryConditionType::Dirichlet,
     211             :             elliptic::BoundaryConditionType::Dirichlet,
     212             :             elliptic::BoundaryConditionType::Dirichlet};
     213             :   }
     214             : 
     215           0 :   using argument_tags = tmpl::flatten<tmpl::list<
     216             :       domain::Tags::FaceNormal<3>,
     217             :       ::Tags::deriv<domain::Tags::UnnormalizedFaceNormal<3>, tmpl::size_t<3>,
     218             :                     Frame::Inertial>,
     219             :       domain::Tags::UnnormalizedFaceNormalMagnitude<3>,
     220             :       domain::Tags::Coordinates<3, Frame::Inertial>,
     221             :       gr::Tags::TraceExtrinsicCurvature<DataVector>,
     222             :       Tags::ShiftBackground<DataVector, 3, Frame::Inertial>,
     223             :       Tags::LongitudinalShiftBackgroundMinusDtConformalMetric<DataVector, 3,
     224             :                                                               Frame::Inertial>,
     225             :       tmpl::conditional_t<ConformalGeometry == Xcts::Geometry::Curved,
     226             :                           tmpl::list<Tags::InverseConformalMetric<
     227             :                                          DataVector, 3, Frame::Inertial>,
     228             :                                      Tags::ConformalChristoffelSecondKind<
     229             :                                          DataVector, 3, Frame::Inertial>>,
     230             :                           tmpl::list<>>>>;
     231           0 :   using volume_tags = tmpl::list<>;
     232             : 
     233           0 :   void apply(
     234             :       gsl::not_null<Scalar<DataVector>*> conformal_factor_minus_one,
     235             :       gsl::not_null<Scalar<DataVector>*> lapse_times_conformal_factor_minus_one,
     236             :       gsl::not_null<tnsr::I<DataVector, 3>*> shift_excess,
     237             :       gsl::not_null<Scalar<DataVector>*> n_dot_conformal_factor_gradient,
     238             :       gsl::not_null<Scalar<DataVector>*>
     239             :           n_dot_lapse_times_conformal_factor_gradient,
     240             :       gsl::not_null<tnsr::I<DataVector, 3>*> n_dot_longitudinal_shift_excess,
     241             :       const tnsr::i<DataVector, 3>& deriv_conformal_factor,
     242             :       const tnsr::i<DataVector, 3>& deriv_lapse_times_conformal_factor,
     243             :       const tnsr::iJ<DataVector, 3>& deriv_shift_excess,
     244             :       const tnsr::i<DataVector, 3>& face_normal,
     245             :       const tnsr::ij<DataVector, 3>& deriv_unnormalized_face_normal,
     246             :       const Scalar<DataVector>& face_normal_magnitude,
     247             :       const tnsr::I<DataVector, 3>& x,
     248             :       const Scalar<DataVector>& extrinsic_curvature_trace,
     249             :       const tnsr::I<DataVector, 3>& shift_background,
     250             :       const tnsr::II<DataVector, 3>& longitudinal_shift_background) const;
     251             : 
     252           0 :   void apply(
     253             :       gsl::not_null<Scalar<DataVector>*> conformal_factor_minus_one,
     254             :       gsl::not_null<Scalar<DataVector>*> lapse_times_conformal_factor_minus_one,
     255             :       gsl::not_null<tnsr::I<DataVector, 3>*> shift_excess,
     256             :       gsl::not_null<Scalar<DataVector>*> n_dot_conformal_factor_gradient,
     257             :       gsl::not_null<Scalar<DataVector>*>
     258             :           n_dot_lapse_times_conformal_factor_gradient,
     259             :       gsl::not_null<tnsr::I<DataVector, 3>*> n_dot_longitudinal_shift_excess,
     260             :       const tnsr::i<DataVector, 3>& deriv_conformal_factor,
     261             :       const tnsr::i<DataVector, 3>& deriv_lapse_times_conformal_factor,
     262             :       const tnsr::iJ<DataVector, 3>& deriv_shift_excess,
     263             :       const tnsr::i<DataVector, 3>& face_normal,
     264             :       const tnsr::ij<DataVector, 3>& deriv_unnormalized_face_normal,
     265             :       const Scalar<DataVector>& face_normal_magnitude,
     266             :       const tnsr::I<DataVector, 3>& x,
     267             :       const Scalar<DataVector>& extrinsic_curvature_trace,
     268             :       const tnsr::I<DataVector, 3>& shift_background,
     269             :       const tnsr::II<DataVector, 3>& longitudinal_shift_background,
     270             :       const tnsr::II<DataVector, 3>& inv_conformal_metric,
     271             :       const tnsr::Ijj<DataVector, 3>& conformal_christoffel_second_kind) const;
     272             : 
     273           0 :   using argument_tags_linearized = tmpl::flatten<tmpl::list<
     274             :       ::Tags::Normalized<
     275             :           domain::Tags::UnnormalizedFaceNormal<3, Frame::Inertial>>,
     276             :       ::Tags::deriv<domain::Tags::UnnormalizedFaceNormal<3, Frame::Inertial>,
     277             :                     tmpl::size_t<3>, Frame::Inertial>,
     278             :       ::Tags::Magnitude<
     279             :           domain::Tags::UnnormalizedFaceNormal<3, Frame::Inertial>>,
     280             :       domain::Tags::Coordinates<3, Frame::Inertial>,
     281             :       gr::Tags::TraceExtrinsicCurvature<DataVector>,
     282             :       Tags::LongitudinalShiftBackgroundMinusDtConformalMetric<DataVector, 3,
     283             :                                                               Frame::Inertial>,
     284             :       Tags::ConformalFactorMinusOne<DataVector>,
     285             :       Tags::LapseTimesConformalFactorMinusOne<DataVector>,
     286             :       ::Tags::NormalDotFlux<Tags::ShiftExcess<DataVector, 3, Frame::Inertial>>,
     287             :       tmpl::conditional_t<ConformalGeometry == Xcts::Geometry::Curved,
     288             :                           tmpl::list<Tags::InverseConformalMetric<
     289             :                                          DataVector, 3, Frame::Inertial>,
     290             :                                      Tags::ConformalChristoffelSecondKind<
     291             :                                          DataVector, 3, Frame::Inertial>>,
     292             :                           tmpl::list<>>>>;
     293           0 :   using volume_tags_linearized = tmpl::list<>;
     294             : 
     295           0 :   void apply_linearized(
     296             :       gsl::not_null<Scalar<DataVector>*> conformal_factor_correction,
     297             :       gsl::not_null<Scalar<DataVector>*>
     298             :           lapse_times_conformal_factor_correction,
     299             :       gsl::not_null<tnsr::I<DataVector, 3>*> shift_excess_correction,
     300             :       gsl::not_null<Scalar<DataVector>*>
     301             :           n_dot_conformal_factor_gradient_correction,
     302             :       gsl::not_null<Scalar<DataVector>*>
     303             :           n_dot_lapse_times_conformal_factor_gradient_correction,
     304             :       gsl::not_null<tnsr::I<DataVector, 3>*>
     305             :           n_dot_longitudinal_shift_excess_correction,
     306             :       const tnsr::i<DataVector, 3>& deriv_conformal_factor_correction,
     307             :       const tnsr::i<DataVector, 3>&
     308             :           deriv_lapse_times_conformal_factor_correction,
     309             :       const tnsr::iJ<DataVector, 3>& deriv_shift_excess_correction,
     310             :       const tnsr::i<DataVector, 3>& face_normal,
     311             :       const tnsr::ij<DataVector, 3>& deriv_unnormalized_face_normal,
     312             :       const Scalar<DataVector>& face_normal_magnitude,
     313             :       const tnsr::I<DataVector, 3>& x,
     314             :       const Scalar<DataVector>& extrinsic_curvature_trace,
     315             :       const tnsr::II<DataVector, 3>& longitudinal_shift_background,
     316             :       const Scalar<DataVector>& conformal_factor_minus_one,
     317             :       const Scalar<DataVector>& lapse_times_conformal_factor_minus_one,
     318             :       const tnsr::I<DataVector, 3>& n_dot_longitudinal_shift_excess) const;
     319             : 
     320           0 :   void apply_linearized(
     321             :       gsl::not_null<Scalar<DataVector>*> conformal_factor_correction,
     322             :       gsl::not_null<Scalar<DataVector>*>
     323             :           lapse_times_conformal_factor_correction,
     324             :       gsl::not_null<tnsr::I<DataVector, 3>*> shift_excess_correction,
     325             :       gsl::not_null<Scalar<DataVector>*>
     326             :           n_dot_conformal_factor_gradient_correction,
     327             :       gsl::not_null<Scalar<DataVector>*>
     328             :           n_dot_lapse_times_conformal_factor_gradient_correction,
     329             :       gsl::not_null<tnsr::I<DataVector, 3>*>
     330             :           n_dot_longitudinal_shift_excess_correction,
     331             :       const tnsr::i<DataVector, 3>& deriv_conformal_factor_correction,
     332             :       const tnsr::i<DataVector, 3>&
     333             :           deriv_lapse_times_conformal_factor_correction,
     334             :       const tnsr::iJ<DataVector, 3>& deriv_shift_excess_correction,
     335             :       const tnsr::i<DataVector, 3>& face_normal,
     336             :       const tnsr::ij<DataVector, 3>& deriv_unnormalized_face_normal,
     337             :       const Scalar<DataVector>& face_normal_magnitude,
     338             :       const tnsr::I<DataVector, 3>& x,
     339             :       const Scalar<DataVector>& extrinsic_curvature_trace,
     340             :       const tnsr::II<DataVector, 3>& longitudinal_shift_background,
     341             :       const Scalar<DataVector>& conformal_factor_minus_one,
     342             :       const Scalar<DataVector>& lapse_times_conformal_factor_minus_one,
     343             :       const tnsr::I<DataVector, 3>& n_dot_longitudinal_shift_excess,
     344             :       const tnsr::II<DataVector, 3>& inv_conformal_metric,
     345             :       const tnsr::Ijj<DataVector, 3>& conformal_christoffel_second_kind) const;
     346             : 
     347             :   // NOLINTNEXTLINE(google-runtime-references)
     348           0 :   void pup(PUP::er& p) override;
     349             : 
     350             :  private:
     351           0 :   std::array<double, 3> center_ =
     352             :       make_array<3>(std::numeric_limits<double>::signaling_NaN());
     353           0 :   std::array<double, 3> rotation_ =
     354             :       make_array<3>(std::numeric_limits<double>::signaling_NaN());
     355             :   std::optional<std::unique_ptr<elliptic::analytic_data::AnalyticSolution>>
     356           0 :       solution_for_lapse_{};
     357             :   std::optional<std::unique_ptr<elliptic::analytic_data::AnalyticSolution>>
     358           0 :       solution_for_negative_expansion_{};
     359             : };
     360             : 
     361             : }  // namespace Xcts::BoundaryConditions

Generated by: LCOV version 1.14