Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <deque>
7 : #include <utility>
8 :
9 : #include "DataStructures/Tensor/IndexType.hpp"
10 : #include "DataStructures/Tensor/TypeAliases.hpp"
11 : #include "NumericalAlgorithms/SphericalHarmonics/TagsTypeAliases.hpp"
12 :
13 : /// \cond
14 : class DataVector;
15 : namespace ylm {
16 : template <typename Fr>
17 : class Strahlkorper;
18 : } // namespace ylm
19 : template <typename X, typename Symm, typename IndexList>
20 : class Tensor;
21 : namespace gsl {
22 : template <typename>
23 : struct not_null;
24 : } // namespace gsl
25 : /// \endcond
26 :
27 : /// \ingroup SurfacesGroup
28 : /// Contains functions that depend on a Strahlkorper but not on a metric.
29 : namespace ylm {
30 : /// @{
31 : /*!
32 : * \f$(\theta,\phi)\f$ on the Strahlkorper surface.
33 : * Doesn't depend on the shape of the surface.
34 : *
35 : * We need to choose upper vs lower indices for theta_phi; it doesn't
36 : * matter because these are coordinates and not geometric objects, so
37 : * we choose lower indices arbitrarily.
38 : */
39 : template <typename Fr>
40 1 : tnsr::i<DataVector, 2, ::Frame::Spherical<Fr>> theta_phi(
41 : const Strahlkorper<Fr>& strahlkorper);
42 :
43 : template <typename Fr>
44 1 : void theta_phi(
45 : const gsl::not_null<tnsr::i<DataVector, 2, ::Frame::Spherical<Fr>>*>
46 : theta_phi,
47 : const Strahlkorper<Fr>& strahlkorper);
48 : /// @}
49 :
50 : /// @{
51 : /*!
52 : * `r_hat(i)` is \f$\hat{r}_i = x_i/\sqrt{x^2+y^2+z^2}\f$ on the
53 : * Strahlkorper surface. Doesn't depend on the shape of the surface.
54 : *
55 : * We need to choose upper vs lower indices for rhat; it doesn't
56 : * matter because rhat is a quantity defined with a Euclidean metric,
57 : * so we choose the lower index arbitrarily.
58 : */
59 : template <typename Fr>
60 1 : tnsr::i<DataVector, 3, Fr> rhat(
61 : const tnsr::i<DataVector, 2, ::Frame::Spherical<Fr>>& theta_phi);
62 :
63 : template <typename Fr>
64 1 : void rhat(const gsl::not_null<tnsr::i<DataVector, 3, Fr>*> r_hat,
65 : const tnsr::i<DataVector, 2, ::Frame::Spherical<Fr>>& theta_phi);
66 : /// @}
67 :
68 : /// @{
69 : /*!
70 : * `Jacobian(i,0)` is \f$\frac{1}{r}\partial x^i/\partial\theta\f$,
71 : * and `Jacobian(i,1)`
72 : * is \f$\frac{1}{r\sin\theta}\partial x^i/\partial\phi\f$.
73 : * Here \f$r\f$ means \f$\sqrt{x^2+y^2+z^2}\f$.
74 : * `Jacobian` doesn't depend on the shape of the surface.
75 : */
76 : template <typename Fr>
77 1 : ylm::Tags::aliases::Jacobian<Fr> jacobian(
78 : const tnsr::i<DataVector, 2, ::Frame::Spherical<Fr>>& theta_phi);
79 :
80 : template <typename Fr>
81 1 : void jacobian(const gsl::not_null<ylm::Tags::aliases::Jacobian<Fr>*> jac,
82 : const tnsr::i<DataVector, 2, ::Frame::Spherical<Fr>>& theta_phi);
83 : /// @}
84 :
85 : /// @{
86 : /*!
87 : * `InvJacobian(0,i)` is \f$r\partial\theta/\partial x^i\f$,
88 : * and `InvJacobian(1,i)` is \f$r\sin\theta\partial\phi/\partial x^i\f$.
89 : * Here \f$r\f$ means \f$\sqrt{x^2+y^2+z^2}\f$.
90 : * `InvJacobian` doesn't depend on the shape of the surface.
91 : */
92 : template <typename Fr>
93 1 : ylm::Tags::aliases::InvJacobian<Fr> inv_jacobian(
94 : const tnsr::i<DataVector, 2, ::Frame::Spherical<Fr>>& theta_phi);
95 :
96 : template <typename Fr>
97 1 : void inv_jacobian(
98 : const gsl::not_null<ylm::Tags::aliases::InvJacobian<Fr>*> inv_jac,
99 : const tnsr::i<DataVector, 2, ::Frame::Spherical<Fr>>& theta_phi);
100 : /// @}
101 :
102 : /// @{
103 : /*!
104 : * `InvHessian(k,i,j)` is \f$r\partial (J^{-1}){}^k_j/\partial x^i\f$,
105 : * where \f$(J^{-1}){}^k_j\f$ is the inverse Jacobian.
106 : * Here \f$r\f$ means \f$\sqrt{x^2+y^2+z^2}\f$.
107 : * `InvHessian` is not symmetric because the Jacobians are Pfaffian.
108 : * `InvHessian` doesn't depend on the shape of the surface.
109 : */
110 : template <typename Fr>
111 1 : ylm::Tags::aliases::InvHessian<Fr> inv_hessian(
112 : const tnsr::i<DataVector, 2, ::Frame::Spherical<Fr>>& theta_phi);
113 :
114 : template <typename Fr>
115 1 : void inv_hessian(
116 : const gsl::not_null<ylm::Tags::aliases::InvHessian<Fr>*> inv_hess,
117 : const tnsr::i<DataVector, 2, ::Frame::Spherical<Fr>>& theta_phi);
118 : /// @}
119 :
120 : /// @{
121 : /*!
122 : * (Euclidean) distance \f$r_{\rm surf}(\theta,\phi)\f$ from the
123 : * expansion center to each point of the Strahlkorper surface.
124 : */
125 : template <typename Fr>
126 1 : Scalar<DataVector> radius(const Strahlkorper<Fr>& strahlkorper);
127 :
128 : template <typename Fr>
129 1 : void radius(const gsl::not_null<Scalar<DataVector>*> result,
130 : const Strahlkorper<Fr>& strahlkorper);
131 : /// @}
132 :
133 : /// @{
134 : /*!
135 : * `cartesian_coords(i)` is \f$x_{\rm surf}^i\f$, the vector of \f$(x,y,z)\f$
136 : * coordinates of each point on the Strahlkorper surface.
137 : *
138 : * \param strahlkorper The Strahlkorper surface.
139 : * \param radius The radius as a function of angle, as returned by
140 : * `ylm::radius`.
141 : * \param r_hat The Euclidean radial unit vector as returned by
142 : * `ylm::rhat`.
143 : */
144 : template <typename Fr>
145 1 : tnsr::I<DataVector, 3, Fr> cartesian_coords(
146 : const Strahlkorper<Fr>& strahlkorper, const Scalar<DataVector>& radius,
147 : const tnsr::i<DataVector, 3, Fr>& r_hat);
148 :
149 : /*!
150 : * \param coords The returned Cartesian coordinates.
151 : * \param strahlkorper The Strahlkorper surface.
152 : * \param radius The radius as a function of angle, as returned by
153 : * `ylm::radius`.
154 : * \param r_hat The Euclidean radial unit vector as returned by
155 : * `ylm::rhat`.
156 : */
157 : template <typename Fr>
158 1 : void cartesian_coords(const gsl::not_null<tnsr::I<DataVector, 3, Fr>*> coords,
159 : const Strahlkorper<Fr>& strahlkorper,
160 : const Scalar<DataVector>& radius,
161 : const tnsr::i<DataVector, 3, Fr>& r_hat);
162 : /// @}
163 :
164 : /// @{
165 : /*!
166 : * `dx_scalar(i)` is \f$\partial f/\partial x^i\f$ evaluated on the
167 : * surface. Here \f$f=f(r,\theta,\phi)=f(\theta,\phi)\f$ is some
168 : * scalar function independent of the radial coordinate. \f$f\f$ is
169 : * considered a function of Cartesian coordinates
170 : * \f$f=f(\theta(x,y,z),\phi(x,y,z))\f$ for this operation.
171 : *
172 : * \param scalar The scalar to be differentiated.
173 : * \param strahlkorper The Strahlkorper surface.
174 : * \param radius_of_strahlkorper The radius of the Strahlkorper at each
175 : * point, as returned by `ylm::radius`.
176 : * \param inv_jac The inverse Jacobian as returned by
177 : * `ylm::inv_jacobian`
178 : */
179 : template <typename Fr>
180 1 : tnsr::i<DataVector, 3, Fr> cartesian_derivs_of_scalar(
181 : const Scalar<DataVector>& scalar, const Strahlkorper<Fr>& strahlkorper,
182 : const Scalar<DataVector>& radius_of_strahlkorper,
183 : const ylm::Tags::aliases::InvJacobian<Fr>& inv_jac);
184 :
185 : /*!
186 : * \param dx_scalar The returned derivatives of the scalar.
187 : * \param scalar The scalar to be differentiated.
188 : * \param strahlkorper The Strahlkorper surface.
189 : * \param radius_of_strahlkorper The radius of the Strahlkorper at each
190 : * point, as returned by `ylm::radius`.
191 : * \param inv_jac The inverse Jacobian as returned by
192 : * `ylm::inv_jacobian`
193 : */
194 : template <typename Fr>
195 1 : void cartesian_derivs_of_scalar(
196 : const gsl::not_null<tnsr::i<DataVector, 3, Fr>*> dx_scalar,
197 : const Scalar<DataVector>& scalar, const Strahlkorper<Fr>& strahlkorper,
198 : const Scalar<DataVector>& radius_of_strahlkorper,
199 : const ylm::Tags::aliases::InvJacobian<Fr>& inv_jac);
200 : /// @}
201 :
202 : /// @{
203 : /*!
204 : * `d2x_scalar(i,j)` is \f$\partial^2 f/\partial x^i\partial x^j\f$
205 : * evaluated on the surface. Here
206 : * \f$f=f(r,\theta,\phi)=f(\theta,\phi)\f$ is some scalar function
207 : * independent of the radial coordinate. \f$f\f$ is considered a
208 : * function of Cartesian coordinates
209 : * \f$f=f(\theta(x,y,z),\phi(x,y,z))\f$ for this operation.
210 : *
211 : * \param scalar The scalar to be differentiated.
212 : * \param strahlkorper The Strahlkorper surface.
213 : * \param radius_of_strahlkorper The radius of the Strahlkorper at each
214 : * point, as returned by `ylm::radius`.
215 : * \param inv_jac The inverse Jacobian as returned by
216 : * `ylm::inv_jacobian`
217 : * \param inv_hess The inverse Hessian as returned by
218 : * `ylm::inv_hessian.
219 : */
220 : template <typename Fr>
221 1 : tnsr::ii<DataVector, 3, Fr> cartesian_second_derivs_of_scalar(
222 : const Scalar<DataVector>& scalar, const Strahlkorper<Fr>& strahlkorper,
223 : const Scalar<DataVector>& radius_of_strahlkorper,
224 : const ylm::Tags::aliases::InvJacobian<Fr>& inv_jac,
225 : const ylm::Tags::aliases::InvHessian<Fr>& inv_hess);
226 :
227 : /*!
228 : * \param d2x_scalar The returned 2nd derivatives of the scalar.
229 : * \param scalar The scalar to be differentiated.
230 : * \param strahlkorper The Strahlkorper surface.
231 : * \param radius_of_strahlkorper The radius of the Strahlkorper at each
232 : * point, as returned by `ylm::radius`.
233 : * \param inv_jac The inverse Jacobian as returned by
234 : * `ylm::inv_jacobian`
235 : * \param inv_hess The inverse Hessian as returned by
236 : * `ylm::inv_hessian.
237 : */
238 : template <typename Fr>
239 1 : void cartesian_second_derivs_of_scalar(
240 : const gsl::not_null<tnsr::ii<DataVector, 3, Fr>*> d2x_scalar,
241 : const Scalar<DataVector>& scalar, const Strahlkorper<Fr>& strahlkorper,
242 : const Scalar<DataVector>& radius_of_strahlkorper,
243 : const ylm::Tags::aliases::InvJacobian<Fr>& inv_jac,
244 : const ylm::Tags::aliases::InvHessian<Fr>& inv_hess);
245 : /// @}
246 :
247 : /// @{
248 : /*!
249 : * \f$\nabla^2 f\f$, the flat Laplacian of a scalar \f$f\f$ on the surface.
250 : * This is \f$\eta^{ij}\partial^2 f/\partial x^i\partial x^j\f$,
251 : * where \f$f=f(r,\theta,\phi)=f(\theta,\phi)\f$ is some scalar function
252 : * independent of the radial coordinate. \f$f\f$ is considered a
253 : * function of Cartesian coordinates
254 : * \f$f=f(\theta(x,y,z),\phi(x,y,z))\f$ for this operation.
255 : *
256 : */
257 : template <typename Fr>
258 1 : Scalar<DataVector> laplacian_of_scalar(
259 : const Scalar<DataVector>& scalar, const Strahlkorper<Fr>& strahlkorper,
260 : const tnsr::i<DataVector, 2, ::Frame::Spherical<Fr>>& theta_phi);
261 :
262 : template <typename Fr>
263 1 : void laplacian_of_scalar(
264 : const gsl::not_null<Scalar<DataVector>*> laplacian,
265 : const Scalar<DataVector>& scalar, const Strahlkorper<Fr>& strahlkorper,
266 : const tnsr::i<DataVector, 2, ::Frame::Spherical<Fr>>& theta_phi);
267 : /// @}
268 :
269 : /// @{
270 : /*!
271 : * `tangents(i,j)` is \f$\partial x_{\rm surf}^i/\partial q^j\f$,
272 : * where \f$x_{\rm surf}^i\f$ are the Cartesian coordinates of the
273 : * surface (i.e. `cartesian_coords`) and are considered functions of
274 : * \f$(\theta,\phi)\f$.
275 : *
276 : * \f$\partial/\partial q^0\f$ means
277 : * \f$\partial/\partial\theta\f$; and \f$\partial/\partial q^1\f$
278 : * means \f$\csc\theta\,\,\partial/\partial\phi\f$. Note that the
279 : * vectors `tangents(i,0)` and `tangents(i,1)` are orthogonal to the
280 : * `normal_one_form` \f$s_i\f$, i.e.
281 : * \f$s_i \partial x_{\rm surf}^i/\partial q^j = 0\f$; this statement
282 : * is independent of a metric. Also, `tangents(i,0)` and
283 : * `tangents(i,1)` are not necessarily orthogonal to each other,
284 : * since orthogonality between 2 vectors (as opposed to a vector and
285 : * a one-form) is metric-dependent.
286 : *
287 : * \param strahlkorper The Strahlkorper surface.
288 : * \param radius The radius of the Strahlkorper at each
289 : * point, as returned by `ylm::radius`.
290 : * \param r_hat The radial unit vector as returned by
291 : * `ylm::rhat`.
292 : * \param jac The jacobian as returned by `ylm::jacobian`.
293 : */
294 : template <typename Fr>
295 1 : ylm::Tags::aliases::Jacobian<Fr> tangents(
296 : const Strahlkorper<Fr>& strahlkorper, const Scalar<DataVector>& radius,
297 : const tnsr::i<DataVector, 3, Fr>& r_hat,
298 : const ylm::Tags::aliases::Jacobian<Fr>& jac);
299 :
300 : /*!
301 : * \param result The computed tangent vectors.
302 : * \param strahlkorper The Strahlkorper surface.
303 : * \param radius The radius of the Strahlkorper at each
304 : * point, as returned by `ylm::radius`.
305 : * \param r_hat The radial unit vector as returned by
306 : * `ylm::rhat`.
307 : * \param jac The jacobian as returned by `ylm::jacobian`.
308 : */
309 : template <typename Fr>
310 1 : void tangents(const gsl::not_null<ylm::Tags::aliases::Jacobian<Fr>*> result,
311 : const Strahlkorper<Fr>& strahlkorper,
312 : const Scalar<DataVector>& radius,
313 : const tnsr::i<DataVector, 3, Fr>& r_hat,
314 : const ylm::Tags::aliases::Jacobian<Fr>& jac);
315 : /// @}
316 :
317 : /// @{
318 : /*!
319 : * `normal_one_form(i)` is \f$s_i\f$, the (unnormalized) normal one-form
320 : * to the surface, expressed in Cartesian components.
321 : * This is computed by \f$x_i/r-\partial r_{\rm surf}/\partial x^i\f$,
322 : * where \f$x_i/r\f$ is `Rhat` and
323 : * \f$\partial r_{\rm surf}/\partial x^i\f$ is `DxRadius`.
324 : * See Eq. (8) of \cite Baumgarte1996hh.
325 : * Note on the word "normal": \f$s_i\f$ points in the correct direction
326 : * (it is "normal" to the surface), but it does not have unit length
327 : * (it is not "normalized"; normalization requires a metric).
328 : *
329 : * \param dx_radius The Cartesian derivatives of the radius, as
330 : * returned by ylm::cartesian_derivs_of_scalar with
331 : * `ylm::radius` passed in as the scalar.
332 : * \param r_hat The radial unit vector as returned by
333 : * `ylm::rhat`.
334 : */
335 : template <typename Fr>
336 1 : tnsr::i<DataVector, 3, Fr> normal_one_form(
337 : const tnsr::i<DataVector, 3, Fr>& dx_radius,
338 : const tnsr::i<DataVector, 3, Fr>& r_hat);
339 :
340 : /*!
341 : * \param one_form The returned normal one form.
342 : * \param dx_radius The Cartesian derivatives of the radius, as
343 : * returned by ylm::cartesian_derivs_of_scalar with
344 : * `ylm::radius` passed in as the scalar.
345 : * \param r_hat The radial unit vector as returned by
346 : * `ylm::rhat`.
347 : */
348 : template <typename Fr>
349 1 : void normal_one_form(const gsl::not_null<tnsr::i<DataVector, 3, Fr>*> one_form,
350 : const tnsr::i<DataVector, 3, Fr>& dx_radius,
351 : const tnsr::i<DataVector, 3, Fr>& r_hat);
352 :
353 : /// @{
354 : /*!
355 : * The linear least squares fit of the polynomial of order 3
356 : * given a `std::vector` of `Strahlkorper`s to their \f$Y_l^m\f$ coefficients.
357 : * Assumes the the \f$l_{\max}\f$ and \f$m_{\max}\f$ of each `Strahlkorper`
358 : * are the same, and the returned vector consists of \f$2l_{\max}m_{\max}\f$
359 : * (the number of \f$Y_l^m\f$ coefficients) `std::array<double, 4>`s, each of
360 : * which consists of the four coefficients that define the best fit cubic to
361 : * each \f$Y_l^m\f$ coefficient of the `Strahlkorper` as a function of time.
362 : *
363 : * \param times The time corresponding to each `Strahlkorper` to be fit to.
364 : * \param strahlkorpers The `Strahlkorper` surfaces which consists of a set
365 : * of \f$Y_l^m\f$ coefficients corresponding to the shape of the `Strahlkorper`
366 : * at a particular time.
367 : */
368 : template <typename Fr>
369 1 : std::vector<std::array<double, 4>> fit_ylm_coeffs(
370 : const DataVector& times,
371 : const std::vector<Strahlkorper<Fr>>& strahlkorpers);
372 :
373 : /*!
374 : * \brief Compute the time derivative of a Strahlkorper from a number of
375 : * previous Strahlkorpers
376 : *
377 : * \details Does simple 1D FD with non-uniform spacing using
378 : * `fd::non_uniform_1d_weights`.
379 : * \param time_deriv Strahlkorper whose coefficients are the time derivative of
380 : * `previous_strahlkorpers`' coefficients.
381 : * \param previous_strahlkorpers All previous Strahlkorpers and the times they
382 : * are at. They are expected to have the most recent Strahlkorper in the front
383 : * and the Strahlkorper furthest in the past in the back of the deque.
384 : */
385 : template <typename Frame>
386 1 : void time_deriv_of_strahlkorper(
387 : gsl::not_null<Strahlkorper<Frame>*> time_deriv,
388 : const std::deque<std::pair<double, Strahlkorper<Frame>>>&
389 : previous_strahlkorpers);
390 : } // namespace ylm
|
