Constraints.hpp
Go to the documentation of this file.
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 
12 
13 /// \cond
14 namespace gsl {
15 template <typename T>
16 class not_null;
17 } // namespace gsl
18 /// \endcond
19 
20 namespace GeneralizedHarmonic {
21 // @{
22 /*!
23  * \brief Computes the generalized-harmonic 3-index constraint.
24  *
25  * \details Computes the generalized-harmonic 3-index constraint,
26  * \f$C_{iab} = \partial_i\psi_{ab} - \Phi_{iab},\f$ which is
27  * given by Eq. (26) of \cite Lindblom2005qh
28  */
29 template <size_t SpatialDim, typename Frame, typename DataType>
30 tnsr::iaa<DataType, SpatialDim, Frame> three_index_constraint(
31  const tnsr::iaa<DataType, SpatialDim, Frame>& d_spacetime_metric,
32  const tnsr::iaa<DataType, SpatialDim, Frame>& phi) noexcept;
33 
34 template <size_t SpatialDim, typename Frame, typename DataType>
36  gsl::not_null<tnsr::iaa<DataType, SpatialDim, Frame>*> constraint,
37  const tnsr::iaa<DataType, SpatialDim, Frame>& d_spacetime_metric,
38  const tnsr::iaa<DataType, SpatialDim, Frame>& phi) noexcept;
39 // @}
40 
41 // @{
42 /*!
43  * \brief Computes the generalized-harmonic gauge constraint.
44  *
45  * \details Computes the generalized-harmonic gauge constraint
46  * [Eq. (40) of \cite Lindblom2005qh],
47  * \f[
48  * C_a = H_a + g^{ij} \Phi_{ija} + t^b \Pi_{ba}
49  * - \frac{1}{2} g^i_a \psi^{bc} \Phi_{ibc}
50  * - \frac{1}{2} t_a \psi^{bc} \Pi_{bc},
51  * \f]
52  * where \f$H_a\f$ is the gauge function,
53  * \f$\psi_{ab}\f$ is the spacetime metric,
54  * \f$\Pi_{ab}=-t^c\partial_c \psi_{ab}\f$, and
55  * \f$\Phi_{iab} = \partial_i\psi_{ab}\f$; \f$t^a\f$ is the timelike unit
56  * normal vector to the spatial slice, \f$g^{ij}\f$ is the inverse spatial
57  * metric, and \f$g^b_c = \delta^b_c + t^b t_c\f$.
58  */
59 template <size_t SpatialDim, typename Frame, typename DataType>
60 tnsr::a<DataType, SpatialDim, Frame> gauge_constraint(
61  const tnsr::a<DataType, SpatialDim, Frame>& gauge_function,
62  const tnsr::a<DataType, SpatialDim, Frame>& spacetime_normal_one_form,
63  const tnsr::A<DataType, SpatialDim, Frame>& spacetime_normal_vector,
64  const tnsr::II<DataType, SpatialDim, Frame>& inverse_spatial_metric,
65  const tnsr::AA<DataType, SpatialDim, Frame>& inverse_spacetime_metric,
66  const tnsr::aa<DataType, SpatialDim, Frame>& pi,
67  const tnsr::iaa<DataType, SpatialDim, Frame>& phi) noexcept;
68 
69 template <size_t SpatialDim, typename Frame, typename DataType>
70 void gauge_constraint(
71  gsl::not_null<tnsr::a<DataType, SpatialDim, Frame>*> constraint,
72  const tnsr::a<DataType, SpatialDim, Frame>& gauge_function,
73  const tnsr::a<DataType, SpatialDim, Frame>& spacetime_normal_one_form,
74  const tnsr::A<DataType, SpatialDim, Frame>& spacetime_normal_vector,
75  const tnsr::II<DataType, SpatialDim, Frame>& inverse_spatial_metric,
76  const tnsr::AA<DataType, SpatialDim, Frame>& inverse_spacetime_metric,
77  const tnsr::aa<DataType, SpatialDim, Frame>& pi,
78  const tnsr::iaa<DataType, SpatialDim, Frame>& phi) noexcept;
79 // @}
80 
81 // @{
82 /*!
83  * \brief Computes the generalized-harmonic 2-index constraint.
84  *
85  * \details Computes the generalized-harmonic 2-index constraint
86  * [Eq. (44) of \cite Lindblom2005qh],
87  * \f{eqnarray}{
88  * C_{ia} &\equiv& g^{jk}\partial_j \Phi_{ika}
89  * - \frac{1}{2} g_a^j\psi^{cd}\partial_j \Phi_{icd}
90  * + t^b \partial_i \Pi_{ba}
91  * - \frac{1}{2} t_a \psi^{cd}\partial_i\Pi_{cd}
92  * \nonumber\\&&
93  * + \partial_i H_a
94  * + \frac{1}{2} g_a^j \Phi_{jcd} \Phi_{ief}
95  * \psi^{ce}\psi^{df}
96  * + \frac{1}{2} g^{jk} \Phi_{jcd} \Phi_{ike}
97  * \psi^{cd}t^e t_a
98  * \nonumber\\&&
99  * - g^{jk}g^{mn}\Phi_{jma}\Phi_{ikn}
100  * + \frac{1}{2} \Phi_{icd} \Pi_{be} t_a
101  * \left(\psi^{cb}\psi^{de}
102  * +\frac{1}{2}\psi^{be} t^c t^d\right)
103  * \nonumber\\&&
104  * - \Phi_{icd} \Pi_{ba} t^c \left(\psi^{bd}
105  * +\frac{1}{2} t^b t^d\right)
106  * + \frac{1}{2} \gamma_2 \left(t_a \psi^{cd}
107  * - 2 \delta^c_a t^d\right) C_{icd}.
108  * \f}
109  * where \f$H_a\f$ is the gauge function,
110  * \f$\psi_{ab}\f$ is the spacetime metric,
111  * \f$\Pi_{ab}=-t^c\partial_c \psi_{ab}\f$, and
112  * \f$\Phi_{iab} = \partial_i\psi_{ab}\f$; \f$t^a\f$ is the timelike unit
113  * normal vector to the spatial slice, \f$g^{ij}\f$ is the inverse spatial
114  * metric, and \f$g^b_c = \delta^b_c + t^b t_c\f$.
115  */
116 template <size_t SpatialDim, typename Frame, typename DataType>
117 tnsr::ia<DataType, SpatialDim, Frame> two_index_constraint(
118  const tnsr::ia<DataType, SpatialDim, Frame>& d_gauge_function,
119  const tnsr::a<DataType, SpatialDim, Frame>& spacetime_normal_one_form,
120  const tnsr::A<DataType, SpatialDim, Frame>& spacetime_normal_vector,
121  const tnsr::II<DataType, SpatialDim, Frame>& inverse_spatial_metric,
122  const tnsr::AA<DataType, SpatialDim, Frame>& inverse_spacetime_metric,
123  const tnsr::aa<DataType, SpatialDim, Frame>& pi,
124  const tnsr::iaa<DataType, SpatialDim, Frame>& phi,
125  const tnsr::iaa<DataType, SpatialDim, Frame>& d_pi,
126  const tnsr::ijaa<DataType, SpatialDim, Frame>& d_phi,
127  const Scalar<DataType>& gamma2,
128  const tnsr::iaa<DataType, SpatialDim, Frame>&
129  three_index_constraint) noexcept;
130 
131 template <size_t SpatialDim, typename Frame, typename DataType>
133  gsl::not_null<tnsr::ia<DataType, SpatialDim, Frame>*> constraint,
134  const tnsr::ia<DataType, SpatialDim, Frame>& d_gauge_function,
135  const tnsr::a<DataType, SpatialDim, Frame>& spacetime_normal_one_form,
136  const tnsr::A<DataType, SpatialDim, Frame>& spacetime_normal_vector,
137  const tnsr::II<DataType, SpatialDim, Frame>& inverse_spatial_metric,
138  const tnsr::AA<DataType, SpatialDim, Frame>& inverse_spacetime_metric,
139  const tnsr::aa<DataType, SpatialDim, Frame>& pi,
140  const tnsr::iaa<DataType, SpatialDim, Frame>& phi,
141  const tnsr::iaa<DataType, SpatialDim, Frame>& d_pi,
142  const tnsr::ijaa<DataType, SpatialDim, Frame>& d_phi,
143  const Scalar<DataType>& gamma2,
144  const tnsr::iaa<DataType, SpatialDim, Frame>&
145  three_index_constraint) noexcept;
146 // @}
147 
148 // @{
149 /*!
150  * \brief Computes the generalized-harmonic 4-index constraint.
151  *
152  * \details Computes the independent components of the generalized-harmonic
153  * 4-index constraint. The constraint itself is given by Eq. (45) of
154  * \cite Lindblom2005qh,
155  * \f{eqnarray}{
156  * C_{ijab} = 2 \partial_{[i}\Phi_{j]ab},
157  * \f}
158  * where \f$\Phi_{iab} = \partial_i\psi_{ab}\f$. Because the constraint is
159  * antisymmetric on the two spatial indices, here we compute and store
160  * only the independent components of \f$C_{ijab}\f$. Specifically, we
161  * compute
162  * \f{eqnarray}{
163  * D_{iab} \equiv \frac{1}{2} \epsilon_{i}{}^{jk} C_{jkab}
164  * = \epsilon_{i}{}^{jk} \partial_j \Phi_{kab},
165  * \f}
166  * where \f$\epsilon_{ijk}\f$ is the flat-space Levi-Civita symbol,
167  * which is raised and lowered with the Kronecker delta.
168  * In terms
169  * of \f$D_{iab}\f$, the full 4-index constraint is
170  * \f{eqnarray}{
171  * C_{jkab} = \epsilon^{i}{}_{jk} D_{iab}.
172  * \f}
173  */
174 template <size_t SpatialDim, typename Frame, typename DataType>
175 tnsr::iaa<DataType, SpatialDim, Frame> four_index_constraint(
176  const tnsr::ijaa<DataType, SpatialDim, Frame>& d_phi) noexcept;
177 
178 template <size_t SpatialDim, typename Frame, typename DataType>
180  gsl::not_null<tnsr::iaa<DataType, SpatialDim, Frame>*> constraint,
181  const tnsr::ijaa<DataType, SpatialDim, Frame>& d_phi) noexcept;
182 // @}
183 
184 // @{
185 /*!
186  * \brief Computes the generalized-harmonic F constraint.
187  *
188  * \details Computes the generalized-harmonic F constraint
189  * [Eq. (43) of \cite Lindblom2005qh],
190  * \f{eqnarray}{
191  * {\cal F}_a &\equiv&
192  * \frac{1}{2} g_a^i \psi^{bc}\partial_i \Pi_{bc}
193  * - g^{ij} \partial_i \Pi_{ja}
194  * - g^{ij} t^b \partial_i \Phi_{jba}
195  * + \frac{1}{2} t_a \psi^{bc} g^{ij} \partial_i \Phi_{jbc}
196  * \nonumber \\ &&
197  * + t_a g^{ij} \partial_i H_j
198  * + g_a^i \Phi_{ijb} g^{jk}\Phi_{kcd} \psi^{bd} t^c
199  * - \frac{1}{2} g_a^i \Phi_{ijb} g^{jk}
200  * \Phi_{kcd} \psi^{cd} t^b
201  * \nonumber \\ &&
202  * - g_a^i t^b \partial_i H_b
203  * + g^{ij} \Phi_{icd} \Phi_{jba} \psi^{bc} t^d
204  * - \frac{1}{2} t_a g^{ij} g^{mn} \Phi_{imc} \Phi_{njd}\psi^{cd}
205  * \nonumber \\ &&
206  * - \frac{1}{4} t_a g^{ij}\Phi_{icd}\Phi_{jbe}
207  * \psi^{cb}\psi^{de}
208  * + \frac{1}{4} t_a \Pi_{cd} \Pi_{be}
209  * \psi^{cb}\psi^{de}
210  * - g^{ij} H_i \Pi_{ja}
211  * \nonumber \\ &&
212  * - t^b g^{ij} \Pi_{b i} \Pi_{ja}
213  * - \frac{1}{4} g_a^i \Phi_{icd} t^c t^d \Pi_{be}
214  * \psi^{be}
215  * + \frac{1}{2} t_a \Pi_{cd} \Pi_{be}\psi^{ce}
216  * t^d t^b
217  * \nonumber \\ &&
218  * + g_a^i \Phi_{icd} \Pi_{be} t^c t^b \psi^{de}
219  * - g^{ij}\Phi_{iba} t^b \Pi_{je} t^e
220  * - \frac{1}{2} g^{ij}\Phi_{icd} t^c t^d \Pi_{ja}
221  * \nonumber \\ &&
222  * - g^{ij} H_i \Phi_{jba} t^b
223  * + g_{a}^i \Phi_{icd} H_b \psi^{bc} t^d
224  * +\gamma_2\bigl(g^{id}{\cal C}_{ida}
225  * -\frac{1}{2} g_a^i\psi^{cd}{\cal C}_{icd}\bigr)
226  * \nonumber \\ &&
227  * + \frac{1}{2} t_a \Pi_{cd}\psi^{cd} H_b t^b
228  * - t_a g^{ij} \Phi_{ijc} H_d \psi^{cd}
229  * +\frac{1}{2} t_a g^{ij} H_i \Phi_{jcd}\psi^{cd},
230  * \f}
231  * where \f$H_a\f$ is the gauge function,
232  * \f$\psi_{ab}\f$ is the spacetime metric,
233  * \f$\Pi_{ab}=-t^c\partial_c \psi_{ab}\f$, and
234  * \f$\Phi_{iab} = \partial_i\psi_{ab}\f$; \f$t^a\f$ is the timelike unit
235  * normal vector to the spatial slice, \f$g^{ij}\f$ is the inverse spatial
236  * metric, and \f$g^b_c = \delta^b_c + t^b t_c\f$.
237  */
238 template <size_t SpatialDim, typename Frame, typename DataType>
239 tnsr::a<DataType, SpatialDim, Frame> f_constraint(
240  const tnsr::a<DataType, SpatialDim, Frame>& gauge_function,
241  const tnsr::ia<DataType, SpatialDim, Frame>& d_gauge_function,
242  const tnsr::a<DataType, SpatialDim, Frame>& spacetime_normal_one_form,
243  const tnsr::A<DataType, SpatialDim, Frame>& spacetime_normal_vector,
244  const tnsr::II<DataType, SpatialDim, Frame>& inverse_spatial_metric,
245  const tnsr::AA<DataType, SpatialDim, Frame>& inverse_spacetime_metric,
246  const tnsr::aa<DataType, SpatialDim, Frame>& pi,
247  const tnsr::iaa<DataType, SpatialDim, Frame>& phi,
248  const tnsr::iaa<DataType, SpatialDim, Frame>& d_pi,
249  const tnsr::ijaa<DataType, SpatialDim, Frame>& d_phi,
250  const Scalar<DataType>& gamma2,
251  const tnsr::iaa<DataType, SpatialDim, Frame>&
252  three_index_constraint) noexcept;
253 
254 template <size_t SpatialDim, typename Frame, typename DataType>
255 void f_constraint(
256  gsl::not_null<tnsr::a<DataType, SpatialDim, Frame>*> constraint,
257  const tnsr::a<DataType, SpatialDim, Frame>& gauge_function,
258  const tnsr::ia<DataType, SpatialDim, Frame>& d_gauge_function,
259  const tnsr::a<DataType, SpatialDim, Frame>& spacetime_normal_one_form,
260  const tnsr::A<DataType, SpatialDim, Frame>& spacetime_normal_vector,
261  const tnsr::II<DataType, SpatialDim, Frame>& inverse_spatial_metric,
262  const tnsr::AA<DataType, SpatialDim, Frame>& inverse_spacetime_metric,
263  const tnsr::aa<DataType, SpatialDim, Frame>& pi,
264  const tnsr::iaa<DataType, SpatialDim, Frame>& phi,
265  const tnsr::iaa<DataType, SpatialDim, Frame>& d_pi,
266  const tnsr::ijaa<DataType, SpatialDim, Frame>& d_phi,
267  const Scalar<DataType>& gamma2,
268  const tnsr::iaa<DataType, SpatialDim, Frame>&
269  three_index_constraint) noexcept;
270 // @}
271 
272 // @{
273 /*!
274  * \brief Computes the generalized-harmonic (unnormalized) constraint energy.
275  *
276  * \details Computes the generalized-harmonic unnormalized constraint energy
277  * [Eq. (53) of \cite Lindblom2005qh with \f$m^{ab}=\delta^{ab}\f$ and with each
278  * term in the sum scaled by an arbitrary coefficient],
279  * \f{eqnarray}{
280  * E & = & K_1 C_a C_a + K_2\left(F_a F_a
281  + C_{ia} C_{ja} g^{ij}\right)\nonumber\\
282  & + & K_3 C_{iab} C_{jab} g^{ij} + K_4 C_{ikab} C_{jlab}g^{ij} g^{kl}.
283  * \f}
284  * Here \f$C_a\f$ is the gauge constraint, \f$F_a\f$ is the f constraint,
285  * \f$C_{ia}\f$ is the two-index constraint, \f$C_{iab}\f$ is the
286  * three-index constraint, \f$C_{ikab}\f$ is the four-index constraint,
287  * \f$g^{ij}\f$ is the inverse spatial metric, and
288  * \f$K_1\f$, \f$K_2\f$, \f$K_3\f$, and \f$K_4\f$ are constant multipliers
289  * for each term that each default to a value of 1.0. Note that in this
290  * equation, spacetime indices \f$a,b\f$ are raised and lowered with
291  * the Kronecker delta.
292  *
293  * Also note that the argument `four_index_constraint` is a rank-3 tensor.
294  * This is because `GeneralizedHarmonic::four_index_constraint()` takes
295  * advantage of the antisymmetry of the four-index constraint's first two
296  * indices to only compute and return the independent
297  * components of \f$C_{ijab}\f$, which can be written as
298  * \f{eqnarray}{
299  * D_{iab} \equiv \frac{1}{2} \epsilon_{i}{}^{jk} C_{jkab},
300  * \f} where \f$\epsilon_{ijk}\f$ is the flat-space Levi-Civita tensor,
301  * whose inidces are raised and lowered with the Kronecker delta.
302  * The result is
303  * \f{eqnarray}{
304  * E & = & K_1 C_a C_a + K_2\left(F_a F_a
305  * + C_{ia} C_{ja} g^{ij}\right)\nonumber\\
306  * & + & K_3 C_{iab} C_{jab} g^{ij}\nonumber\\
307  & + & 2 K_4 g D_{iab} D_{jab} g^{ij},
308  * \f} where \f$g\f$ is the determinant of the spatial metric.
309  *
310  * To derive this expression for the constraint energy implemented here,
311  * Eq.~(53) of \cite Lindblom2005qh is
312  * \f{eqnarray}{
313  * S_{AB} dc^A dc^B &=&
314  * m^{ab}\Bigl[d F_ad F_b
315  * +g^{ij}\bigl(d C_{ia}d C_{jb}
316  * +g^{kl}m^{cd}d C_{ikac}d C_{jlbd}\bigr)
317  * \nonumber\\
318  * & + & \Lambda^2\bigl(d C_ad C_b
319  * +g^{ij}m^{cd}d C_{iac}d C_{jbd}\bigr)
320  * \Bigr].
321  * \f} Replace \f$dc^A\rightarrow c^A\f$ to get
322  * \f{eqnarray}{
323  * E&=&
324  * m^{ab}\Bigl[ F_a F_b
325  * +g^{ij}\bigl( C_{ia} C_{jb}
326  * +g^{kl}m^{cd} C_{ikac} C_{jlbd}\bigr)
327  * \nonumber\\
328  * & + & \Lambda^2\bigl( C_a C_b
329  * +g^{ij}m^{cd} C_{iac} C_{jbd}\bigr)
330  * \Bigr]\nonumber\\
331  * &=&
332  * m^{ab} F_a F_b
333  * +m^{ab}g^{ij} C_{ia} C_{jb}
334  * +m^{ab}g^{ij} g^{kl}m^{cd} C_{ikac} C_{jlbd}
335  * \nonumber\\
336  * & + & m^{ab}\Lambda^2 C_a C_b
337  * +m^{ab}\Lambda^2 g^{ij}m^{cd} C_{iac} C_{jbd}.
338  * \f} Here \f$m^{ab}\f$ is an arbitrary positive-definite matrix, and
339  * \f$\Lambda\f$ is an arbitrary real scalar.
340  * Choose \f$m^{ab} = \delta^{ab}\f$ but allow an arbitrary coefficient to be
341  * placed in front of each term. Then, absorb \f$\Lambda^2\f$ into one of
342  * these coefficients, to get
343  * \f{eqnarray}{ E &=& K_
344  * F\delta^{ab} F_a F_b +K_2\delta^{ab}g^{ij} C_{ia} C_{jb}
345  +K_4\delta^{ab}g^{ij}
346  * g^{kl}\delta^{cd} C_{ikac} C_{jlbd}
347  * \nonumber\\
348  * & + & K_1\delta^{ab} C_a C_b
349  * +K_3\delta^{ab} g^{ij}\delta^{cd} C_{iac} C_{jbd}.
350  * \f}
351  * Adopting a Euclidean norm for the constraint space (i.e., choosing to raise
352  and
353  * lower spacetime indices with Kronecker deltas) gives
354  * \f{eqnarray}{ E &=& K_ F
355  * F_a F_a +K_2g^{ij} C_{ia} C_{ja} +K_4 g^{ij} g^{kl} C_{ikac} C_{jlac}
356  * \nonumber\\
357  * & + & K_1 C_a C_a
358  * +K_3g^{ij} C_{iac} C_{jac}.
359  * \f} The two-index constraint and f constraint can be viewed as the
360  * time and space components of a combined spacetime constraint. So next
361  * choose
362  * \f$K_ F=K_2\f$, giving \f{eqnarray}{ E&=& K_1 C_a C_a + K_2\left(F_a F_a
363  * + C_{ia} C_{ja} g^{ij}\right)\nonumber\\
364  * & + & K_3 C_{iab} C_{jab} g^{ij}
365  * + K_4 C_{ikab} C_{jlab}g^{ij} g^{kl}.
366  * \f}
367  *
368  * Note that \f$C_{ikab}\f$ is antisymmetric on the first two indices. Next,
369  * replace the four-index constraint \f$C_{ijab}\f$ with \f$D_{iab}\f$, which
370  * contains the independent components of \f$C_{ijab}\f$. Specifically,
371  * \f{eqnarray}{
372  * D_{iab} \equiv \frac{1}{2} \epsilon_{i}{}^{jk} C_{jkab}.
373  * \f} The inverse relationship is
374  * \f{eqnarray}{
375  * C_{jkab} = \epsilon^{i}{}_{jk} D_{iab},
376  * \f} where \f$\epsilon_{ijk}\f$ is the flat-space Levi-Civita tensor, whose
377  * indices are raised and lowered with the Kronecker delta. Inserting this
378  relation
379  * to replace \f$C_{jkab}\f$ with \f$D_{iab}\f$ gives \f{eqnarray}{ E &=& K_1
380  C_a
381  * C_a + K_2\left(F_a F_a
382  * + C_{ia} C_{ja} g^{ij}\right)\nonumber\\
383  * & + & K_3 C_{iab} C_{jab} g^{ij}\nonumber\\
384  * & + & K_4 D_{mab} D_{nab} \epsilon^{m}{}_{ik}
385  * \epsilon^{n}{}_{jl} g^{ij} g^{kl}. \f}
386  *
387  * There's a subtle point here: \f$g^{ij}\f$ is the inverse spatial metric,
388  which
389  * is not necessarily flat. But \f$\epsilon^{i}{}_{jk}\f$ is the flat space
390  * Levi-Civita tensor. In order to raise and lower indices of the Levi-Civita
391  * tensor with the inverse spatial metrics, put in the appropriate factors of
392  * \f$\sqrt{g}\f$, where \f$g\f$ is the metric determinant, to make the
393  * curved-space Levi-Civita tensor compatible with \f$g^{ij}\f$. Let
394  * \f$\varepsilon^{ijk}\f$ represent the curved space Levi-Civita tensor
395  compatible
396  * with \f$g^{ij}\f$: \f{eqnarray}{
397  * \varepsilon^{mik} = g^{-1/2} \epsilon^{mik}\\
398  * \varepsilon_{mik} = g^{1/2} \epsilon_{mik}.
399  * \f} Then we can write the constraint energy as
400  * \f{eqnarray}{
401  * E &=& K_1 C_a C_a + K_2\left(F_a F_a
402  * + C_{ia} C_{ja} g^{ij}\right)\nonumber\\
403  * & + & K_3 C_{iab} C_{jab} g^{ij}\nonumber\\
404  * & + & K_4 D_{mab} D_{nab} g g^{-1/2}\epsilon^{m}{}_{ik}
405  * g^{-1/2}\epsilon^{n}{}_{jl} g^{ij} g^{kl}. \f} The factors of
406  * \f$g^{-1/2}\f$ make the Levi-Civita tensor compatible with \f$g^{ij}\f$.
407  * Swapping which summed indices are raised and which are lowered gives
408  * \f{eqnarray}{ E &=& K_1 C_a
409  C_a +
410  * K_2\left(F_a F_a
411  * + C_{ia} C_{ja} g^{ij}\right)\nonumber\\
412  * & + & K_3 C_{iab} C_{jab} g^{ij}\nonumber\\
413  * & + & K_4 D_{mab} D_{nab} g g^{-1/2}\epsilon^{mik}
414  g^{-1/2}\epsilon^{njl}
415  * g_{ij} g_{kl}, \f} or \f{eqnarray}{ E &=& K_1 C_a C_a + K_2\left(F_a F_a
416  * + C_{ia} C_{ja} g^{ij}\right)\nonumber\\
417  * & + & K_3 C_{iab} C_{jab} g^{ij}\nonumber\\
418  * & + & K_4 D_{mab} D_{nab} g \varepsilon^{mik} \varepsilon^{njl} g_{ij}
419  * g_{kl}, \f} or, reversing up and down repeated indices again,
420  * \f{eqnarray}{ E
421  * &=& K_1 C_a C_a + K_2\left(F_a F_a
422  * + C_{ia} C_{ja} g^{ij}\right)\nonumber\\
423  * & + & K_3 C_{iab} C_{jab} g^{ij}\nonumber\\
424  * & + & K_4 D_{mab} D_{nab} g \varepsilon^{m}{}_{ik}
425  \varepsilon^{n}{}_{jl}
426  * g^{ij} g^{kl}. \f}
427  *
428  * The metric raises and lowers the indices of \f$\varepsilon^{ijk}\f$,
429  * so this can
430  * be written as \f{eqnarray}{ E &=& K_1 C_a C_a + K_2\left(F_a F_a
431  * + C_{ia} C_{ja} g^{ij}\right)\nonumber\\
432  * & + & K_3 C_{iab} C_{jab} g^{ij}\nonumber\\
433  * & + & K_4 g D_{mab} D^{n}{}_{ab} \varepsilon^{mjl} \varepsilon_{njl}.
434  * \f}
435  *
436  * Now, in flat space (Eq. (1.23) of \cite ThorneBlandford2017),
437  * \f{eqnarray}{
438  * \epsilon^{mjl} \epsilon_{njl} = \delta^{mj}_{nj} = \delta^m_n \delta^j_j -
439  * \delta^m_j \delta^j_n = 2 \delta^m_n. \f} But this holds for curved space
440  * as well: multiply the left hand side by \f$1 = g^{1/2} g^{-1/2}\f$ to get
441  * \f{eqnarray}{
442  * g^{-1/2}\epsilon^{mjl} g^{1/2}\epsilon_{njl} = \varepsilon^{mjl}
443  * \varepsilon_{njl} = \delta^{mj}_{nj} = \delta^m_n \delta^j_j - \delta^m_j
444  * \delta^j_n = 2 \delta^m_n. \f} So the constraint energy is \f{eqnarray}{ E
445  &=&
446  * K_1 C_a 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  * & + & 2 K_4 D_{mab} D^{n}{}_{ab} \delta^m_n.
450  * \f}
451  * Simplifying gives the formula implemented here:
452  * \f{eqnarray}{
453  * E &=& K_1 C_a C_a + K_2\left(F_a F_a
454  * + C_{ia} C_{ja} g^{ij}\right)\nonumber\\
455  * & + & K_3 C_{iab} C_{jab} g^{ij}\nonumber\\
456  * & + & 2 K_4 g D_{iab} D_{jab} g^{ij}.
457  * \f}
458  */
459 template <size_t SpatialDim, typename Frame, typename DataType>
461  const tnsr::a<DataType, SpatialDim, Frame>& gauge_constraint,
462  const tnsr::a<DataType, SpatialDim, Frame>& f_constraint,
463  const tnsr::ia<DataType, SpatialDim, Frame>& two_index_constraint,
464  const tnsr::iaa<DataType, SpatialDim, Frame>& three_index_constraint,
465  const tnsr::iaa<DataType, SpatialDim, Frame>& four_index_constraint,
466  const tnsr::II<DataType, SpatialDim, Frame>& inverse_spatial_metric,
467  const Scalar<DataType>& spatial_metric_determinant,
468  double gauge_constraint_multiplier = 1.0,
469  double two_index_constraint_multiplier = 1.0,
470  double three_index_constraint_multiplier = 1.0,
471  double four_index_constraint_multiplier = 1.0) noexcept;
472 
473 template <size_t SpatialDim, typename Frame, typename DataType>
474 void constraint_energy(
475  gsl::not_null<Scalar<DataType>*> energy,
476  const tnsr::a<DataType, SpatialDim, Frame>& gauge_constraint,
477  const tnsr::a<DataType, SpatialDim, Frame>& f_constraint,
478  const tnsr::ia<DataType, SpatialDim, Frame>& two_index_constraint,
479  const tnsr::iaa<DataType, SpatialDim, Frame>& three_index_constraint,
480  const tnsr::iaa<DataType, SpatialDim, Frame>& four_index_constraint,
481  const tnsr::II<DataType, SpatialDim, Frame>& inverse_spatial_metric,
482  const Scalar<DataType>& spatial_metric_determinant,
483  double gauge_constraint_multiplier = 1.0,
484  double two_index_constraint_multiplier = 1.0,
485  double three_index_constraint_multiplier = 1.0,
486  double four_index_constraint_multiplier = 1.0) noexcept;
487 // @}
488 } // namespace GeneralizedHarmonic
Implementations from the Guideline Support Library.
Definition: ConservativeFromPrimitive.hpp:10
void three_index_constraint(gsl::not_null< tnsr::iaa< DataType, SpatialDim, Frame > *> constraint, const tnsr::iaa< DataType, SpatialDim, Frame > &d_spacetime_metric, const tnsr::iaa< DataType, SpatialDim, Frame > &phi) noexcept
Computes the generalized-harmonic 3-index constraint.
tnsr::a< DataType, SpatialDim, Frame > f_constraint(const tnsr::a< DataType, SpatialDim, Frame > &gauge_function, const tnsr::ia< DataType, SpatialDim, Frame > &d_gauge_function, const tnsr::a< DataType, SpatialDim, Frame > &spacetime_normal_one_form, const tnsr::A< DataType, SpatialDim, Frame > &spacetime_normal_vector, const tnsr::II< DataType, SpatialDim, Frame > &inverse_spatial_metric, const tnsr::AA< DataType, SpatialDim, Frame > &inverse_spacetime_metric, const tnsr::aa< DataType, SpatialDim, Frame > &pi, const tnsr::iaa< DataType, SpatialDim, Frame > &phi, const tnsr::iaa< DataType, SpatialDim, Frame > &d_pi, const tnsr::ijaa< DataType, SpatialDim, Frame > &d_phi, const Scalar< DataType > &gamma2, const tnsr::iaa< DataType, SpatialDim, Frame > &three_index_constraint) noexcept
Computes the generalized-harmonic F constraint.
Definition: Constraints.cpp:1117
tnsr::A< DataType, SpatialDim, Frame > spacetime_normal_vector(const Scalar< DataType > &lapse, const tnsr::I< DataType, SpatialDim, Frame > &shift) noexcept
Computes spacetime normal vector from lapse and shift.
Definition: ComputeSpacetimeQuantities.cpp:185
void two_index_constraint(gsl::not_null< tnsr::ia< DataType, SpatialDim, Frame > *> constraint, const tnsr::ia< DataType, SpatialDim, Frame > &d_gauge_function, const tnsr::a< DataType, SpatialDim, Frame > &spacetime_normal_one_form, const tnsr::A< DataType, SpatialDim, Frame > &spacetime_normal_vector, const tnsr::II< DataType, SpatialDim, Frame > &inverse_spatial_metric, const tnsr::AA< DataType, SpatialDim, Frame > &inverse_spacetime_metric, const tnsr::aa< DataType, SpatialDim, Frame > &pi, const tnsr::iaa< DataType, SpatialDim, Frame > &phi, const tnsr::iaa< DataType, SpatialDim, Frame > &d_pi, const tnsr::ijaa< DataType, SpatialDim, Frame > &d_phi, const Scalar< DataType > &gamma2, const tnsr::iaa< DataType, SpatialDim, Frame > &three_index_constraint) noexcept
Computes the generalized-harmonic 2-index constraint.
tnsr::AA< DataType, SpatialDim, Frame > inverse_spacetime_metric(const Scalar< DataType > &lapse, const tnsr::I< DataType, SpatialDim, Frame > &shift, const tnsr::II< DataType, SpatialDim, Frame > &inverse_spatial_metric) noexcept
Compute inverse spacetime metric from inverse spatial metric, lapse and shift.
void gauge_constraint(gsl::not_null< tnsr::a< DataType, SpatialDim, Frame > *> constraint, const tnsr::a< DataType, SpatialDim, Frame > &gauge_function, const tnsr::a< DataType, SpatialDim, Frame > &spacetime_normal_one_form, const tnsr::A< DataType, SpatialDim, Frame > &spacetime_normal_vector, const tnsr::II< DataType, SpatialDim, Frame > &inverse_spatial_metric, const tnsr::AA< DataType, SpatialDim, Frame > &inverse_spacetime_metric, const tnsr::aa< DataType, SpatialDim, Frame > &pi, const tnsr::iaa< DataType, SpatialDim, Frame > &phi) noexcept
Computes the generalized-harmonic gauge constraint.
tnsr::aa< DataType, SpatialDim, Frame > pi(const Scalar< DataType > &lapse, const Scalar< DataType > &dt_lapse, const tnsr::I< DataType, SpatialDim, Frame > &shift, const tnsr::I< DataType, SpatialDim, Frame > &dt_shift, const tnsr::ii< DataType, SpatialDim, Frame > &spatial_metric, const tnsr::ii< DataType, SpatialDim, Frame > &dt_spatial_metric, const tnsr::iaa< DataType, SpatialDim, Frame > &phi) noexcept
Computes the conjugate momentum of the spacetime metric .
Definition: ComputeGhQuantities.cpp:58
Type aliases to construct common Tensors.
Definition: TypeAliases.hpp:31
void four_index_constraint(gsl::not_null< tnsr::iaa< DataType, SpatialDim, Frame > *> constraint, const tnsr::ijaa< DataType, SpatialDim, Frame > &d_phi) noexcept
Computes the generalized-harmonic 4-index constraint.
tnsr::iaa< DataType, SpatialDim, Frame > phi(const Scalar< DataType > &lapse, const tnsr::i< DataType, SpatialDim, Frame > &deriv_lapse, const tnsr::I< DataType, SpatialDim, Frame > &shift, const tnsr::iJ< DataType, SpatialDim, Frame > &deriv_shift, const tnsr::ii< DataType, SpatialDim, Frame > &spatial_metric, const tnsr::ijj< DataType, SpatialDim, Frame > &deriv_spatial_metric) noexcept
Computes the auxiliary variable used by the generalized harmonic formulation of Einstein&#39;s equations...
Definition: ComputeGhQuantities.cpp:22
Items related to evolving the first-order generalized harmonic system.
Definition: Characteristics.cpp:25
Indicates the Frame that a TensorIndexType is in.
Definition: IndexType.hpp:36
Defines a list of useful type aliases for tensors.
Scalar< DataType > constraint_energy(const tnsr::a< DataType, SpatialDim, Frame > &gauge_constraint, const tnsr::a< DataType, SpatialDim, Frame > &f_constraint, const tnsr::ia< DataType, SpatialDim, Frame > &two_index_constraint, const tnsr::iaa< DataType, SpatialDim, Frame > &three_index_constraint, const tnsr::iaa< DataType, SpatialDim, Frame > &four_index_constraint, const tnsr::II< DataType, SpatialDim, Frame > &inverse_spatial_metric, const Scalar< DataType > &spatial_metric_determinant, double gauge_constraint_multiplier, double two_index_constraint_multiplier, double three_index_constraint_multiplier, double four_index_constraint_multiplier) noexcept
Computes the generalized-harmonic (unnormalized) constraint energy.
Definition: Constraints.cpp:1228
Tensor< T, Symmetry<>, index_list<> > Scalar
Scalar type.
Definition: TypeAliases.hpp:21
Require a pointer to not be a nullptr
Definition: ConservativeFromPrimitive.hpp:12