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