Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <array>
7 : #include <cmath>
8 : #include <cstddef>
9 : #include <limits>
10 : #include <optional>
11 :
12 : #include "DataStructures/Tensor/Tensor.hpp"
13 : #include "DataStructures/Tensor/TypeAliases.hpp"
14 : #include "Options/Context.hpp"
15 : #include "Options/String.hpp"
16 :
17 : /// \cond
18 : namespace PUP {
19 : class er;
20 : } // namespace PUP
21 : namespace gsl {
22 : template <typename T>
23 : class not_null;
24 : } // namespace gsl
25 : /// \endcond
26 :
27 : namespace grmhd::AnalyticData {
28 : /*!
29 : * \brief Torus made by removing two polar cones from a spherical shell
30 : *
31 : * Maps source coordinates \f$(\xi, \eta, \zeta)\f$ to
32 : * \f{align}
33 : * \vec{x}(\xi, \eta, \zeta) =
34 : * \begin{bmatrix}
35 : * r \sin\theta\cos\phi \\
36 : * r \sin\theta\sin\phi \\
37 : * r \cos\theta
38 : * \end{bmatrix}
39 : * \f}
40 : *
41 : * where
42 : * \f{align}
43 : * r & = r_\mathrm{min}\frac{1-\xi}{2} + r_\mathrm{max}\frac{1+\xi}{2}, \\
44 : * \eta_\mathrm{new} & = a_\mathrm{comp} \eta^3 + (1-a_\mathrm{comp}) \eta, \\
45 : * \theta & = \pi/2 - (\pi/2 - \theta_\mathrm{min}) \eta_\mathrm{new}, \\
46 : * \phi & = f_\mathrm{torus} \pi \zeta.
47 : * \f}
48 : *
49 : * - $r_\mathrm{min}$ and $r_\mathrm{max}$ are inner and outer radius of torus.
50 : * - $\theta_\mathrm{min}\in(0,\pi/2)$ is the minimum polar angle (measured
51 : * from +z axis) of torus, which is equal to the half of the apex angle of
52 : * the removed polar cones.
53 : * - $f_\mathrm{torus}\in(0, 1)$ is azimuthal fraction that the torus covers.
54 : * - $a_\mathrm{comp}\in[0,1)$ sets the level of equatorial compression
55 : * for theta, with zero being none.
56 : *
57 : */
58 1 : class SphericalTorus {
59 : public:
60 0 : static constexpr size_t dim = 3;
61 :
62 0 : struct RadialRange {
63 0 : using type = std::array<double, 2>;
64 0 : static constexpr Options::String help =
65 : "Radial extent of the torus, "
66 : "[min_radius, max_radius] ";
67 : };
68 :
69 0 : struct MinPolarAngle {
70 0 : using type = double;
71 0 : static constexpr Options::String help =
72 : "Half of the apex angle of excised polar cones. "
73 : "Polar angle (measured from +z axis) of torus has range "
74 : "[MinPolarAngle, pi - MinPolarAngle]";
75 0 : static type lower_bound() { return 0.0; }
76 0 : static type upper_bound() { return 0.5 * M_PI; }
77 : };
78 :
79 0 : struct FractionOfTorus {
80 0 : using type = double;
81 0 : static constexpr Options::String help =
82 : "Fraction of (azimuthal) orbit covered. Azimuthal angle has range "
83 : "[- pi * FractionOfTorus, pi * FractionOfTorus].";
84 0 : static type lower_bound() { return 0.0; }
85 0 : static type upper_bound() { return 1.0; }
86 : };
87 :
88 0 : struct CompressionLevel {
89 0 : using type = double;
90 0 : static constexpr Options::String help =
91 : "Level of Equatorial Compression for the polar angle.";
92 0 : static type lower_bound() { return 0.0; }
93 0 : static type upper_bound() { return 1.0; }
94 : };
95 :
96 0 : static constexpr Options::String help =
97 : "Torus made by removing polar cones from a spherical shell";
98 :
99 0 : using options =
100 : tmpl::list<RadialRange, MinPolarAngle, FractionOfTorus, CompressionLevel>;
101 :
102 0 : SphericalTorus(const std::array<double, 2>& radial_range,
103 : double min_polar_angle, double fraction_of_torus,
104 : double compression_level,
105 : const Options::Context& context = {});
106 :
107 0 : SphericalTorus(double r_min, double r_max, double min_polar_angle,
108 : double fraction_of_torus = 1.0, double compression_level = 0.0,
109 : const Options::Context& context = {});
110 :
111 0 : SphericalTorus() = default;
112 :
113 : template <typename T>
114 0 : tnsr::I<T, 3> operator()(const tnsr::I<T, 3>& source_coords) const;
115 :
116 0 : tnsr::I<double, 3> inverse(const tnsr::I<double, 3>& target_coords) const;
117 :
118 : template <typename T>
119 0 : Jacobian<T, 3, Frame::BlockLogical, Frame::Inertial> jacobian(
120 : const tnsr::I<T, 3>& source_coords) const;
121 :
122 : template <typename T>
123 0 : InverseJacobian<T, 3, Frame::BlockLogical, Frame::Inertial> inv_jacobian(
124 : const tnsr::I<T, 3>& source_coords) const;
125 :
126 : template <typename T>
127 0 : tnsr::Ijj<T, 3, Frame::NoFrame> hessian(
128 : const tnsr::I<T, 3>& source_coords) const;
129 :
130 : template <typename T>
131 0 : tnsr::Ijk<T, 3, Frame::NoFrame> derivative_of_inv_jacobian(
132 : const tnsr::I<T, 3>& source_coords) const;
133 :
134 : // NOLINTNEXTLINE(google-runtime-references)
135 0 : void pup(PUP::er& p);
136 :
137 0 : bool is_identity() const { return false; }
138 :
139 : private:
140 : template <typename T>
141 0 : T radius(const T& x) const;
142 :
143 : template <typename T>
144 0 : void radius(gsl::not_null<T*> r, const T& x) const;
145 :
146 : template <typename T>
147 0 : T radius_inverse(const T& r) const;
148 :
149 : template <typename T>
150 0 : T cubic_compression(const T& x) const;
151 :
152 0 : double cubic_inversion(double x) const;
153 :
154 0 : friend bool operator==(const SphericalTorus& lhs, const SphericalTorus& rhs);
155 :
156 0 : double r_min_ = std::numeric_limits<double>::signaling_NaN();
157 0 : double r_max_ = std::numeric_limits<double>::signaling_NaN();
158 0 : double pi_over_2_minus_theta_min_ =
159 : std::numeric_limits<double>::signaling_NaN();
160 0 : double fraction_of_torus_ = std::numeric_limits<double>::signaling_NaN();
161 0 : double compression_level_ = std::numeric_limits<double>::signaling_NaN();
162 : };
163 :
164 0 : bool operator!=(const SphericalTorus& lhs, const SphericalTorus& rhs);
165 : } // namespace grmhd::AnalyticData
|
