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/ComplexDataVector.hpp"
9 : #include "DataStructures/DataBox/DataBox.hpp"
10 : #include "DataStructures/SpinWeighted.hpp"
11 : #include "Evolution/Systems/Cce/IntegrandInputSteps.hpp"
12 : #include "Evolution/Systems/Cce/OptionTags.hpp"
13 : #include "Evolution/Systems/Cce/Tags.hpp"
14 : #include "NumericalAlgorithms/Spectral/Basis.hpp"
15 : #include "NumericalAlgorithms/Spectral/Quadrature.hpp"
16 : #include "NumericalAlgorithms/Spectral/Spectral.hpp"
17 : #include "NumericalAlgorithms/Spectral/SwshCollocation.hpp"
18 : #include "NumericalAlgorithms/Spectral/SwshTags.hpp"
19 : #include "Utilities/Gsl.hpp"
20 : #include "Utilities/TMPL.hpp"
21 : #include "Utilities/VectorAlgebra.hpp"
22 :
23 : namespace Cce {
24 :
25 : namespace detail {
26 : // A convenience function for computing the spin-weighted derivatives of \f$R\f$
27 : // divided by \f$R\f$, which appears often in Jacobians to transform between
28 : // Bondi coordinates and the numerical coordinates used in CCE.
29 : template <typename DerivKind>
30 : void angular_derivative_of_r_divided_by_r_impl(
31 : gsl::not_null<
32 : SpinWeighted<ComplexDataVector,
33 : Spectral::Swsh::Tags::derivative_spin_weight<DerivKind>>*>
34 : d_r_divided_by_r,
35 : const SpinWeighted<ComplexDataVector, 0>& boundary_r, size_t l_max,
36 : size_t number_of_radial_points);
37 :
38 : } // namespace detail
39 :
40 : /*!
41 : * \brief A set of procedures for computing the set of inputs to the CCE
42 : * integrand computations that can be computed before any of the intermediate
43 : * integrands are evaluated.
44 : *
45 : * \details The template specializations of this template are
46 : * compatible with acting as a the mutator in a \ref DataBoxGroup
47 : * `db::mutate_apply` operation. For flexibility in defining the \ref
48 : * DataBoxGroup structure, the tags for `Tensor`s used in these functions are
49 : * also organized into type lists:
50 : * - type alias `integration_independent_tags`: with a subset of
51 : * `Cce::pre_computation_tags`, used for both input and output.
52 : * - type alias `boundary_values`: with a subset of
53 : * `Cce::pre_computation_boundary_tags`, used only for input.
54 : * - type alias `pre_swsh_derivatives` containing hypersurface quantities. For
55 : * this struct, it will only ever contain `Cce::Tags::BondiJ`, and is used as
56 : * input.
57 : *
58 : * The `BoundaryPrefix` tag allows easy switching between the
59 : * regularity-preserving version and standard CCE
60 : *
61 : */
62 : template <template <typename> class BoundaryPrefix, typename Tag>
63 1 : struct PrecomputeCceDependencies;
64 :
65 : /// Computes \f$1 - y\f$ for the CCE system.
66 : template <template <typename> class BoundaryPrefix>
67 1 : struct PrecomputeCceDependencies<BoundaryPrefix, Tags::OneMinusY> {
68 0 : using boundary_tags = tmpl::list<>;
69 0 : using pre_swsh_derivative_tags = tmpl::list<>;
70 0 : using integration_independent_tags = tmpl::list<>;
71 :
72 0 : using return_tags = tmpl::list<Tags::OneMinusY>;
73 0 : using argument_tags = tmpl::list<Tags::LMax, Tags::NumberOfRadialPoints>;
74 :
75 0 : static void apply(
76 : const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 0>>*>
77 : one_minus_y,
78 : const size_t l_max, const size_t number_of_radial_points) {
79 : const size_t number_of_angular_points =
80 : Spectral::Swsh::number_of_swsh_collocation_points(l_max);
81 : const DataVector one_minus_y_collocation =
82 : 1.0 - Spectral::collocation_points<Spectral::Basis::Legendre,
83 : Spectral::Quadrature::GaussLobatto>(
84 : number_of_radial_points);
85 : // iterate through the angular 'chunks' and set them to their 1-y value
86 : for (size_t i = 0; i < number_of_radial_points; ++i) {
87 : ComplexDataVector angular_view{
88 : get(*one_minus_y).data().data() + number_of_angular_points * i,
89 : number_of_angular_points};
90 : angular_view = one_minus_y_collocation[i];
91 : }
92 : }
93 : };
94 :
95 : /// Computes a volume version of Bondi radius of the worldtube \f$R\f$ from its
96 : /// boundary value (by repeating it over the radial dimension)
97 : template <template <typename> class BoundaryPrefix>
98 1 : struct PrecomputeCceDependencies<BoundaryPrefix, Tags::BondiR> {
99 0 : using boundary_tags = tmpl::list<BoundaryPrefix<Tags::BondiR>>;
100 0 : using pre_swsh_derivative_tags = tmpl::list<>;
101 0 : using integration_independent_tags = tmpl::list<>;
102 :
103 0 : using return_tags = tmpl::list<Tags::BondiR>;
104 0 : using argument_tags =
105 : tmpl::append<boundary_tags, tmpl::list<Tags::NumberOfRadialPoints>>;
106 :
107 0 : static void apply(
108 : const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 0>>*> r,
109 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& boundary_r,
110 : const size_t number_of_radial_points) {
111 : fill_with_n_copies(make_not_null(&get(*r).data()), get(boundary_r).data(),
112 : number_of_radial_points);
113 : }
114 : };
115 :
116 : /// Computes \f$\partial_u R / R\f$ from its boundary value (by repeating it
117 : /// over the radial dimension).
118 : template <template <typename> class BoundaryPrefix>
119 1 : struct PrecomputeCceDependencies<BoundaryPrefix, Tags::DuRDividedByR> {
120 0 : using boundary_tags = tmpl::list<BoundaryPrefix<Tags::DuRDividedByR>>;
121 0 : using pre_swsh_derivative_tags = tmpl::list<>;
122 0 : using integration_independent_tags = tmpl::list<>;
123 :
124 0 : using return_tags = tmpl::list<Tags::DuRDividedByR>;
125 0 : using argument_tags =
126 : tmpl::append<boundary_tags, tmpl::list<Tags::NumberOfRadialPoints>>;
127 :
128 0 : static void apply(
129 : const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 0>>*>
130 : du_r_divided_by_r,
131 : const Scalar<SpinWeighted<ComplexDataVector, 0>>&
132 : boundary_du_r_divided_by_r,
133 : const size_t number_of_radial_points) {
134 : fill_with_n_copies(make_not_null(&get(*du_r_divided_by_r).data()),
135 : get(boundary_du_r_divided_by_r).data(),
136 : number_of_radial_points);
137 : }
138 : };
139 :
140 : /// Computes \f$\eth R / R\f$ by differentiating and repeating the boundary
141 : /// value of \f$R\f$.
142 : template <template <typename> class BoundaryPrefix>
143 1 : struct PrecomputeCceDependencies<BoundaryPrefix, Tags::EthRDividedByR> {
144 0 : using boundary_tags = tmpl::list<BoundaryPrefix<Tags::BondiR>>;
145 0 : using pre_swsh_derivative_tags = tmpl::list<>;
146 0 : using integration_independent_tags = tmpl::list<>;
147 :
148 0 : using return_tags = tmpl::list<Tags::EthRDividedByR>;
149 0 : using argument_tags =
150 : tmpl::append<boundary_tags,
151 : tmpl::list<Tags::LMax, Tags::NumberOfRadialPoints>>;
152 :
153 0 : static void apply(
154 : const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 1>>*>
155 : eth_r_divided_by_r,
156 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& boundary_r,
157 : const size_t l_max, const size_t number_of_radial_points) {
158 : detail::angular_derivative_of_r_divided_by_r_impl<
159 : Spectral::Swsh::Tags::Eth>(make_not_null(&get(*eth_r_divided_by_r)),
160 : get(boundary_r), l_max,
161 : number_of_radial_points);
162 : }
163 : };
164 :
165 : /// Computes \f$\eth \eth R / R\f$ by differentiating and repeating the boundary
166 : /// value of \f$R\f$.
167 : template <template <typename> class BoundaryPrefix>
168 1 : struct PrecomputeCceDependencies<BoundaryPrefix, Tags::EthEthRDividedByR> {
169 0 : using boundary_tags = tmpl::list<BoundaryPrefix<Tags::BondiR>>;
170 0 : using pre_swsh_derivative_tags = tmpl::list<>;
171 0 : using integration_independent_tags = tmpl::list<>;
172 :
173 0 : using return_tags = tmpl::list<Tags::EthEthRDividedByR>;
174 0 : using argument_tags =
175 : tmpl::append<boundary_tags,
176 : tmpl::list<Tags::LMax, Tags::NumberOfRadialPoints>>;
177 :
178 0 : static void apply(
179 : const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 2>>*>
180 : eth_eth_r_divided_by_r,
181 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& boundary_r,
182 : const size_t l_max, const size_t number_of_radial_points) {
183 : detail::angular_derivative_of_r_divided_by_r_impl<
184 : Spectral::Swsh::Tags::EthEth>(
185 : make_not_null(&get(*eth_eth_r_divided_by_r)), get(boundary_r), l_max,
186 : number_of_radial_points);
187 : }
188 : };
189 :
190 : /// Computes \f$\eth \bar{\eth} R / R\f$ by differentiating and repeating the
191 : /// boundary value of \f$R\f$.
192 : template <template <typename> class BoundaryPrefix>
193 1 : struct PrecomputeCceDependencies<BoundaryPrefix, Tags::EthEthbarRDividedByR> {
194 0 : using boundary_tags = tmpl::list<BoundaryPrefix<Tags::BondiR>>;
195 0 : using pre_swsh_derivative_tags = tmpl::list<>;
196 0 : using integration_independent_tags = tmpl::list<>;
197 :
198 0 : using return_tags = tmpl::list<Tags::EthEthbarRDividedByR>;
199 0 : using argument_tags =
200 : tmpl::append<boundary_tags,
201 : tmpl::list<Tags::LMax, Tags::NumberOfRadialPoints>>;
202 :
203 0 : static void apply(
204 : const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 0>>*>
205 : eth_ethbar_r_divided_by_r,
206 : const Scalar<SpinWeighted<ComplexDataVector, 0>>& boundary_r,
207 : const size_t l_max, const size_t number_of_radial_points) {
208 : detail::angular_derivative_of_r_divided_by_r_impl<
209 : Spectral::Swsh::Tags::EthEthbar>(
210 : make_not_null(&get(*eth_ethbar_r_divided_by_r)), get(boundary_r), l_max,
211 : number_of_radial_points);
212 : }
213 : };
214 :
215 : /// Computes \f$K = \sqrt{1 + J \bar{J}}\f$.
216 : template <template <typename> class BoundaryPrefix>
217 1 : struct PrecomputeCceDependencies<BoundaryPrefix, Tags::BondiK> {
218 0 : using boundary_tags = tmpl::list<>;
219 0 : using pre_swsh_derivative_tags = tmpl::list<Tags::BondiJ>;
220 0 : using integration_independent_tags = tmpl::list<>;
221 :
222 0 : using return_tags = tmpl::list<Tags::BondiK>;
223 0 : using argument_tags = tmpl::push_front<pre_swsh_derivative_tags, Tags::LMax>;
224 :
225 0 : static void apply(
226 : const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 0>>*> k,
227 : const size_t /*l_max*/,
228 : const Scalar<SpinWeighted<ComplexDataVector, 2>>& j) {
229 : get(*k).data() = sqrt(1.0 + get(j).data() * conj(get(j)).data());
230 : }
231 : };
232 :
233 : /*!
234 : * \brief Convenience routine for computing all of the CCE inputs to integrand
235 : * computation that do not depend on intermediate integrand results. It should
236 : * be executed before moving through the hierarchy of integrands.
237 : *
238 : * \details Provided a \ref DataBoxGroup with the appropriate tags (including
239 : * `Cce::pre_computation_boundary_tags`, `Cce::pre_computation_tags`,
240 : * `Cce::Tags::BondiJ` and `Tags::LMax`), this function will
241 : * apply all of the necessary mutations to update the
242 : * `Cce::pre_computation_tags` to their correct values for the current values
243 : * for the remaining (input) tags.
244 : *
245 : * The `BoundaryPrefix` template template parameter is to be passed a prefix
246 : * tag associated with the boundary value prefix used in the computation (e.g.
247 : * `Cce::Tags::BoundaryValue`), and allows easy switching between the
248 : * regularity-preserving version and standard CCE.
249 : */
250 : template <template <typename> class BoundaryPrefix, typename DataBoxType>
251 1 : void mutate_all_precompute_cce_dependencies(
252 : const gsl::not_null<DataBoxType*> box) {
253 : tmpl::for_each<pre_computation_tags>([&box](auto x) {
254 : using integration_independent_tag = typename decltype(x)::type;
255 : using mutation =
256 : PrecomputeCceDependencies<BoundaryPrefix, integration_independent_tag>;
257 : db::mutate_apply<mutation>(box);
258 : });
259 : }
260 : } // namespace Cce
|
