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