Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cmath>
7 : #include <cstddef>
8 : #include <type_traits>
9 :
10 : #include "DataStructures/ComplexDataVector.hpp"
11 : #include "DataStructures/DataBox/DataBox.hpp"
12 : #include "DataStructures/Tags.hpp"
13 : #include "DataStructures/Tensor/Tensor.hpp"
14 : #include "DataStructures/Tensor/TypeAliases.hpp"
15 : #include "Evolution/Systems/Cce/IntegrandInputSteps.hpp"
16 : #include "Evolution/Systems/Cce/OptionTags.hpp"
17 : #include "Evolution/Systems/Cce/Tags.hpp"
18 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
19 : #include "NumericalAlgorithms/SpinWeightedSphericalHarmonics/SwshDerivatives.hpp"
20 : #include "NumericalAlgorithms/SpinWeightedSphericalHarmonics/SwshTags.hpp"
21 : #include "Utilities/ForceInline.hpp"
22 : #include "Utilities/Gsl.hpp"
23 : #include "Utilities/TypeTraits.hpp"
24 : #include "Utilities/TypeTraits/IsA.hpp"
25 :
26 : namespace Cce {
27 : namespace detail {
28 : // Precomputation routines for supplying the additional quantities necessary
29 : // to correct the output of the angular derivative routines from the angular
30 : // derivatives evaluated at constant numerical coordinates (which is what is
31 : // returned after the libsharp evaluation) to the angular derivatives at
32 : // constant Bondi radius (which is what appears in the literature equations
33 : // and is simple to combine to obtain the hypersurface integrands).
34 : //
35 : // Warning: this 'on demand' template is a way of taking advantage of the blaze
36 : // expression templates in a generic, modular way. However, this can be
37 : // dangerous. The returned value MUST be a blaze expression template directly,
38 : // and not a wrapper type (like `SpinWeighted`). Otherwise, some information is
39 : // lost on the stack and the expression template is corrupted. So, in these 'on
40 : // demand' returns, the arguments must be fully unpacked to vector types.
41 : //
42 : // The `Select` operand is a `std::bool_constant` that ensures mutual
43 : // exclusivity of the template specializations.
44 : template <typename Tag, typename SpinConstant, typename Select>
45 : struct OnDemandInputsForSwshJacobianImpl;
46 :
47 : // default to just retrieving it from the box if not providing an expression
48 : // template shortcut
49 : template <typename Tag>
50 : struct OnDemandInputsForSwshJacobianImpl<
51 : Tags::Dy<Tag>, std::integral_constant<int, Tag::type::type::spin>,
52 : std::bool_constant<not tt::is_a_v<::Tags::Multiplies, Tag> and
53 : not tt::is_a_v<Tags::Dy, Tag>>> {
54 : template <typename DataBoxTagList>
55 : SPECTRE_ALWAYS_INLINE decltype(auto) operator()(
56 : const db::DataBox<DataBoxTagList>& box) {
57 : return get(db::get<Tags::Dy<Tag>>(box)).data();
58 : }
59 : };
60 :
61 : // default to retrieving from the box if the requested tag is a second
62 : // derivative without an additional evaluation channel
63 : template <typename Tag>
64 : struct OnDemandInputsForSwshJacobianImpl<
65 : Tags::Dy<Tags::Dy<Tag>>, std::integral_constant<int, Tag::type::type::spin>,
66 : std::bool_constant<not tt::is_a_v<::Tags::Multiplies, Tag>>> {
67 : template <typename DataBoxTagList>
68 : SPECTRE_ALWAYS_INLINE decltype(auto) operator()(
69 : const db::DataBox<DataBoxTagList>& box) {
70 : return get(db::get<Tags::Dy<Tags::Dy<Tag>>>(box)).data();
71 : }
72 : };
73 :
74 : // use the product rule to provide an expression template for derivatives of
75 : // products
76 : template <typename LhsTag, typename RhsTag>
77 : struct OnDemandInputsForSwshJacobianImpl<
78 : Tags::Dy<::Tags::Multiplies<LhsTag, RhsTag>>,
79 : std::integral_constant<int,
80 : LhsTag::type::type::spin + RhsTag::type::type::spin>,
81 : std::bool_constant<not std::is_same_v<LhsTag, Tags::BondiJbar> and
82 : not std::is_same_v<LhsTag, Tags::BondiUbar> and
83 : not std::is_same_v<RhsTag, Tags::BondiJbar>>> {
84 : template <typename DataBoxTagList>
85 : SPECTRE_ALWAYS_INLINE decltype(auto) operator()(
86 : const db::DataBox<DataBoxTagList>& box) {
87 : decltype(auto) lhs = get(db::get<LhsTag>(box)).data();
88 : decltype(auto) dy_lhs = get(db::get<Tags::Dy<LhsTag>>(box)).data();
89 : decltype(auto) rhs = get(db::get<RhsTag>(box)).data();
90 : decltype(auto) dy_rhs = get(db::get<Tags::Dy<RhsTag>>(box)).data();
91 : return lhs * dy_rhs + dy_lhs * rhs;
92 : }
93 : };
94 :
95 : // use the product rule and an explicit conjugate to provide an expression
96 : // template for derivatives of products with `Tags::BondiJbar` as the right-hand
97 : // operand
98 : template <typename LhsTag>
99 : struct OnDemandInputsForSwshJacobianImpl<
100 : Tags::Dy<::Tags::Multiplies<LhsTag, Tags::BondiJbar>>,
101 : std::integral_constant<int, LhsTag::type::type::spin - 2>, std::true_type> {
102 : template <typename DataBoxTagList>
103 : SPECTRE_ALWAYS_INLINE decltype(auto) operator()(
104 : const db::DataBox<DataBoxTagList>& box) {
105 : decltype(auto) lhs = get(get<LhsTag>(box)).data();
106 : decltype(auto) dy_lhs = get(get<Tags::Dy<LhsTag>>(box)).data();
107 : decltype(auto) jbar = conj(get(get<Tags::BondiJ>(box)).data());
108 : decltype(auto) dy_jbar = conj(get(get<Tags::Dy<Tags::BondiJ>>(box)).data());
109 : return lhs * dy_jbar + dy_lhs * jbar;
110 : }
111 : };
112 :
113 : // use the product rule and an explicit conjugate to provide an expression
114 : // template for derivatives of products with `Tags::BondiJbar` as the left-hand
115 : // operand.
116 : template <typename RhsTag>
117 : struct OnDemandInputsForSwshJacobianImpl<
118 : Tags::Dy<::Tags::Multiplies<Tags::BondiJbar, RhsTag>>,
119 : std::integral_constant<int, RhsTag::type::type::spin - 2>, std::true_type> {
120 : template <typename DataBoxTagList>
121 : SPECTRE_ALWAYS_INLINE decltype(auto) operator()(
122 : const db::DataBox<DataBoxTagList>& box) {
123 : decltype(auto) rhs = get(get<RhsTag>(box)).data();
124 : decltype(auto) dy_rhs = get(get<Tags::Dy<RhsTag>>(box)).data();
125 : decltype(auto) jbar = conj(get(get<Tags::BondiJ>(box)).data());
126 : decltype(auto) dy_jbar = conj(get(get<Tags::Dy<Tags::BondiJ>>(box)).data());
127 : return dy_jbar * rhs + jbar * dy_rhs;
128 : }
129 : };
130 :
131 : // use the product rule and an explicit conjugate to provide an expression
132 : // template for derivatives of products with `Tags::BondiUbar` as the left-hand
133 : // operand.
134 : template <typename RhsTag>
135 : struct OnDemandInputsForSwshJacobianImpl<
136 : Tags::Dy<::Tags::Multiplies<Tags::BondiUbar, RhsTag>>,
137 : std::integral_constant<int, RhsTag::type::type::spin - 1>, std::true_type> {
138 : template <typename DataBoxTagList>
139 : SPECTRE_ALWAYS_INLINE decltype(auto) operator()(
140 : const db::DataBox<DataBoxTagList>& box) {
141 : decltype(auto) ubar = conj(get(get<Tags::BondiU>(box)).data());
142 : decltype(auto) dy_ubar = conj(get(get<Tags::Dy<Tags::BondiU>>(box)).data());
143 : decltype(auto) rhs = get(get<RhsTag>(box)).data();
144 : decltype(auto) dy_rhs = get(get<Tags::Dy<RhsTag>>(box)).data();
145 : return ubar * dy_rhs + dy_ubar * rhs;
146 : }
147 : };
148 :
149 : // use the product rule and an explicit conjugate to provide an expression
150 : // template for second derivatives of products with `Tags::BondiJbar` as the
151 : // right-hand operand.
152 : template <typename LhsTag>
153 : struct OnDemandInputsForSwshJacobianImpl<
154 : Tags::Dy<Tags::Dy<::Tags::Multiplies<LhsTag, Tags::BondiJbar>>>,
155 : std::integral_constant<int, LhsTag::type::type::spin - 2>, std::true_type> {
156 : template <typename DataBoxTagList>
157 : SPECTRE_ALWAYS_INLINE decltype(auto) operator()(
158 : const db::DataBox<DataBoxTagList>& box) {
159 : decltype(auto) lhs = get(get<LhsTag>(box)).data();
160 : decltype(auto) dy_lhs = get(get<Tags::Dy<LhsTag>>(box)).data();
161 : decltype(auto) dy_dy_lhs = get(get<Tags::Dy<Tags::Dy<LhsTag>>>(box)).data();
162 : decltype(auto) jbar = conj(get(get<Tags::BondiJ>(box)).data());
163 : decltype(auto) dy_jbar = conj(get(get<Tags::Dy<Tags::BondiJ>>(box)).data());
164 : decltype(auto) dy_dy_jbar =
165 : conj(get(get<Tags::Dy<Tags::Dy<Tags::BondiJ>>>(box)).data());
166 : return lhs * dy_dy_jbar + 2.0 * dy_lhs * dy_jbar + dy_dy_lhs * jbar;
167 : }
168 : };
169 :
170 : // default to extracting directly from the box for spin-weighted derivatives of
171 : // radial (y) derivatives
172 : template <typename Tag, typename DerivKind>
173 : struct OnDemandInputsForSwshJacobianImpl<
174 : Spectral::Swsh::Tags::Derivative<Tags::Dy<Tag>, DerivKind>,
175 : std::integral_constant<
176 : int, Spectral::Swsh::Tags::Derivative<Tags::Dy<Tag>, DerivKind>::spin>,
177 : std::true_type> {
178 : template <typename DataBoxTagList>
179 : SPECTRE_ALWAYS_INLINE decltype(auto) operator()(
180 : const db::DataBox<DataBoxTagList>& box) {
181 : return get(get<Spectral::Swsh::Tags::Derivative<Tags::Dy<Tag>, DerivKind>>(
182 : box))
183 : .data();
184 : }
185 : };
186 :
187 : // compute the derivative of the `jbar * (q - 2 eth_beta)` using the product
188 : // rule and the commutation rule for the partials with respect to y and the
189 : // spin_weighted derivatives
190 : template <>
191 : struct OnDemandInputsForSwshJacobianImpl<Tags::Dy<Tags::JbarQMinus2EthBeta>,
192 : std::integral_constant<int, -1>,
193 : std::true_type> {
194 : template <typename DataBoxTagList>
195 : SPECTRE_ALWAYS_INLINE decltype(auto) operator()(
196 : const db::DataBox<DataBoxTagList>& box) {
197 : decltype(auto) dy_beta = get(get<Tags::Dy<Tags::BondiBeta>>(box)).data();
198 : decltype(auto) dy_j = get(get<Tags::Dy<Tags::BondiJ>>(box)).data();
199 : decltype(auto) dy_q = get(get<Tags::Dy<Tags::BondiQ>>(box)).data();
200 : decltype(auto) eth_beta =
201 : get(get<Spectral::Swsh::Tags::Derivative<Tags::BondiBeta,
202 : Spectral::Swsh::Tags::Eth>>(
203 : box))
204 : .data();
205 : decltype(auto) eth_dy_beta =
206 : get(get<Spectral::Swsh::Tags::Derivative<Tags::Dy<Tags::BondiBeta>,
207 : Spectral::Swsh::Tags::Eth>>(
208 : box))
209 : .data();
210 : decltype(auto) eth_r_divided_by_r =
211 : get(get<Tags::EthRDividedByR>(box)).data();
212 : decltype(auto) j = get(get<Tags::BondiJ>(box)).data();
213 : decltype(auto) q = get(get<Tags::BondiQ>(box)).data();
214 : return conj(j) * dy_q + conj(dy_j) * q - 2.0 * conj(j) * eth_dy_beta -
215 : 2.0 * eth_beta * conj(dy_j) -
216 : 2.0 * conj(j) * eth_r_divided_by_r * dy_beta;
217 : }
218 : };
219 : } // namespace detail
220 :
221 : /// Provide an expression template or reference to `Tag`, intended for
222 : /// situations for which a repeated computation is more
223 : /// desirable than storing a value in the \ref DataBoxGroup (e.g. for
224 : /// conjugation and simple product rule expansion).
225 : template <typename Tag>
226 1 : using OnDemandInputsForSwshJacobian = detail::OnDemandInputsForSwshJacobianImpl<
227 : Tag, std::integral_constant<int, Tag::type::type::spin>, std::true_type>;
228 :
229 : /*!
230 : * \brief Performs a mutation to a spin-weighted spherical harmonic derivative
231 : * value from the numerical coordinate (the spin-weighted derivative at
232 : * fixed \f$y\f$) to the Bondi coordinates (the spin-weighted derivative at
233 : * fixed \f$r\f$), inplace to the requested tag.
234 : *
235 : * \details This should be performed only once for each derivative evaluation
236 : * for each tag, as a repeated inplace evaluation will compound and result in
237 : * incorrect values in the \ref DataBoxGroup. This is compatible with acting as
238 : * a mutation in `db::mutate_apply`.
239 : * \note In each specialization, there is an additional type alias
240 : * `on_demand_argument_tags` that contains tags that represent additional
241 : * quantities to be passed as arguments that need not be in the \ref
242 : * DataBoxGroup. These quantities are suggested to be evaluated by the 'on
243 : * demand' mechanism provided by `Cce::OnDemandInputsForSwshJacobian`, which
244 : * provides the additional quantities as blaze expression templates rather than
245 : * unnecessarily caching intermediate results that aren't re-used.
246 : */
247 : template <typename DerivativeTag>
248 1 : struct ApplySwshJacobianInplace;
249 :
250 : /*!
251 : * \brief Specialization for the spin-weighted derivative \f$\eth\f$.
252 : *
253 : * \details The implemented equation is:
254 : *
255 : * \f[ \eth F = \eth^\prime F - (1 - y) \frac{\eth R}{R} \partial_y F,
256 : * \f]
257 : *
258 : * where \f$\eth\f$ is the derivative at constant Bondi radius \f$r\f$ and
259 : * \f$\eth^\prime\f$ is the derivative at constant numerical radius \f$y\f$.
260 : */
261 : template <typename ArgumentTag>
262 1 : struct ApplySwshJacobianInplace<
263 : Spectral::Swsh::Tags::Derivative<ArgumentTag, Spectral::Swsh::Tags::Eth>> {
264 0 : using pre_swsh_derivative_tags = tmpl::list<>;
265 0 : using swsh_derivative_tags = tmpl::list<>;
266 0 : using integration_independent_tags =
267 : tmpl::list<Tags::OneMinusY, Tags::EthRDividedByR>;
268 :
269 0 : using return_tags = tmpl::list<
270 : Spectral::Swsh::Tags::Derivative<ArgumentTag, Spectral::Swsh::Tags::Eth>>;
271 0 : using argument_tags = tmpl::append<integration_independent_tags>;
272 0 : using on_demand_argument_tags = tmpl::list<Tags::Dy<ArgumentTag>>;
273 :
274 0 : static constexpr int spin =
275 : Spectral::Swsh::Tags::Derivative<ArgumentTag,
276 : Spectral::Swsh::Tags::Eth>::spin;
277 : template <typename DyArgumentType>
278 0 : static void apply(
279 : const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, spin>>*>
280 : eth_argument,
281 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& one_minus_y,
282 : const Scalar<SpinWeighted<ComplexDataVector, 1>>& eth_r_divided_by_r,
283 : const DyArgumentType& dy_argument) {
284 : get(*eth_argument) -= get(one_minus_y) * get(eth_r_divided_by_r) *
285 : SpinWeighted<DyArgumentType, spin - 1>{dy_argument};
286 : }
287 : };
288 :
289 : /*!
290 : * \brief Specialization for the spin-weighted derivative \f$\bar{\eth}\f$.
291 : *
292 : * \details The implemented equation is:
293 : *
294 : * \f[
295 : * \bar{\eth} F = \bar{\eth}^\prime F
296 : * - (1 - y) \frac{\bar{\eth} R}{R} \partial_y F,
297 : *\f]
298 : *
299 : * where \f$\bar{\eth}\f$ is the derivative at constant Bondi radius \f$r\f$ and
300 : * \f$\bar{\eth}^\prime\f$ is the derivative at constant numerical radius
301 : * \f$y\f$.
302 : */
303 : template <typename ArgumentTag>
304 1 : struct ApplySwshJacobianInplace<Spectral::Swsh::Tags::Derivative<
305 : ArgumentTag, Spectral::Swsh::Tags::Ethbar>> {
306 0 : using pre_swsh_derivative_tags = tmpl::list<>;
307 0 : using swsh_derivative_tags = tmpl::list<>;
308 0 : using integration_independent_tags =
309 : tmpl::list<Tags::OneMinusY, Tags::EthRDividedByR>;
310 :
311 0 : using return_tags = tmpl::list<Spectral::Swsh::Tags::Derivative<
312 : ArgumentTag, Spectral::Swsh::Tags::Ethbar>>;
313 0 : using argument_tags = tmpl::append<integration_independent_tags>;
314 0 : using on_demand_argument_tags = tmpl::list<Tags::Dy<ArgumentTag>>;
315 :
316 0 : static constexpr int spin =
317 : Spectral::Swsh::Tags::Derivative<ArgumentTag,
318 : Spectral::Swsh::Tags::Ethbar>::spin;
319 : template <typename DyArgumentType>
320 0 : static void apply(
321 : const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, spin>>*>
322 : ethbar_argument,
323 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& one_minus_y,
324 : const Scalar<SpinWeighted<ComplexDataVector, 1>>& eth_r_divided_by_r,
325 : const DyArgumentType& dy_argument) {
326 : get(*ethbar_argument) -=
327 : get(one_minus_y) * conj(get(eth_r_divided_by_r)) *
328 : SpinWeighted<DyArgumentType, spin + 1>{dy_argument};
329 : }
330 : };
331 :
332 : /*!
333 : * \brief Specialization for the spin-weighted derivative \f$\eth \bar{\eth}\f$.
334 : *
335 : * \details The implemented equation is:
336 : *
337 : * \f[
338 : * \eth \bar{\eth} F = \eth^\prime \bar{\eth}^\prime F
339 : * - \frac{\eth R \bar{\eth} R}{R^2} (1 - y)^2 \partial_y^2 F
340 : * - (1 - y)\left(\frac{\eth R}{R} \bar{\eth} \partial_y F
341 : * + \frac{\bar{\eth} R}{R} \eth \partial_y F
342 : * + \frac{\eth \bar\eth R}{R} \partial_y F\right),
343 : * \f]
344 : *
345 : * where \f$\eth \bar{\eth}\f$ is the derivative at constant Bondi radius
346 : * \f$r\f$ and \f$\eth^\prime \bar{\eth}^\prime\f$ is the derivative at constant
347 : * numerical radius \f$y\f$.
348 : */
349 : template <typename ArgumentTag>
350 1 : struct ApplySwshJacobianInplace<Spectral::Swsh::Tags::Derivative<
351 : ArgumentTag, Spectral::Swsh::Tags::EthEthbar>> {
352 0 : using pre_swsh_derivative_tags = tmpl::list<>;
353 0 : using swsh_derivative_tags = tmpl::list<>;
354 0 : using integration_independent_tags =
355 : tmpl::list<Tags::OneMinusY, Tags::EthRDividedByR,
356 : Tags::EthEthbarRDividedByR>;
357 :
358 0 : using return_tags = tmpl::list<Spectral::Swsh::Tags::Derivative<
359 : ArgumentTag, Spectral::Swsh::Tags::EthEthbar>>;
360 0 : using argument_tags = tmpl::append<integration_independent_tags>;
361 0 : using on_demand_argument_tags =
362 : tmpl::list<Tags::Dy<ArgumentTag>, Tags::Dy<Tags::Dy<ArgumentTag>>,
363 : Spectral::Swsh::Tags::Derivative<Tags::Dy<ArgumentTag>,
364 : Spectral::Swsh::Tags::Eth>,
365 : Spectral::Swsh::Tags::Derivative<
366 : Tags::Dy<ArgumentTag>, Spectral::Swsh::Tags::Ethbar>>;
367 :
368 0 : static constexpr int spin =
369 : Spectral::Swsh::Tags::Derivative<ArgumentTag,
370 : Spectral::Swsh::Tags::EthEthbar>::spin;
371 : template <typename DyArgumentType, typename DyDyArgumentType,
372 : typename EthDyArgumentType, typename EthbarDyArgumentType>
373 0 : static void apply(
374 : const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, spin>>*>
375 : eth_ethbar_argument,
376 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& one_minus_y,
377 : const Scalar<SpinWeighted<ComplexDataVector, 1>>& eth_r_divided_by_r,
378 : const Scalar<SpinWeighted<ComplexDataVector, 0>>&
379 : eth_ethbar_r_divided_by_r,
380 : const DyArgumentType& dy_argument, const DyDyArgumentType& dy_dy_argument,
381 : const EthDyArgumentType& eth_dy_argument,
382 : const EthbarDyArgumentType ethbar_dy_argument) {
383 : get(*eth_ethbar_argument) -=
384 : get(eth_r_divided_by_r) * conj(get(eth_r_divided_by_r)) *
385 : (square(get(one_minus_y)) *
386 : SpinWeighted<DyDyArgumentType, spin>{dy_dy_argument}) +
387 : get(one_minus_y) *
388 : (get(eth_r_divided_by_r) *
389 : SpinWeighted<EthbarDyArgumentType, spin - 1>{
390 : ethbar_dy_argument} +
391 : conj(get(eth_r_divided_by_r)) *
392 : SpinWeighted<EthDyArgumentType, spin + 1>{eth_dy_argument} +
393 : get(eth_ethbar_r_divided_by_r) *
394 : SpinWeighted<DyArgumentType, spin>{dy_argument});
395 : }
396 : };
397 :
398 : /*!
399 : * \brief Specialization for the spin-weighted derivative \f$\bar{\eth} \eth\f$.
400 : *
401 : * \details The implemented equation is:
402 : *
403 : * \f[
404 : * \bar{\eth} \eth F = \bar{\eth}^\prime \eth^\prime F
405 : * - \frac{\eth R \bar{\eth} R}{R^2} (1 - y)^2 \partial_y^2 F
406 : * - (1 - y)\left(\frac{\eth R}{R} \bar{\eth} \partial_y F
407 : * + \frac{\bar{\eth} R}{R} \eth \partial_y F
408 : * + \frac{\eth \bar\eth R}{R} \partial_y F\right),
409 : * \f]
410 : *
411 : * where \f$\bar{\eth} \eth\f$ is the derivative at constant Bondi radius
412 : * \f$r\f$ and \f$\bar{\eth}^\prime \eth^\prime\f$ is the derivative at constant
413 : * numerical radius \f$y\f$.
414 : */
415 : template <typename ArgumentTag>
416 1 : struct ApplySwshJacobianInplace<Spectral::Swsh::Tags::Derivative<
417 : ArgumentTag, Spectral::Swsh::Tags::EthbarEth>> {
418 0 : using pre_swsh_derivative_tags = tmpl::list<>;
419 0 : using swsh_derivative_tags = tmpl::list<>;
420 0 : using integration_independent_tags =
421 : tmpl::list<Tags::OneMinusY, Tags::EthRDividedByR,
422 : Tags::EthEthbarRDividedByR>;
423 :
424 0 : using return_tags = tmpl::list<Spectral::Swsh::Tags::Derivative<
425 : ArgumentTag, Spectral::Swsh::Tags::EthbarEth>>;
426 0 : using argument_tags = tmpl::append<integration_independent_tags>;
427 0 : using on_demand_argument_tags =
428 : tmpl::list<Tags::Dy<ArgumentTag>, Tags::Dy<Tags::Dy<ArgumentTag>>,
429 : Spectral::Swsh::Tags::Derivative<Tags::Dy<ArgumentTag>,
430 : Spectral::Swsh::Tags::Eth>,
431 : Spectral::Swsh::Tags::Derivative<
432 : Tags::Dy<ArgumentTag>, Spectral::Swsh::Tags::Ethbar>>;
433 :
434 0 : static constexpr int spin =
435 : Spectral::Swsh::Tags::Derivative<ArgumentTag,
436 : Spectral::Swsh::Tags::EthbarEth>::spin;
437 : template <typename DyArgumentType, typename DyDyArgumentType,
438 : typename EthDyArgumentType, typename EthbarDyArgumentType>
439 0 : static void apply(
440 : const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, spin>>*>
441 : ethbar_eth_argument,
442 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& one_minus_y,
443 : const Scalar<SpinWeighted<ComplexDataVector, 1>>& eth_r_divided_by_r,
444 : const Scalar<SpinWeighted<ComplexDataVector, 0>>&
445 : eth_ethbar_r_divided_by_r,
446 : const DyArgumentType& dy_argument, const DyDyArgumentType& dy_dy_argument,
447 : const EthDyArgumentType& eth_dy_argument,
448 : const EthbarDyArgumentType ethbar_dy_argument) {
449 : get(*ethbar_eth_argument) -=
450 : get(eth_r_divided_by_r) * conj(get(eth_r_divided_by_r)) *
451 : (square(get(one_minus_y)) *
452 : SpinWeighted<DyDyArgumentType, spin>{dy_dy_argument}) +
453 : get(one_minus_y) *
454 : (get(eth_r_divided_by_r) *
455 : SpinWeighted<EthbarDyArgumentType, spin - 1>{
456 : ethbar_dy_argument} +
457 : conj(get(eth_r_divided_by_r)) *
458 : SpinWeighted<EthDyArgumentType, spin + 1>{eth_dy_argument} +
459 : get(eth_ethbar_r_divided_by_r) *
460 : SpinWeighted<DyArgumentType, spin>{dy_argument});
461 : }
462 : };
463 :
464 : /*!
465 : * \brief Specialization for the spin-weighted derivative \f$\eth \eth\f$.
466 : *
467 : * \details The implemented equation is:
468 : *
469 : * \f[
470 : * \eth \eth F = \eth^\prime \eth^\prime F
471 : * - (1 - y)^2 \frac{(\eth R)^2}{R^2} \partial_y^2 F
472 : * - (1 - y) \left( 2 \frac{\eth R}{R} \eth \partial_y F
473 : * + \frac{\eth \eth R}{R} \partial_y F\right),
474 : * \f]
475 : *
476 : * where \f$\eth \eth\f$ is the derivative at constant Bondi radius \f$r\f$ and
477 : * \f$\eth^\prime \eth^\prime\f$ is the derivative at constant numerical radius
478 : * \f$y\f$.
479 : */
480 : template <typename ArgumentTag>
481 1 : struct ApplySwshJacobianInplace<Spectral::Swsh::Tags::Derivative<
482 : ArgumentTag, Spectral::Swsh::Tags::EthEth>> {
483 0 : using pre_swsh_derivative_tags = tmpl::list<>;
484 0 : using swsh_derivative_tags = tmpl::list<>;
485 0 : using integration_independent_tags =
486 : tmpl::list<Tags::OneMinusY, Tags::EthRDividedByR,
487 : Tags::EthEthRDividedByR>;
488 :
489 0 : using return_tags = tmpl::list<Spectral::Swsh::Tags::Derivative<
490 : ArgumentTag, Spectral::Swsh::Tags::EthEth>>;
491 0 : using argument_tags = tmpl::append<integration_independent_tags>;
492 0 : using on_demand_argument_tags =
493 : tmpl::list<Tags::Dy<ArgumentTag>, Tags::Dy<Tags::Dy<ArgumentTag>>,
494 : Spectral::Swsh::Tags::Derivative<Tags::Dy<ArgumentTag>,
495 : Spectral::Swsh::Tags::Eth>>;
496 :
497 0 : static constexpr int spin =
498 : Spectral::Swsh::Tags::Derivative<ArgumentTag,
499 : Spectral::Swsh::Tags::EthEth>::spin;
500 : template <typename DyArgumentType, typename DyDyArgumentType,
501 : typename EthDyArgumentType>
502 0 : static void apply(
503 : const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, spin>>*>
504 : eth_eth_argument,
505 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& one_minus_y,
506 : const Scalar<SpinWeighted<ComplexDataVector, 1>>& eth_r_divided_by_r,
507 : const Scalar<SpinWeighted<ComplexDataVector, 2>>& eth_eth_r_divided_by_r,
508 : const DyArgumentType& dy_argument, const DyDyArgumentType& dy_dy_argument,
509 : const EthDyArgumentType& eth_dy_argument) {
510 : get(*eth_eth_argument) -=
511 : square(get(eth_r_divided_by_r)) *
512 : (square(get(one_minus_y)) *
513 : SpinWeighted<DyDyArgumentType, spin - 2>{dy_dy_argument}) +
514 : get(one_minus_y) *
515 : (2.0 * get(eth_r_divided_by_r) *
516 : SpinWeighted<EthDyArgumentType, spin - 1>{eth_dy_argument} +
517 : get(eth_eth_r_divided_by_r) *
518 : SpinWeighted<DyArgumentType, spin - 2>{dy_argument});
519 : }
520 : };
521 :
522 : /*!
523 : * \brief Specialization for the spin-weighted derivative \f$\bar{\eth}
524 : * \bar{\eth}\f$.
525 : *
526 : * \details The implemented equation is:
527 : *
528 : * \f[
529 : * \bar{\eth} \bar{\eth} F = \bar{\eth}^\prime \bar{\eth}^\prime F
530 : * - (1 - y)^2 \frac{(\bar{\eth} R)^2}{R^2} \partial_y^2 F
531 : * - (1 - y) \left( 2 \frac{\bar{\eth} R}{R} \bar{\eth} \partial_y F
532 : * + \frac{\bar{\eth} \bar{\eth} R}{R} \partial_y F\right),
533 : * \f]
534 : *
535 : * where \f$\bar{\eth} \bar{\eth}\f$ is the derivative at constant Bondi radius
536 : * \f$r\f$ and \f$\bar{\eth}^\prime \bar{\eth}^\prime\f$ is the derivative at
537 : * constant numerical radius \f$y\f$.
538 : */
539 : template <typename ArgumentTag>
540 1 : struct ApplySwshJacobianInplace<Spectral::Swsh::Tags::Derivative<
541 : ArgumentTag, Spectral::Swsh::Tags::EthbarEthbar>> {
542 0 : using pre_swsh_derivative_tags = tmpl::list<>;
543 0 : using swsh_derivative_tags = tmpl::list<>;
544 0 : using integration_independent_tags =
545 : tmpl::list<Tags::OneMinusY, Tags::EthRDividedByR,
546 : Tags::EthEthRDividedByR>;
547 :
548 0 : using return_tags = tmpl::list<Spectral::Swsh::Tags::Derivative<
549 : ArgumentTag, Spectral::Swsh::Tags::EthbarEthbar>>;
550 0 : using argument_tags = tmpl::append<integration_independent_tags>;
551 0 : using on_demand_argument_tags =
552 : tmpl::list<Tags::Dy<ArgumentTag>, Tags::Dy<Tags::Dy<ArgumentTag>>,
553 : Spectral::Swsh::Tags::Derivative<
554 : Tags::Dy<ArgumentTag>, Spectral::Swsh::Tags::Ethbar>>;
555 :
556 0 : static constexpr int spin = Spectral::Swsh::Tags::Derivative<
557 : ArgumentTag, Spectral::Swsh::Tags::EthbarEthbar>::spin;
558 : template <typename DyArgumentType, typename DyDyArgumentType,
559 : typename EthbarDyArgumentType>
560 0 : static void apply(
561 : const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, spin>>*>
562 : ethbar_ethbar_argument,
563 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& one_minus_y,
564 : const Scalar<SpinWeighted<ComplexDataVector, 1>>& eth_r_divided_by_r,
565 : const Scalar<SpinWeighted<ComplexDataVector, 2>>& eth_eth_r_divided_by_r,
566 : const DyArgumentType& dy_argument, const DyDyArgumentType& dy_dy_argument,
567 : const EthbarDyArgumentType& ethbar_dy_argument) {
568 : get(*ethbar_ethbar_argument) -=
569 : square(conj(get(eth_r_divided_by_r))) *
570 : (square(get(one_minus_y)) *
571 : SpinWeighted<DyDyArgumentType, spin + 2>{dy_dy_argument}) +
572 : get(one_minus_y) *
573 : (2.0 * conj(get(eth_r_divided_by_r)) *
574 : SpinWeighted<EthbarDyArgumentType, spin + 1>{
575 : ethbar_dy_argument} +
576 : conj(get(eth_eth_r_divided_by_r)) *
577 : SpinWeighted<DyArgumentType, spin + 2>{dy_argument});
578 : }
579 : };
580 :
581 : namespace detail {
582 : // A helper to forward to the `ApplySwshJacobianInplace` mutators that takes
583 : // advantage of the `OnDemandInputsForSwshJacobian`, which computes blaze
584 : // template expressions for those quantities for which it is anticipated to be
585 : // sufficiently cheap to repeatedly compute that it is worth saving the cost of
586 : // additional storage (e.g. conjugates and derivatives for which the product
587 : // rule applies)
588 : template <typename DerivativeTag, typename DataBoxTagList,
589 : typename... OnDemandTags>
590 : void apply_swsh_jacobian_helper(
591 : const gsl::not_null<db::DataBox<DataBoxTagList>*> box,
592 : tmpl::list<OnDemandTags...> /*meta*/) {
593 : db::mutate_apply<ApplySwshJacobianInplace<DerivativeTag>>(
594 : box, OnDemandInputsForSwshJacobian<OnDemandTags>{}(*box)...);
595 : }
596 : } // namespace detail
597 :
598 : /*!
599 : * \brief This routine evaluates the set of inputs to the CCE integrand for
600 : * `BondiValueTag` which are spin-weighted angular derivatives.
601 :
602 : * \details This function is called on the \ref DataBoxGroup holding the
603 : * relevant CCE data during each hypersurface integration step, after evaluating
604 : * `mutate_all_pre_swsh_derivatives_for_tag()` with template argument
605 : * `BondiValueTag` and before evaluating `ComputeBondiIntegrand<BondiValueTag>`.
606 : * Provided a \ref DataBoxGroup with the appropriate tags (including
607 : * `Cce::all_pre_swsh_derivative_tags`, `Cce::all_swsh_derivative_tags`,
608 : * `Cce::all_transform_buffer_tags`, `Cce::pre_computation_tags`, and
609 : * `Cce::Tags::LMax`), this function will apply all of the necessary
610 : * mutations to update
611 : * `Cce::single_swsh_derivative_tags_to_compute_for<BondiValueTag>` and
612 : * `Cce::second_swsh_derivative_tags_to_compute_for<BondiValueTag>` to their
613 : * correct values for the current values of the remaining (input) tags.
614 : */
615 : template <typename BondiValueTag, typename DataBoxTagList>
616 1 : void mutate_all_swsh_derivatives_for_tag(
617 : const gsl::not_null<db::DataBox<DataBoxTagList>*> box) {
618 : // The collection of spin-weighted derivatives cannot be applied as individual
619 : // compute items, because it is better to aggregate similar spins and dispatch
620 : // to libsharp in groups. So, we supply a bulk mutate operation which takes in
621 : // multiple Variables from the presumed DataBox, and alters their values as
622 : // necessary.
623 : db::mutate_apply<Spectral::Swsh::AngularDerivatives<
624 : single_swsh_derivative_tags_to_compute_for_t<BondiValueTag>>>(box);
625 : tmpl::for_each<single_swsh_derivative_tags_to_compute_for_t<BondiValueTag>>(
626 : [&box](auto derivative_tag_v) {
627 : using derivative_tag = typename decltype(derivative_tag_v)::type;
628 : detail::apply_swsh_jacobian_helper<derivative_tag>(
629 : box, typename ApplySwshJacobianInplace<
630 : derivative_tag>::on_demand_argument_tags{});
631 : });
632 :
633 : db::mutate_apply<Spectral::Swsh::AngularDerivatives<
634 : second_swsh_derivative_tags_to_compute_for_t<BondiValueTag>>>(box);
635 : tmpl::for_each<second_swsh_derivative_tags_to_compute_for_t<BondiValueTag>>(
636 : [&box](auto derivative_tag_v) {
637 : using derivative_tag = typename decltype(derivative_tag_v)::type;
638 : detail::apply_swsh_jacobian_helper<derivative_tag>(
639 : box, typename ApplySwshJacobianInplace<
640 : derivative_tag>::on_demand_argument_tags{});
641 : });
642 : }
643 :
644 : namespace Tags {
645 :
646 : /*!
647 : * \brief Compute tag for a manually-handled Swsh derivative in the volume.
648 : *
649 : * \details Numerical swsh derivatives are $\breve\eth$ (at fixed $y$), but we
650 : * prefer to write our calculations here in terms of $\eth$ (at fixed $r$). This
651 : * involves a correction term that depends on the $\partial_{\breve{y}}$
652 : * derivative of the field being differentiated. Most of these are handled in
653 : * place by `ApplySwshJacobianInplace`. However, that infrastructure is specific
654 : * to the Cce evolution system itself, so it won't work for other uses
655 : * (e.g. observers that are invoked periodically). This compute tag instead
656 : * manually computes a certain swsh derivative (at fixed $r$). Right now this
657 : * will only work when `DerivativeKind` is either `Spectral::Swsh::Tags::Eth` or
658 : * `Spectral::Swsh::Tags::Ethbar`. If you want to extend it to second
659 : * derivatives, there will have to be more template specializations.
660 : */
661 : template <typename ArgumentTag, typename DerivativeKind>
662 1 : struct SwshDerivativeCompute
663 : : Spectral::Swsh::Tags::Derivative<ArgumentTag, DerivativeKind>,
664 : db::ComputeTag {
665 0 : using base = Spectral::Swsh::Tags::Derivative<ArgumentTag, DerivativeKind>;
666 0 : using return_type = typename base::type;
667 :
668 0 : using ApplySwshJac = ApplySwshJacobianInplace<base>;
669 :
670 0 : using argument_tags = tmpl::append<
671 : tmpl::list<Tags::LMax, Tags::NumberOfRadialPoints, ArgumentTag>,
672 : typename ApplySwshJac::argument_tags,
673 : typename ApplySwshJac::on_demand_argument_tags>;
674 :
675 0 : static constexpr auto function(
676 : const gsl::not_null<return_type*> eth_field, size_t l_max,
677 : size_t number_of_radial_points, const typename ArgumentTag::type& field,
678 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& one_minus_y,
679 : const Scalar<SpinWeighted<ComplexDataVector, 1>>& eth_r_divided_by_r,
680 : const typename Tags::Dy<ArgumentTag>::type& dy_field) {
681 : get(*eth_field) = Spectral::Swsh::angular_derivative<DerivativeKind>(
682 : l_max, number_of_radial_points, get(field));
683 :
684 : ApplySwshJac::apply(eth_field, one_minus_y, eth_r_divided_by_r,
685 : get(dy_field).data());
686 : }
687 : };
688 :
689 : } // namespace Tags
690 :
691 : } // namespace Cce
|
