Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 : #include <limits>
8 :
9 : #include "DataStructures/TaggedTuple.hpp"
10 : #include "DataStructures/Tensor/TypeAliases.hpp"
11 : #include "Options/String.hpp"
12 : #include "PointwiseFunctions/AnalyticData/AnalyticData.hpp"
13 : #include "PointwiseFunctions/AnalyticData/GrMhd/AnalyticData.hpp"
14 : #include "PointwiseFunctions/AnalyticSolutions/RelativisticEuler/FishboneMoncriefDisk.hpp"
15 : #include "PointwiseFunctions/GeneralRelativity/KerrSchildCoords.hpp"
16 : #include "PointwiseFunctions/Hydro/EquationsOfState/PolytropicFluid.hpp"
17 : #include "PointwiseFunctions/Hydro/TagsDeclarations.hpp"
18 : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
19 : #include "Utilities/Serialization/CharmPupable.hpp"
20 : #include "Utilities/TMPL.hpp"
21 :
22 : /// \cond
23 : namespace PUP {
24 : class er;
25 : } // namespace PUP
26 : /// \endcond
27 :
28 : namespace grmhd::AnalyticData {
29 :
30 : /*!
31 : * \brief Magnetized fluid disk orbiting a Kerr black hole.
32 : *
33 : * In the context of simulating accretion disks, this class implements a widely
34 : * used (e.g. \cite Gammie2003, \cite Porth2016rfi, \cite White2015omx)
35 : * initial setup for the GRMHD variables, consisting of a Fishbone-Moncrief disk
36 : * \cite Fishbone1976apj (see also
37 : * RelativisticEuler::Solutions::FishboneMoncriefDisk),
38 : * threaded by a weak poloidal magnetic field. The magnetic field is constructed
39 : * from an axially symmetric toroidal magnetic potential which, in Kerr
40 : * ("spherical Kerr-Schild") coordinates, has the form
41 : *
42 : * \f{align*}
43 : * A_\phi(r,\theta) \propto \text{max}(\rho(r,\theta) - \rho_\text{thresh}, 0),
44 : * \f}
45 : *
46 : * where \f$\rho_\text{thresh}\f$ is a user-specified threshold density that
47 : * confines the magnetic flux to exist inside of the fluid disk only. A commonly
48 : * used value for this parameter is
49 : * \f$\rho_\text{thresh} = 0.2\rho_\text{max}\f$, where \f$\rho_\text{max}\f$
50 : * is the maximum value of
51 : * the rest mass density in the disk. Using this potential, the Eulerian
52 : * magnetic field takes the form
53 : *
54 : * \f{align*}
55 : * B^r = \frac{F_{\theta\phi}}{\sqrt{\gamma}},\quad
56 : * B^\theta = \frac{F_{\phi r}}{\sqrt{\gamma}},\quad B^\phi = 0,
57 : * \f}
58 : *
59 : * where \f$F_{ij} = \partial_i A_j - \partial_j A_i\f$ are the spatial
60 : * components of the Faraday tensor, and \f$\gamma\f$ is the determinant of the
61 : * spatial metric. The magnetic field is then normalized so that the
62 : * plasma-\f$\beta\f$ parameter, \f$\beta = 2p/b^2\f$, equals some value
63 : * specified by the user. Here, \f$p\f$ is the fluid pressure, and
64 : *
65 : * \f{align*}
66 : * b^2 = b^\mu b_\mu = \frac{B_iB^i}{W^2} + (B^iv_i)^2
67 : * \f}
68 : *
69 : * is the norm of the magnetic field in the fluid frame, with \f$v_i\f$ being
70 : * the spatial velocity, and \f$W\f$ the Lorentz factor.
71 : *
72 : * \note When using Kerr-Schild coordinates, the horizon that is at
73 : * constant \f$r\f$ is not spherical, but instead spheroidal. This could make
74 : * application of boundary condition and computing various fluxes
75 : * across the horizon more complicated than they need to be.
76 : * Thus, similar to RelativisticEuler::Solutions::FishboneMoncriefDisk
77 : * we use Spherical Kerr-Schild coordinates,
78 : * see gr::Solutions::SphericalKerrSchild, in which constant \f$r\f$
79 : * is spherical. Because we compute variables in Kerr-Schild coordinates,
80 : * there is a necessary extra step of transforming them back to
81 : * Spherical Kerr-Schild coordinates.
82 : *
83 : * \warning Spherical Kerr-Schild coordinates and "spherical Kerr-Schild"
84 : * coordinates are not same.
85 : *
86 : */
87 1 : class MagnetizedFmDisk : public virtual evolution::initial_data::InitialData,
88 : public MarkAsAnalyticData {
89 : private:
90 0 : using FmDisk = RelativisticEuler::Solutions::FishboneMoncriefDisk;
91 :
92 : public:
93 : /// The rest mass density (in units of the maximum rest mass density in the
94 : /// disk) below which the matter in the disk is initially unmagetized.
95 1 : struct ThresholdDensity {
96 0 : using type = double;
97 0 : static constexpr Options::String help = {
98 : "Frac. rest mass density below which B-field vanishes."};
99 0 : static type lower_bound() { return 0.0; }
100 0 : static type upper_bound() { return 1.0; }
101 : };
102 : /// The maximum-magnetic-pressure-to-maximum-fluid-pressure ratio.
103 1 : struct InversePlasmaBeta {
104 0 : using type = double;
105 0 : static constexpr Options::String help = {
106 : "Ratio of max magnetic pressure to max fluid pressure."};
107 0 : static type lower_bound() { return 0.0; }
108 : };
109 : /// Grid resolution used in magnetic field normalization.
110 1 : struct BFieldNormGridRes {
111 0 : using type = size_t;
112 0 : static constexpr Options::String help = {
113 : "Grid Resolution for b-field normalization."};
114 0 : static type suggested_value() { return 255; }
115 0 : static type lower_bound() { return 4; }
116 : };
117 :
118 : // Unlike the other analytic data classes, we cannot get these from the
119 : // `AnalyticDataBase` because this case causes clang-tidy to believe that
120 : // there is an ambiguous inheritance problem
121 0 : static constexpr size_t volume_dim = 3_st;
122 :
123 : template <typename DataType>
124 0 : using tags =
125 : tmpl::push_back<typename gr::AnalyticSolution<3>::template tags<DataType>,
126 : hydro::Tags::RestMassDensity<DataType>,
127 : hydro::Tags::ElectronFraction<DataType>,
128 : hydro::Tags::SpecificInternalEnergy<DataType>,
129 : hydro::Tags::Temperature<DataType>,
130 : hydro::Tags::Pressure<DataType>,
131 : hydro::Tags::SpatialVelocity<DataType, 3>,
132 : hydro::Tags::MagneticField<DataType, 3>,
133 : hydro::Tags::DivergenceCleaningField<DataType>,
134 : hydro::Tags::LorentzFactor<DataType>,
135 : hydro::Tags::SpecificEnthalpy<DataType>>;
136 :
137 0 : using options = tmpl::push_back<FmDisk::options, ThresholdDensity,
138 : InversePlasmaBeta, BFieldNormGridRes>;
139 :
140 0 : static constexpr Options::String help = {
141 : "Magnetized Fishbone-Moncrief disk."};
142 :
143 0 : MagnetizedFmDisk() = default;
144 0 : MagnetizedFmDisk(const MagnetizedFmDisk& /*rhs*/) = default;
145 0 : MagnetizedFmDisk& operator=(const MagnetizedFmDisk& /*rhs*/) = default;
146 0 : MagnetizedFmDisk(MagnetizedFmDisk&& /*rhs*/) = default;
147 0 : MagnetizedFmDisk& operator=(MagnetizedFmDisk&& /*rhs*/) = default;
148 0 : ~MagnetizedFmDisk() override = default;
149 :
150 0 : MagnetizedFmDisk(
151 : double bh_mass, double bh_dimless_spin, double inner_edge_radius,
152 : double max_pressure_radius, double polytropic_constant,
153 : double polytropic_exponent, double noise, double threshold_density,
154 : double inverse_plasma_beta,
155 : size_t normalization_grid_res = BFieldNormGridRes::suggested_value());
156 :
157 0 : auto get_clone() const
158 : -> std::unique_ptr<evolution::initial_data::InitialData> override;
159 :
160 : /// \cond
161 : explicit MagnetizedFmDisk(CkMigrateMessage* msg);
162 : using PUP::able::register_constructor;
163 : WRAPPED_PUPable_decl_template(MagnetizedFmDisk);
164 : /// \endcond
165 :
166 : // Overload the variables function from the base class.
167 0 : using equation_of_state_type = typename FmDisk::equation_of_state_type;
168 :
169 : /// @{
170 : /// The variables in Cartesian Kerr-Schild coordinates at `(x, t)`.
171 : template <typename DataType, typename... Tags>
172 1 : tuples::TaggedTuple<Tags...> variables(const tnsr::I<DataType, 3>& x,
173 : tmpl::list<Tags...> /*meta*/) const {
174 : // Can't store IntermediateVariables as member variable because we
175 : // need to be threadsafe.
176 : typename FmDisk::IntermediateVariables<DataType> vars(x);
177 : return {std::move(get<Tags>([this, &x, &vars]() {
178 : if constexpr (std::is_same_v<hydro::Tags::MagneticField<DataType, 3>,
179 : Tags>) {
180 : return variables(x, tmpl::list<Tags>{}, make_not_null(&vars));
181 : } else {
182 : return fm_disk_.variables(x, tmpl::list<Tags>{}, make_not_null(&vars));
183 : }
184 : }()))...};
185 : }
186 :
187 : template <typename DataType, typename Tag>
188 1 : tuples::TaggedTuple<Tag> variables(const tnsr::I<DataType, 3>& x,
189 : tmpl::list<Tag> /*meta*/) const {
190 : // Can't store IntermediateVariables as member variable because we need to
191 : // be threadsafe.
192 : typename FmDisk::IntermediateVariables<DataType> vars(x);
193 : if constexpr (std::is_same_v<hydro::Tags::MagneticField<DataType, 3>,
194 : Tag>) {
195 : return variables(x, tmpl::list<Tag>{}, make_not_null(&vars));
196 : } else {
197 : return fm_disk_.variables(x, tmpl::list<Tag>{}, make_not_null(&vars));
198 : }
199 : }
200 : /// @}
201 :
202 : // NOLINTNEXTLINE(google-runtime-references)
203 0 : void pup(PUP::er& /*p*/) override;
204 :
205 0 : const equation_of_state_type& equation_of_state() const {
206 : return fm_disk_.equation_of_state();
207 : }
208 :
209 : private:
210 : template <typename DataType>
211 0 : auto variables(const tnsr::I<DataType, 3>& x,
212 : tmpl::list<hydro::Tags::MagneticField<DataType, 3>> /*meta*/,
213 : gsl::not_null<FmDisk::IntermediateVariables<DataType>*> vars)
214 : const -> tuples::TaggedTuple<hydro::Tags::MagneticField<DataType, 3>>;
215 :
216 : template <typename DataType>
217 0 : tnsr::I<DataType, 3> unnormalized_magnetic_field(
218 : const tnsr::I<DataType, 3>& x) const;
219 :
220 0 : friend bool operator==(const MagnetizedFmDisk& lhs,
221 : const MagnetizedFmDisk& rhs);
222 :
223 0 : RelativisticEuler::Solutions::FishboneMoncriefDisk fm_disk_{};
224 0 : double threshold_density_ = std::numeric_limits<double>::signaling_NaN();
225 0 : double inverse_plasma_beta_ = std::numeric_limits<double>::signaling_NaN();
226 0 : double b_field_normalization_ = std::numeric_limits<double>::signaling_NaN();
227 0 : size_t normalization_grid_res_ = 255;
228 0 : gr::KerrSchildCoords kerr_schild_coords_{};
229 : };
230 :
231 0 : bool operator!=(const MagnetizedFmDisk& lhs, const MagnetizedFmDisk& rhs);
232 :
233 : } // namespace grmhd::AnalyticData
|
