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
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>& face_normal,
242 : const tnsr::ij<DataVector, 3>& deriv_unnormalized_face_normal,
243 : const Scalar<DataVector>& face_normal_magnitude,
244 : const tnsr::I<DataVector, 3>& x,
245 : const Scalar<DataVector>& extrinsic_curvature_trace,
246 : const tnsr::I<DataVector, 3>& shift_background,
247 : const tnsr::II<DataVector, 3>& longitudinal_shift_background) const;
248 :
249 0 : void apply(
250 : gsl::not_null<Scalar<DataVector>*> conformal_factor_minus_one,
251 : gsl::not_null<Scalar<DataVector>*> lapse_times_conformal_factor_minus_one,
252 : gsl::not_null<tnsr::I<DataVector, 3>*> shift_excess,
253 : gsl::not_null<Scalar<DataVector>*> n_dot_conformal_factor_gradient,
254 : gsl::not_null<Scalar<DataVector>*>
255 : n_dot_lapse_times_conformal_factor_gradient,
256 : gsl::not_null<tnsr::I<DataVector, 3>*> n_dot_longitudinal_shift_excess,
257 : const tnsr::i<DataVector, 3>& face_normal,
258 : const tnsr::ij<DataVector, 3>& deriv_unnormalized_face_normal,
259 : const Scalar<DataVector>& face_normal_magnitude,
260 : const tnsr::I<DataVector, 3>& x,
261 : const Scalar<DataVector>& extrinsic_curvature_trace,
262 : const tnsr::I<DataVector, 3>& shift_background,
263 : const tnsr::II<DataVector, 3>& longitudinal_shift_background,
264 : const tnsr::II<DataVector, 3>& inv_conformal_metric,
265 : const tnsr::Ijj<DataVector, 3>& conformal_christoffel_second_kind) const;
266 :
267 0 : using argument_tags_linearized = tmpl::flatten<tmpl::list<
268 : ::Tags::Normalized<
269 : domain::Tags::UnnormalizedFaceNormal<3, Frame::Inertial>>,
270 : ::Tags::deriv<domain::Tags::UnnormalizedFaceNormal<3, Frame::Inertial>,
271 : tmpl::size_t<3>, Frame::Inertial>,
272 : ::Tags::Magnitude<
273 : domain::Tags::UnnormalizedFaceNormal<3, Frame::Inertial>>,
274 : domain::Tags::Coordinates<3, Frame::Inertial>,
275 : gr::Tags::TraceExtrinsicCurvature<DataVector>,
276 : Tags::LongitudinalShiftBackgroundMinusDtConformalMetric<DataVector, 3,
277 : Frame::Inertial>,
278 : Tags::ConformalFactorMinusOne<DataVector>,
279 : Tags::LapseTimesConformalFactorMinusOne<DataVector>,
280 : ::Tags::NormalDotFlux<Tags::ShiftExcess<DataVector, 3, Frame::Inertial>>,
281 : tmpl::conditional_t<ConformalGeometry == Xcts::Geometry::Curved,
282 : tmpl::list<Tags::InverseConformalMetric<
283 : DataVector, 3, Frame::Inertial>,
284 : Tags::ConformalChristoffelSecondKind<
285 : DataVector, 3, Frame::Inertial>>,
286 : tmpl::list<>>>>;
287 0 : using volume_tags_linearized = tmpl::list<>;
288 :
289 0 : void apply_linearized(
290 : gsl::not_null<Scalar<DataVector>*> conformal_factor_correction,
291 : gsl::not_null<Scalar<DataVector>*>
292 : lapse_times_conformal_factor_correction,
293 : gsl::not_null<tnsr::I<DataVector, 3>*> shift_excess_correction,
294 : gsl::not_null<Scalar<DataVector>*>
295 : n_dot_conformal_factor_gradient_correction,
296 : gsl::not_null<Scalar<DataVector>*>
297 : n_dot_lapse_times_conformal_factor_gradient_correction,
298 : gsl::not_null<tnsr::I<DataVector, 3>*>
299 : n_dot_longitudinal_shift_excess_correction,
300 : const tnsr::i<DataVector, 3>& face_normal,
301 : const tnsr::ij<DataVector, 3>& deriv_unnormalized_face_normal,
302 : const Scalar<DataVector>& face_normal_magnitude,
303 : const tnsr::I<DataVector, 3>& x,
304 : const Scalar<DataVector>& extrinsic_curvature_trace,
305 : const tnsr::II<DataVector, 3>& longitudinal_shift_background,
306 : const Scalar<DataVector>& conformal_factor_minus_one,
307 : const Scalar<DataVector>& lapse_times_conformal_factor_minus_one,
308 : const tnsr::I<DataVector, 3>& n_dot_longitudinal_shift_excess) const;
309 :
310 0 : void apply_linearized(
311 : gsl::not_null<Scalar<DataVector>*> conformal_factor_correction,
312 : gsl::not_null<Scalar<DataVector>*>
313 : lapse_times_conformal_factor_correction,
314 : gsl::not_null<tnsr::I<DataVector, 3>*> shift_excess_correction,
315 : gsl::not_null<Scalar<DataVector>*>
316 : n_dot_conformal_factor_gradient_correction,
317 : gsl::not_null<Scalar<DataVector>*>
318 : n_dot_lapse_times_conformal_factor_gradient_correction,
319 : gsl::not_null<tnsr::I<DataVector, 3>*>
320 : n_dot_longitudinal_shift_excess_correction,
321 : const tnsr::i<DataVector, 3>& face_normal,
322 : const tnsr::ij<DataVector, 3>& deriv_unnormalized_face_normal,
323 : const Scalar<DataVector>& face_normal_magnitude,
324 : const tnsr::I<DataVector, 3>& x,
325 : const Scalar<DataVector>& extrinsic_curvature_trace,
326 : const tnsr::II<DataVector, 3>& longitudinal_shift_background,
327 : const Scalar<DataVector>& conformal_factor_minus_one,
328 : const Scalar<DataVector>& lapse_times_conformal_factor_minus_one,
329 : const tnsr::I<DataVector, 3>& n_dot_longitudinal_shift_excess,
330 : const tnsr::II<DataVector, 3>& inv_conformal_metric,
331 : const tnsr::Ijj<DataVector, 3>& conformal_christoffel_second_kind) const;
332 :
333 : // NOLINTNEXTLINE(google-runtime-references)
334 0 : void pup(PUP::er& p) override;
335 :
336 : private:
337 0 : std::array<double, 3> center_ =
338 : make_array<3>(std::numeric_limits<double>::signaling_NaN());
339 0 : std::array<double, 3> rotation_ =
340 : make_array<3>(std::numeric_limits<double>::signaling_NaN());
341 : std::optional<std::unique_ptr<elliptic::analytic_data::AnalyticSolution>>
342 0 : solution_for_lapse_{};
343 : std::optional<std::unique_ptr<elliptic::analytic_data::AnalyticSolution>>
344 0 : solution_for_negative_expansion_{};
345 : };
346 :
347 : } // namespace Xcts::BoundaryConditions
|