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 <limits>
8 :
9 : #include "DataStructures/Tensor/TypeAliases.hpp"
10 : #include "Options/Context.hpp"
11 : #include "Options/String.hpp"
12 : #include "PointwiseFunctions/AnalyticData/AnalyticData.hpp"
13 : #include "PointwiseFunctions/AnalyticData/GrMhd/AnalyticData.hpp"
14 : #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/Minkowski.hpp"
15 : #include "PointwiseFunctions/Hydro/EquationsOfState/IdealFluid.hpp"
16 : #include "PointwiseFunctions/Hydro/TagsDeclarations.hpp"
17 : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
18 : #include "Utilities/Serialization/CharmPupable.hpp"
19 : #include "Utilities/TMPL.hpp"
20 : #include "Utilities/TaggedTuple.hpp"
21 :
22 : /// \cond
23 : namespace PUP {
24 : class er;
25 : } // namespace PUP
26 : /// \endcond
27 :
28 : namespace grmhd::AnalyticData {
29 :
30 : /*!
31 : * \brief Analytic initial data for an advecting magnetic field loop.
32 : *
33 : * This test, originally proposed in \cite Gardiner2005hy and presented in a
34 : * slightly modified form by \cite Mignone2010br, a region with annular
35 : * cross section with the specified `InnerRadius` and `OuterRadius` is given a
36 : * non-zero azimuthal magnetic field of constant magnitude `MagFieldStrength`
37 : * with zero magnetic field outside the loop. Inside the `InnerRadius` the
38 : * magnetic field strength falls to zero quadratically. The loop is embedded in
39 : * an ideal fluid with the given `AdiabaticIndex`, `RestMassDensity` and
40 : * `Pressure` with a uniform `AdvectionVelocity`. The magnetic field loop
41 : * should advect across the grid, maintaining its shape and strength, as long
42 : * as the magnetic pressure is negligible compared to the thermal pressure.
43 : *
44 : * This test diagnoses how well the evolution scheme preserves the no-monopole
45 : * condition, as well as the diffusivity of the scheme.
46 : *
47 : * The standard test setup is done on \f$x \in [-1,1]\f$, \f$y \in [-0.5,
48 : * 0.5]\f$, with periodic boundary conditions and with the following values
49 : * given for the options:
50 : * - InnerRadius: 0.06
51 : * - OuterRadius: 0.3
52 : * - RestMassDensity: 1.0
53 : * - Pressure: 1.0
54 : * - AdvectionVelocity: [0.08164965809277261, 0.040824829046386304, 0.0]
55 : * - MagFieldStrength: 0.001
56 : * - AdiabaticIndex: 1.66666666666666667
57 : *
58 : */
59 1 : class MagneticFieldLoop
60 : : public evolution::initial_data::InitialData,
61 : public MarkAsAnalyticData,
62 : public AnalyticDataBase,
63 : public hydro::TemperatureInitialization<MagneticFieldLoop> {
64 : public:
65 0 : using equation_of_state_type = EquationsOfState::IdealFluid<true>;
66 :
67 : /// The pressure throughout the fluid.
68 1 : struct Pressure {
69 0 : using type = double;
70 0 : static constexpr Options::String help = {
71 : "The constant pressure throughout the fluid."};
72 0 : static type lower_bound() { return 0.0; }
73 : };
74 :
75 : /// The rest mass density throughout the fluid.
76 1 : struct RestMassDensity {
77 0 : using type = double;
78 0 : static constexpr Options::String help = {
79 : "The constant density throughout the fluid."};
80 0 : static type lower_bound() { return 0.0; }
81 : };
82 :
83 : /// The adiabatic index for the ideal fluid.
84 1 : struct AdiabaticIndex {
85 0 : using type = double;
86 0 : static constexpr Options::String help = {
87 : "The adiabatic index for the ideal fluid."};
88 0 : static type lower_bound() { return 1.0; }
89 : };
90 :
91 : /// The fluid velocity.
92 1 : struct AdvectionVelocity {
93 0 : using type = std::array<double, 3>;
94 0 : static constexpr Options::String help = {"The advection velocity."};
95 0 : static type lower_bound() { return {{-1.0, -1.0, -1.0}}; }
96 0 : static type upper_bound() { return {{1.0, 1.0, 1.0}}; }
97 : };
98 :
99 : /// The strength of the magnetic field.
100 1 : struct MagFieldStrength {
101 0 : using type = double;
102 0 : static constexpr Options::String help = {
103 : "The magnitude of the magnetic field."};
104 0 : static type lower_bound() { return 0.0; }
105 : };
106 :
107 : /// The inner radius of the magnetic loop.
108 1 : struct InnerRadius {
109 0 : using type = double;
110 0 : static constexpr Options::String help = {
111 : "The inner radius of the magnetic loop."};
112 0 : static type lower_bound() { return 0.0; }
113 : };
114 :
115 : /// The outer radius of the magnetic loop.
116 1 : struct OuterRadius {
117 0 : using type = double;
118 0 : static constexpr Options::String help = {
119 : "The outer radius of the magnetic loop."};
120 0 : static type lower_bound() { return 0.0; }
121 : };
122 :
123 0 : using options =
124 : tmpl::list<Pressure, RestMassDensity, AdiabaticIndex, AdvectionVelocity,
125 : MagFieldStrength, InnerRadius, OuterRadius>;
126 0 : static constexpr Options::String help = {
127 : "Periodic advection of a magnetic field loop in Minkowski."};
128 :
129 0 : MagneticFieldLoop() = default;
130 0 : MagneticFieldLoop(const MagneticFieldLoop& /*rhs*/) = default;
131 0 : MagneticFieldLoop& operator=(const MagneticFieldLoop& /*rhs*/) = default;
132 0 : MagneticFieldLoop(MagneticFieldLoop&& /*rhs*/) = default;
133 0 : MagneticFieldLoop& operator=(MagneticFieldLoop&& /*rhs*/) = default;
134 0 : ~MagneticFieldLoop() override = default;
135 :
136 0 : MagneticFieldLoop(double pressure, double rest_mass_density,
137 : double adiabatic_index,
138 : const std::array<double, 3>& advection_velocity,
139 : double magnetic_field_magnitude, double inner_radius,
140 : double outer_radius, const Options::Context& context = {});
141 :
142 0 : auto get_clone() const
143 : -> std::unique_ptr<evolution::initial_data::InitialData> override;
144 :
145 : /// \cond
146 : explicit MagneticFieldLoop(CkMigrateMessage* msg);
147 : using PUP::able::register_constructor;
148 : WRAPPED_PUPable_decl_template(MagneticFieldLoop);
149 : /// \endcond
150 :
151 : /// @{
152 : /// Retrieve the GRMHD variables at a given position.
153 : template <typename DataType>
154 1 : auto variables(const tnsr::I<DataType, 3>& x,
155 : tmpl::list<hydro::Tags::RestMassDensity<DataType>> /*meta*/)
156 : const -> tuples::TaggedTuple<hydro::Tags::RestMassDensity<DataType>>;
157 :
158 : template <typename DataType>
159 1 : auto variables(const tnsr::I<DataType, 3>& x,
160 : tmpl::list<hydro::Tags::ElectronFraction<DataType>> /*meta*/)
161 : const -> tuples::TaggedTuple<hydro::Tags::ElectronFraction<DataType>>;
162 :
163 : template <typename DataType>
164 1 : auto variables(
165 : const tnsr::I<DataType, 3>& x,
166 : tmpl::list<hydro::Tags::SpecificInternalEnergy<DataType>> /*meta*/) const
167 : -> tuples::TaggedTuple<hydro::Tags::SpecificInternalEnergy<DataType>>;
168 :
169 : template <typename DataType>
170 1 : auto variables(const tnsr::I<DataType, 3>& x,
171 : tmpl::list<hydro::Tags::Pressure<DataType>> /*meta*/) const
172 : -> tuples::TaggedTuple<hydro::Tags::Pressure<DataType>>;
173 :
174 : template <typename DataType>
175 1 : auto variables(const tnsr::I<DataType, 3>& x,
176 : tmpl::list<hydro::Tags::SpatialVelocity<DataType, 3>> /*meta*/)
177 : const -> tuples::TaggedTuple<hydro::Tags::SpatialVelocity<DataType, 3>>;
178 :
179 : template <typename DataType>
180 1 : auto variables(const tnsr::I<DataType, 3>& x,
181 : tmpl::list<hydro::Tags::MagneticField<DataType, 3>> /*meta*/)
182 : const -> tuples::TaggedTuple<hydro::Tags::MagneticField<DataType, 3>>;
183 :
184 : template <typename DataType>
185 1 : auto variables(
186 : const tnsr::I<DataType, 3>& x,
187 : tmpl::list<hydro::Tags::DivergenceCleaningField<DataType>> /*meta*/) const
188 : -> tuples::TaggedTuple<hydro::Tags::DivergenceCleaningField<DataType>>;
189 :
190 : template <typename DataType>
191 1 : auto variables(const tnsr::I<DataType, 3>& x,
192 : tmpl::list<hydro::Tags::LorentzFactor<DataType>> /*meta*/)
193 : const -> tuples::TaggedTuple<hydro::Tags::LorentzFactor<DataType>>;
194 :
195 : template <typename DataType>
196 1 : auto variables(const tnsr::I<DataType, 3>& x,
197 : tmpl::list<hydro::Tags::SpecificEnthalpy<DataType>> /*meta*/)
198 : const -> tuples::TaggedTuple<hydro::Tags::SpecificEnthalpy<DataType>>;
199 :
200 : template <typename DataType>
201 1 : auto variables(const tnsr::I<DataType, 3>& x,
202 : tmpl::list<hydro::Tags::Temperature<DataType>> /*meta*/) const
203 : -> tuples::TaggedTuple<hydro::Tags::Temperature<DataType>> {
204 : return TemperatureInitialization::variables(
205 : x, tmpl::list<hydro::Tags::Temperature<DataType>>{});
206 : }
207 : /// @}
208 :
209 : /// Retrieve a collection of hydrodynamic variables at position x
210 : template <typename DataType, typename Tag1, typename Tag2, typename... Tags>
211 1 : tuples::TaggedTuple<Tag1, Tag2, Tags...> variables(
212 : const tnsr::I<DataType, 3>& x,
213 : tmpl::list<Tag1, Tag2, Tags...> /*meta*/) const {
214 : return {tuples::get<Tag1>(variables(x, tmpl::list<Tag1>{})),
215 : tuples::get<Tag2>(variables(x, tmpl::list<Tag2>{})),
216 : tuples::get<Tags>(variables(x, tmpl::list<Tags>{}))...};
217 : }
218 :
219 : /// Retrieve the metric variables
220 : template <typename DataType, typename Tag,
221 : Requires<tmpl::list_contains_v<
222 : gr::analytic_solution_tags<3, DataType>, Tag>> = nullptr>
223 1 : tuples::TaggedTuple<Tag> variables(const tnsr::I<DataType, 3>& x,
224 : tmpl::list<Tag> /*meta*/) const {
225 : constexpr double dummy_time = 0.0;
226 : return background_spacetime_.variables(x, dummy_time, tmpl::list<Tag>{});
227 : }
228 :
229 0 : const EquationsOfState::IdealFluid<true>& equation_of_state() const {
230 : return equation_of_state_;
231 : }
232 :
233 : // NOLINTNEXTLINE(google-runtime-references)
234 0 : void pup(PUP::er& /*p*/) override;
235 :
236 : private:
237 0 : double pressure_ = std::numeric_limits<double>::signaling_NaN();
238 0 : double rest_mass_density_ = std::numeric_limits<double>::signaling_NaN();
239 0 : double adiabatic_index_ = std::numeric_limits<double>::signaling_NaN();
240 0 : std::array<double, 3> advection_velocity_{
241 : {std::numeric_limits<double>::signaling_NaN(),
242 : std::numeric_limits<double>::signaling_NaN(),
243 : std::numeric_limits<double>::signaling_NaN()}};
244 0 : double magnetic_field_magnitude_ =
245 : std::numeric_limits<double>::signaling_NaN();
246 0 : double inner_radius_ = std::numeric_limits<double>::signaling_NaN();
247 0 : double outer_radius_ = std::numeric_limits<double>::signaling_NaN();
248 :
249 0 : EquationsOfState::IdealFluid<true> equation_of_state_{};
250 0 : gr::Solutions::Minkowski<3> background_spacetime_{};
251 :
252 0 : friend bool operator==(const MagneticFieldLoop& lhs,
253 : const MagneticFieldLoop& rhs);
254 :
255 0 : friend bool operator!=(const MagneticFieldLoop& lhs,
256 : const MagneticFieldLoop& rhs);
257 : };
258 :
259 : } // namespace grmhd::AnalyticData
|
