Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <algorithm>
7 : #include <array>
8 : #include <cmath>
9 : #include <cstddef>
10 : #include <iterator>
11 : #include <type_traits>
12 :
13 : #include "DataStructures/TaggedTuple.hpp"
14 : #include "DataStructures/Tensor/EagerMath/Determinant.hpp"
15 : #include "DataStructures/Tensor/EagerMath/FrameTransform.hpp"
16 : #include "DataStructures/Tensor/Tensor.hpp"
17 : #include "Options/Options.hpp"
18 : #include "PointwiseFunctions/AnalyticData/AnalyticData.hpp"
19 : #include "PointwiseFunctions/AnalyticData/GrMhd/AnalyticData.hpp"
20 : #include "PointwiseFunctions/AnalyticData/GrMhd/MagnetizedFmDisk.hpp"
21 : #include "PointwiseFunctions/AnalyticData/GrMhd/SphericalTorus.hpp"
22 : #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/Solutions.hpp"
23 : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
24 : #include "Utilities/Gsl.hpp"
25 : #include "Utilities/Serialization/CharmPupable.hpp"
26 : #include "Utilities/TMPL.hpp"
27 :
28 : /// \cond
29 : namespace PUP {
30 : class er;
31 : } // namespace PUP
32 : /// \endcond
33 :
34 : namespace grmhd::AnalyticData {
35 : /*!
36 : * \brief Magnetized fluid disk orbiting a Kerr black hole in the Kerr-Schild
37 : * Cartesian coordinates, but in a toroidal mesh defined from a torus map
38 : * (see grmhd::AnalyticData::SphericalTorus).
39 : */
40 1 : class PolarMagnetizedFmDisk
41 : : public virtual evolution::initial_data::InitialData,
42 : public MarkAsAnalyticData,
43 : public AnalyticDataBase {
44 : public:
45 0 : struct DiskParameters {
46 0 : using type = MagnetizedFmDisk;
47 0 : static constexpr Options::String help = "Parameters for the disk.";
48 : };
49 :
50 0 : struct TorusParameters {
51 0 : using type = grmhd::AnalyticData::SphericalTorus;
52 0 : static constexpr Options::String help =
53 : "Parameters for the evolution region.";
54 : };
55 :
56 0 : using options = tmpl::list<DiskParameters, TorusParameters>;
57 :
58 0 : static constexpr Options::String help =
59 : "Magnetized Fishbone-Moncrief disk in polar coordinates.";
60 :
61 0 : PolarMagnetizedFmDisk() = default;
62 0 : PolarMagnetizedFmDisk(const PolarMagnetizedFmDisk& /*rhs*/) = default;
63 0 : PolarMagnetizedFmDisk& operator=(const PolarMagnetizedFmDisk& /*rhs*/) =
64 : default;
65 0 : PolarMagnetizedFmDisk(PolarMagnetizedFmDisk&& /*rhs*/) = default;
66 0 : PolarMagnetizedFmDisk& operator=(PolarMagnetizedFmDisk&& /*rhs*/) = default;
67 0 : ~PolarMagnetizedFmDisk() override = default;
68 :
69 0 : PolarMagnetizedFmDisk(MagnetizedFmDisk fm_disk,
70 : grmhd::AnalyticData::SphericalTorus torus_map);
71 :
72 0 : auto get_clone() const
73 : -> std::unique_ptr<evolution::initial_data::InitialData> override;
74 :
75 : /// \cond
76 : explicit PolarMagnetizedFmDisk(CkMigrateMessage* msg);
77 : using PUP::able::register_constructor;
78 : WRAPPED_PUPable_decl_template(PolarMagnetizedFmDisk);
79 : /// \endcond
80 :
81 : /// The grmhd variables.
82 : ///
83 : /// \note The functions are optimized for retrieving the hydro variables
84 : /// before the metric variables.
85 : template <typename DataType, typename... Tags>
86 1 : tuples::TaggedTuple<Tags...> variables(const tnsr::I<DataType, 3>& x,
87 : tmpl::list<Tags...> /*meta*/) const {
88 : // In this function, we label the coordinates this solution works
89 : // in with Frame::BlockLogical, and the coordinates the wrapped
90 : // solution uses Frame::Inertial. This means the input and output
91 : // have to be converted to the correct label.
92 :
93 : const tnsr::I<DataType, 3> observation_coordinates(torus_map_(x));
94 :
95 : using dependencies = tmpl::map<
96 : tmpl::pair<gr::AnalyticSolution<3>::DerivShift<DataType>,
97 : gr::Tags::Shift<DataType, 3, Frame::Inertial>>,
98 : tmpl::pair<gr::AnalyticSolution<3>::DerivSpatialMetric<DataType>,
99 : gr::Tags::SpatialMetric<DataType, 3, Frame::Inertial>>>;
100 : using required_tags = tmpl::remove_duplicates<
101 : tmpl::remove<tmpl::list<Tags..., tmpl::at<dependencies, Tags>...>,
102 : tmpl::no_such_type_>>;
103 :
104 : auto observation_data =
105 : fm_disk_.variables(observation_coordinates, required_tags{});
106 :
107 : const auto jacobian = torus_map_.jacobian(x);
108 : const auto inv_jacobian = torus_map_.inv_jacobian(x);
109 :
110 : const auto change_frame = [this, &inv_jacobian, &jacobian, &x](
111 : const auto& data, auto tag) {
112 : using Tag = decltype(tag);
113 : auto result =
114 : transform::to_different_frame(get<Tag>(data), jacobian, inv_jacobian);
115 :
116 : if constexpr (std::is_same_v<
117 : Tag, gr::AnalyticSolution<3>::DerivShift<DataType>>) {
118 : const auto deriv_inv_jacobian =
119 : torus_map_.derivative_of_inv_jacobian(x);
120 : const auto& shift =
121 : get<gr::Tags::Shift<DataType, 3, Frame::Inertial>>(data);
122 : for (size_t i = 0; i < 3; ++i) {
123 : for (size_t j = 0; j < 3; ++j) {
124 : for (size_t k = 0; k < 3; ++k) {
125 : result.get(i, j) +=
126 : deriv_inv_jacobian.get(j, k, i) * shift.get(k);
127 : }
128 : }
129 : }
130 : } else if constexpr (std::is_same_v<
131 : Tag, gr::AnalyticSolution<3>::DerivSpatialMetric<
132 : DataType>>) {
133 : const auto hessian = torus_map_.hessian(x);
134 : const auto& spatial_metric =
135 : get<gr::Tags::SpatialMetric<DataType, 3, Frame::Inertial>>(data);
136 : for (size_t i = 0; i < 3; ++i) {
137 : for (size_t j = 0; j < 3; ++j) {
138 : for (size_t k = j; k < 3; ++k) {
139 : for (size_t l = 0; l < 3; ++l) {
140 : for (size_t m = 0; m < 3; ++m) {
141 : result.get(i, j, k) +=
142 : (hessian.get(l, j, i) * jacobian.get(m, k) +
143 : hessian.get(l, k, i) * jacobian.get(m, j)) *
144 : spatial_metric.get(l, m);
145 : }
146 : }
147 : }
148 : }
149 : }
150 : } else if constexpr (std::is_same_v<
151 : Tag, gr::Tags::SqrtDetSpatialMetric<DataType>>) {
152 : get(result) *= abs(get(determinant(jacobian)));
153 : }
154 :
155 : typename Tag::type result_with_replaced_frame{};
156 : std::copy(std::move_iterator(result.begin()),
157 : std::move_iterator(result.end()),
158 : result_with_replaced_frame.begin());
159 : return result_with_replaced_frame;
160 : };
161 :
162 : return {change_frame(observation_data, Tags{})...};
163 : }
164 :
165 0 : using equation_of_state_type = MagnetizedFmDisk::equation_of_state_type;
166 0 : const equation_of_state_type& equation_of_state() const {
167 : return fm_disk_.equation_of_state();
168 : }
169 :
170 : // NOLINTNEXTLINE(google-runtime-references)
171 0 : void pup(PUP::er& p) override;
172 :
173 : private:
174 0 : friend bool operator==(const PolarMagnetizedFmDisk& lhs,
175 : const PolarMagnetizedFmDisk& rhs);
176 :
177 0 : MagnetizedFmDisk fm_disk_;
178 0 : grmhd::AnalyticData::SphericalTorus torus_map_;
179 : };
180 :
181 0 : bool operator!=(const PolarMagnetizedFmDisk& lhs,
182 : const PolarMagnetizedFmDisk& rhs);
183 : } // namespace grmhd::AnalyticData
|
