Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 :
8 : #include "DataStructures/DataBox/Tag.hpp"
9 : #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
10 : #include "DataStructures/Tensor/Tensor.hpp"
11 : #include "DataStructures/Variables.hpp"
12 : #include "Domain/FaceNormal.hpp"
13 : #include "Evolution/Systems/CurvedScalarWave/Constraints.hpp"
14 : #include "Evolution/Systems/CurvedScalarWave/TagsDeclarations.hpp"
15 : #include "PointwiseFunctions/GeneralRelativity/TagsDeclarations.hpp"
16 : #include "Utilities/TMPL.hpp"
17 :
18 : namespace CurvedScalarWave {
19 : /// @{
20 : /*!
21 : * \brief Compute the characteristic speeds for the scalar wave system in curved
22 : * spacetime.
23 : *
24 : * Computes the speeds as described in "Optimal constraint projection for
25 : * hyperbolic evolution systems" by Holst et. al \cite Holst2004wt
26 : * [see text following Eq. (32)]. The characteristic fields' names used here
27 : * are similar to the paper:
28 : *
29 : * \f{align*}
30 : * \mathrm{SpECTRE} && \mathrm{Holst} \\
31 : * v^{\hat \psi} && Z^1 \\
32 : * v^{\hat 0}_{i} && Z^{2}_{i} \\
33 : * v^{\hat \pm} && u^{1\pm}
34 : * \f}
35 : *
36 : * The corresponding characteristic speeds \f$\lambda\f$ are given in the text
37 : * following Eq. (38) of \cite Holst2004wt :
38 : *
39 : * \f{align*}
40 : * \lambda_{\hat \psi} =& -(1 + \gamma_1) n_k N^k \\
41 : * \lambda_{\hat 0} =& -n_k N^k \\
42 : * \lambda_{\hat \pm} =& -n_k N^k \pm N
43 : * \f}
44 : *
45 : * where \f$n_k\f$ is the unit normal to the surface.
46 : */
47 : template <size_t SpatialDim>
48 1 : std::array<DataVector, 4> characteristic_speeds(
49 : const Scalar<DataVector>& gamma_1, const Scalar<DataVector>& lapse,
50 : const tnsr::I<DataVector, SpatialDim, Frame::Inertial>& shift,
51 : const tnsr::i<DataVector, SpatialDim, Frame::Inertial>&
52 : unit_normal_one_form);
53 :
54 : template <size_t SpatialDim>
55 1 : void characteristic_speeds(
56 : gsl::not_null<std::array<DataVector, 4>*> char_speeds,
57 : const Scalar<DataVector>& gamma_1, const Scalar<DataVector>& lapse,
58 : const tnsr::I<DataVector, SpatialDim, Frame::Inertial>& shift,
59 : const tnsr::i<DataVector, SpatialDim, Frame::Inertial>&
60 : unit_normal_one_form);
61 :
62 : template <size_t SpatialDim>
63 1 : void characteristic_speeds(
64 : gsl::not_null<tnsr::a<DataVector, 3, Frame::Inertial>*> char_speeds,
65 : const Scalar<DataVector>& gamma_1, const Scalar<DataVector>& lapse,
66 : const tnsr::I<DataVector, SpatialDim, Frame::Inertial>& shift,
67 : const tnsr::i<DataVector, SpatialDim, Frame::Inertial>&
68 : unit_normal_one_form);
69 :
70 : template <size_t SpatialDim>
71 0 : struct CharacteristicSpeedsCompute : Tags::CharacteristicSpeeds<SpatialDim>,
72 : db::ComputeTag {
73 0 : using base = Tags::CharacteristicSpeeds<SpatialDim>;
74 0 : using return_type = typename base::type;
75 0 : using argument_tags = tmpl::list<
76 : Tags::ConstraintGamma1, gr::Tags::Lapse<DataVector>,
77 : gr::Tags::Shift<DataVector, SpatialDim>,
78 : ::Tags::Normalized<domain::Tags::UnnormalizedFaceNormal<SpatialDim>>>;
79 :
80 0 : static constexpr void function(
81 : gsl::not_null<return_type*> result, const Scalar<DataVector>& gamma_1,
82 : const Scalar<DataVector>& lapse,
83 : const tnsr::I<DataVector, SpatialDim, Frame::Inertial>& shift,
84 : const tnsr::i<DataVector, SpatialDim, Frame::Inertial>&
85 : unit_normal_one_form) {
86 : characteristic_speeds<SpatialDim>(result, gamma_1, lapse, shift,
87 : unit_normal_one_form);
88 : }
89 : };
90 : /// @}
91 :
92 : /// @{
93 : /*!
94 : * \brief Computes characteristic fields from evolved fields
95 : *
96 : * \ref CharacteristicFieldsCompute and
97 : * \ref EvolvedFieldsFromCharacteristicFieldsCompute convert between
98 : * characteristic and evolved fields for the scalar-wave system in curved
99 : * spacetime.
100 : *
101 : * \ref CharacteristicFieldsCompute computes
102 : * characteristic fields as described in "Optimal constraint projection for
103 : * hyperbolic evolution systems" by Holst et. al \cite Holst2004wt .
104 : * Their names used here differ from this paper:
105 : *
106 : * \f{align*}
107 : * \mathrm{SpECTRE} && \mathrm{Holst} \\
108 : * v^{\hat \psi} && Z^1 \\
109 : * v^{\hat 0}_{i} && Z^{2}_{i} \\
110 : * v^{\hat \pm} && u^{1\pm}
111 : * \f}
112 : *
113 : * The characteristic fields \f$u\f$ are given in terms of the evolved fields by
114 : * Eq. (33) - (35) of \cite Holst2004wt, respectively:
115 : *
116 : * \f{align*}
117 : * v^{\hat \psi} =& \psi \\
118 : * v^{\hat 0}_{i} =& (\delta^k_i - n_i n^k) \Phi_{k} := P^k_i \Phi_{k} \\
119 : * v^{\hat \pm} =& \Pi \pm n^i \Phi_{i} - \gamma_2\psi
120 : * \f}
121 : *
122 : * where \f$\psi\f$ is the scalar field, \f$\Pi\f$ and \f$\Phi_{i}\f$ are
123 : * evolved fields introduced by first derivatives of \f$\psi\f$, \f$\gamma_2\f$
124 : * is a constraint damping parameter, and \f$n_k\f$ is the unit normal to the
125 : * surface.
126 : *
127 : * \ref EvolvedFieldsFromCharacteristicFieldsCompute computes evolved fields
128 : * \f$w\f$ in terms of the characteristic fields. This uses the inverse of
129 : * above relations (c.f. Eq. (36) - (38) of \cite Holst2004wt ):
130 : *
131 : * \f{align*}
132 : * \psi =& v^{\hat \psi}, \\
133 : * \Pi =& \frac{1}{2}(v^{\hat +} + v^{\hat -}) + \gamma_2 v^{\hat \psi}, \\
134 : * \Phi_{i} =& \frac{1}{2}(v^{\hat +} - v^{\hat -}) n_i + v^{\hat 0}_{i}.
135 : * \f}
136 : *
137 : * The corresponding characteristic speeds \f$\lambda\f$ are computed by
138 : * \ref CharacteristicSpeedsCompute .
139 : */
140 : template <size_t SpatialDim>
141 : Variables<
142 : tmpl::list<Tags::VPsi, Tags::VZero<SpatialDim>, Tags::VPlus, Tags::VMinus>>
143 1 : characteristic_fields(
144 : const Scalar<DataVector>& gamma_2, const Scalar<DataVector>& psi,
145 : const Scalar<DataVector>& pi,
146 : const tnsr::i<DataVector, SpatialDim, Frame::Inertial>& phi,
147 : const tnsr::i<DataVector, SpatialDim, Frame::Inertial>&
148 : unit_normal_one_form,
149 : const tnsr::I<DataVector, SpatialDim, Frame::Inertial>& unit_normal_vector);
150 :
151 : template <size_t SpatialDim>
152 1 : void characteristic_fields(
153 : gsl::not_null<Variables<tmpl::list<Tags::VPsi, Tags::VZero<SpatialDim>,
154 : Tags::VPlus, Tags::VMinus>>*>
155 : char_fields,
156 : const Scalar<DataVector>& gamma_2, const Scalar<DataVector>& psi,
157 : const Scalar<DataVector>& pi,
158 : const tnsr::i<DataVector, SpatialDim, Frame::Inertial>& phi,
159 : const tnsr::i<DataVector, SpatialDim, Frame::Inertial>&
160 : unit_normal_one_form,
161 : const tnsr::I<DataVector, SpatialDim, Frame::Inertial>& unit_normal_vector);
162 :
163 : template <size_t SpatialDim>
164 1 : void characteristic_fields(
165 : const gsl::not_null<Scalar<DataVector>*>& v_psi,
166 : const gsl::not_null<tnsr::i<DataVector, SpatialDim, Frame::Inertial>*>&
167 : v_zero,
168 : const gsl::not_null<Scalar<DataVector>*>& v_plus,
169 : const gsl::not_null<Scalar<DataVector>*>& v_minus,
170 : const Scalar<DataVector>& gamma_2, const Scalar<DataVector>& psi,
171 : const Scalar<DataVector>& pi,
172 : const tnsr::i<DataVector, SpatialDim, Frame::Inertial>& phi,
173 : const tnsr::i<DataVector, SpatialDim, Frame::Inertial>&
174 : unit_normal_one_form,
175 : const tnsr::I<DataVector, SpatialDim, Frame::Inertial>& unit_normal_vector);
176 :
177 : template <size_t SpatialDim>
178 0 : struct CharacteristicFieldsCompute : Tags::CharacteristicFields<SpatialDim>,
179 : db::ComputeTag {
180 0 : using base = Tags::CharacteristicFields<SpatialDim>;
181 0 : using return_type = typename base::type;
182 0 : using argument_tags = tmpl::list<
183 : Tags::ConstraintGamma2,
184 : gr::Tags::InverseSpatialMetric<DataVector, SpatialDim>, Tags::Psi,
185 : Tags::Pi, Tags::Phi<SpatialDim>,
186 : ::Tags::Normalized<domain::Tags::UnnormalizedFaceNormal<SpatialDim>>>;
187 :
188 0 : static constexpr void function(
189 : gsl::not_null<return_type*> result, const Scalar<DataVector>& gamma_2,
190 : const tnsr::II<DataVector, SpatialDim, Frame::Inertial>&
191 : inverse_spatial_metric,
192 : const Scalar<DataVector>& psi, const Scalar<DataVector>& pi,
193 : const tnsr::i<DataVector, SpatialDim, Frame::Inertial>& phi,
194 : const tnsr::i<DataVector, SpatialDim, Frame::Inertial>&
195 : unit_normal_one_form) {
196 : const auto unit_normal_vector = tenex::evaluate<ti::I>(
197 : inverse_spatial_metric(ti::I, ti::J) * unit_normal_one_form(ti::j));
198 : characteristic_fields<SpatialDim>(result, gamma_2, psi, pi, phi,
199 : unit_normal_one_form, unit_normal_vector);
200 : }
201 : };
202 : /// @}
203 :
204 : /// @{
205 : /*!
206 : * \brief For expressions used here to compute evolved fields from
207 : * characteristic ones, see \ref CharacteristicFieldsCompute.
208 : */
209 : template <size_t SpatialDim>
210 : Variables<tmpl::list<Tags::Psi, Tags::Pi, Tags::Phi<SpatialDim>>>
211 1 : evolved_fields_from_characteristic_fields(
212 : const Scalar<DataVector>& gamma_2, const Scalar<DataVector>& v_psi,
213 : const tnsr::i<DataVector, SpatialDim, Frame::Inertial>& v_zero,
214 : const Scalar<DataVector>& v_plus, const Scalar<DataVector>& v_minus,
215 : const tnsr::i<DataVector, SpatialDim, Frame::Inertial>&
216 : unit_normal_one_form);
217 :
218 : template <size_t SpatialDim>
219 1 : void evolved_fields_from_characteristic_fields(
220 : gsl::not_null<
221 : Variables<tmpl::list<Tags::Psi, Tags::Pi, Tags::Phi<SpatialDim>>>*>
222 : evolved_fields,
223 : const Scalar<DataVector>& gamma_2, const Scalar<DataVector>& v_psi,
224 : const tnsr::i<DataVector, SpatialDim, Frame::Inertial>& v_zero,
225 : const Scalar<DataVector>& v_plus, const Scalar<DataVector>& v_minus,
226 : const tnsr::i<DataVector, SpatialDim, Frame::Inertial>&
227 : unit_normal_one_form);
228 :
229 : template <size_t SpatialDim>
230 1 : void evolved_fields_from_characteristic_fields(
231 : gsl::not_null<Scalar<DataVector>*> psi,
232 : gsl::not_null<Scalar<DataVector>*> pi,
233 : gsl::not_null<tnsr::i<DataVector, SpatialDim, Frame::Inertial>*> phi,
234 : const Scalar<DataVector>& gamma_2, const Scalar<DataVector>& v_psi,
235 : const tnsr::i<DataVector, SpatialDim, Frame::Inertial>& v_zero,
236 : const Scalar<DataVector>& v_plus, const Scalar<DataVector>& v_minus,
237 : const tnsr::i<DataVector, SpatialDim, Frame::Inertial>&
238 : unit_normal_one_form);
239 :
240 : template <size_t SpatialDim>
241 0 : struct EvolvedFieldsFromCharacteristicFieldsCompute
242 : : Tags::EvolvedFieldsFromCharacteristicFields<SpatialDim>,
243 : db::ComputeTag {
244 0 : using base = Tags::EvolvedFieldsFromCharacteristicFields<SpatialDim>;
245 0 : using return_type = typename base::type;
246 0 : using argument_tags = tmpl::list<
247 : Tags::ConstraintGamma2, Tags::VPsi, Tags::VZero<SpatialDim>, Tags::VPlus,
248 : Tags::VMinus,
249 : ::Tags::Normalized<domain::Tags::UnnormalizedFaceNormal<SpatialDim>>>;
250 :
251 0 : static constexpr void function(
252 : gsl::not_null<return_type*> result, const Scalar<DataVector>& gamma_2,
253 : const Scalar<DataVector>& v_psi,
254 : const tnsr::i<DataVector, SpatialDim, Frame::Inertial>& v_zero,
255 : const Scalar<DataVector>& v_plus, const Scalar<DataVector>& v_minus,
256 : const tnsr::i<DataVector, SpatialDim, Frame::Inertial>&
257 : unit_normal_one_form) {
258 : evolved_fields_from_characteristic_fields<SpatialDim>(
259 : result, gamma_2, v_psi, v_zero, v_plus, v_minus, unit_normal_one_form);
260 : }
261 : };
262 : /// @}
263 :
264 : namespace Tags {
265 : /*!
266 : * \brief Computes the largest magnitude of the characteristic speeds.
267 : *
268 : * \details Returns the magnitude of the largest characteristic speed
269 : * along any direction at a given point in space, considering all
270 : * characteristic fields. This is useful, for e.g., in computing the
271 : * Courant factor. The coordinate characteristic speeds for this system are
272 : * \f$\{-(1+\gamma_1)n_k N^k, -n_k N^k, -n_k N^k \pm N\}\f$. At any point
273 : * in space, these are maximized when the normal vector is parallel to the
274 : * shift vector, i.e. \f$ n^j = N^j / \sqrt{N^i N_i}\f$, and \f$ n_k N^k
275 : * = g_{jk} N^j N^k / \sqrt{N^i N_i} = \sqrt{N^i N_i} =\f$
276 : * `magnitude(shift, spatial_metric)`. The maximum characteristic speed
277 : * is therefore calculated as \f$ \rm{max}(\vert 1+\gamma_1\vert\sqrt{N^i
278 : * N_i},\, \sqrt{N^i N_i}+\vert N\vert) \f$.
279 : */
280 : template <size_t SpatialDim>
281 1 : struct ComputeLargestCharacteristicSpeed : LargestCharacteristicSpeed,
282 : db::ComputeTag {
283 0 : using argument_tags =
284 : tmpl::list<Tags::ConstraintGamma1, gr::Tags::Lapse<DataVector>,
285 : gr::Tags::Shift<DataVector, SpatialDim>,
286 : gr::Tags::SpatialMetric<DataVector, SpatialDim>>;
287 0 : using return_type = double;
288 0 : using base = LargestCharacteristicSpeed;
289 0 : static void function(
290 : const gsl::not_null<double*> max_speed, const Scalar<DataVector>& gamma_1,
291 : const Scalar<DataVector>& lapse,
292 : const tnsr::I<DataVector, SpatialDim, Frame::Inertial>& shift,
293 : const tnsr::ii<DataVector, SpatialDim, Frame::Inertial>& spatial_metric);
294 : };
295 : } // namespace Tags
296 : } // namespace CurvedScalarWave
|
