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 : deserialize<
175 : std::unique_ptr<elliptic::analytic_data::InitialGuess>>(
176 : serialize(solution_for_lapse_.value()).data()))
177 : : std::nullopt,
178 : solution_for_negative_expansion_.has_value()
179 : ? std::make_optional(deserialize<std::unique_ptr<
180 : elliptic::analytic_data::InitialGuess>>(
181 : serialize(solution_for_negative_expansion_.value()).data()))
182 : : std::nullopt);
183 : }
184 :
185 0 : ApparentHorizon(
186 : std::array<double, 3> center, std::array<double, 3> rotation,
187 : std::optional<std::unique_ptr<elliptic::analytic_data::InitialGuess>>
188 : solution_for_lapse,
189 : std::optional<std::unique_ptr<elliptic::analytic_data::InitialGuess>>
190 : solution_for_negative_expansion,
191 : const Options::Context& context = {});
192 :
193 0 : const std::array<double, 3>& center() const { return center_; }
194 0 : const std::array<double, 3>& rotation() const { return rotation_; }
195 : const std::optional<std::unique_ptr<elliptic::analytic_data::InitialGuess>>&
196 0 : solution_for_lapse() const {
197 : return solution_for_lapse_;
198 : }
199 : const std::optional<std::unique_ptr<elliptic::analytic_data::InitialGuess>>&
200 0 : solution_for_negative_expansion() const {
201 : return solution_for_negative_expansion_;
202 : }
203 :
204 0 : std::vector<elliptic::BoundaryConditionType> boundary_condition_types()
205 : const override {
206 : return {// Conformal factor
207 : elliptic::BoundaryConditionType::Neumann,
208 : // Lapse times conformal factor
209 : this->solution_for_lapse_.has_value()
210 : ? elliptic::BoundaryConditionType::Dirichlet
211 : : elliptic::BoundaryConditionType::Neumann,
212 : // Shift
213 : elliptic::BoundaryConditionType::Dirichlet,
214 : elliptic::BoundaryConditionType::Dirichlet,
215 : elliptic::BoundaryConditionType::Dirichlet};
216 : }
217 :
218 0 : using argument_tags = tmpl::flatten<tmpl::list<
219 : domain::Tags::FaceNormal<3>,
220 : ::Tags::deriv<domain::Tags::UnnormalizedFaceNormal<3>, tmpl::size_t<3>,
221 : Frame::Inertial>,
222 : domain::Tags::UnnormalizedFaceNormalMagnitude<3>,
223 : domain::Tags::Coordinates<3, Frame::Inertial>,
224 : gr::Tags::TraceExtrinsicCurvature<DataVector>,
225 : Tags::ShiftBackground<DataVector, 3, Frame::Inertial>,
226 : Tags::LongitudinalShiftBackgroundMinusDtConformalMetric<DataVector, 3,
227 : Frame::Inertial>,
228 : tmpl::conditional_t<ConformalGeometry == Xcts::Geometry::Curved,
229 : tmpl::list<Tags::InverseConformalMetric<
230 : DataVector, 3, Frame::Inertial>,
231 : Tags::ConformalChristoffelSecondKind<
232 : DataVector, 3, Frame::Inertial>>,
233 : tmpl::list<>>>>;
234 0 : using volume_tags = tmpl::list<>;
235 :
236 0 : void apply(
237 : gsl::not_null<Scalar<DataVector>*> conformal_factor_minus_one,
238 : gsl::not_null<Scalar<DataVector>*> lapse_times_conformal_factor_minus_one,
239 : gsl::not_null<tnsr::I<DataVector, 3>*> shift_excess,
240 : gsl::not_null<Scalar<DataVector>*> n_dot_conformal_factor_gradient,
241 : gsl::not_null<Scalar<DataVector>*>
242 : n_dot_lapse_times_conformal_factor_gradient,
243 : gsl::not_null<tnsr::I<DataVector, 3>*> n_dot_longitudinal_shift_excess,
244 : const tnsr::i<DataVector, 3>& deriv_conformal_factor,
245 : const tnsr::i<DataVector, 3>& deriv_lapse_times_conformal_factor,
246 : const tnsr::iJ<DataVector, 3>& deriv_shift_excess,
247 : const tnsr::i<DataVector, 3>& face_normal,
248 : const tnsr::ij<DataVector, 3>& deriv_unnormalized_face_normal,
249 : const Scalar<DataVector>& face_normal_magnitude,
250 : const tnsr::I<DataVector, 3>& x,
251 : const Scalar<DataVector>& extrinsic_curvature_trace,
252 : const tnsr::I<DataVector, 3>& shift_background,
253 : const tnsr::II<DataVector, 3>& longitudinal_shift_background) const;
254 :
255 0 : void apply(
256 : gsl::not_null<Scalar<DataVector>*> conformal_factor_minus_one,
257 : gsl::not_null<Scalar<DataVector>*> lapse_times_conformal_factor_minus_one,
258 : gsl::not_null<tnsr::I<DataVector, 3>*> shift_excess,
259 : gsl::not_null<Scalar<DataVector>*> n_dot_conformal_factor_gradient,
260 : gsl::not_null<Scalar<DataVector>*>
261 : n_dot_lapse_times_conformal_factor_gradient,
262 : gsl::not_null<tnsr::I<DataVector, 3>*> n_dot_longitudinal_shift_excess,
263 : const tnsr::i<DataVector, 3>& deriv_conformal_factor,
264 : const tnsr::i<DataVector, 3>& deriv_lapse_times_conformal_factor,
265 : const tnsr::iJ<DataVector, 3>& deriv_shift_excess,
266 : const tnsr::i<DataVector, 3>& face_normal,
267 : const tnsr::ij<DataVector, 3>& deriv_unnormalized_face_normal,
268 : const Scalar<DataVector>& face_normal_magnitude,
269 : const tnsr::I<DataVector, 3>& x,
270 : const Scalar<DataVector>& extrinsic_curvature_trace,
271 : const tnsr::I<DataVector, 3>& shift_background,
272 : const tnsr::II<DataVector, 3>& longitudinal_shift_background,
273 : const tnsr::II<DataVector, 3>& inv_conformal_metric,
274 : const tnsr::Ijj<DataVector, 3>& conformal_christoffel_second_kind) const;
275 :
276 0 : using argument_tags_linearized = tmpl::flatten<tmpl::list<
277 : ::Tags::Normalized<
278 : domain::Tags::UnnormalizedFaceNormal<3, Frame::Inertial>>,
279 : ::Tags::deriv<domain::Tags::UnnormalizedFaceNormal<3, Frame::Inertial>,
280 : tmpl::size_t<3>, Frame::Inertial>,
281 : ::Tags::Magnitude<
282 : domain::Tags::UnnormalizedFaceNormal<3, Frame::Inertial>>,
283 : domain::Tags::Coordinates<3, Frame::Inertial>,
284 : gr::Tags::TraceExtrinsicCurvature<DataVector>,
285 : Tags::LongitudinalShiftBackgroundMinusDtConformalMetric<DataVector, 3,
286 : Frame::Inertial>,
287 : Tags::ConformalFactorMinusOne<DataVector>,
288 : Tags::LapseTimesConformalFactorMinusOne<DataVector>,
289 : ::Tags::NormalDotFlux<Tags::ShiftExcess<DataVector, 3, Frame::Inertial>>,
290 : tmpl::conditional_t<ConformalGeometry == Xcts::Geometry::Curved,
291 : tmpl::list<Tags::InverseConformalMetric<
292 : DataVector, 3, Frame::Inertial>,
293 : Tags::ConformalChristoffelSecondKind<
294 : DataVector, 3, Frame::Inertial>>,
295 : tmpl::list<>>>>;
296 0 : using volume_tags_linearized = tmpl::list<>;
297 :
298 0 : void apply_linearized(
299 : gsl::not_null<Scalar<DataVector>*> conformal_factor_correction,
300 : gsl::not_null<Scalar<DataVector>*>
301 : lapse_times_conformal_factor_correction,
302 : gsl::not_null<tnsr::I<DataVector, 3>*> shift_excess_correction,
303 : gsl::not_null<Scalar<DataVector>*>
304 : n_dot_conformal_factor_gradient_correction,
305 : gsl::not_null<Scalar<DataVector>*>
306 : n_dot_lapse_times_conformal_factor_gradient_correction,
307 : gsl::not_null<tnsr::I<DataVector, 3>*>
308 : n_dot_longitudinal_shift_excess_correction,
309 : const tnsr::i<DataVector, 3>& deriv_conformal_factor_correction,
310 : const tnsr::i<DataVector, 3>&
311 : deriv_lapse_times_conformal_factor_correction,
312 : const tnsr::iJ<DataVector, 3>& deriv_shift_excess_correction,
313 : const tnsr::i<DataVector, 3>& face_normal,
314 : const tnsr::ij<DataVector, 3>& deriv_unnormalized_face_normal,
315 : const Scalar<DataVector>& face_normal_magnitude,
316 : const tnsr::I<DataVector, 3>& x,
317 : const Scalar<DataVector>& extrinsic_curvature_trace,
318 : const tnsr::II<DataVector, 3>& longitudinal_shift_background,
319 : const Scalar<DataVector>& conformal_factor_minus_one,
320 : const Scalar<DataVector>& lapse_times_conformal_factor_minus_one,
321 : const tnsr::I<DataVector, 3>& n_dot_longitudinal_shift_excess) const;
322 :
323 0 : void apply_linearized(
324 : gsl::not_null<Scalar<DataVector>*> conformal_factor_correction,
325 : gsl::not_null<Scalar<DataVector>*>
326 : lapse_times_conformal_factor_correction,
327 : gsl::not_null<tnsr::I<DataVector, 3>*> shift_excess_correction,
328 : gsl::not_null<Scalar<DataVector>*>
329 : n_dot_conformal_factor_gradient_correction,
330 : gsl::not_null<Scalar<DataVector>*>
331 : n_dot_lapse_times_conformal_factor_gradient_correction,
332 : gsl::not_null<tnsr::I<DataVector, 3>*>
333 : n_dot_longitudinal_shift_excess_correction,
334 : const tnsr::i<DataVector, 3>& deriv_conformal_factor_correction,
335 : const tnsr::i<DataVector, 3>&
336 : deriv_lapse_times_conformal_factor_correction,
337 : const tnsr::iJ<DataVector, 3>& deriv_shift_excess_correction,
338 : const tnsr::i<DataVector, 3>& face_normal,
339 : const tnsr::ij<DataVector, 3>& deriv_unnormalized_face_normal,
340 : const Scalar<DataVector>& face_normal_magnitude,
341 : const tnsr::I<DataVector, 3>& x,
342 : const Scalar<DataVector>& extrinsic_curvature_trace,
343 : const tnsr::II<DataVector, 3>& longitudinal_shift_background,
344 : const Scalar<DataVector>& conformal_factor_minus_one,
345 : const Scalar<DataVector>& lapse_times_conformal_factor_minus_one,
346 : const tnsr::I<DataVector, 3>& n_dot_longitudinal_shift_excess,
347 : const tnsr::II<DataVector, 3>& inv_conformal_metric,
348 : const tnsr::Ijj<DataVector, 3>& conformal_christoffel_second_kind) const;
349 :
350 : // NOLINTNEXTLINE(google-runtime-references)
351 0 : void pup(PUP::er& p) override;
352 :
353 : private:
354 0 : std::array<double, 3> center_ =
355 : make_array<3>(std::numeric_limits<double>::signaling_NaN());
356 0 : std::array<double, 3> rotation_ =
357 : make_array<3>(std::numeric_limits<double>::signaling_NaN());
358 : std::optional<std::unique_ptr<elliptic::analytic_data::InitialGuess>>
359 0 : solution_for_lapse_{};
360 : std::optional<std::unique_ptr<elliptic::analytic_data::InitialGuess>>
361 0 : solution_for_negative_expansion_{};
362 : };
363 :
364 : } // namespace Xcts::BoundaryConditions
|
