Line data Source code
1 1 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : ///\file
5 : /// Defines functions to calculate the generalized harmonic constraints
6 :
7 : #pragma once
8 :
9 : #include <cstddef>
10 :
11 : #include "DataStructures/DataBox/Prefixes.hpp"
12 : #include "DataStructures/DataBox/Tag.hpp"
13 : #include "DataStructures/DataVector.hpp"
14 : #include "DataStructures/Tensor/TypeAliases.hpp"
15 : #include "Evolution/Systems/GeneralizedHarmonic/Tags.hpp"
16 : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
17 : #include "PointwiseFunctions/GeneralRelativity/GeneralizedHarmonic/ConstraintDampingTags.hpp"
18 : #include "Utilities/SetNumberOfGridPoints.hpp"
19 : #include "Utilities/TMPL.hpp"
20 :
21 : /// \cond
22 : namespace gsl {
23 : template <typename T>
24 : class not_null;
25 : } // namespace gsl
26 : /// \endcond
27 :
28 : namespace gh {
29 : /// @{
30 : /*!
31 : * \brief Computes the generalized-harmonic 3-index constraint.
32 : *
33 : * \details Computes the generalized-harmonic 3-index constraint,
34 : * \f$C_{iab} = \partial_i g_{ab} - \Phi_{iab},\f$ which is
35 : * given by Eq. (26) of \cite Lindblom2005qh
36 : */
37 : template <typename DataType, size_t SpatialDim, typename Frame>
38 1 : tnsr::iaa<DataType, SpatialDim, Frame> three_index_constraint(
39 : const tnsr::iaa<DataType, SpatialDim, Frame>& d_spacetime_metric,
40 : const tnsr::iaa<DataType, SpatialDim, Frame>& phi);
41 :
42 : template <typename DataType, size_t SpatialDim, typename Frame>
43 1 : void three_index_constraint(
44 : gsl::not_null<tnsr::iaa<DataType, SpatialDim, Frame>*> constraint,
45 : const tnsr::iaa<DataType, SpatialDim, Frame>& d_spacetime_metric,
46 : const tnsr::iaa<DataType, SpatialDim, Frame>& phi);
47 : /// @}
48 :
49 : /// @{
50 : /*!
51 : * \brief Computes the generalized-harmonic gauge constraint.
52 : *
53 : * \details Computes the generalized-harmonic gauge constraint
54 : * [Eq. (40) of \cite Lindblom2005qh],
55 : * \f[
56 : * C_a = H_a + \gamma^{ij} \Phi_{ija} + n^b \Pi_{ba}
57 : * - \frac{1}{2} \gamma^i_a g^{bc} \Phi_{ibc}
58 : * - \frac{1}{2} n_a g^{bc} \Pi_{bc},
59 : * \f]
60 : * where \f$H_a\f$ is the gauge function,
61 : * \f$g_{ab}\f$ is the spacetime metric,
62 : * \f$\Pi_{ab}=-n^c\partial_c g_{ab}\f$, and
63 : * \f$\Phi_{iab} = \partial_i g_{ab}\f$; \f$n^a\f$ is the timelike unit
64 : * normal vector to the spatial slice, \f$\gamma^{ij}\f$ is the inverse spatial
65 : * metric, and \f$\gamma^b_c = \delta^b_c + n^b n_c\f$.
66 : */
67 : template <typename DataType, size_t SpatialDim, typename Frame>
68 1 : tnsr::a<DataType, SpatialDim, Frame> gauge_constraint(
69 : const tnsr::a<DataType, SpatialDim, Frame>& gauge_function,
70 : const tnsr::a<DataType, SpatialDim, Frame>& spacetime_normal_one_form,
71 : const tnsr::A<DataType, SpatialDim, Frame>& spacetime_normal_vector,
72 : const tnsr::II<DataType, SpatialDim, Frame>& inverse_spatial_metric,
73 : const tnsr::AA<DataType, SpatialDim, Frame>& inverse_spacetime_metric,
74 : const tnsr::aa<DataType, SpatialDim, Frame>& pi,
75 : const tnsr::iaa<DataType, SpatialDim, Frame>& phi);
76 :
77 : template <typename DataType, size_t SpatialDim, typename Frame>
78 1 : void gauge_constraint(
79 : gsl::not_null<tnsr::a<DataType, SpatialDim, Frame>*> constraint,
80 : const tnsr::a<DataType, SpatialDim, Frame>& gauge_function,
81 : const tnsr::a<DataType, SpatialDim, Frame>& spacetime_normal_one_form,
82 : const tnsr::A<DataType, SpatialDim, Frame>& spacetime_normal_vector,
83 : const tnsr::II<DataType, SpatialDim, Frame>& inverse_spatial_metric,
84 : const tnsr::AA<DataType, SpatialDim, Frame>& inverse_spacetime_metric,
85 : const tnsr::aa<DataType, SpatialDim, Frame>& pi,
86 : const tnsr::iaa<DataType, SpatialDim, Frame>& phi);
87 : /// @}
88 :
89 : /// @{
90 : /*!
91 : * \brief Computes the generalized-harmonic 2-index constraint.
92 : *
93 : * \details Computes the generalized-harmonic 2-index constraint
94 : * [Eq. (44) of \cite Lindblom2005qh],
95 : * \f{eqnarray}{
96 : * C_{ia} &\equiv& \gamma^{jk}\partial_j \Phi_{ika}
97 : * - \frac{1}{2} \gamma_a^j g^{cd}\partial_j \Phi_{icd}
98 : * + n^b \partial_i \Pi_{ba}
99 : * - \frac{1}{2} n_a g^{cd}\partial_i\Pi_{cd}
100 : * \nonumber\\&&
101 : * + \partial_i H_a
102 : * + \frac{1}{2} \gamma_a^j \Phi_{jcd} \Phi_{ief}
103 : * g^{ce}g^{df}
104 : * + \frac{1}{2} \gamma^{jk} \Phi_{jcd} \Phi_{ike}
105 : * g^{cd}n^e n_a
106 : * \nonumber\\&&
107 : * - \gamma^{jk}\gamma^{mn}\Phi_{jma}\Phi_{ikn}
108 : * + \frac{1}{2} \Phi_{icd} \Pi_{be} n_a
109 : * \left(g^{cb}g^{de}
110 : * +\frac{1}{2}g^{be} n^c n^d\right)
111 : * \nonumber\\&&
112 : * - \Phi_{icd} \Pi_{ba} n^c \left(g^{bd}
113 : * +\frac{1}{2} n^b n^d\right)
114 : * + \frac{1}{2} \gamma_2 \left(n_a g^{cd}
115 : * - 2 \delta^c_a n^d\right) C_{icd}.
116 : * \f}
117 : * where \f$H_a\f$ is the gauge function,
118 : * \f$g_{ab}\f$ is the spacetime metric,
119 : * \f$\Pi_{ab}=-n^c\partial_c g_{ab}\f$, and
120 : * \f$\Phi_{iab} = \partial_i g_{ab}\f$; \f$n^a\f$ is the timelike unit
121 : * normal vector to the spatial slice, \f$\gamma^{ij}\f$ is the inverse spatial
122 : * metric, and \f$\gamma^b_c = \delta^b_c + n^b n_c\f$.
123 : */
124 : template <typename DataType, size_t SpatialDim, typename Frame>
125 1 : tnsr::ia<DataType, SpatialDim, Frame> two_index_constraint(
126 : const tnsr::ab<DataType, SpatialDim, Frame>& spacetime_d_gauge_function,
127 : const tnsr::a<DataType, SpatialDim, Frame>& spacetime_normal_one_form,
128 : const tnsr::A<DataType, SpatialDim, Frame>& spacetime_normal_vector,
129 : const tnsr::II<DataType, SpatialDim, Frame>& inverse_spatial_metric,
130 : const tnsr::AA<DataType, SpatialDim, Frame>& inverse_spacetime_metric,
131 : const tnsr::aa<DataType, SpatialDim, Frame>& pi,
132 : const tnsr::iaa<DataType, SpatialDim, Frame>& phi,
133 : const tnsr::iaa<DataType, SpatialDim, Frame>& d_pi,
134 : const tnsr::ijaa<DataType, SpatialDim, Frame>& d_phi,
135 : const Scalar<DataType>& gamma2,
136 : const tnsr::iaa<DataType, SpatialDim, Frame>& three_index_constraint);
137 :
138 : template <typename DataType, size_t SpatialDim, typename Frame>
139 1 : void two_index_constraint(
140 : gsl::not_null<tnsr::ia<DataType, SpatialDim, Frame>*> constraint,
141 : const tnsr::ab<DataType, SpatialDim, Frame>& spacetime_d_gauge_function,
142 : const tnsr::a<DataType, SpatialDim, Frame>& spacetime_normal_one_form,
143 : const tnsr::A<DataType, SpatialDim, Frame>& spacetime_normal_vector,
144 : const tnsr::II<DataType, SpatialDim, Frame>& inverse_spatial_metric,
145 : const tnsr::AA<DataType, SpatialDim, Frame>& inverse_spacetime_metric,
146 : const tnsr::aa<DataType, SpatialDim, Frame>& pi,
147 : const tnsr::iaa<DataType, SpatialDim, Frame>& phi,
148 : const tnsr::iaa<DataType, SpatialDim, Frame>& d_pi,
149 : const tnsr::ijaa<DataType, SpatialDim, Frame>& d_phi,
150 : const Scalar<DataType>& gamma2,
151 : const tnsr::iaa<DataType, SpatialDim, Frame>& three_index_constraint);
152 : /// @}
153 :
154 : /// @{
155 : /*!
156 : * \brief Computes the generalized-harmonic 4-index constraint.
157 : *
158 : * \details Computes the independent components of the generalized-harmonic
159 : * 4-index constraint. The constraint itself is given by Eq. (45) of
160 : * \cite Lindblom2005qh,
161 : * \f{eqnarray}{
162 : * C_{ijab} = 2 \partial_{[i}\Phi_{j]ab},
163 : * \f}
164 : * where \f$\Phi_{iab} = \partial_i g_{ab}\f$. Because the constraint is
165 : * antisymmetric on the two spatial indices, here we compute and store
166 : * only the independent components of \f$C_{ijab}\f$. Specifically, we
167 : * compute
168 : * \f{eqnarray}{
169 : * D_{iab} \equiv \frac{1}{2} \epsilon_{i}{}^{jk} C_{jkab}
170 : * = \epsilon_{i}{}^{jk} \partial_j \Phi_{kab},
171 : * \f}
172 : * where \f$\epsilon_{ijk}\f$ is the flat-space Levi-Civita symbol,
173 : * which is raised and lowered with the Kronecker delta.
174 : * In terms
175 : * of \f$D_{iab}\f$, the full 4-index constraint is
176 : * \f{eqnarray}{
177 : * C_{jkab} = \epsilon^{i}{}_{jk} D_{iab}.
178 : * \f}
179 : */
180 : template <typename DataType, size_t SpatialDim, typename Frame>
181 1 : tnsr::iaa<DataType, SpatialDim, Frame> four_index_constraint(
182 : const tnsr::ijaa<DataType, SpatialDim, Frame>& d_phi);
183 :
184 : template <typename DataType, size_t SpatialDim, typename Frame>
185 1 : void four_index_constraint(
186 : gsl::not_null<tnsr::iaa<DataType, SpatialDim, Frame>*> constraint,
187 : const tnsr::ijaa<DataType, SpatialDim, Frame>& d_phi);
188 : /// @}
189 :
190 : /// @{
191 : /*!
192 : * \brief Computes the generalized-harmonic F constraint.
193 : *
194 : * \details Computes the generalized-harmonic F constraint
195 : * [Eq. (43) of \cite Lindblom2005qh],
196 : * \f{eqnarray}{
197 : * {\cal F}_a &\equiv&
198 : * \frac{1}{2} \gamma_a^i g^{bc}\partial_i \Pi_{bc}
199 : * - \gamma^{ij} \partial_i \Pi_{ja}
200 : * - \gamma^{ij} n^b \partial_i \Phi_{jba}
201 : * + \frac{1}{2} n_a g^{bc} \gamma^{ij} \partial_i \Phi_{jbc}
202 : * \nonumber \\ &&
203 : * + n_a \gamma^{ij} \partial_i H_j
204 : * + \gamma_a^i \Phi_{ijb} \gamma^{jk}\Phi_{kcd} g^{bd} n^c
205 : * - \frac{1}{2} \gamma_a^i \Phi_{ijb} \gamma^{jk}
206 : * \Phi_{kcd} g^{cd} n^b
207 : * \nonumber \\ &&
208 : * - \gamma_a^i n^b \partial_i H_b
209 : * + \gamma^{ij} \Phi_{icd} \Phi_{jba} g^{bc} n^d
210 : * - \frac{1}{2} n_a \gamma^{ij} \gamma^{mn} \Phi_{imc} \Phi_{njd}g^{cd}
211 : * \nonumber \\ &&
212 : * - \frac{1}{4} n_a \gamma^{ij}\Phi_{icd}\Phi_{jbe}
213 : * g^{cb}g^{de}
214 : * + \frac{1}{4} n_a \Pi_{cd} \Pi_{be}
215 : * g^{cb}g^{de}
216 : * - \gamma^{ij} H_i \Pi_{ja}
217 : * \nonumber \\ &&
218 : * - n^b \gamma^{ij} \Pi_{b i} \Pi_{ja}
219 : * - \frac{1}{4} \gamma_a^i \Phi_{icd} n^c n^d \Pi_{be}
220 : * g^{be}
221 : * + \frac{1}{2} n_a \Pi_{cd} \Pi_{be}g^{ce}
222 : * n^d n^b
223 : * \nonumber \\ &&
224 : * + \gamma_a^i \Phi_{icd} \Pi_{be} n^c n^b g^{de}
225 : * - \gamma^{ij}\Phi_{iba} n^b \Pi_{je} n^e
226 : * - \frac{1}{2} \gamma^{ij}\Phi_{icd} n^c n^d \Pi_{ja}
227 : * \nonumber \\ &&
228 : * - \gamma^{ij} H_i \Phi_{jba} n^b
229 : * + \gamma_{a}^i \Phi_{icd} H_b g^{bc} n^d
230 : * +\gamma_2\bigl(\gamma^{id}{\cal C}_{ida}
231 : * -\frac{1}{2} \gamma_a^i g^{cd}{\cal C}_{icd}\bigr)
232 : * \nonumber \\ &&
233 : * + \frac{1}{2} n_a \Pi_{cd}g^{cd} H_b n^b
234 : * - n_a \gamma^{ij} \Phi_{ijc} H_d g^{cd}
235 : * +\frac{1}{2} n_a \gamma^{ij} H_i \Phi_{jcd}g^{cd}
236 : * \nonumber \\ &&
237 : * - 16 \pi n^a T_{a b}
238 : * \f}
239 : * where \f$H_a\f$ is the gauge function,
240 : * \f$g_{ab}\f$ is the spacetime metric,
241 : * \f$\Pi_{ab}=-n^c\partial_c g_{ab}\f$, and
242 : * \f$\Phi_{iab} = \partial_i g_{ab}\f$; \f$n^a\f$ is the timelike unit
243 : * normal vector to the spatial slice, \f$\gamma^{ij}\f$ is the inverse spatial
244 : * metric, \f$\gamma^b_c = \delta^b_c + n^b n_c\f$, and \f$T_{a b}\f$ is the
245 : * stress-energy tensor if nonzero (if using the overload with no stress-energy
246 : * tensor provided, the stress energy term is omitted).
247 : *
248 : * To justify the stress-energy contribution to the F constraint, note that
249 : * the stress-energy tensor appears in the dynamics of the Generalized
250 : * Harmonic system only through \f$\partial_t \Pi_{a b}\f$.
251 : * That dependence arises from (using \f$\dots\f$ to indicate collections of
252 : * terms that are known to be independent of the stress-energy tensor):
253 : *
254 : * \f[
255 : * {\cal F}_a = \dots \alpha^{-1}(\partial_t {\cal C}_a),
256 : * \f]
257 : *
258 : * where
259 : *
260 : * \f[
261 : * {\cal C}_a = H_a + \gamma^{i j} \Phi_{ij a} + n^b \Pi_{ba}
262 : * - \frac{1}{2} \gamma_a{}^i g^{bc} \Phi_{i b c}
263 : * - \frac{1}{2} n_a g^{bc} \Pi_{b c}.
264 : * \f].
265 : *
266 : * Therefore, the Stress-energy contribution can be calculated from the
267 : * trace-reversed contribution appearing in
268 : * `grmhd::GhValenciaDivClean::add_stress_energy_term_to_dt_pi` -- the
269 : * trace reversal in that function and the trace-reversal that appears
270 : * explicitly in \f${\cal C}_a\f$ cancel.
271 : */
272 : template <typename DataType, size_t SpatialDim, typename Frame>
273 1 : tnsr::a<DataType, SpatialDim, Frame> f_constraint(
274 : const tnsr::a<DataType, SpatialDim, Frame>& gauge_function,
275 : const tnsr::ab<DataType, SpatialDim, Frame>& spacetime_d_gauge_function,
276 : const tnsr::a<DataType, SpatialDim, Frame>& spacetime_normal_one_form,
277 : const tnsr::A<DataType, SpatialDim, Frame>& spacetime_normal_vector,
278 : const tnsr::II<DataType, SpatialDim, Frame>& inverse_spatial_metric,
279 : const tnsr::AA<DataType, SpatialDim, Frame>& inverse_spacetime_metric,
280 : const tnsr::aa<DataType, SpatialDim, Frame>& pi,
281 : const tnsr::iaa<DataType, SpatialDim, Frame>& phi,
282 : const tnsr::iaa<DataType, SpatialDim, Frame>& d_pi,
283 : const tnsr::ijaa<DataType, SpatialDim, Frame>& d_phi,
284 : const Scalar<DataType>& gamma2,
285 : const tnsr::iaa<DataType, SpatialDim, Frame>& three_index_constraint);
286 :
287 : template <typename DataType, size_t SpatialDim, typename Frame>
288 1 : void f_constraint(
289 : gsl::not_null<tnsr::a<DataType, SpatialDim, Frame>*> constraint,
290 : const tnsr::a<DataType, SpatialDim, Frame>& gauge_function,
291 : const tnsr::ab<DataType, SpatialDim, Frame>& spacetime_d_gauge_function,
292 : const tnsr::a<DataType, SpatialDim, Frame>& spacetime_normal_one_form,
293 : const tnsr::A<DataType, SpatialDim, Frame>& spacetime_normal_vector,
294 : const tnsr::II<DataType, SpatialDim, Frame>& inverse_spatial_metric,
295 : const tnsr::AA<DataType, SpatialDim, Frame>& inverse_spacetime_metric,
296 : const tnsr::aa<DataType, SpatialDim, Frame>& pi,
297 : const tnsr::iaa<DataType, SpatialDim, Frame>& phi,
298 : const tnsr::iaa<DataType, SpatialDim, Frame>& d_pi,
299 : const tnsr::ijaa<DataType, SpatialDim, Frame>& d_phi,
300 : const Scalar<DataType>& gamma2,
301 : const tnsr::iaa<DataType, SpatialDim, Frame>& three_index_constraint);
302 :
303 : template <typename DataType, size_t SpatialDim, typename Frame>
304 1 : tnsr::a<DataType, SpatialDim, Frame> f_constraint(
305 : const tnsr::a<DataType, SpatialDim, Frame>& gauge_function,
306 : const tnsr::ab<DataType, SpatialDim, Frame>& spacetime_d_gauge_function,
307 : const tnsr::a<DataType, SpatialDim, Frame>& spacetime_normal_one_form,
308 : const tnsr::A<DataType, SpatialDim, Frame>& spacetime_normal_vector,
309 : const tnsr::II<DataType, SpatialDim, Frame>& inverse_spatial_metric,
310 : const tnsr::AA<DataType, SpatialDim, Frame>& inverse_spacetime_metric,
311 : const tnsr::aa<DataType, SpatialDim, Frame>& pi,
312 : const tnsr::iaa<DataType, SpatialDim, Frame>& phi,
313 : const tnsr::iaa<DataType, SpatialDim, Frame>& d_pi,
314 : const tnsr::ijaa<DataType, SpatialDim, Frame>& d_phi,
315 : const Scalar<DataType>& gamma2,
316 : const tnsr::iaa<DataType, SpatialDim, Frame>& three_index_constraint,
317 : const tnsr::aa<DataType, SpatialDim, Frame>& trace_reversed_stress_energy);
318 :
319 : template <typename DataType, size_t SpatialDim, typename Frame>
320 1 : void f_constraint(
321 : gsl::not_null<tnsr::a<DataType, SpatialDim, Frame>*> constraint,
322 : const tnsr::a<DataType, SpatialDim, Frame>& gauge_function,
323 : const tnsr::ab<DataType, SpatialDim, Frame>& spacetime_d_gauge_function,
324 : const tnsr::a<DataType, SpatialDim, Frame>& spacetime_normal_one_form,
325 : const tnsr::A<DataType, SpatialDim, Frame>& spacetime_normal_vector,
326 : const tnsr::II<DataType, SpatialDim, Frame>& inverse_spatial_metric,
327 : const tnsr::AA<DataType, SpatialDim, Frame>& inverse_spacetime_metric,
328 : const tnsr::aa<DataType, SpatialDim, Frame>& pi,
329 : const tnsr::iaa<DataType, SpatialDim, Frame>& phi,
330 : const tnsr::iaa<DataType, SpatialDim, Frame>& d_pi,
331 : const tnsr::ijaa<DataType, SpatialDim, Frame>& d_phi,
332 : const Scalar<DataType>& gamma2,
333 : const tnsr::iaa<DataType, SpatialDim, Frame>& three_index_constraint,
334 : const tnsr::aa<DataType, SpatialDim, Frame>& trace_reversed_stress_energy);
335 : /// @}
336 :
337 : /// @{
338 : /*!
339 : * \brief Computes the generalized-harmonic (unnormalized) constraint energy.
340 : *
341 : * \details Computes the generalized-harmonic unnormalized constraint energy
342 : * [Eq. (53) of \cite Lindblom2005qh with \f$m^{ab}=\delta^{ab}\f$ and with each
343 : * term in the sum scaled by an arbitrary coefficient],
344 : * \f{eqnarray}{
345 : * E & = & K_1 C_a C_a + K_2\left(F_a F_a
346 : + C_{ia} C_{ja} \gamma^{ij}\right)\nonumber\\
347 : & + & K_3 C_{iab} C_{jab} \gamma^{ij} + K_4 C_{ikab} C_{jlab}\gamma^{ij}
348 : \gamma^{kl}.
349 : * \f}
350 : * Here \f$C_a\f$ is the gauge constraint, \f$F_a\f$ is the f constraint,
351 : * \f$C_{ia}\f$ is the two-index constraint, \f$C_{iab}\f$ is the
352 : * three-index constraint, \f$C_{ikab}\f$ is the four-index constraint,
353 : * \f$\gamma^{ij}\f$ is the inverse spatial metric, and
354 : * \f$K_1\f$, \f$K_2\f$, \f$K_3\f$, and \f$K_4\f$ are constant multipliers
355 : * for each term that each default to a value of 1.0. Note that in this
356 : * equation, spacetime indices \f$a,b\f$ are raised and lowered with
357 : * the Kronecker delta.
358 : *
359 : * Also note that the argument `four_index_constraint` is a rank-3 tensor.
360 : * This is because `gh::four_index_constraint()` takes
361 : * advantage of the antisymmetry of the four-index constraint's first two
362 : * indices to only compute and return the independent
363 : * components of \f$C_{ijab}\f$, which can be written as
364 : * \f{eqnarray}{
365 : * D_{iab} \equiv \frac{1}{2} \epsilon_{i}{}^{jk} C_{jkab},
366 : * \f} where \f$\epsilon_{ijk}\f$ is the flat-space Levi-Civita tensor,
367 : * whose inidces are raised and lowered with the Kronecker delta.
368 : * The result is
369 : * \f{eqnarray}{
370 : * E & = & K_1 C_a C_a + K_2\left(F_a F_a
371 : * + C_{ia} C_{ja} \gamma^{ij}\right)\nonumber\\
372 : * & + & K_3 C_{iab} C_{jab} \gamma^{ij}\nonumber\\
373 : & + & 2 K_4 \gamma D_{iab} D_{jab} \gamma^{ij},
374 : * \f} where \f$\gamma\f$ is the determinant of the spatial metric.
375 : *
376 : * To derive this expression for the constraint energy implemented here,
377 : * Eq.~(53) of \cite Lindblom2005qh is
378 : * \f{eqnarray}{
379 : * S_{AB} dc^A dc^B &=&
380 : * m^{ab}\Bigl[d F_ad F_b
381 : * +\gamma^{ij}\bigl(d C_{ia}d C_{jb}
382 : * +\gamma^{kl}m^{cd}d C_{ikac}d C_{jlbd}\bigr)
383 : * \nonumber\\
384 : * & + & \Lambda^2\bigl(d C_ad C_b
385 : * +\gamma^{ij}m^{cd}d C_{iac}d C_{jbd}\bigr)
386 : * \Bigr].
387 : * \f} Replace \f$dc^A\rightarrow c^A\f$ to get
388 : * \f{eqnarray}{
389 : * E&=&
390 : * m^{ab}\Bigl[ F_a F_b
391 : * +\gamma^{ij}\bigl( C_{ia} C_{jb}
392 : * +\gamma^{kl}m^{cd} C_{ikac} C_{jlbd}\bigr)
393 : * \nonumber\\
394 : * & + & \Lambda^2\bigl( C_a C_b
395 : * +\gamma^{ij}m^{cd} C_{iac} C_{jbd}\bigr)
396 : * \Bigr]\nonumber\\
397 : * &=&
398 : * m^{ab} F_a F_b
399 : * +m^{ab}\gamma^{ij} C_{ia} C_{jb}
400 : * +m^{ab}\gamma^{ij} \gamma^{kl}m^{cd} C_{ikac} C_{jlbd}
401 : * \nonumber\\
402 : * & + & m^{ab}\Lambda^2 C_a C_b
403 : * +m^{ab}\Lambda^2 \gamma^{ij}m^{cd} C_{iac} C_{jbd}.
404 : * \f} Here \f$m^{ab}\f$ is an arbitrary positive-definite matrix, and
405 : * \f$\Lambda\f$ is an arbitrary real scalar.
406 : * Choose \f$m^{ab} = \delta^{ab}\f$ but allow an arbitrary coefficient to be
407 : * placed in front of each term. Then, absorb \f$\Lambda^2\f$ into one of
408 : * these coefficients, to get
409 : * \f{eqnarray}{ E &=& K_
410 : * F\delta^{ab} F_a F_b +K_2\delta^{ab}\gamma^{ij} C_{ia} C_{jb}
411 : +K_4\delta^{ab}\gamma^{ij}
412 : * \gamma^{kl}\delta^{cd} C_{ikac} C_{jlbd}
413 : * \nonumber\\
414 : * & + & K_1\delta^{ab} C_a C_b
415 : * +K_3\delta^{ab} \gamma^{ij}\delta^{cd} C_{iac} C_{jbd}.
416 : * \f}
417 : * Adopting a Euclidean norm for the constraint space (i.e., choosing to raise
418 : and
419 : * lower spacetime indices with Kronecker deltas) gives
420 : * \f{eqnarray}{ E &=& K_ F
421 : * F_a F_a +K_2\gamma^{ij} C_{ia} C_{ja} +K_4 \gamma^{ij} \gamma^{kl} C_{ikac}
422 : C_{jlac}
423 : * \nonumber\\
424 : * & + & K_1 C_a C_a
425 : * +K_3\gamma^{ij} C_{iac} C_{jac}.
426 : * \f} The two-index constraint and f constraint can be viewed as the
427 : * time and space components of a combined spacetime constraint. So next
428 : * choose
429 : * \f$K_ F=K_2\f$, giving \f{eqnarray}{ E&=& K_1 C_a C_a + K_2\left(F_a F_a
430 : * + C_{ia} C_{ja} \gamma^{ij}\right)\nonumber\\
431 : * & + & K_3 C_{iab} C_{jab} \gamma^{ij}
432 : * + K_4 C_{ikab} C_{jlab}\gamma^{ij} \gamma^{kl}.
433 : * \f}
434 : *
435 : * Note that \f$C_{ikab}\f$ is antisymmetric on the first two indices. Next,
436 : * replace the four-index constraint \f$C_{ijab}\f$ with \f$D_{iab}\f$, which
437 : * contains the independent components of \f$C_{ijab}\f$. Specifically,
438 : * \f{eqnarray}{
439 : * D_{iab} \equiv \frac{1}{2} \epsilon_{i}{}^{jk} C_{jkab}.
440 : * \f} The inverse relationship is
441 : * \f{eqnarray}{
442 : * C_{jkab} = \epsilon^{i}{}_{jk} D_{iab},
443 : * \f} where \f$\epsilon_{ijk}\f$ is the flat-space Levi-Civita tensor, whose
444 : * indices are raised and lowered with the Kronecker delta. Inserting this
445 : relation
446 : * to replace \f$C_{jkab}\f$ with \f$D_{iab}\f$ gives \f{eqnarray}{ E &=& K_1
447 : C_a
448 : * C_a + K_2\left(F_a F_a
449 : * + C_{ia} C_{ja} \gamma^{ij}\right)\nonumber\\
450 : * & + & K_3 C_{iab} C_{jab} \gamma^{ij}\nonumber\\
451 : * & + & K_4 D_{mab} D_{nab} \epsilon^{m}{}_{ik}
452 : * \epsilon^{n}{}_{jl} \gamma^{ij} \gamma^{kl}. \f}
453 : *
454 : * There's a subtle point here: \f$\gamma^{ij}\f$ is the inverse spatial metric,
455 : which
456 : * is not necessarily flat. But \f$\epsilon^{i}{}_{jk}\f$ is the flat space
457 : * Levi-Civita tensor. In order to raise and lower indices of the Levi-Civita
458 : * tensor with the inverse spatial metrics, put in the appropriate factors of
459 : * \f$\sqrt{\gamma}\f$, where \f$\gamma\f$ is the metric determinant, to make
460 : the
461 : * curved-space Levi-Civita tensor compatible with \f$\gamma^{ij}\f$. Let
462 : * \f$\varepsilon^{ijk}\f$ represent the curved space Levi-Civita tensor
463 : compatible
464 : * with \f$\gamma^{ij}\f$: \f{eqnarray}{
465 : * \varepsilon^{mik} = \gamma^{-1/2} \epsilon^{mik}\\
466 : * \varepsilon_{mik} = \gamma^{1/2} \epsilon_{mik}.
467 : * \f} Then we can write the constraint energy as
468 : * \f{eqnarray}{
469 : * E &=& K_1 C_a C_a + K_2\left(F_a F_a
470 : * + C_{ia} C_{ja} \gamma^{ij}\right)\nonumber\\
471 : * & + & K_3 C_{iab} C_{jab} \gamma^{ij}\nonumber\\
472 : * & + & K_4 D_{mab} D_{nab} \gamma \gamma^{-1/2}\epsilon^{m}{}_{ik}
473 : * \gamma^{-1/2}\epsilon^{n}{}_{jl} \gamma^{ij} \gamma^{kl}. \f} The factors of
474 : * \f$\gamma^{-1/2}\f$ make the Levi-Civita tensor compatible with
475 : \f$\gamma^{ij}\f$.
476 : * Swapping which summed indices are raised and which are lowered gives
477 : * \f{eqnarray}{ E &=& K_1 C_a
478 : C_a +
479 : * K_2\left(F_a F_a
480 : * + C_{ia} C_{ja} \gamma^{ij}\right)\nonumber\\
481 : * & + & K_3 C_{iab} C_{jab} \gamma^{ij}\nonumber\\
482 : * & + & K_4 D_{mab} D_{nab} \gamma \gamma^{-1/2}\epsilon^{mik}
483 : \gamma^{-1/2}\epsilon^{njl}
484 : * \gamma_{ij} \gamma_{kl}, \f} or \f{eqnarray}{ E &=& K_1 C_a C_a +
485 : K_2\left(F_a F_a
486 : * + C_{ia} C_{ja} \gamma^{ij}\right)\nonumber\\
487 : * & + & K_3 C_{iab} C_{jab} \gamma^{ij}\nonumber\\
488 : * & + & K_4 D_{mab} D_{nab} \gamma \varepsilon^{mik} \varepsilon^{njl}
489 : \gamma_{ij}
490 : * \gamma_{kl}, \f} or, reversing up and down repeated indices again,
491 : * \f{eqnarray}{ E
492 : * &=& K_1 C_a C_a + K_2\left(F_a F_a
493 : * + C_{ia} C_{ja} \gamma^{ij}\right)\nonumber\\
494 : * & + & K_3 C_{iab} C_{jab} \gamma^{ij}\nonumber\\
495 : * & + & K_4 D_{mab} D_{nab} \gamma \varepsilon^{m}{}_{ik}
496 : \varepsilon^{n}{}_{jl}
497 : * \gamma^{ij} \gamma^{kl}. \f}
498 : *
499 : * The metric raises and lowers the indices of \f$\varepsilon^{ijk}\f$,
500 : * so this can
501 : * be written as \f{eqnarray}{ E &=& K_1 C_a C_a + K_2\left(F_a F_a
502 : * + C_{ia} C_{ja} \gamma^{ij}\right)\nonumber\\
503 : * & + & K_3 C_{iab} C_{jab} \gamma^{ij}\nonumber\\
504 : * & + & K_4 \gamma D_{mab} D^{n}{}_{ab} \varepsilon^{mjl}
505 : \varepsilon_{njl}.
506 : * \f}
507 : *
508 : * Now, in flat space (Eq. (1.23) of \cite ThorneBlandford2017),
509 : * \f{eqnarray}{
510 : * \epsilon^{mjl} \epsilon_{njl} = \delta^{mj}_{nj} = \delta^m_n \delta^j_j -
511 : * \delta^m_j \delta^j_n = 2 \delta^m_n. \f} But this holds for curved space
512 : * as well: multiply the left hand side by \f$1 = \gamma^{1/2} \gamma^{-1/2}\f$
513 : to get
514 : * \f{eqnarray}{
515 : * \gamma^{-1/2}\epsilon^{mjl} \gamma^{1/2}\epsilon_{njl} = \varepsilon^{mjl}
516 : * \varepsilon_{njl} = \delta^{mj}_{nj} = \delta^m_n \delta^j_j - \delta^m_j
517 : * \delta^j_n = 2 \delta^m_n. \f} So the constraint energy is \f{eqnarray}{ E
518 : &=&
519 : * K_1 C_a C_a + K_2\left(F_a F_a
520 : * + C_{ia} C_{ja} \gamma^{ij}\right)\nonumber\\
521 : * & + & K_3 C_{iab} C_{jab} \gamma^{ij}\nonumber\\
522 : * & + & 2 K_4 D_{mab} D^{n}{}_{ab} \delta^m_n.
523 : * \f}
524 : * Simplifying gives the formula implemented here:
525 : * \f{eqnarray}{
526 : * E &=& K_1 C_a C_a + K_2\left(F_a F_a
527 : * + C_{ia} C_{ja} \gamma^{ij}\right)\nonumber\\
528 : * & + & K_3 C_{iab} C_{jab} \gamma^{ij}\nonumber\\
529 : * & + & 2 K_4 \gamma D_{iab} D_{jab} \gamma^{ij}.
530 : * \f}
531 : */
532 : template <typename DataType, size_t SpatialDim, typename Frame>
533 1 : Scalar<DataType> constraint_energy(
534 : const tnsr::a<DataType, SpatialDim, Frame>& gauge_constraint,
535 : const tnsr::a<DataType, SpatialDim, Frame>& f_constraint,
536 : const tnsr::ia<DataType, SpatialDim, Frame>& two_index_constraint,
537 : const tnsr::iaa<DataType, SpatialDim, Frame>& three_index_constraint,
538 : const tnsr::iaa<DataType, SpatialDim, Frame>& four_index_constraint,
539 : const tnsr::II<DataType, SpatialDim, Frame>& inverse_spatial_metric,
540 : const Scalar<DataType>& spatial_metric_determinant,
541 : double gauge_constraint_multiplier = 1.0,
542 : double two_index_constraint_multiplier = 1.0,
543 : double three_index_constraint_multiplier = 1.0,
544 : double four_index_constraint_multiplier = 1.0);
545 :
546 : template <typename DataType, size_t SpatialDim, typename Frame>
547 1 : void constraint_energy(
548 : gsl::not_null<Scalar<DataType>*> energy,
549 : const tnsr::a<DataType, SpatialDim, Frame>& gauge_constraint,
550 : const tnsr::a<DataType, SpatialDim, Frame>& f_constraint,
551 : const tnsr::ia<DataType, SpatialDim, Frame>& two_index_constraint,
552 : const tnsr::iaa<DataType, SpatialDim, Frame>& three_index_constraint,
553 : const tnsr::iaa<DataType, SpatialDim, Frame>& four_index_constraint,
554 : const tnsr::II<DataType, SpatialDim, Frame>& inverse_spatial_metric,
555 : const Scalar<DataType>& spatial_metric_determinant,
556 : double gauge_constraint_multiplier = 1.0,
557 : double two_index_constraint_multiplier = 1.0,
558 : double three_index_constraint_multiplier = 1.0,
559 : double four_index_constraint_multiplier = 1.0);
560 : /// @}
561 :
562 : /// @{
563 : /*!
564 : * \brief Computes the generalized-harmonic constraint energy normalization.
565 : *
566 : * \details Computes the generalized-harmonic normalized constraint energy
567 : * integrand [Eq. (70) of \cite Lindblom2005qh with \f$m^{ab}=\delta^{ab}\f$],
568 : * \f{eqnarray}{
569 : * ||E||=\gamma^{ij}\delta^{ab}\delta^{cd}
570 : * (\Lambda^{2}\partial_{i}g_{ac}\partial_{j}g_{bd}
571 : * + \partial_{i}\Pi_{ac}\partial_{j}\Pi_{bd} +
572 : * \gamma^{kl}\partial_{i}\Phi_{kac}\partial_{j}\Phi_{lbd})\gamma^{1/2}
573 : * \f}
574 : *
575 : * Here \f$\gamma^{ij}\f$ is the inverse spatial metric, and \f$\Lambda^{2}\f$
576 : * is the squared dimensional constant. Note
577 : * that in this equation, spacetime indices \f$a,b\f$ are raised and lowered
578 : * with the Kronecker delta.
579 : */
580 : template <typename DataType, size_t SpatialDim, typename Frame>
581 1 : Scalar<DataType> constraint_energy_normalization(
582 : const tnsr::iaa<DataType, SpatialDim, Frame>& d_spacetime_metric,
583 : const tnsr::iaa<DataType, SpatialDim, Frame>& d_pi,
584 : const tnsr::ijaa<DataType, SpatialDim, Frame>& d_phi,
585 : const tnsr::II<DataType, SpatialDim, Frame>& inverse_spatial_metric,
586 : const Scalar<DataType>& sqrt_spatial_metric_determinant,
587 : double dimensional_constant);
588 :
589 : template <typename DataType, size_t SpatialDim, typename Frame>
590 1 : void constraint_energy_normalization(
591 : gsl::not_null<Scalar<DataType>*> energy_norm,
592 : const tnsr::iaa<DataType, SpatialDim, Frame>& d_spacetime_metric,
593 : const tnsr::iaa<DataType, SpatialDim, Frame>& d_pi,
594 : const tnsr::ijaa<DataType, SpatialDim, Frame>& d_phi,
595 : const tnsr::II<DataType, SpatialDim, Frame>& inverse_spatial_metric,
596 : const Scalar<DataType>& sqrt_spatial_metric_determinant,
597 : double dimensional_constant);
598 : /// @}
599 :
600 : /// @{
601 : /*!
602 : * \brief Computes the generalized-harmonic normalized constraint energy.
603 : *
604 : * \details Computes `constraint_energy()` and then normalizes by dividing
605 : * the unnormalized constraint energy. by `(eps +
606 : * constraint_energy_normalization())`, where eps (following SpEC) has a value
607 : * of 1.e-10, The `eps` term is present so that situations like Minkowski space,
608 : * in which `constraint_energy_normalization()` vanishes, do not cause
609 : * floating-point exceptions when computing the normalized constraint energy.
610 : */
611 : template <typename DataType>
612 1 : Scalar<DataType> normalized_constraint_energy(
613 : const Scalar<DataType>& constraint_energy,
614 : const Scalar<DataType>& constraint_energy_normalization);
615 :
616 : template <typename DataType>
617 1 : void normalized_constraint_energy(
618 : gsl::not_null<Scalar<DataType>*> normalized_energy,
619 : const Scalar<DataType>& constraint_energy,
620 : const Scalar<DataType>& constraint_energy_normalization);
621 : /// @}
622 :
623 : namespace Tags {
624 : /*!
625 : * \brief Compute item to get the gauge constraint for the
626 : * generalized harmonic evolution system.
627 : *
628 : * \details See `gauge_constraint()`. Can be retrieved using
629 : * `gh::Tags::GaugeConstraint`.
630 : */
631 : template <size_t SpatialDim, typename Frame>
632 1 : struct GaugeConstraintCompute : GaugeConstraint<DataVector, SpatialDim, Frame>,
633 : db::ComputeTag {
634 0 : using argument_tags = tmpl::list<
635 : GaugeH<DataVector, SpatialDim, Frame>,
636 : gr::Tags::SpacetimeNormalOneForm<DataVector, SpatialDim, Frame>,
637 : gr::Tags::SpacetimeNormalVector<DataVector, SpatialDim, Frame>,
638 : gr::Tags::InverseSpatialMetric<DataVector, SpatialDim, Frame>,
639 : gr::Tags::InverseSpacetimeMetric<DataVector, SpatialDim, Frame>,
640 : Pi<DataVector, SpatialDim, Frame>, Phi<DataVector, SpatialDim, Frame>>;
641 :
642 0 : using return_type = tnsr::a<DataVector, SpatialDim, Frame>;
643 :
644 0 : static constexpr auto function = static_cast<void (*)(
645 : gsl::not_null<tnsr::a<DataVector, SpatialDim, Frame>*>,
646 : const tnsr::a<DataVector, SpatialDim, Frame>&,
647 : const tnsr::a<DataVector, SpatialDim, Frame>&,
648 : const tnsr::A<DataVector, SpatialDim, Frame>&,
649 : const tnsr::II<DataVector, SpatialDim, Frame>&,
650 : const tnsr::AA<DataVector, SpatialDim, Frame>&,
651 : const tnsr::aa<DataVector, SpatialDim, Frame>&,
652 : const tnsr::iaa<DataVector, SpatialDim, Frame>&)>(
653 : &gauge_constraint<DataVector, SpatialDim, Frame>);
654 :
655 0 : using base = GaugeConstraint<DataVector, SpatialDim, Frame>;
656 : };
657 :
658 : /*!
659 : * \brief Compute item to get the F-constraint for the
660 : * generalized harmonic evolution system.
661 : *
662 : * \details See `f_constraint()`. Can be retrieved using
663 : * `gh::Tags::FConstraint`.
664 : *
665 : * \note If the system contains matter, you will need to use a system-specific
666 : * version of this compute tag that passes the appropriate stress-energy tensor
667 : * to the F-constraint calculation.
668 : */
669 : template <size_t SpatialDim, typename Frame>
670 1 : struct FConstraintCompute : FConstraint<DataVector, SpatialDim, Frame>,
671 : db::ComputeTag {
672 0 : using argument_tags = tmpl::list<
673 : GaugeH<DataVector, SpatialDim, Frame>,
674 : SpacetimeDerivGaugeH<DataVector, SpatialDim, Frame>,
675 : gr::Tags::SpacetimeNormalOneForm<DataVector, SpatialDim, Frame>,
676 : gr::Tags::SpacetimeNormalVector<DataVector, SpatialDim, Frame>,
677 : gr::Tags::InverseSpatialMetric<DataVector, SpatialDim, Frame>,
678 : gr::Tags::InverseSpacetimeMetric<DataVector, SpatialDim, Frame>,
679 : Pi<DataVector, SpatialDim, Frame>, Phi<DataVector, SpatialDim, Frame>,
680 : ::Tags::deriv<Pi<DataVector, SpatialDim, Frame>, tmpl::size_t<SpatialDim>,
681 : Frame>,
682 : ::Tags::deriv<Phi<DataVector, SpatialDim, Frame>,
683 : tmpl::size_t<SpatialDim>, Frame>,
684 : ::gh::Tags::ConstraintGamma2,
685 : ThreeIndexConstraint<DataVector, SpatialDim, Frame>>;
686 :
687 0 : using return_type = tnsr::a<DataVector, SpatialDim, Frame>;
688 :
689 0 : static constexpr auto function = static_cast<void (*)(
690 : gsl::not_null<tnsr::a<DataVector, SpatialDim, Frame>*>,
691 : const tnsr::a<DataVector, SpatialDim, Frame>&,
692 : const tnsr::ab<DataVector, SpatialDim, Frame>&,
693 : const tnsr::a<DataVector, SpatialDim, Frame>&,
694 : const tnsr::A<DataVector, SpatialDim, Frame>&,
695 : const tnsr::II<DataVector, SpatialDim, Frame>&,
696 : const tnsr::AA<DataVector, SpatialDim, Frame>&,
697 : const tnsr::aa<DataVector, SpatialDim, Frame>&,
698 : const tnsr::iaa<DataVector, SpatialDim, Frame>&,
699 : const tnsr::iaa<DataVector, SpatialDim, Frame>&,
700 : const tnsr::ijaa<DataVector, SpatialDim, Frame>&,
701 : const Scalar<DataVector>&,
702 : const tnsr::iaa<DataVector, SpatialDim, Frame>&)>(
703 : &f_constraint<DataVector, SpatialDim, Frame>);
704 :
705 0 : using base = FConstraint<DataVector, SpatialDim, Frame>;
706 : };
707 :
708 : /*!
709 : * \brief Compute item to get the two-index constraint for the
710 : * generalized harmonic evolution system.
711 : *
712 : * \details See `two_index_constraint()`. Can be retrieved using
713 : * `gh::Tags::TwoIndexConstraint`.
714 : */
715 : template <size_t SpatialDim, typename Frame>
716 1 : struct TwoIndexConstraintCompute
717 : : TwoIndexConstraint<DataVector, SpatialDim, Frame>,
718 : db::ComputeTag {
719 0 : using argument_tags = tmpl::list<
720 : SpacetimeDerivGaugeH<DataVector, SpatialDim, Frame>,
721 : gr::Tags::SpacetimeNormalOneForm<DataVector, SpatialDim, Frame>,
722 : gr::Tags::SpacetimeNormalVector<DataVector, SpatialDim, Frame>,
723 : gr::Tags::InverseSpatialMetric<DataVector, SpatialDim, Frame>,
724 : gr::Tags::InverseSpacetimeMetric<DataVector, SpatialDim, Frame>,
725 : Pi<DataVector, SpatialDim, Frame>, Phi<DataVector, SpatialDim, Frame>,
726 : ::Tags::deriv<Pi<DataVector, SpatialDim, Frame>, tmpl::size_t<SpatialDim>,
727 : Frame>,
728 : ::Tags::deriv<Phi<DataVector, SpatialDim, Frame>,
729 : tmpl::size_t<SpatialDim>, Frame>,
730 : ::gh::Tags::ConstraintGamma2,
731 : ThreeIndexConstraint<DataVector, SpatialDim, Frame>>;
732 :
733 0 : using return_type = tnsr::ia<DataVector, SpatialDim, Frame>;
734 :
735 0 : static constexpr auto function = static_cast<void (*)(
736 : gsl::not_null<tnsr::ia<DataVector, SpatialDim, Frame>*>,
737 : const tnsr::ab<DataVector, SpatialDim, Frame>&,
738 : const tnsr::a<DataVector, SpatialDim, Frame>&,
739 : const tnsr::A<DataVector, SpatialDim, Frame>&,
740 : const tnsr::II<DataVector, SpatialDim, Frame>&,
741 : const tnsr::AA<DataVector, SpatialDim, Frame>&,
742 : const tnsr::aa<DataVector, SpatialDim, Frame>&,
743 : const tnsr::iaa<DataVector, SpatialDim, Frame>&,
744 : const tnsr::iaa<DataVector, SpatialDim, Frame>&,
745 : const tnsr::ijaa<DataVector, SpatialDim, Frame>&,
746 : const Scalar<DataVector>&,
747 : const tnsr::iaa<DataVector, SpatialDim, Frame>&)>(
748 : &two_index_constraint<DataVector, SpatialDim, Frame>);
749 :
750 0 : using base = TwoIndexConstraint<DataVector, SpatialDim, Frame>;
751 : };
752 :
753 : /*!
754 : * \brief Compute item to get the three-index constraint for the
755 : * generalized harmonic evolution system.
756 : *
757 : * \details See `three_index_constraint()`. Can be retrieved using
758 : * `gh::Tags::ThreeIndexConstraint`.
759 : */
760 : template <size_t SpatialDim, typename Frame>
761 1 : struct ThreeIndexConstraintCompute
762 : : ThreeIndexConstraint<DataVector, SpatialDim, Frame>,
763 : db::ComputeTag {
764 0 : using argument_tags = tmpl::list<
765 : ::Tags::deriv<gr::Tags::SpacetimeMetric<DataVector, SpatialDim, Frame>,
766 : tmpl::size_t<SpatialDim>, Frame>,
767 : Phi<DataVector, SpatialDim, Frame>>;
768 :
769 0 : using return_type = tnsr::iaa<DataVector, SpatialDim, Frame>;
770 :
771 0 : static constexpr auto function = static_cast<void (*)(
772 : gsl::not_null<tnsr::iaa<DataVector, SpatialDim, Frame>*>,
773 : const tnsr::iaa<DataVector, SpatialDim, Frame>&,
774 : const tnsr::iaa<DataVector, SpatialDim, Frame>&)>(
775 : &three_index_constraint<DataVector, SpatialDim, Frame>);
776 :
777 0 : using base = ThreeIndexConstraint<DataVector, SpatialDim, Frame>;
778 : };
779 :
780 : /*!
781 : * \brief Compute item to get the four-index constraint for the
782 : * generalized harmonic evolution system.
783 : *
784 : * \details See `four_index_constraint()`. Can be retrieved using
785 : * `gh::Tags::FourIndexConstraint`.
786 : */
787 : template <size_t SpatialDim, typename Frame>
788 1 : struct FourIndexConstraintCompute
789 : : FourIndexConstraint<DataVector, SpatialDim, Frame>,
790 : db::ComputeTag {
791 0 : using argument_tags =
792 : tmpl::list<::Tags::deriv<Phi<DataVector, SpatialDim, Frame>,
793 : tmpl::size_t<SpatialDim>, Frame>>;
794 :
795 0 : using return_type = tnsr::iaa<DataVector, SpatialDim, Frame>;
796 :
797 0 : static constexpr auto function = static_cast<void (*)(
798 : gsl::not_null<tnsr::iaa<DataVector, SpatialDim, Frame>*>,
799 : const tnsr::ijaa<DataVector, SpatialDim, Frame>&)>(
800 : &four_index_constraint<DataVector, SpatialDim, Frame>);
801 :
802 0 : using base = FourIndexConstraint<DataVector, SpatialDim, Frame>;
803 : };
804 :
805 : /*!
806 : * \brief Compute item to get combined energy in all constraints for the
807 : * generalized harmonic evolution system.
808 : *
809 : * \details See `constraint_energy()`. Can be retrieved using
810 : * `gh::Tags::ConstraintEnergy`.
811 : */
812 : template <size_t SpatialDim, typename Frame>
813 1 : struct ConstraintEnergyCompute
814 : : ConstraintEnergy<DataVector, SpatialDim, Frame>,
815 : db::ComputeTag {
816 0 : using argument_tags =
817 : tmpl::list<GaugeConstraint<DataVector, SpatialDim, Frame>,
818 : FConstraint<DataVector, SpatialDim, Frame>,
819 : TwoIndexConstraint<DataVector, SpatialDim, Frame>,
820 : ThreeIndexConstraint<DataVector, SpatialDim, Frame>,
821 : FourIndexConstraint<DataVector, SpatialDim, Frame>,
822 : gr::Tags::InverseSpatialMetric<DataVector, SpatialDim, Frame>,
823 : gr::Tags::DetSpatialMetric<DataVector>>;
824 :
825 0 : using return_type = Scalar<DataVector>;
826 :
827 0 : static constexpr auto function(
828 : const gsl::not_null<Scalar<DataVector>*> energy,
829 : const tnsr::a<DataVector, SpatialDim, Frame>& gauge_constraint,
830 : const tnsr::a<DataVector, SpatialDim, Frame>& f_constraint,
831 : const tnsr::ia<DataVector, SpatialDim, Frame>& two_index_constraint,
832 : const tnsr::iaa<DataVector, SpatialDim, Frame>& three_index_constraint,
833 : const tnsr::iaa<DataVector, SpatialDim, Frame>& four_index_constraint,
834 : const tnsr::II<DataVector, SpatialDim, Frame>& inverse_spatial_metric,
835 : const Scalar<DataVector>& spatial_metric_determinant) {
836 : set_number_of_grid_points(energy, spatial_metric_determinant);
837 : constraint_energy<DataVector, SpatialDim, Frame>(
838 : energy, gauge_constraint, f_constraint, two_index_constraint,
839 : three_index_constraint, four_index_constraint, inverse_spatial_metric,
840 : spatial_metric_determinant);
841 : }
842 :
843 0 : using base = ConstraintEnergy<DataVector, SpatialDim, Frame>;
844 : };
845 :
846 : /*!
847 : * \brief Compute item to get the normalized combined energy in all constraints
848 : * for the generalized harmonic evolution system.
849 : *
850 : * \details See `normalized_constraint_energy()`. Can be retrieved using
851 : * `gh::Tags::NormalizedConstraintEnergy`.
852 : */
853 : template <size_t SpatialDim, typename Frame>
854 1 : struct NormalizedConstraintEnergyCompute
855 : : NormalizedConstraintEnergy<DataVector, SpatialDim, Frame>,
856 : db::ComputeTag {
857 0 : using argument_tags = tmpl::list<
858 : ConstraintEnergy<DataVector, SpatialDim, Frame>,
859 : ::Tags::deriv<gr::Tags::SpacetimeMetric<DataVector, SpatialDim, Frame>,
860 : tmpl::size_t<SpatialDim>, Frame>,
861 : ::Tags::deriv<Pi<DataVector, SpatialDim, Frame>, tmpl::size_t<SpatialDim>,
862 : Frame>,
863 : ::Tags::deriv<Phi<DataVector, SpatialDim, Frame>,
864 : tmpl::size_t<SpatialDim>, Frame>,
865 : gr::Tags::InverseSpatialMetric<DataVector, SpatialDim, Frame>,
866 : gr::Tags::SqrtDetSpatialMetric<DataVector>>;
867 :
868 0 : using return_type = Scalar<DataVector>;
869 :
870 0 : static constexpr auto function(
871 : const gsl::not_null<Scalar<DataVector>*> normalized_energy,
872 : const Scalar<DataVector>& constraint_energy,
873 : const tnsr::iaa<DataVector, SpatialDim, Frame>& d_spacetime_metric,
874 : const tnsr::iaa<DataVector, SpatialDim, Frame>& d_pi,
875 : const tnsr::ijaa<DataVector, SpatialDim, Frame>& d_phi,
876 : const tnsr::II<DataVector, SpatialDim, Frame>& inverse_spatial_metric,
877 : const Scalar<DataVector>& sqrt_spatial_metric_determinant) {
878 : set_number_of_grid_points(normalized_energy, constraint_energy);
879 : const auto normalization =
880 : constraint_energy_normalization<DataVector, SpatialDim, Frame>(
881 : d_spacetime_metric, d_pi, d_phi, inverse_spatial_metric,
882 : sqrt_spatial_metric_determinant, 1.0);
883 : normalized_constraint_energy<DataVector>(normalized_energy,
884 : constraint_energy, normalization);
885 : }
886 :
887 0 : using base = NormalizedConstraintEnergy<DataVector, SpatialDim, Frame>;
888 : };
889 : } // namespace Tags
890 : } // namespace gh
|
