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/ConstraintDamping/Tags.hpp"
16 : #include "Evolution/Systems/GeneralizedHarmonic/Tags.hpp"
17 : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.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\psi_{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 + g^{ij} \Phi_{ija} + t^b \Pi_{ba}
57 : * - \frac{1}{2} g^i_a \psi^{bc} \Phi_{ibc}
58 : * - \frac{1}{2} t_a \psi^{bc} \Pi_{bc},
59 : * \f]
60 : * where \f$H_a\f$ is the gauge function,
61 : * \f$\psi_{ab}\f$ is the spacetime metric,
62 : * \f$\Pi_{ab}=-t^c\partial_c \psi_{ab}\f$, and
63 : * \f$\Phi_{iab} = \partial_i\psi_{ab}\f$; \f$t^a\f$ is the timelike unit
64 : * normal vector to the spatial slice, \f$g^{ij}\f$ is the inverse spatial
65 : * metric, and \f$g^b_c = \delta^b_c + t^b t_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& g^{jk}\partial_j \Phi_{ika}
97 : * - \frac{1}{2} g_a^j\psi^{cd}\partial_j \Phi_{icd}
98 : * + t^b \partial_i \Pi_{ba}
99 : * - \frac{1}{2} t_a \psi^{cd}\partial_i\Pi_{cd}
100 : * \nonumber\\&&
101 : * + \partial_i H_a
102 : * + \frac{1}{2} g_a^j \Phi_{jcd} \Phi_{ief}
103 : * \psi^{ce}\psi^{df}
104 : * + \frac{1}{2} g^{jk} \Phi_{jcd} \Phi_{ike}
105 : * \psi^{cd}t^e t_a
106 : * \nonumber\\&&
107 : * - g^{jk}g^{mn}\Phi_{jma}\Phi_{ikn}
108 : * + \frac{1}{2} \Phi_{icd} \Pi_{be} t_a
109 : * \left(\psi^{cb}\psi^{de}
110 : * +\frac{1}{2}\psi^{be} t^c t^d\right)
111 : * \nonumber\\&&
112 : * - \Phi_{icd} \Pi_{ba} t^c \left(\psi^{bd}
113 : * +\frac{1}{2} t^b t^d\right)
114 : * + \frac{1}{2} \gamma_2 \left(t_a \psi^{cd}
115 : * - 2 \delta^c_a t^d\right) C_{icd}.
116 : * \f}
117 : * where \f$H_a\f$ is the gauge function,
118 : * \f$\psi_{ab}\f$ is the spacetime metric,
119 : * \f$\Pi_{ab}=-t^c\partial_c \psi_{ab}\f$, and
120 : * \f$\Phi_{iab} = \partial_i\psi_{ab}\f$; \f$t^a\f$ is the timelike unit
121 : * normal vector to the spatial slice, \f$g^{ij}\f$ is the inverse spatial
122 : * metric, and \f$g^b_c = \delta^b_c + t^b t_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\psi_{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} g_a^i \psi^{bc}\partial_i \Pi_{bc}
199 : * - g^{ij} \partial_i \Pi_{ja}
200 : * - g^{ij} t^b \partial_i \Phi_{jba}
201 : * + \frac{1}{2} t_a \psi^{bc} g^{ij} \partial_i \Phi_{jbc}
202 : * \nonumber \\ &&
203 : * + t_a g^{ij} \partial_i H_j
204 : * + g_a^i \Phi_{ijb} g^{jk}\Phi_{kcd} \psi^{bd} t^c
205 : * - \frac{1}{2} g_a^i \Phi_{ijb} g^{jk}
206 : * \Phi_{kcd} \psi^{cd} t^b
207 : * \nonumber \\ &&
208 : * - g_a^i t^b \partial_i H_b
209 : * + g^{ij} \Phi_{icd} \Phi_{jba} \psi^{bc} t^d
210 : * - \frac{1}{2} t_a g^{ij} g^{mn} \Phi_{imc} \Phi_{njd}\psi^{cd}
211 : * \nonumber \\ &&
212 : * - \frac{1}{4} t_a g^{ij}\Phi_{icd}\Phi_{jbe}
213 : * \psi^{cb}\psi^{de}
214 : * + \frac{1}{4} t_a \Pi_{cd} \Pi_{be}
215 : * \psi^{cb}\psi^{de}
216 : * - g^{ij} H_i \Pi_{ja}
217 : * \nonumber \\ &&
218 : * - t^b g^{ij} \Pi_{b i} \Pi_{ja}
219 : * - \frac{1}{4} g_a^i \Phi_{icd} t^c t^d \Pi_{be}
220 : * \psi^{be}
221 : * + \frac{1}{2} t_a \Pi_{cd} \Pi_{be}\psi^{ce}
222 : * t^d t^b
223 : * \nonumber \\ &&
224 : * + g_a^i \Phi_{icd} \Pi_{be} t^c t^b \psi^{de}
225 : * - g^{ij}\Phi_{iba} t^b \Pi_{je} t^e
226 : * - \frac{1}{2} g^{ij}\Phi_{icd} t^c t^d \Pi_{ja}
227 : * \nonumber \\ &&
228 : * - g^{ij} H_i \Phi_{jba} t^b
229 : * + g_{a}^i \Phi_{icd} H_b \psi^{bc} t^d
230 : * +\gamma_2\bigl(g^{id}{\cal C}_{ida}
231 : * -\frac{1}{2} g_a^i\psi^{cd}{\cal C}_{icd}\bigr)
232 : * \nonumber \\ &&
233 : * + \frac{1}{2} t_a \Pi_{cd}\psi^{cd} H_b t^b
234 : * - t_a g^{ij} \Phi_{ijc} H_d \psi^{cd}
235 : * +\frac{1}{2} t_a g^{ij} H_i \Phi_{jcd}\psi^{cd}
236 : * \nonumber \\ &&
237 : * - 16 \pi t^a T_{a b}
238 : * \f}
239 : * where \f$H_a\f$ is the gauge function,
240 : * \f$\psi_{ab}\f$ is the spacetime metric,
241 : * \f$\Pi_{ab}=-t^c\partial_c \psi_{ab}\f$, and
242 : * \f$\Phi_{iab} = \partial_i\psi_{ab}\f$; \f$t^a\f$ is the timelike unit
243 : * normal vector to the spatial slice, \f$g^{ij}\f$ is the inverse spatial
244 : * metric, \f$g^b_c = \delta^b_c + t^b t_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 + g^{i j} \Phi_{ij a} + t^b \Pi_{ba}
262 : * - \frac{1}{2} g_a{}^i \psi^{bc} \Phi_{i b c}
263 : * - \frac{1}{2} t_a \psi^{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} g^{ij}\right)\nonumber\\
347 : & + & K_3 C_{iab} C_{jab} g^{ij} + K_4 C_{ikab} C_{jlab}g^{ij} g^{kl}.
348 : * \f}
349 : * Here \f$C_a\f$ is the gauge constraint, \f$F_a\f$ is the f constraint,
350 : * \f$C_{ia}\f$ is the two-index constraint, \f$C_{iab}\f$ is the
351 : * three-index constraint, \f$C_{ikab}\f$ is the four-index constraint,
352 : * \f$g^{ij}\f$ is the inverse spatial metric, and
353 : * \f$K_1\f$, \f$K_2\f$, \f$K_3\f$, and \f$K_4\f$ are constant multipliers
354 : * for each term that each default to a value of 1.0. Note that in this
355 : * equation, spacetime indices \f$a,b\f$ are raised and lowered with
356 : * the Kronecker delta.
357 : *
358 : * Also note that the argument `four_index_constraint` is a rank-3 tensor.
359 : * This is because `gh::four_index_constraint()` takes
360 : * advantage of the antisymmetry of the four-index constraint's first two
361 : * indices to only compute and return the independent
362 : * components of \f$C_{ijab}\f$, which can be written as
363 : * \f{eqnarray}{
364 : * D_{iab} \equiv \frac{1}{2} \epsilon_{i}{}^{jk} C_{jkab},
365 : * \f} where \f$\epsilon_{ijk}\f$ is the flat-space Levi-Civita tensor,
366 : * whose inidces are raised and lowered with the Kronecker delta.
367 : * The result is
368 : * \f{eqnarray}{
369 : * E & = & K_1 C_a C_a + K_2\left(F_a F_a
370 : * + C_{ia} C_{ja} g^{ij}\right)\nonumber\\
371 : * & + & K_3 C_{iab} C_{jab} g^{ij}\nonumber\\
372 : & + & 2 K_4 g D_{iab} D_{jab} g^{ij},
373 : * \f} where \f$g\f$ is the determinant of the spatial metric.
374 : *
375 : * To derive this expression for the constraint energy implemented here,
376 : * Eq.~(53) of \cite Lindblom2005qh is
377 : * \f{eqnarray}{
378 : * S_{AB} dc^A dc^B &=&
379 : * m^{ab}\Bigl[d F_ad F_b
380 : * +g^{ij}\bigl(d C_{ia}d C_{jb}
381 : * +g^{kl}m^{cd}d C_{ikac}d C_{jlbd}\bigr)
382 : * \nonumber\\
383 : * & + & \Lambda^2\bigl(d C_ad C_b
384 : * +g^{ij}m^{cd}d C_{iac}d C_{jbd}\bigr)
385 : * \Bigr].
386 : * \f} Replace \f$dc^A\rightarrow c^A\f$ to get
387 : * \f{eqnarray}{
388 : * E&=&
389 : * m^{ab}\Bigl[ F_a F_b
390 : * +g^{ij}\bigl( C_{ia} C_{jb}
391 : * +g^{kl}m^{cd} C_{ikac} C_{jlbd}\bigr)
392 : * \nonumber\\
393 : * & + & \Lambda^2\bigl( C_a C_b
394 : * +g^{ij}m^{cd} C_{iac} C_{jbd}\bigr)
395 : * \Bigr]\nonumber\\
396 : * &=&
397 : * m^{ab} F_a F_b
398 : * +m^{ab}g^{ij} C_{ia} C_{jb}
399 : * +m^{ab}g^{ij} g^{kl}m^{cd} C_{ikac} C_{jlbd}
400 : * \nonumber\\
401 : * & + & m^{ab}\Lambda^2 C_a C_b
402 : * +m^{ab}\Lambda^2 g^{ij}m^{cd} C_{iac} C_{jbd}.
403 : * \f} Here \f$m^{ab}\f$ is an arbitrary positive-definite matrix, and
404 : * \f$\Lambda\f$ is an arbitrary real scalar.
405 : * Choose \f$m^{ab} = \delta^{ab}\f$ but allow an arbitrary coefficient to be
406 : * placed in front of each term. Then, absorb \f$\Lambda^2\f$ into one of
407 : * these coefficients, to get
408 : * \f{eqnarray}{ E &=& K_
409 : * F\delta^{ab} F_a F_b +K_2\delta^{ab}g^{ij} C_{ia} C_{jb}
410 : +K_4\delta^{ab}g^{ij}
411 : * g^{kl}\delta^{cd} C_{ikac} C_{jlbd}
412 : * \nonumber\\
413 : * & + & K_1\delta^{ab} C_a C_b
414 : * +K_3\delta^{ab} g^{ij}\delta^{cd} C_{iac} C_{jbd}.
415 : * \f}
416 : * Adopting a Euclidean norm for the constraint space (i.e., choosing to raise
417 : and
418 : * lower spacetime indices with Kronecker deltas) gives
419 : * \f{eqnarray}{ E &=& K_ F
420 : * F_a F_a +K_2g^{ij} C_{ia} C_{ja} +K_4 g^{ij} g^{kl} C_{ikac} C_{jlac}
421 : * \nonumber\\
422 : * & + & K_1 C_a C_a
423 : * +K_3g^{ij} C_{iac} C_{jac}.
424 : * \f} The two-index constraint and f constraint can be viewed as the
425 : * time and space components of a combined spacetime constraint. So next
426 : * choose
427 : * \f$K_ F=K_2\f$, giving \f{eqnarray}{ E&=& K_1 C_a C_a + K_2\left(F_a F_a
428 : * + C_{ia} C_{ja} g^{ij}\right)\nonumber\\
429 : * & + & K_3 C_{iab} C_{jab} g^{ij}
430 : * + K_4 C_{ikab} C_{jlab}g^{ij} g^{kl}.
431 : * \f}
432 : *
433 : * Note that \f$C_{ikab}\f$ is antisymmetric on the first two indices. Next,
434 : * replace the four-index constraint \f$C_{ijab}\f$ with \f$D_{iab}\f$, which
435 : * contains the independent components of \f$C_{ijab}\f$. Specifically,
436 : * \f{eqnarray}{
437 : * D_{iab} \equiv \frac{1}{2} \epsilon_{i}{}^{jk} C_{jkab}.
438 : * \f} The inverse relationship is
439 : * \f{eqnarray}{
440 : * C_{jkab} = \epsilon^{i}{}_{jk} D_{iab},
441 : * \f} where \f$\epsilon_{ijk}\f$ is the flat-space Levi-Civita tensor, whose
442 : * indices are raised and lowered with the Kronecker delta. Inserting this
443 : relation
444 : * to replace \f$C_{jkab}\f$ with \f$D_{iab}\f$ gives \f{eqnarray}{ E &=& K_1
445 : C_a
446 : * C_a + K_2\left(F_a F_a
447 : * + C_{ia} C_{ja} g^{ij}\right)\nonumber\\
448 : * & + & K_3 C_{iab} C_{jab} g^{ij}\nonumber\\
449 : * & + & K_4 D_{mab} D_{nab} \epsilon^{m}{}_{ik}
450 : * \epsilon^{n}{}_{jl} g^{ij} g^{kl}. \f}
451 : *
452 : * There's a subtle point here: \f$g^{ij}\f$ is the inverse spatial metric,
453 : which
454 : * is not necessarily flat. But \f$\epsilon^{i}{}_{jk}\f$ is the flat space
455 : * Levi-Civita tensor. In order to raise and lower indices of the Levi-Civita
456 : * tensor with the inverse spatial metrics, put in the appropriate factors of
457 : * \f$\sqrt{g}\f$, where \f$g\f$ is the metric determinant, to make the
458 : * curved-space Levi-Civita tensor compatible with \f$g^{ij}\f$. Let
459 : * \f$\varepsilon^{ijk}\f$ represent the curved space Levi-Civita tensor
460 : compatible
461 : * with \f$g^{ij}\f$: \f{eqnarray}{
462 : * \varepsilon^{mik} = g^{-1/2} \epsilon^{mik}\\
463 : * \varepsilon_{mik} = g^{1/2} \epsilon_{mik}.
464 : * \f} Then we can write the constraint energy as
465 : * \f{eqnarray}{
466 : * E &=& K_1 C_a C_a + K_2\left(F_a F_a
467 : * + C_{ia} C_{ja} g^{ij}\right)\nonumber\\
468 : * & + & K_3 C_{iab} C_{jab} g^{ij}\nonumber\\
469 : * & + & K_4 D_{mab} D_{nab} g g^{-1/2}\epsilon^{m}{}_{ik}
470 : * g^{-1/2}\epsilon^{n}{}_{jl} g^{ij} g^{kl}. \f} The factors of
471 : * \f$g^{-1/2}\f$ make the Levi-Civita tensor compatible with \f$g^{ij}\f$.
472 : * Swapping which summed indices are raised and which are lowered gives
473 : * \f{eqnarray}{ E &=& K_1 C_a
474 : C_a +
475 : * K_2\left(F_a F_a
476 : * + C_{ia} C_{ja} g^{ij}\right)\nonumber\\
477 : * & + & K_3 C_{iab} C_{jab} g^{ij}\nonumber\\
478 : * & + & K_4 D_{mab} D_{nab} g g^{-1/2}\epsilon^{mik}
479 : g^{-1/2}\epsilon^{njl}
480 : * g_{ij} g_{kl}, \f} or \f{eqnarray}{ E &=& K_1 C_a C_a + K_2\left(F_a F_a
481 : * + C_{ia} C_{ja} g^{ij}\right)\nonumber\\
482 : * & + & K_3 C_{iab} C_{jab} g^{ij}\nonumber\\
483 : * & + & K_4 D_{mab} D_{nab} g \varepsilon^{mik} \varepsilon^{njl} g_{ij}
484 : * g_{kl}, \f} or, reversing up and down repeated indices again,
485 : * \f{eqnarray}{ E
486 : * &=& K_1 C_a C_a + K_2\left(F_a F_a
487 : * + C_{ia} C_{ja} g^{ij}\right)\nonumber\\
488 : * & + & K_3 C_{iab} C_{jab} g^{ij}\nonumber\\
489 : * & + & K_4 D_{mab} D_{nab} g \varepsilon^{m}{}_{ik}
490 : \varepsilon^{n}{}_{jl}
491 : * g^{ij} g^{kl}. \f}
492 : *
493 : * The metric raises and lowers the indices of \f$\varepsilon^{ijk}\f$,
494 : * so this can
495 : * be written as \f{eqnarray}{ E &=& K_1 C_a C_a + K_2\left(F_a F_a
496 : * + C_{ia} C_{ja} g^{ij}\right)\nonumber\\
497 : * & + & K_3 C_{iab} C_{jab} g^{ij}\nonumber\\
498 : * & + & K_4 g D_{mab} D^{n}{}_{ab} \varepsilon^{mjl} \varepsilon_{njl}.
499 : * \f}
500 : *
501 : * Now, in flat space (Eq. (1.23) of \cite ThorneBlandford2017),
502 : * \f{eqnarray}{
503 : * \epsilon^{mjl} \epsilon_{njl} = \delta^{mj}_{nj} = \delta^m_n \delta^j_j -
504 : * \delta^m_j \delta^j_n = 2 \delta^m_n. \f} But this holds for curved space
505 : * as well: multiply the left hand side by \f$1 = g^{1/2} g^{-1/2}\f$ to get
506 : * \f{eqnarray}{
507 : * g^{-1/2}\epsilon^{mjl} g^{1/2}\epsilon_{njl} = \varepsilon^{mjl}
508 : * \varepsilon_{njl} = \delta^{mj}_{nj} = \delta^m_n \delta^j_j - \delta^m_j
509 : * \delta^j_n = 2 \delta^m_n. \f} So the constraint energy is \f{eqnarray}{ E
510 : &=&
511 : * K_1 C_a C_a + K_2\left(F_a F_a
512 : * + C_{ia} C_{ja} g^{ij}\right)\nonumber\\
513 : * & + & K_3 C_{iab} C_{jab} g^{ij}\nonumber\\
514 : * & + & 2 K_4 D_{mab} D^{n}{}_{ab} \delta^m_n.
515 : * \f}
516 : * Simplifying gives the formula implemented here:
517 : * \f{eqnarray}{
518 : * E &=& K_1 C_a C_a + K_2\left(F_a F_a
519 : * + C_{ia} C_{ja} g^{ij}\right)\nonumber\\
520 : * & + & K_3 C_{iab} C_{jab} g^{ij}\nonumber\\
521 : * & + & 2 K_4 g D_{iab} D_{jab} g^{ij}.
522 : * \f}
523 : */
524 : template <typename DataType, size_t SpatialDim, typename Frame>
525 1 : Scalar<DataType> constraint_energy(
526 : const tnsr::a<DataType, SpatialDim, Frame>& gauge_constraint,
527 : const tnsr::a<DataType, SpatialDim, Frame>& f_constraint,
528 : const tnsr::ia<DataType, SpatialDim, Frame>& two_index_constraint,
529 : const tnsr::iaa<DataType, SpatialDim, Frame>& three_index_constraint,
530 : const tnsr::iaa<DataType, SpatialDim, Frame>& four_index_constraint,
531 : const tnsr::II<DataType, SpatialDim, Frame>& inverse_spatial_metric,
532 : const Scalar<DataType>& spatial_metric_determinant,
533 : double gauge_constraint_multiplier = 1.0,
534 : double two_index_constraint_multiplier = 1.0,
535 : double three_index_constraint_multiplier = 1.0,
536 : double four_index_constraint_multiplier = 1.0);
537 :
538 : template <typename DataType, size_t SpatialDim, typename Frame>
539 1 : void constraint_energy(
540 : gsl::not_null<Scalar<DataType>*> energy,
541 : const tnsr::a<DataType, SpatialDim, Frame>& gauge_constraint,
542 : const tnsr::a<DataType, SpatialDim, Frame>& f_constraint,
543 : const tnsr::ia<DataType, SpatialDim, Frame>& two_index_constraint,
544 : const tnsr::iaa<DataType, SpatialDim, Frame>& three_index_constraint,
545 : const tnsr::iaa<DataType, SpatialDim, Frame>& four_index_constraint,
546 : const tnsr::II<DataType, SpatialDim, Frame>& inverse_spatial_metric,
547 : const Scalar<DataType>& spatial_metric_determinant,
548 : double gauge_constraint_multiplier = 1.0,
549 : double two_index_constraint_multiplier = 1.0,
550 : double three_index_constraint_multiplier = 1.0,
551 : double four_index_constraint_multiplier = 1.0);
552 : /// @}
553 :
554 : /// @{
555 : /*!
556 : * \brief Computes the generalized-harmonic normalized constraint energy.
557 : *
558 : * \details Computes the generalized-harmonic normalized constraint energy
559 : * integrand [Eq. (70) of \cite Lindblom2005qh with \f$m^{ab}=\delta^{ab}\f$],
560 : * \f{eqnarray}{
561 : * ||E||=\gamma^{ij}\delta^{ab}\delta^{cd}
562 : * (\Lambda^{2}\partial_{i}g_{ac}\partial_{j}g_{bd}
563 : * + \partial_{i}\Pi_{ac}\partial_{j}\Pi_{bd} +
564 : * \gamma^{kl}\partial_{i}\Phi_{kac}\partial_{j}\Phi_{lbd})\gamma^{1/2}
565 : * \f}
566 : *
567 : * Here \f$\gamma^{ij}\f$ is the inverse spatial metric, and \f$\Lambda^{2}\f$
568 : * is the squared dimensional constant. Note
569 : * that in this equation, spacetime indices \f$a,b\f$ are raised and lowered
570 : * with the Kronecker delta.
571 : */
572 : template <typename DataType, size_t SpatialDim, typename Frame>
573 1 : Scalar<DataType> constraint_energy_normalization(
574 : const tnsr::iaa<DataType, SpatialDim, Frame>& d_spacetime_metric,
575 : const tnsr::iaa<DataType, SpatialDim, Frame>& d_pi,
576 : const tnsr::ijaa<DataType, SpatialDim, Frame>& d_phi,
577 : const tnsr::II<DataType, SpatialDim, Frame>& inverse_spatial_metric,
578 : const Scalar<DataType>& sqrt_spatial_metric_determinant,
579 : double dimensional_constant);
580 :
581 : template <typename DataType, size_t SpatialDim, typename Frame>
582 1 : void constraint_energy_normalization(
583 : gsl::not_null<Scalar<DataType>*> energy_norm,
584 : const tnsr::iaa<DataType, SpatialDim, Frame>& d_spacetime_metric,
585 : const tnsr::iaa<DataType, SpatialDim, Frame>& d_pi,
586 : const tnsr::ijaa<DataType, SpatialDim, Frame>& d_phi,
587 : const tnsr::II<DataType, SpatialDim, Frame>& inverse_spatial_metric,
588 : const Scalar<DataType>& sqrt_spatial_metric_determinant,
589 : double dimensional_constant);
590 : /// @}
591 :
592 : namespace Tags {
593 : /*!
594 : * \brief Compute item to get the gauge constraint for the
595 : * generalized harmonic evolution system.
596 : *
597 : * \details See `gauge_constraint()`. Can be retrieved using
598 : * `gh::Tags::GaugeConstraint`.
599 : */
600 : template <size_t SpatialDim, typename Frame>
601 1 : struct GaugeConstraintCompute : GaugeConstraint<DataVector, SpatialDim, Frame>,
602 : db::ComputeTag {
603 0 : using argument_tags = tmpl::list<
604 : GaugeH<DataVector, SpatialDim, Frame>,
605 : gr::Tags::SpacetimeNormalOneForm<DataVector, SpatialDim, Frame>,
606 : gr::Tags::SpacetimeNormalVector<DataVector, SpatialDim, Frame>,
607 : gr::Tags::InverseSpatialMetric<DataVector, SpatialDim, Frame>,
608 : gr::Tags::InverseSpacetimeMetric<DataVector, SpatialDim, Frame>,
609 : Pi<DataVector, SpatialDim, Frame>, Phi<DataVector, SpatialDim, Frame>>;
610 :
611 0 : using return_type = tnsr::a<DataVector, SpatialDim, Frame>;
612 :
613 0 : static constexpr auto function = static_cast<void (*)(
614 : gsl::not_null<tnsr::a<DataVector, SpatialDim, Frame>*>,
615 : const tnsr::a<DataVector, SpatialDim, Frame>&,
616 : const tnsr::a<DataVector, SpatialDim, Frame>&,
617 : const tnsr::A<DataVector, SpatialDim, Frame>&,
618 : const tnsr::II<DataVector, SpatialDim, Frame>&,
619 : const tnsr::AA<DataVector, SpatialDim, Frame>&,
620 : const tnsr::aa<DataVector, SpatialDim, Frame>&,
621 : const tnsr::iaa<DataVector, SpatialDim, Frame>&)>(
622 : &gauge_constraint<DataVector, SpatialDim, Frame>);
623 :
624 0 : using base = GaugeConstraint<DataVector, SpatialDim, Frame>;
625 : };
626 :
627 : /*!
628 : * \brief Compute item to get the F-constraint for the
629 : * generalized harmonic evolution system.
630 : *
631 : * \details See `f_constraint()`. Can be retrieved using
632 : * `gh::Tags::FConstraint`.
633 : *
634 : * \note If the system contains matter, you will need to use a system-specific
635 : * version of this compute tag that passes the appropriate stress-energy tensor
636 : * to the F-constraint calculation.
637 : */
638 : template <size_t SpatialDim, typename Frame>
639 1 : struct FConstraintCompute : FConstraint<DataVector, SpatialDim, Frame>,
640 : db::ComputeTag {
641 0 : using argument_tags = tmpl::list<
642 : GaugeH<DataVector, SpatialDim, Frame>,
643 : SpacetimeDerivGaugeH<DataVector, SpatialDim, Frame>,
644 : gr::Tags::SpacetimeNormalOneForm<DataVector, SpatialDim, Frame>,
645 : gr::Tags::SpacetimeNormalVector<DataVector, SpatialDim, Frame>,
646 : gr::Tags::InverseSpatialMetric<DataVector, SpatialDim, Frame>,
647 : gr::Tags::InverseSpacetimeMetric<DataVector, SpatialDim, Frame>,
648 : Pi<DataVector, SpatialDim, Frame>, Phi<DataVector, SpatialDim, Frame>,
649 : ::Tags::deriv<Pi<DataVector, SpatialDim, Frame>, tmpl::size_t<SpatialDim>,
650 : Frame>,
651 : ::Tags::deriv<Phi<DataVector, SpatialDim, Frame>,
652 : tmpl::size_t<SpatialDim>, Frame>,
653 : ::gh::ConstraintDamping::Tags::ConstraintGamma2,
654 : ThreeIndexConstraint<DataVector, SpatialDim, Frame>>;
655 :
656 0 : using return_type = tnsr::a<DataVector, SpatialDim, Frame>;
657 :
658 0 : static constexpr auto function = static_cast<void (*)(
659 : gsl::not_null<tnsr::a<DataVector, SpatialDim, Frame>*>,
660 : const tnsr::a<DataVector, SpatialDim, Frame>&,
661 : const tnsr::ab<DataVector, SpatialDim, Frame>&,
662 : const tnsr::a<DataVector, SpatialDim, Frame>&,
663 : const tnsr::A<DataVector, SpatialDim, Frame>&,
664 : const tnsr::II<DataVector, SpatialDim, Frame>&,
665 : const tnsr::AA<DataVector, SpatialDim, Frame>&,
666 : const tnsr::aa<DataVector, SpatialDim, Frame>&,
667 : const tnsr::iaa<DataVector, SpatialDim, Frame>&,
668 : const tnsr::iaa<DataVector, SpatialDim, Frame>&,
669 : const tnsr::ijaa<DataVector, SpatialDim, Frame>&,
670 : const Scalar<DataVector>&,
671 : const tnsr::iaa<DataVector, SpatialDim, Frame>&)>(
672 : &f_constraint<DataVector, SpatialDim, Frame>);
673 :
674 0 : using base = FConstraint<DataVector, SpatialDim, Frame>;
675 : };
676 :
677 : /*!
678 : * \brief Compute item to get the two-index constraint for the
679 : * generalized harmonic evolution system.
680 : *
681 : * \details See `two_index_constraint()`. Can be retrieved using
682 : * `gh::Tags::TwoIndexConstraint`.
683 : */
684 : template <size_t SpatialDim, typename Frame>
685 1 : struct TwoIndexConstraintCompute
686 : : TwoIndexConstraint<DataVector, SpatialDim, Frame>,
687 : db::ComputeTag {
688 0 : using argument_tags = tmpl::list<
689 : SpacetimeDerivGaugeH<DataVector, SpatialDim, Frame>,
690 : gr::Tags::SpacetimeNormalOneForm<DataVector, SpatialDim, Frame>,
691 : gr::Tags::SpacetimeNormalVector<DataVector, SpatialDim, Frame>,
692 : gr::Tags::InverseSpatialMetric<DataVector, SpatialDim, Frame>,
693 : gr::Tags::InverseSpacetimeMetric<DataVector, SpatialDim, Frame>,
694 : Pi<DataVector, SpatialDim, Frame>, Phi<DataVector, SpatialDim, Frame>,
695 : ::Tags::deriv<Pi<DataVector, SpatialDim, Frame>, tmpl::size_t<SpatialDim>,
696 : Frame>,
697 : ::Tags::deriv<Phi<DataVector, SpatialDim, Frame>,
698 : tmpl::size_t<SpatialDim>, Frame>,
699 : ::gh::ConstraintDamping::Tags::ConstraintGamma2,
700 : ThreeIndexConstraint<DataVector, SpatialDim, Frame>>;
701 :
702 0 : using return_type = tnsr::ia<DataVector, SpatialDim, Frame>;
703 :
704 0 : static constexpr auto function = static_cast<void (*)(
705 : gsl::not_null<tnsr::ia<DataVector, SpatialDim, Frame>*>,
706 : const tnsr::ab<DataVector, SpatialDim, Frame>&,
707 : const tnsr::a<DataVector, SpatialDim, Frame>&,
708 : const tnsr::A<DataVector, SpatialDim, Frame>&,
709 : const tnsr::II<DataVector, SpatialDim, Frame>&,
710 : const tnsr::AA<DataVector, SpatialDim, Frame>&,
711 : const tnsr::aa<DataVector, SpatialDim, Frame>&,
712 : const tnsr::iaa<DataVector, SpatialDim, Frame>&,
713 : const tnsr::iaa<DataVector, SpatialDim, Frame>&,
714 : const tnsr::ijaa<DataVector, SpatialDim, Frame>&,
715 : const Scalar<DataVector>&,
716 : const tnsr::iaa<DataVector, SpatialDim, Frame>&)>(
717 : &two_index_constraint<DataVector, SpatialDim, Frame>);
718 :
719 0 : using base = TwoIndexConstraint<DataVector, SpatialDim, Frame>;
720 : };
721 :
722 : /*!
723 : * \brief Compute item to get the three-index constraint for the
724 : * generalized harmonic evolution system.
725 : *
726 : * \details See `three_index_constraint()`. Can be retrieved using
727 : * `gh::Tags::ThreeIndexConstraint`.
728 : */
729 : template <size_t SpatialDim, typename Frame>
730 1 : struct ThreeIndexConstraintCompute
731 : : ThreeIndexConstraint<DataVector, SpatialDim, Frame>,
732 : db::ComputeTag {
733 0 : using argument_tags = tmpl::list<
734 : ::Tags::deriv<gr::Tags::SpacetimeMetric<DataVector, SpatialDim, Frame>,
735 : tmpl::size_t<SpatialDim>, Frame>,
736 : Phi<DataVector, SpatialDim, Frame>>;
737 :
738 0 : using return_type = tnsr::iaa<DataVector, SpatialDim, Frame>;
739 :
740 0 : static constexpr auto function = static_cast<void (*)(
741 : gsl::not_null<tnsr::iaa<DataVector, SpatialDim, Frame>*>,
742 : const tnsr::iaa<DataVector, SpatialDim, Frame>&,
743 : const tnsr::iaa<DataVector, SpatialDim, Frame>&)>(
744 : &three_index_constraint<DataVector, SpatialDim, Frame>);
745 :
746 0 : using base = ThreeIndexConstraint<DataVector, SpatialDim, Frame>;
747 : };
748 :
749 : /*!
750 : * \brief Compute item to get the four-index constraint for the
751 : * generalized harmonic evolution system.
752 : *
753 : * \details See `four_index_constraint()`. Can be retrieved using
754 : * `gh::Tags::FourIndexConstraint`.
755 : */
756 : template <size_t SpatialDim, typename Frame>
757 1 : struct FourIndexConstraintCompute
758 : : FourIndexConstraint<DataVector, SpatialDim, Frame>,
759 : db::ComputeTag {
760 0 : using argument_tags =
761 : tmpl::list<::Tags::deriv<Phi<DataVector, SpatialDim, Frame>,
762 : tmpl::size_t<SpatialDim>, Frame>>;
763 :
764 0 : using return_type = tnsr::iaa<DataVector, SpatialDim, Frame>;
765 :
766 0 : static constexpr auto function = static_cast<void (*)(
767 : gsl::not_null<tnsr::iaa<DataVector, SpatialDim, Frame>*>,
768 : const tnsr::ijaa<DataVector, SpatialDim, Frame>&)>(
769 : &four_index_constraint<DataVector, SpatialDim, Frame>);
770 :
771 0 : using base = FourIndexConstraint<DataVector, SpatialDim, Frame>;
772 : };
773 :
774 : /*!
775 : * \brief Compute item to get combined energy in all constraints for the
776 : * generalized harmonic evolution system.
777 : *
778 : * \details See `constraint_energy()`. Can be retrieved using
779 : * `gh::Tags::ConstraintEnergy`.
780 : */
781 : template <size_t SpatialDim, typename Frame>
782 1 : struct ConstraintEnergyCompute
783 : : ConstraintEnergy<DataVector, SpatialDim, Frame>,
784 : db::ComputeTag {
785 0 : using argument_tags =
786 : tmpl::list<GaugeConstraint<DataVector, SpatialDim, Frame>,
787 : FConstraint<DataVector, SpatialDim, Frame>,
788 : TwoIndexConstraint<DataVector, SpatialDim, Frame>,
789 : ThreeIndexConstraint<DataVector, SpatialDim, Frame>,
790 : FourIndexConstraint<DataVector, SpatialDim, Frame>,
791 : gr::Tags::InverseSpatialMetric<DataVector, SpatialDim, Frame>,
792 : gr::Tags::DetSpatialMetric<DataVector>>;
793 :
794 0 : using return_type = Scalar<DataVector>;
795 :
796 0 : static constexpr auto function(
797 : const gsl::not_null<Scalar<DataVector>*> energy,
798 : const tnsr::a<DataVector, SpatialDim, Frame>& gauge_constraint,
799 : const tnsr::a<DataVector, SpatialDim, Frame>& f_constraint,
800 : const tnsr::ia<DataVector, SpatialDim, Frame>& two_index_constraint,
801 : const tnsr::iaa<DataVector, SpatialDim, Frame>& three_index_constraint,
802 : const tnsr::iaa<DataVector, SpatialDim, Frame>& four_index_constraint,
803 : const tnsr::II<DataVector, SpatialDim, Frame>& inverse_spatial_metric,
804 : const Scalar<DataVector>& spatial_metric_determinant) {
805 : set_number_of_grid_points(energy, spatial_metric_determinant);
806 : constraint_energy<DataVector, SpatialDim, Frame>(
807 : energy, gauge_constraint, f_constraint, two_index_constraint,
808 : three_index_constraint, four_index_constraint, inverse_spatial_metric,
809 : spatial_metric_determinant);
810 : }
811 :
812 0 : using base = ConstraintEnergy<DataVector, SpatialDim, Frame>;
813 : };
814 : } // namespace Tags
815 : } // namespace gh
|
