Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <array>
7 : #include <cstddef>
8 : #include <limits>
9 :
10 : #include "DataStructures/Tensor/TypeAliases.hpp"
11 : #include "Utilities/ContainerHelpers.hpp"
12 :
13 : /// \cond
14 : class DataVector;
15 : template <size_t VolumeDim, typename Fr>
16 : class MathFunction;
17 : namespace gsl {
18 : template <class T>
19 : class not_null;
20 : } // namespace gsl
21 : /// \endcond
22 :
23 : namespace gh::BoundaryConditions {
24 : /// \brief Detailed implementation of Bjorhus-type boundary corrections
25 1 : namespace Bjorhus {
26 : /*!
27 : * \brief Computes the expression needed to set boundary conditions on the time
28 : * derivative of the characteristic field \f$v^{g}_{ab}\f$
29 : *
30 : * \details In the Bjorhus scheme, the time derivatives of evolved variables are
31 : * characteristic projected. A constraint-preserving correction term is added
32 : * here to the resulting characteristic (time-derivative) field:
33 : *
34 : * \f{align}
35 : * \Delta \partial_t v^{g}_{ab} = \lambda_{g} n^i C_{iab}
36 : * \f}
37 : *
38 : * where \f$n^i\f$ is the local unit normal to the external boundary,
39 : * \f$C_{iab} = \partial_i g_{ab} - \Phi_{iab}\f$ is the three-index
40 : * constraint, and \f$\lambda_{g}\f$ is the characteristic speed of the field
41 : * \f$v^{g}_{ab}\f$.
42 : */
43 : template <size_t VolumeDim, typename DataType>
44 1 : void constraint_preserving_corrections_dt_v_psi(
45 : gsl::not_null<tnsr::aa<DataType, VolumeDim, Frame::Inertial>*> bc_dt_v_psi,
46 : const tnsr::I<DataType, VolumeDim, Frame::Inertial>&
47 : unit_interface_normal_vector,
48 : const tnsr::iaa<DataType, VolumeDim, Frame::Inertial>&
49 : three_index_constraint,
50 : const std::array<DataType, 4>& char_speeds);
51 :
52 : /*!
53 : * \brief Computes the expression needed to set boundary conditions on the time
54 : * derivative of the characteristic field \f$v^{0}_{iab}\f$
55 : *
56 : * \details In the Bjorhus scheme, the time derivatives of evolved variables are
57 : * characteristic projected. A constraint-preserving correction term is added
58 : * here to the resulting characteristic (time-derivative) field \f$v^0_{iab}\f$:
59 : *
60 : * \f{align}
61 : * \Delta \partial_t v^{0}_{iab} = \lambda_{0} n^j C_{jiab}
62 : * \f}
63 : *
64 : * where \f$n^i\f$ is the local unit normal to the external boundary,
65 : * \f$C_{ijab} = \partial_i\Phi_{jab} - \partial_j\Phi_{iab}\f$ is the
66 : * four-index constraint, and \f$\lambda_{0}\f$ is the characteristic speed of
67 : * the field \f$v^0_{iab}\f$.
68 : *
69 : * \note In 3D, only the non-zero and unique components of the four-index
70 : * constraint are stored, as \f$\hat{C}_{iab} = \epsilon_i^{jk} C_{jkab}\f$. In
71 : * 2D the input expected here is \f$\hat{C}_{0ab} = C_{01ab}, \hat{C}_{1ab} =
72 : * C_{10ab}\f$, and in 1D \f$\hat{C}_{0ab} = 0\f$.
73 : */
74 : template <size_t VolumeDim, typename DataType>
75 1 : void constraint_preserving_corrections_dt_v_zero(
76 : gsl::not_null<tnsr::iaa<DataType, VolumeDim, Frame::Inertial>*>
77 : bc_dt_v_zero,
78 : const tnsr::I<DataType, VolumeDim, Frame::Inertial>&
79 : unit_interface_normal_vector,
80 : const tnsr::iaa<DataType, VolumeDim, Frame::Inertial>&
81 : four_index_constraint,
82 : const std::array<DataType, 4>& char_speeds);
83 :
84 : /// @{
85 : /*!
86 : * \brief Computes the expression needed to set boundary conditions on the time
87 : * derivative of the characteristic field \f$v^{-}_{ab}\f$
88 : *
89 : * \details In the Bjorhus scheme, the time derivatives of evolved variables are
90 : * characteristic projected. Constraint-preserving correction terms
91 : * \f$T^{\mathrm C}_{ab}\f$, Sommerfeld condition terms for gauge degrees of
92 : * freedom \f$T^{\mathrm G}_{ab}\f$, and terms constraining physical degrees of
93 : * freedom \f$T^{\mathrm P}_{ab}\f$ are added here to the resulting
94 : * characteristic (time-derivative) field:
95 : *
96 : * \f{align}
97 : * \Delta \partial_t v^{-}_{ab} = T^{\mathrm C}_{ab} + T^{\mathrm G}_{ab}
98 : * + T^{\mathrm P}_{ab}.
99 : * \f}
100 : *
101 : * These terms are given by Eq. (64) of \cite Lindblom2005qh :
102 : *
103 : * \f{eqnarray}
104 : * T^{\mathrm C}_{ab} &=& \frac{1}{2}
105 : * \left(2 k^c k^d l_a l_b - k^c l_b P^d_a - k^c l_a P^d_b - k^d l_b P^c_a
106 : * - k^d l_a P^c_b + P^{cd} P_{ab}\right) \partial_t v^{-}_{cd}\\
107 : * &&+ \frac{1}{\sqrt{2}} \lambda_{-}c^{\hat{0}-}_c \left(l_a l_b k^c +
108 : * P_{ab} l^c - P^c_b l_a - P^c_a l_b \right) \nonumber
109 : * \f}
110 : *
111 : * where \f$l^a (l_a)\f$ is the outgoing null vector (one-form), \f$k^a (k_a)\f$
112 : * is the incoming null vector (one-form), \f$P_{ab}, P^{ab}, P^a_b\f$ are
113 : * spacetime projection operators defined in `transverse_projection_operator()`,
114 : * \f$\partial_t v^{-}_{ab}\f$ is the characteristic projected time derivative
115 : * of evolved variables (corresponding to the \f$v^{-}\f$ field), and
116 : * \f$c^{\hat{0}\pm}_a\f$ are characteristic modes of the constraint evolution
117 : * system:
118 : *
119 : * \f{align}\nonumber
120 : * c^{\hat{0}\pm}_a = F_a \mp n^k C_{ka},
121 : * \f}
122 : *
123 : * where \f$F_a\f$ is the generalized-harmonic (GH) F constraint [Eq. (43) of
124 : * \cite Lindblom2005qh], \f$C_{ka}\f$ is the GH 2-index constraint [Eq. (44)
125 : * of \cite Lindblom2005qh], and \f$n^k\f$ is the unit spatial normal to the
126 : * outer boundary. Boundary correction terms that prevent strong reflections of
127 : * gauge perturbations are given by Eq. (25) of \cite Rinne2007ui :
128 : *
129 : * \f{align}
130 : * T^{\mathrm G}_{ab} =
131 : * \left(k_a P^c_b l^d + k_b P^c_a l^d -
132 : * \left(k_a l_b k^c l^d + k_b l_a k_c l^d + k_a k_b l^c l^d
133 : * \right) \right) \left(\gamma_2 - \frac{1}{r}
134 : * \right) \partial_t v^{g}_{cd}
135 : * \f}
136 : *
137 : * where \f$r\f$ is the radial coordinate at the outer boundary, which is
138 : * assumed to be spherical, \f$\gamma_2\f$ is a GH constraint damping parameter,
139 : * and \f$\partial_t v^{g}_{ab}\f$ is the characteristic projected time
140 : * derivative of evolved variables (corresponding to the \f$v^{g}\f$ field).
141 : * Finally, we constrain physical degrees of freedom using corrections from
142 : * Eq. (68) of \cite Lindblom2005qh :
143 : *
144 : * \f{align}
145 : * T^{\mathrm P}_{ab} = \left( P^c_a P^d_b - \frac{1}{2} P_{ab} P^{cd} \right)
146 : * \left(\partial_t v^{-}_{cd} + \lambda_{-}
147 : * \left(U^{3-}_{cd} - \gamma_2 n^i C_{icd}\right)\right),
148 : * \f}
149 : *
150 : * where \f$C_{icd}\f$ is the GH 3-index constraint [c.f. Eq. (26) of
151 : * \cite Lindblom2005qh, see also `three_index_constraint()`], and
152 : *
153 : * \f{align}\nonumber
154 : * U^{3-}_{ab} = 2 P^i_a P^j_b U^{8-}_{ij}
155 : * \f}
156 : *
157 : * is the inward propagating characteristic mode of the Weyl tensor evolution
158 : * \f$U^{8-}_{ab}\f$ [c.f. Eq. 75 of \cite Kidder2004rw ] projected onto the
159 : * outer boundary using the spatial-spacetime projection operators \f$P^i_a\f$.
160 : * Note that (A) the covariant derivative of extrinsic curvature needed to
161 : * get the Weyl propagating modes is calculated substituting the evolved
162 : * variable \f$\Phi_{iab}\f$ for spatial derivatives of the spacetime metric;
163 : * and (B) the spatial Ricci tensor used in the same calculation is also
164 : * calculated using the same substituion, and also includes corrections
165 : * proportional to the GH 4-index constraint \f$C_{ijab}\f$ [c.f. Eq. (45) of
166 : * \cite Lindblom2005qh , see also `four_index_constraint()`]:
167 : *
168 : * \f{align}\nonumber
169 : * R_{ij} \rightarrow R_{ij} + C_{(iklj)} + \frac{1}{2} n^k q^a C_{(ikj)a},
170 : * \f}
171 : *
172 : * where \f$q^a\f$ is the future-directed spacetime normal vector.
173 : */
174 : template <size_t VolumeDim, typename DataType>
175 1 : void constraint_preserving_corrections_dt_v_minus(
176 : gsl::not_null<tnsr::aa<DataType, VolumeDim, Frame::Inertial>*>
177 : bc_dt_v_minus,
178 : const tnsr::a<DataType, VolumeDim, Frame::Inertial>& outgoing_null_one_form,
179 : const tnsr::A<DataType, VolumeDim, Frame::Inertial>& incoming_null_vector,
180 : const tnsr::A<DataType, VolumeDim, Frame::Inertial>& outgoing_null_vector,
181 : const tnsr::aa<DataType, VolumeDim, Frame::Inertial>& projection_ab,
182 : const tnsr::Ab<DataType, VolumeDim, Frame::Inertial>& projection_Ab,
183 : const tnsr::AA<DataType, VolumeDim, Frame::Inertial>& projection_AB,
184 : const tnsr::aa<DataType, VolumeDim, Frame::Inertial>&
185 : char_projected_rhs_dt_v_minus,
186 : const tnsr::a<DataType, VolumeDim, Frame::Inertial>&
187 : constraint_char_zero_plus,
188 : const tnsr::a<DataType, VolumeDim, Frame::Inertial>&
189 : constraint_char_zero_minus,
190 : const std::array<DataType, 4>& char_speeds);
191 :
192 : template <size_t VolumeDim, typename DataType>
193 1 : void constraint_preserving_gauge_corrections_dt_v_minus(
194 : gsl::not_null<tnsr::aa<DataType, VolumeDim, Frame::Inertial>*>
195 : bc_dt_v_minus,
196 : const Scalar<DataType>& gamma2,
197 : const tnsr::I<DataType, VolumeDim, Frame::Inertial>& inertial_coords,
198 : const tnsr::a<DataType, VolumeDim, Frame::Inertial>& incoming_null_one_form,
199 : const tnsr::a<DataType, VolumeDim, Frame::Inertial>& outgoing_null_one_form,
200 : const tnsr::A<DataType, VolumeDim, Frame::Inertial>& incoming_null_vector,
201 : const tnsr::A<DataType, VolumeDim, Frame::Inertial>& outgoing_null_vector,
202 : const tnsr::aa<DataType, VolumeDim, Frame::Inertial>& projection_ab,
203 : const tnsr::Ab<DataType, VolumeDim, Frame::Inertial>& projection_Ab,
204 : const tnsr::AA<DataType, VolumeDim, Frame::Inertial>& projection_AB,
205 : const tnsr::aa<DataType, VolumeDim, Frame::Inertial>&
206 : char_projected_rhs_dt_v_psi,
207 : const tnsr::aa<DataType, VolumeDim, Frame::Inertial>&
208 : char_projected_rhs_dt_v_minus,
209 : const tnsr::a<DataType, VolumeDim, Frame::Inertial>&
210 : constraint_char_zero_plus,
211 : const tnsr::a<DataType, VolumeDim, Frame::Inertial>&
212 : constraint_char_zero_minus,
213 : const std::array<DataType, 4>& char_speeds);
214 :
215 : template <size_t VolumeDim, typename DataType>
216 1 : void constraint_preserving_gauge_physical_corrections_dt_v_minus(
217 : gsl::not_null<tnsr::aa<DataType, VolumeDim, Frame::Inertial>*>
218 : bc_dt_v_minus,
219 : const Scalar<DataType>& gamma2,
220 : const tnsr::I<DataType, VolumeDim, Frame::Inertial>& inertial_coords,
221 : double time,
222 : const tnsr::i<DataType, VolumeDim, Frame::Inertial>&
223 : unit_interface_normal_one_form,
224 : const tnsr::I<DataType, VolumeDim, Frame::Inertial>&
225 : unit_interface_normal_vector,
226 : const tnsr::A<DataType, VolumeDim, Frame::Inertial>&
227 : spacetime_unit_normal_vector,
228 : const tnsr::a<DataType, VolumeDim, Frame::Inertial>& incoming_null_one_form,
229 : const tnsr::a<DataType, VolumeDim, Frame::Inertial>& outgoing_null_one_form,
230 : const tnsr::A<DataType, VolumeDim, Frame::Inertial>& incoming_null_vector,
231 : const tnsr::A<DataType, VolumeDim, Frame::Inertial>& outgoing_null_vector,
232 : const tnsr::aa<DataType, VolumeDim, Frame::Inertial>& projection_ab,
233 : const tnsr::Ab<DataType, VolumeDim, Frame::Inertial>& projection_Ab,
234 : const tnsr::AA<DataType, VolumeDim, Frame::Inertial>& projection_AB,
235 : const tnsr::II<DataType, VolumeDim, Frame::Inertial>&
236 : inverse_spatial_metric,
237 : const tnsr::ii<DataType, VolumeDim, Frame::Inertial>& extrinsic_curvature,
238 : const tnsr::aa<DataType, VolumeDim, Frame::Inertial>& spacetime_metric,
239 : const tnsr::AA<DataType, VolumeDim, Frame::Inertial>&
240 : inverse_spacetime_metric,
241 : const tnsr::iaa<DataType, VolumeDim, Frame::Inertial>&
242 : three_index_constraint,
243 : const tnsr::aa<DataType, VolumeDim, Frame::Inertial>&
244 : char_projected_rhs_dt_v_psi,
245 : const tnsr::aa<DataType, VolumeDim, Frame::Inertial>&
246 : char_projected_rhs_dt_v_minus,
247 : const tnsr::a<DataType, VolumeDim, Frame::Inertial>&
248 : constraint_char_zero_plus,
249 : const tnsr::a<DataType, VolumeDim, Frame::Inertial>&
250 : constraint_char_zero_minus,
251 : const tnsr::iaa<DataType, VolumeDim, Frame::Inertial>& phi,
252 : const tnsr::ijaa<DataType, VolumeDim, Frame::Inertial>& d_phi,
253 : const tnsr::iaa<DataType, VolumeDim, Frame::Inertial>& d_pi,
254 : const std::array<DataType, 4>& char_speeds,
255 : const MathFunction<1, Frame::Inertial>* incoming_wave_profile = nullptr);
256 : /// @}
257 :
258 : namespace detail {
259 : template <size_t VolumeDim, typename DataType>
260 : void add_gauge_sommerfeld_terms_to_dt_v_minus(
261 : const gsl::not_null<tnsr::aa<DataType, VolumeDim, Frame::Inertial>*>
262 : bc_dt_v_minus,
263 : const Scalar<DataType>& gamma2,
264 : const tnsr::I<DataType, VolumeDim, Frame::Inertial>& inertial_coords,
265 : const tnsr::a<DataType, VolumeDim, Frame::Inertial>& incoming_null_one_form,
266 : const tnsr::a<DataType, VolumeDim, Frame::Inertial>& outgoing_null_one_form,
267 : const tnsr::A<DataType, VolumeDim, Frame::Inertial>& incoming_null_vector,
268 : const tnsr::A<DataType, VolumeDim, Frame::Inertial>& outgoing_null_vector,
269 : const tnsr::Ab<DataType, VolumeDim, Frame::Inertial>& projection_Ab,
270 : const tnsr::aa<DataType, VolumeDim, Frame::Inertial>&
271 : char_projected_rhs_dt_v_psi);
272 :
273 : template <size_t VolumeDim, typename DataType>
274 : void add_constraint_dependent_terms_to_dt_v_minus(
275 : const gsl::not_null<tnsr::aa<DataType, VolumeDim, Frame::Inertial>*>
276 : bc_dt_v_minus,
277 : const tnsr::a<DataType, VolumeDim, Frame::Inertial>& outgoing_null_one_form,
278 : const tnsr::A<DataType, VolumeDim, Frame::Inertial>& incoming_null_vector,
279 : const tnsr::A<DataType, VolumeDim, Frame::Inertial>& outgoing_null_vector,
280 : const tnsr::aa<DataType, VolumeDim, Frame::Inertial>& projection_ab,
281 : const tnsr::Ab<DataType, VolumeDim, Frame::Inertial>& projection_Ab,
282 : const tnsr::AA<DataType, VolumeDim, Frame::Inertial>& projection_AB,
283 : const tnsr::a<DataType, VolumeDim, Frame::Inertial>&
284 : constraint_char_zero_plus,
285 : const tnsr::a<DataType, VolumeDim, Frame::Inertial>&
286 : constraint_char_zero_minus,
287 : const tnsr::aa<DataType, VolumeDim, Frame::Inertial>&
288 : char_projected_rhs_dt_v_minus,
289 : const std::array<DataType, 4>& char_speeds);
290 :
291 : template <size_t VolumeDim, typename DataType>
292 : void add_physical_terms_to_dt_v_minus(
293 : const gsl::not_null<tnsr::aa<DataType, VolumeDim, Frame::Inertial>*>
294 : bc_dt_v_minus,
295 : const Scalar<DataType>& gamma2,
296 : const tnsr::i<DataType, VolumeDim, Frame::Inertial>&
297 : unit_interface_normal_one_form,
298 : const tnsr::I<DataType, VolumeDim, Frame::Inertial>&
299 : unit_interface_normal_vector,
300 : const tnsr::A<DataType, VolumeDim, Frame::Inertial>&
301 : spacetime_unit_normal_vector,
302 : const tnsr::aa<DataType, VolumeDim, Frame::Inertial>& projection_ab,
303 : const tnsr::Ab<DataType, VolumeDim, Frame::Inertial>& projection_Ab,
304 : const tnsr::AA<DataType, VolumeDim, Frame::Inertial>& projection_AB,
305 : const tnsr::II<DataType, VolumeDim, Frame::Inertial>&
306 : inverse_spatial_metric,
307 : const tnsr::ii<DataType, VolumeDim, Frame::Inertial>& extrinsic_curvature,
308 : const tnsr::aa<DataType, VolumeDim, Frame::Inertial>& spacetime_metric,
309 : const tnsr::AA<DataType, VolumeDim, Frame::Inertial>&
310 : inverse_spacetime_metric,
311 : const tnsr::iaa<DataType, VolumeDim, Frame::Inertial>&
312 : three_index_constraint,
313 : const tnsr::aa<DataType, VolumeDim, Frame::Inertial>&
314 : char_projected_rhs_dt_v_minus,
315 : const tnsr::iaa<DataType, VolumeDim, Frame::Inertial>& phi,
316 : const tnsr::ijaa<DataType, VolumeDim, Frame::Inertial>& d_phi,
317 : const tnsr::iaa<DataType, VolumeDim, Frame::Inertial>& d_pi,
318 : const std::array<DataType, 4>& char_speeds,
319 : double time = std::numeric_limits<double>::signaling_NaN(),
320 : const MathFunction<1, Frame::Inertial>* incoming_wave_profile = nullptr);
321 : } // namespace detail
322 : } // namespace Bjorhus
323 : } // namespace gh::BoundaryConditions
|
