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/Strahlkorper/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 : * This overload computes `radius`, `theta_phi`, and `r_hat` internally.
165 : * Use the other overloads if you already have these quantities.
166 : *
167 : * \param strahlkorper The Strahlkorper surface.
168 : */
169 : template <typename Fr>
170 1 : tnsr::I<DataVector, 3, Fr> cartesian_coords(
171 : const Strahlkorper<Fr>& strahlkorper);
172 : /// @}
173 :
174 : /// @{
175 : /*!
176 : * `dx_scalar(i)` is \f$\partial f/\partial x^i\f$ evaluated on the
177 : * surface. Here \f$f=f(r,\theta,\phi)=f(\theta,\phi)\f$ is some
178 : * scalar function independent of the radial coordinate. \f$f\f$ is
179 : * considered a function of Cartesian coordinates
180 : * \f$f=f(\theta(x,y,z),\phi(x,y,z))\f$ for this operation.
181 : *
182 : * \param scalar The scalar to be differentiated.
183 : * \param strahlkorper The Strahlkorper surface.
184 : * \param radius_of_strahlkorper The radius of the Strahlkorper at each
185 : * point, as returned by `ylm::radius`.
186 : * \param inv_jac The inverse Jacobian as returned by
187 : * `ylm::inv_jacobian`
188 : */
189 : template <typename Fr>
190 1 : tnsr::i<DataVector, 3, Fr> cartesian_derivs_of_scalar(
191 : const Scalar<DataVector>& scalar, const Strahlkorper<Fr>& strahlkorper,
192 : const Scalar<DataVector>& radius_of_strahlkorper,
193 : const ylm::Tags::aliases::InvJacobian<Fr>& inv_jac);
194 :
195 : /*!
196 : * \param dx_scalar The returned derivatives of the scalar.
197 : * \param scalar The scalar to be differentiated.
198 : * \param strahlkorper The Strahlkorper surface.
199 : * \param radius_of_strahlkorper The radius of the Strahlkorper at each
200 : * point, as returned by `ylm::radius`.
201 : * \param inv_jac The inverse Jacobian as returned by
202 : * `ylm::inv_jacobian`
203 : */
204 : template <typename Fr>
205 1 : void cartesian_derivs_of_scalar(
206 : const gsl::not_null<tnsr::i<DataVector, 3, Fr>*> dx_scalar,
207 : const Scalar<DataVector>& scalar, const Strahlkorper<Fr>& strahlkorper,
208 : const Scalar<DataVector>& radius_of_strahlkorper,
209 : const ylm::Tags::aliases::InvJacobian<Fr>& inv_jac);
210 : /// @}
211 :
212 : /// @{
213 : /*!
214 : * `d2x_scalar(i,j)` is \f$\partial^2 f/\partial x^i\partial x^j\f$
215 : * evaluated on the surface. Here
216 : * \f$f=f(r,\theta,\phi)=f(\theta,\phi)\f$ is some scalar function
217 : * independent of the radial coordinate. \f$f\f$ is considered a
218 : * function of Cartesian coordinates
219 : * \f$f=f(\theta(x,y,z),\phi(x,y,z))\f$ for this operation.
220 : *
221 : * \param scalar The scalar to be differentiated.
222 : * \param strahlkorper The Strahlkorper surface.
223 : * \param radius_of_strahlkorper The radius of the Strahlkorper at each
224 : * point, as returned by `ylm::radius`.
225 : * \param inv_jac The inverse Jacobian as returned by
226 : * `ylm::inv_jacobian`
227 : * \param inv_hess The inverse Hessian as returned by
228 : * `ylm::inv_hessian`.
229 : */
230 : template <typename Fr>
231 1 : tnsr::ii<DataVector, 3, Fr> cartesian_second_derivs_of_scalar(
232 : const Scalar<DataVector>& scalar, const Strahlkorper<Fr>& strahlkorper,
233 : const Scalar<DataVector>& radius_of_strahlkorper,
234 : const ylm::Tags::aliases::InvJacobian<Fr>& inv_jac,
235 : const ylm::Tags::aliases::InvHessian<Fr>& inv_hess);
236 :
237 : /*!
238 : * \param d2x_scalar The returned 2nd derivatives of the scalar.
239 : * \param scalar The scalar to be differentiated.
240 : * \param strahlkorper The Strahlkorper surface.
241 : * \param radius_of_strahlkorper The radius of the Strahlkorper at each
242 : * point, as returned by `ylm::radius`.
243 : * \param inv_jac The inverse Jacobian as returned by
244 : * `ylm::inv_jacobian`
245 : * \param inv_hess The inverse Hessian as returned by
246 : * `ylm::inv_hessian`.
247 : */
248 : template <typename Fr>
249 1 : void cartesian_second_derivs_of_scalar(
250 : const gsl::not_null<tnsr::ii<DataVector, 3, Fr>*> d2x_scalar,
251 : const Scalar<DataVector>& scalar, const Strahlkorper<Fr>& strahlkorper,
252 : const Scalar<DataVector>& radius_of_strahlkorper,
253 : const ylm::Tags::aliases::InvJacobian<Fr>& inv_jac,
254 : const ylm::Tags::aliases::InvHessian<Fr>& inv_hess);
255 : /// @}
256 :
257 : /// @{
258 : /*!
259 : * \f$\nabla^2 f\f$, the flat Laplacian of a scalar \f$f\f$ on the surface.
260 : * This is \f$\eta^{ij}\partial^2 f/\partial x^i\partial x^j\f$,
261 : * where \f$f=f(r,\theta,\phi)=f(\theta,\phi)\f$ is some scalar function
262 : * independent of the radial coordinate. \f$f\f$ is considered a
263 : * function of Cartesian coordinates
264 : * \f$f=f(\theta(x,y,z),\phi(x,y,z))\f$ for this operation.
265 : *
266 : */
267 : template <typename Fr>
268 1 : Scalar<DataVector> laplacian_of_scalar(
269 : const Scalar<DataVector>& scalar, const Strahlkorper<Fr>& strahlkorper,
270 : const tnsr::i<DataVector, 2, ::Frame::Spherical<Fr>>& theta_phi);
271 :
272 : template <typename Fr>
273 1 : void laplacian_of_scalar(
274 : const gsl::not_null<Scalar<DataVector>*> laplacian,
275 : const Scalar<DataVector>& scalar, const Strahlkorper<Fr>& strahlkorper,
276 : const tnsr::i<DataVector, 2, ::Frame::Spherical<Fr>>& theta_phi);
277 : /// @}
278 :
279 : /// @{
280 : /*!
281 : * `tangents(i,j)` is \f$\partial x_{\rm surf}^i/\partial q^j\f$,
282 : * where \f$x_{\rm surf}^i\f$ are the Cartesian coordinates of the
283 : * surface (i.e. `cartesian_coords`) and are considered functions of
284 : * \f$(\theta,\phi)\f$.
285 : *
286 : * \f$\partial/\partial q^0\f$ means
287 : * \f$\partial/\partial\theta\f$; and \f$\partial/\partial q^1\f$
288 : * means \f$\csc\theta\,\,\partial/\partial\phi\f$. Note that the
289 : * vectors `tangents(i,0)` and `tangents(i,1)` are orthogonal to the
290 : * `normal_one_form` \f$s_i\f$, i.e.
291 : * \f$s_i \partial x_{\rm surf}^i/\partial q^j = 0\f$; this statement
292 : * is independent of a metric. Also, `tangents(i,0)` and
293 : * `tangents(i,1)` are not necessarily orthogonal to each other,
294 : * since orthogonality between 2 vectors (as opposed to a vector and
295 : * a one-form) is metric-dependent.
296 : *
297 : * \param strahlkorper The Strahlkorper surface.
298 : * \param radius The radius of the Strahlkorper at each
299 : * point, as returned by `ylm::radius`.
300 : * \param r_hat The radial unit vector as returned by
301 : * `ylm::rhat`.
302 : * \param jac The jacobian as returned by `ylm::jacobian`.
303 : */
304 : template <typename Fr>
305 1 : ylm::Tags::aliases::Jacobian<Fr> tangents(
306 : const Strahlkorper<Fr>& strahlkorper, const Scalar<DataVector>& radius,
307 : const tnsr::i<DataVector, 3, Fr>& r_hat,
308 : const ylm::Tags::aliases::Jacobian<Fr>& jac);
309 :
310 : /*!
311 : * \param result The computed tangent vectors.
312 : * \param strahlkorper The Strahlkorper surface.
313 : * \param radius The radius of the Strahlkorper at each
314 : * point, as returned by `ylm::radius`.
315 : * \param r_hat The radial unit vector as returned by
316 : * `ylm::rhat`.
317 : * \param jac The jacobian as returned by `ylm::jacobian`.
318 : */
319 : template <typename Fr>
320 1 : void tangents(const gsl::not_null<ylm::Tags::aliases::Jacobian<Fr>*> result,
321 : const Strahlkorper<Fr>& strahlkorper,
322 : const Scalar<DataVector>& radius,
323 : const tnsr::i<DataVector, 3, Fr>& r_hat,
324 : const ylm::Tags::aliases::Jacobian<Fr>& jac);
325 : /// @}
326 :
327 : /// @{
328 : /*!
329 : * `normal_one_form(i)` is \f$s_i\f$, the (unnormalized) normal one-form
330 : * to the surface, expressed in Cartesian components.
331 : * This is computed by \f$x_i/r-\partial r_{\rm surf}/\partial x^i\f$,
332 : * where \f$x_i/r\f$ is `Rhat` and
333 : * \f$\partial r_{\rm surf}/\partial x^i\f$ is `DxRadius`.
334 : * See Eq. (8) of \cite Baumgarte1996hh.
335 : * Note on the word "normal": \f$s_i\f$ points in the correct direction
336 : * (it is "normal" to the surface), but it does not have unit length
337 : * (it is not "normalized"; normalization requires a metric).
338 : *
339 : * \param dx_radius The Cartesian derivatives of the radius, as
340 : * returned by ylm::cartesian_derivs_of_scalar with
341 : * `ylm::radius` passed in as the scalar.
342 : * \param r_hat The radial unit vector as returned by
343 : * `ylm::rhat`.
344 : */
345 : template <typename Fr>
346 1 : tnsr::i<DataVector, 3, Fr> normal_one_form(
347 : const tnsr::i<DataVector, 3, Fr>& dx_radius,
348 : const tnsr::i<DataVector, 3, Fr>& r_hat);
349 :
350 : /*!
351 : * \param one_form The returned normal one form.
352 : * \param dx_radius The Cartesian derivatives of the radius, as
353 : * returned by ylm::cartesian_derivs_of_scalar with
354 : * `ylm::radius` passed in as the scalar.
355 : * \param r_hat The radial unit vector as returned by
356 : * `ylm::rhat`.
357 : */
358 : template <typename Fr>
359 1 : void normal_one_form(const gsl::not_null<tnsr::i<DataVector, 3, Fr>*> one_form,
360 : const tnsr::i<DataVector, 3, Fr>& dx_radius,
361 : const tnsr::i<DataVector, 3, Fr>& r_hat);
362 :
363 : /// @{
364 : /*!
365 : * The linear least squares fit of the polynomial of order 3
366 : * given a `std::vector` of `Strahlkorper`s to their \f$Y_l^m\f$ coefficients.
367 : * Assumes the the \f$l_{\max}\f$ and \f$m_{\max}\f$ of each `Strahlkorper`
368 : * are the same, and the returned vector consists of \f$2l_{\max}m_{\max}\f$
369 : * (the number of \f$Y_l^m\f$ coefficients) `std::array<double, 4>`s, each of
370 : * which consists of the four coefficients that define the best fit cubic to
371 : * each \f$Y_l^m\f$ coefficient of the `Strahlkorper` as a function of time.
372 : *
373 : * \param times The time corresponding to each `Strahlkorper` to be fit to.
374 : * \param strahlkorpers The `Strahlkorper` surfaces which consists of a set
375 : * of \f$Y_l^m\f$ coefficients corresponding to the shape of the `Strahlkorper`
376 : * at a particular time.
377 : */
378 : template <typename Fr>
379 1 : std::vector<std::array<double, 4>> fit_ylm_coeffs(
380 : const DataVector& times,
381 : const std::vector<Strahlkorper<Fr>>& strahlkorpers);
382 :
383 : /*!
384 : * \brief Compute the time derivative of a Strahlkorper from a number of
385 : * previous Strahlkorpers
386 : *
387 : * \details Does simple 1D FD with non-uniform spacing using
388 : * `fd::non_uniform_1d_weights`.
389 : * \param time_deriv Strahlkorper whose coefficients are the time derivative of
390 : * `previous_strahlkorpers`' coefficients.
391 : * \param previous_strahlkorpers All previous Strahlkorpers and the times they
392 : * are at. They are expected to have the most recent Strahlkorper in the front
393 : * and the Strahlkorper furthest in the past in the back of the deque.
394 : */
395 : template <typename Frame>
396 1 : void time_deriv_of_strahlkorper(
397 : gsl::not_null<Strahlkorper<Frame>*> time_deriv,
398 : const std::deque<std::pair<double, Strahlkorper<Frame>>>&
399 : previous_strahlkorpers);
400 :
401 : /// @{
402 : /*!
403 : * \brief Compute the power monitor of a Strahlkorper
404 : *
405 : * \details Computes a Strahlkorper represented as a spectral expansion
406 : * in scalar spherical harmonics $Y_{lm}(\theta,\phi)$ up to
407 : * $l=l_{\rm max}$ and $m=m_{\rm max}$. Then this function
408 : * computes the power in each $l$. If
409 : * a Strahlkorper's radius $r(\theta,\phi)$ is expanded as
410 : * \begin{equation}
411 : * r(\theta,\phi) = \sum_{l=0}^{l_{\rm max}}
412 : * \sum_{m=-\mathcal{M}}^{\mathcal{M}} C_{lm} Y_{lm}(\theta,\phi),
413 : * \end{equation}
414 : * where $\mathcal{M}$ is $l$ or $m_{\rm max}$, whichever is smaller,
415 : * then the power monitor $P_l$ is the square root of a normalized sum of
416 : * the magnitudes squared of the coefficients $C_{lm}$:
417 : * \begin{equation}
418 : * P_l = \sqrt{\frac{1}{2 \mathcal{M}+1}\sum_{m=-\mathcal{M}}^\mathcal{M}
419 : * \left|C_{lm}\right|^2}.
420 : * \end{equation}
421 : * See Eq. (51) and the surrounding discussion in \cite Szilagyi2014fna
422 : * for the motivation of this power monitor, which is useful in adaptively
423 : * selecting the angular resolution $l_{\rm max}$ when finding
424 : * apparent horizons. But note that because internally the Strahlkorper is
425 : * represented as a Spherepack expansion, the implementation of this function
426 : * uses that expansion rather than a literal implementation of the
427 : * above equations. As a result, it's the Spherepack spectral coefficients
428 : * that get squared and (for each $l$) summed over.
429 : * \param strahlkorper The Strahlkorper surface.
430 : */
431 : template <typename Frame>
432 1 : DataVector power_monitor(const Strahlkorper<Frame>& strahlkorper);
433 :
434 : /*!
435 : * \copydoc power_monitor
436 : * \param result A DataVector containing the power monitor $P_l$ for each
437 : * $l$ up to the Strahlkorper's `l_max`.
438 : */
439 : template <typename Frame>
440 1 : void power_monitor(gsl::not_null<DataVector*> result,
441 : const Strahlkorper<Frame>& strahlkorper);
442 : /// @}
443 : } // namespace ylm
|
