Line data Source code
1 1 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : /// \file
5 : /// Defines functions computing partial derivatives.
6 :
7 : #pragma once
8 :
9 : #include <array>
10 : #include <cstddef>
11 : #include <string>
12 :
13 : #include "DataStructures/DataBox/PrefixHelpers.hpp"
14 : #include "DataStructures/DataBox/Prefixes.hpp"
15 : #include "DataStructures/DataBox/Tag.hpp"
16 : #include "DataStructures/DataBox/TagName.hpp"
17 : #include "DataStructures/Variables.hpp"
18 : #include "Utilities/Requires.hpp"
19 : #include "Utilities/TMPL.hpp"
20 : #include "Utilities/TypeTraits/IsA.hpp"
21 :
22 : /// \cond
23 : class DataVector;
24 : class Matrix;
25 : template <size_t Dim>
26 : class Mesh;
27 :
28 : namespace domain {
29 : namespace Tags {
30 : template <size_t Dim>
31 : struct Mesh;
32 : } // namespace Tags
33 : } // namespace domain
34 : namespace Tags {
35 : template <class TagList>
36 : struct Variables;
37 : } // namespace Tags
38 : /// \endcond
39 :
40 : namespace partial_derivatives_detail {
41 : void apply_matrix_in_first_dim(double* result, const double* input,
42 : const Matrix& matrix, size_t size,
43 : bool add_to_result = false);
44 : void apply_matrix_in_first_dim(std::complex<double>* result,
45 : const std::complex<double>* input,
46 : const Matrix& matrix, size_t size,
47 : bool add_to_result = false);
48 : } // namespace partial_derivatives_detail
49 :
50 : /// @{
51 : /// \ingroup NumericalAlgorithmsGroup
52 : /// \brief Compute the partial derivatives of each variable with respect to
53 : /// the element logical coordinate.
54 : ///
55 : /// \requires `DerivativeTags` to be the head of `VariableTags`
56 : ///
57 : /// Returns a `Variables` with a spatial tensor index appended to the front
58 : /// of each tensor within `u` and each `Tag` wrapped with a `Tags::deriv`.
59 : ///
60 : /// \tparam DerivativeTags the subset of `VariableTags` for which derivatives
61 : /// are computed.
62 : template <typename DerivativeTags, typename VariableTags, size_t Dim>
63 1 : void logical_partial_derivatives(
64 : gsl::not_null<std::array<Variables<DerivativeTags>, Dim>*>
65 : logical_partial_derivatives_of_u,
66 : const Variables<VariableTags>& u, const Mesh<Dim>& mesh);
67 :
68 : template <typename DerivativeTags, typename VariableTags, size_t Dim>
69 1 : auto logical_partial_derivatives(const Variables<VariableTags>& u,
70 : const Mesh<Dim>& mesh)
71 : -> std::array<Variables<DerivativeTags>, Dim>;
72 : /// @}
73 :
74 : /// @{
75 : /*!
76 : * \ingroup NumericalAlgorithmsGroup
77 : * \brief Computes the logical partial derivative of a tensor, prepending the
78 : * spatial derivative index, e.g. for \f$\partial_i T_{a}{}^{b}\f$ the C++ call
79 : * is `get(i, a, b)`.
80 : *
81 : * There is an overload that accepts a `buffer` of size
82 : * `mesh.number_of_grid_points()` or larger. When passed this function performs
83 : * no memory allocations, which helps improve performance.
84 : *
85 : * If you have a `Variables` with several tensors you need to differentiate you
86 : * should use the `logical_partial_derivatives` function that operates on
87 : * `Variables` since that'll be more efficient.
88 : */
89 : template <typename DataType, typename SymmList, typename IndexList, size_t Dim>
90 1 : void logical_partial_derivative(
91 : gsl::not_null<TensorMetafunctions::prepend_spatial_index<
92 : Tensor<DataType, SymmList, IndexList>, Dim, UpLo::Lo,
93 : Frame::ElementLogical>*>
94 : logical_derivative_of_u,
95 : gsl::not_null<gsl::span<typename DataType::value_type>*> buffer,
96 : const Tensor<DataType, SymmList, IndexList>& u, const Mesh<Dim>& mesh);
97 :
98 : template <typename DataType, typename SymmList, typename IndexList, size_t Dim>
99 1 : void logical_partial_derivative(
100 : gsl::not_null<TensorMetafunctions::prepend_spatial_index<
101 : Tensor<DataType, SymmList, IndexList>, Dim, UpLo::Lo,
102 : Frame::ElementLogical>*>
103 : logical_derivative_of_u,
104 : const Tensor<DataType, SymmList, IndexList>& u, const Mesh<Dim>& mesh);
105 :
106 : template <typename DataType, typename SymmList, typename IndexList, size_t Dim>
107 1 : auto logical_partial_derivative(const Tensor<DataType, SymmList, IndexList>& u,
108 : const Mesh<Dim>& mesh)
109 : -> TensorMetafunctions::prepend_spatial_index<
110 : Tensor<DataType, SymmList, IndexList>, Dim, UpLo::Lo,
111 : Frame::ElementLogical>;
112 : /// @}
113 :
114 : /// @{
115 : /*!
116 : * \ingroup NumericalAlgorithmsGroup
117 : * \brief Compute the partial derivatives of each variable, using the Cartoon
118 : * method for directions not in the computational domain.
119 : *
120 : * For a spacetime with some Killing vector \f$\xi^a\f$, not only does the
121 : * metric respect the isometry via \f$\mathcal{L}_\xi g_{ab} = 0\f$, but the
122 : * fields must as well.
123 : *
124 : * In the case of axial symmetry about the \f$y\f$-axis, the Killing vector is
125 : * \f$\xi^a = (0, -z, 0, x)\f$, so we have that for a vector \f$V^a\f$,
126 : *
127 : * \f{align*}{
128 : * \mathcal{L}_\xi V^a &= \xi^b \partial_b V^a - V^b \partial_b \xi^a
129 : * = -z \partial_x V^a + x \partial_z V^a - V^b \partial_b \xi^a = 0.
130 : * \f}
131 : *
132 : * Choosing the computational domain to be the \f$z=0\f$ plane, get the result:
133 : *
134 : * \f{align*}{
135 : * \partial_z V^a &= \frac{V^b \partial_b \xi^a}{x}.
136 : * \f}
137 : *
138 : * In the case that $x=0$ is in the computational domain, L'Hôpital's
139 : * rule is used. For a general tensor,
140 : * \f$T^{a_0,\ldots,a_n}{}_{b_0,\ldots,b_m}\f$ we have
141 : *
142 : * \f{align*}{
143 : * \partial_z T^{a_0,\ldots,a_n}{}_{b_0,\ldots,b_m} &=
144 : * \left\{
145 : * \begin{array}{ll}
146 : * \partial_x(\sum_k
147 : * T^{a_0,\ldots,c_k,\ldots,a_n}{}_{b_0,\ldots,b_m}\partial_{c_k}
148 : * \xi^{a_k} & \\
149 : * \;\;\;\;\;-\sum_k
150 : * T^{a_0,\ldots,a_n}{}_{b_0,\ldots,c_k,\ldots,b_m}\partial_{b_k}
151 : * \xi^{c_k}
152 : * ) & \mathrm{if} \; x=0 \\
153 : * (\sum_k
154 : * T^{a_0,\ldots,c_k,\ldots,a_n}{}_{b_0,\ldots,b_m}\partial_{c_k}
155 : * \xi^{a_k} & \\
156 : * \;-\sum_k
157 : * T^{a_0,\ldots,a_n}{}_{b_0,\ldots,c_k,\ldots,b_m}\partial_{b_k}
158 : * \xi^{c_k})\frac{1}{x} & \mathrm{otherwise}
159 : * \end{array}\right.
160 : * \f}
161 : *
162 : * In the case of spherical symmetry, another Killing vector
163 : * \f$\xi'^a = (0, -y, x, 0) \f$ is used as well, with the above equation
164 : * easily generalizing. In that case, the computational domain is only the
165 : * \f$x\f$-axis.
166 : *
167 : * This function computes the partial derivatives of _all_ variables in
168 : * `VariableTags`. The additional parameter `inertial_coords` is used for
169 : * division by the \f$x\f$ coordinates. If \f$x=0\f$ is included in the domain,
170 : * it is assumed to be present only at the first index and is handled by
171 : * L'Hôpital's rule.
172 : *
173 : * The mesh is required to have the Cartoon basis in the last and potentially
174 : * second-to-last coordinates and the inverse Jacobian is accordingly used only
175 : * in the first and potentially second dimensions.
176 : */
177 : template <typename ResultTags, typename VariableTags, size_t Dim,
178 : typename DerivativeFrame, Requires<Dim == 3> = nullptr>
179 1 : void cartoon_partial_derivatives(
180 : gsl::not_null<Variables<ResultTags>*> du, const Variables<VariableTags>& u,
181 : const Mesh<Dim>& mesh,
182 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
183 : DerivativeFrame>& inverse_jacobian_3d,
184 : const tnsr::I<DataVector, Dim, Frame::Inertial>& inertial_coords);
185 : /// @}
186 :
187 : /// @{
188 : /// \ingroup NumericalAlgorithmsGroup
189 : /// \brief Compute the partial derivatives of each variable with respect to
190 : /// the coordinates of `DerivativeFrame`.
191 : ///
192 : /// Either compute partial derivatives of _all_ variables in `VariableTags`, or
193 : /// of a subset of the `VariablesTags`. The subset of tags (`DerivativeTags`)
194 : /// must be the head of `VariablesTags`.
195 : ///
196 : /// The return-by-reference overload infers all template parameters from the
197 : /// arguments. The tensor types in the output buffer must have a spatial index
198 : /// appended to the front.
199 : ///
200 : /// The return-by-value overload requires that the `DerivativeTags` are
201 : /// specified explicitly as the first template parameter. It returns a
202 : /// `Variables` with the `DerivativeTags` wrapped in `Tags::deriv`.
203 : template <typename ResultTags, typename DerivativeTags, size_t Dim,
204 : typename DerivativeFrame>
205 1 : void partial_derivatives(
206 : gsl::not_null<Variables<ResultTags>*> du,
207 : const std::array<Variables<DerivativeTags>, Dim>&
208 : logical_partial_derivatives_of_u,
209 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
210 : DerivativeFrame>& inverse_jacobian);
211 :
212 : template <typename ResultTags, typename VariableTags, size_t Dim,
213 : typename DerivativeFrame>
214 1 : void partial_derivatives(
215 : gsl::not_null<Variables<ResultTags>*> du, const Variables<VariableTags>& u,
216 : const Mesh<Dim>& mesh,
217 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
218 : DerivativeFrame>& inverse_jacobian);
219 :
220 : template <typename DerivativeTags, typename VariableTags, size_t Dim,
221 : typename DerivativeFrame>
222 1 : auto partial_derivatives(
223 : const Variables<VariableTags>& u, const Mesh<Dim>& mesh,
224 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
225 : DerivativeFrame>& inverse_jacobian)
226 : -> Variables<db::wrap_tags_in<Tags::deriv, DerivativeTags,
227 : tmpl::size_t<Dim>, DerivativeFrame>>;
228 : /// @}
229 :
230 : /// @{
231 : /// \ingroup NumericalAlgorithmsGroup
232 : /// \brief Compute the partial derivative of a `Tensor` with respect to
233 : /// the coordinates of `DerivativeFrame`.
234 : ///
235 : /// Returns a `Tensor` with a spatial tensor index appended to the front
236 : /// of the input `Tensor`.
237 : ///
238 : /// If you have a `Variables` with several tensors you need to differentiate,
239 : /// you should use the `partial_derivatives` function that operates on
240 : /// `Variables` since that'll be more efficient.
241 : template <typename DataType, typename SymmList, typename IndexList, size_t Dim,
242 : typename DerivativeFrame>
243 1 : void partial_derivative(
244 : const gsl::not_null<TensorMetafunctions::prepend_spatial_index<
245 : Tensor<DataType, SymmList, IndexList>, Dim, UpLo::Lo, DerivativeFrame>*>
246 : du,
247 : const TensorMetafunctions::prepend_spatial_index<
248 : Tensor<DataType, SymmList, IndexList>, Dim, UpLo::Lo,
249 : Frame::ElementLogical>& logical_partial_derivative_of_u,
250 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
251 : DerivativeFrame>& inverse_jacobian);
252 :
253 : template <typename DataType, typename SymmList, typename IndexList, size_t Dim,
254 : typename DerivativeFrame>
255 1 : void partial_derivative(
256 : const gsl::not_null<TensorMetafunctions::prepend_spatial_index<
257 : Tensor<DataType, SymmList, IndexList>, Dim, UpLo::Lo, DerivativeFrame>*>
258 : du,
259 : const Tensor<DataType, SymmList, IndexList>& u, const Mesh<Dim>& mesh,
260 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
261 : DerivativeFrame>& inverse_jacobian);
262 :
263 : template <typename DataType, typename SymmList, typename IndexList, size_t Dim,
264 : typename DerivativeFrame>
265 1 : auto partial_derivative(
266 : const Tensor<DataType, SymmList, IndexList>& u, const Mesh<Dim>& mesh,
267 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
268 : DerivativeFrame>& inverse_jacobian)
269 : -> TensorMetafunctions::prepend_spatial_index<
270 : Tensor<DataType, SymmList, IndexList>, Dim, UpLo::Lo, DerivativeFrame>;
271 : /// @}
272 :
273 : /// @{
274 : /// \ingroup NumericalAlgorithmsGroup
275 : /// \brief Compute the partial derivative of a `Tensor` with respect to
276 : /// the coordinates of `DerivativeFrame` when a Cartoon basis is being used.
277 : ///
278 : /// Returns a `Tensor` with a spatial tensor index appended to the front
279 : /// of the input `Tensor`.
280 : ///
281 : /// If you have a `Variables` with several tensors with Cartoon bases you need
282 : /// to differentiate, you should use the `cartoon_partial_derivatives` function
283 : /// that operates on `Variables` since that'll be more efficient.
284 : template <typename SymmList, typename IndexList, size_t Dim,
285 : typename DerivativeFrame, Requires<Dim == 3> = nullptr>
286 1 : void cartoon_partial_derivative(
287 : gsl::not_null<TensorMetafunctions::prepend_spatial_index<
288 : Tensor<DataVector, SymmList, IndexList>, Dim, UpLo::Lo,
289 : DerivativeFrame>*>
290 : du,
291 : const Tensor<DataVector, SymmList, IndexList>& u, const Mesh<Dim>& mesh,
292 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
293 : DerivativeFrame>& inverse_jacobian,
294 : const tnsr::I<DataVector, Dim, Frame::Inertial>& inertial_coords);
295 :
296 : template <typename SymmList, typename IndexList, size_t Dim,
297 : typename DerivativeFrame, Requires<Dim == 3> = nullptr>
298 1 : auto cartoon_partial_derivative(
299 : const Tensor<DataVector, SymmList, IndexList>& u, const Mesh<Dim>& mesh,
300 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
301 : DerivativeFrame>& inverse_jacobian,
302 : const tnsr::I<DataVector, Dim, Frame::Inertial>& inertial_coords)
303 : -> TensorMetafunctions::prepend_spatial_index<
304 : Tensor<DataVector, SymmList, IndexList>, Dim, UpLo::Lo,
305 : DerivativeFrame>;
306 : /// @}
307 :
308 : /// @{
309 : /// \ingroup NumericalAlgorithmsGroup
310 : /// \brief Calls the correct partial derivatives function, either normal
311 : /// partials or cartoon partials, as determined by mesh basis.
312 : ///
313 : /// To be used in executables that allow a Cartoon basis, a choice that is
314 : /// only known at runtime.
315 : template <typename ResultTags, typename VariableTags, size_t Dim,
316 : typename DerivativeFrame>
317 1 : void partial_derivatives(
318 : gsl::not_null<Variables<ResultTags>*> du,
319 : const Variables<VariableTags>& u, const Mesh<Dim>& mesh,
320 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
321 : DerivativeFrame>& inverse_jacobian,
322 : const tnsr::I<DataVector, Dim, Frame::Inertial>& inertial_coords);
323 : /// @}
324 :
325 : /// @{
326 : /// \ingroup NumericalAlgorithmsGroup
327 : /// \brief Calls the correct partial derivative function, either normal
328 : /// partial or cartoon partial, as determined by mesh basis.
329 : ///
330 : /// To be used in executables that allow a Cartoon basis, a choice that is
331 : /// only known at runtime.
332 : template <typename SymmList, typename IndexList, size_t Dim,
333 : typename DerivativeFrame>
334 1 : void partial_derivative(
335 : gsl::not_null<TensorMetafunctions::prepend_spatial_index<
336 : Tensor<DataVector, SymmList, IndexList>, Dim, UpLo::Lo,
337 : DerivativeFrame>*>
338 : du,
339 : const Tensor<DataVector, SymmList, IndexList>& u, const Mesh<Dim>& mesh,
340 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
341 : DerivativeFrame>& inverse_jacobian,
342 : const tnsr::I<DataVector, Dim, Frame::Inertial>& inertial_coords);
343 :
344 : template <typename SymmList, typename IndexList, size_t Dim,
345 : typename DerivativeFrame>
346 1 : auto partial_derivative(
347 : const Tensor<DataVector, SymmList, IndexList>& u, const Mesh<Dim>& mesh,
348 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
349 : DerivativeFrame>& inverse_jacobian,
350 : const tnsr::I<DataVector, Dim, Frame::Inertial>& inertial_coords)
351 : -> TensorMetafunctions::prepend_spatial_index<
352 : Tensor<DataVector, SymmList, IndexList>, Dim, UpLo::Lo,
353 : DerivativeFrame>;
354 : /// @}
355 :
356 : namespace Tags {
357 :
358 : /*!
359 : * \ingroup DataBoxTagsGroup
360 : * \brief Compute the spatial derivatives of tags in a Variables
361 : *
362 : * Computes the spatial derivatives of the Tensors in the Variables represented
363 : * by `VariablesTag` in the frame mapped to by `InverseJacobianTag`. To only
364 : * take the derivatives of a subset of these Tensors you can set the
365 : * `DerivTags` template parameter. It takes a `tmpl::list` of the desired
366 : * tags and defaults to the full `tags_list` of the Variables.
367 : *
368 : * For an executable that does not allow a Cartoon basis, the last parameter,
369 : * `InertialCoordsTag`, should not be passed.
370 : *
371 : * This tag may be retrieved via `::Tags::Variables<db::wrap_tags_in<deriv,
372 : * DerivTags, Dim, deriv_frame>`.
373 : */
374 : template <typename VariablesTag, typename MeshTag, typename InverseJacobianTag,
375 : typename DerivTags = typename VariablesTag::type::tags_list,
376 : typename InertialCoordsTag = void>
377 1 : struct DerivCompute
378 : : db::add_tag_prefix<
379 : deriv, ::Tags::Variables<DerivTags>,
380 : tmpl::size_t<
381 : tmpl::back<typename InverseJacobianTag::type::index_list>::dim>,
382 : typename tmpl::back<
383 : typename InverseJacobianTag::type::index_list>::Frame>,
384 : db::ComputeTag {
385 : private:
386 0 : using inv_jac_indices = typename InverseJacobianTag::type::index_list;
387 0 : static constexpr auto Dim = tmpl::back<inv_jac_indices>::dim;
388 0 : using deriv_frame = typename tmpl::back<inv_jac_indices>::Frame;
389 :
390 : public:
391 0 : using base = db::add_tag_prefix<
392 : deriv, ::Tags::Variables<DerivTags>,
393 : tmpl::size_t<
394 : tmpl::back<typename InverseJacobianTag::type::index_list>::dim>,
395 : typename tmpl::back<
396 : typename InverseJacobianTag::type::index_list>::Frame>;
397 0 : using return_type = typename base::type;
398 0 : static constexpr void function(
399 : gsl::not_null<return_type*> du, const typename VariablesTag::type& u,
400 : const Mesh<Dim>& mesh,
401 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
402 : deriv_frame>& inverse_jacobian) {
403 : partial_derivatives(du, u, mesh, inverse_jacobian);
404 : }
405 0 : static constexpr void function(
406 : gsl::not_null<return_type*> du, const typename VariablesTag::type& u,
407 : const Mesh<Dim>& mesh,
408 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
409 : deriv_frame>& inverse_jacobian,
410 : const tnsr::I<DataVector, Dim, Frame::Inertial>& inertial_coords) {
411 : partial_derivatives(du, u, mesh, inverse_jacobian, inertial_coords);
412 : }
413 0 : using argument_tags = tmpl::conditional_t<
414 : std::is_same_v<void, InertialCoordsTag>,
415 : tmpl::list<VariablesTag, MeshTag, InverseJacobianTag>,
416 : tmpl::list<VariablesTag, MeshTag, InverseJacobianTag, InertialCoordsTag>>;
417 : };
418 :
419 : /*!
420 : * \ingroup DataBoxTagsGroup
421 : * \brief Computes the spatial derivative of a single tensor tag not in a
422 : * Variables.
423 : *
424 : * Computes the spatial derivative for a single tensor represented by
425 : * 'TensorTag' in the frame mapped to by 'InverseJacobianTag'. It takes a
426 : * single Tensor designated by 'TensorTag', the inverse Jacobian, and a mesh.
427 : */
428 : template <typename TensorTag, typename InverseJacobianTag, typename MeshTag>
429 1 : struct DerivTensorCompute
430 : : ::Tags::deriv<TensorTag,
431 : tmpl::size_t<tmpl::back<
432 : typename InverseJacobianTag::type::index_list>::dim>,
433 : typename tmpl::back<
434 : typename InverseJacobianTag::type::index_list>::Frame>,
435 : db::ComputeTag {
436 : private:
437 0 : using inv_jac_indices = typename InverseJacobianTag::type::index_list;
438 0 : static constexpr auto Dim = tmpl::back<inv_jac_indices>::dim;
439 0 : using deriv_frame = typename tmpl::back<inv_jac_indices>::Frame;
440 :
441 : public:
442 0 : using base = ::Tags::deriv<TensorTag, tmpl::size_t<Dim>, deriv_frame>;
443 0 : using return_type = typename base::type;
444 0 : static constexpr void (*function)(
445 : gsl::not_null<return_type*>, const typename TensorTag::type&,
446 : const Mesh<Dim>&,
447 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
448 : deriv_frame>&) =
449 : partial_derivative<
450 : typename TensorTag::type::type, typename TensorTag::type::symmetry,
451 : typename TensorTag::type::index_list, Dim, deriv_frame>;
452 0 : using argument_tags =
453 : tmpl::list<TensorTag, domain::Tags::Mesh<Dim>, InverseJacobianTag>;
454 : };
455 : } // namespace Tags
|
