Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <limits>
7 : #include <memory>
8 : #include <optional>
9 :
10 : #include "DataStructures/DataBox/Prefixes.hpp"
11 : #include "DataStructures/Tensor/TypeAliases.hpp"
12 : #include "Evolution/BoundaryCorrection.hpp"
13 : #include "Evolution/Systems/GrMhd/ValenciaDivClean/Tags.hpp"
14 : #include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
15 : #include "Options/String.hpp"
16 : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
17 : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
18 : #include "PointwiseFunctions/Hydro/Tags.hpp"
19 : #include "Utilities/Gsl.hpp"
20 : #include "Utilities/Serialization/CharmPupable.hpp"
21 : #include "Utilities/TMPL.hpp"
22 :
23 : /// \cond
24 : class DataVector;
25 : namespace gsl {
26 : template <typename T>
27 : class not_null;
28 : } // namespace gsl
29 : namespace PUP {
30 : class er;
31 : } // namespace PUP
32 : /// \endcond
33 :
34 : namespace grmhd::ValenciaDivClean::BoundaryCorrections {
35 : /*!
36 : * \brief An HLL Riemann solver
37 : *
38 : * Let \f$U\f$ be the evolved variable, \f$F^i\f$ the flux, and \f$n_i\f$ be the
39 : * outward directed unit normal to the interface. Denoting \f$F := n_i F^i\f$,
40 : * the HLL boundary correction is \cite Harten1983
41 : *
42 : * \f{align*}
43 : * G_\text{HLL} = \frac{\lambda_\text{max} F_\text{int} +
44 : * \lambda_\text{min} F_\text{ext}}{\lambda_\text{max} - \lambda_\text{min}}
45 : * - \frac{\lambda_\text{min}\lambda_\text{max}}{\lambda_\text{max} -
46 : * \lambda_\text{min}} \left(U_\text{int} - U_\text{ext}\right)
47 : * \f}
48 : *
49 : * where "int" and "ext" stand for interior and exterior.
50 : * \f$\lambda_\text{min}\f$ and \f$\lambda_\text{max}\f$ are defined as
51 : *
52 : * \f{align*}
53 : * \lambda_\text{min} &=
54 : * \text{min}\left(\lambda^{-}_\text{int},-\lambda^{+}_\text{ext}, 0\right) \\
55 : * \lambda_\text{max} &=
56 : * \text{max}\left(\lambda^{+}_\text{int},-\lambda^{-}_\text{ext}, 0\right)
57 : * \f}
58 : *
59 : * where \f$\lambda^{+}\f$ (\f$\lambda^{-}\f$) is the largest characteristic
60 : * speed in the outgoing (ingoing) direction. Note the minus signs in front of
61 : * \f$\lambda^{\pm}_\text{ext}\f$, which is because an outgoing speed w.r.t. the
62 : * neighboring element is an ingoing speed w.r.t. the local element, and vice
63 : * versa. Similarly, the \f$F_{\text{ext}}\f$ term in \f$G_\text{HLL}\f$ has a
64 : * positive sign because the outward directed normal of the neighboring element
65 : * has the opposite sign, i.e. \f$n_i^{\text{ext}}=-n_i^{\text{int}}\f$.
66 : *
67 : * The characteristic/signal speeds are given in the documentation for
68 : * `grmhd::ValenciaDivClean::characteristic_speeds()`. Since the fluid is
69 : * travelling slower than the speed of light, the speeds we are interested in
70 : * are
71 : *
72 : * \f{align*}{
73 : * \lambda^{\pm}&=\pm\alpha-\beta^i n_i,
74 : * \f}
75 : *
76 : * which correspond to the divergence cleaning field.
77 : *
78 : * \note
79 : * - In the strong form the `dg_boundary_terms` function returns
80 : * \f$G - F_\text{int}\f$
81 : * - For either \f$\lambda_\text{min} = 0\f$ or \f$\lambda_\text{max} = 0\f$
82 : * (i.e. all characteristics move in the same direction) the HLL boundary
83 : * correction reduces to pure upwinding.
84 : * - Some references use \f$S\f$ instead of \f$\lambda\f$ for the
85 : * signal/characteristic speeds
86 : * - It may be possible to use the slower speeds for the magnetic field and
87 : * fluid part of the system in order to make the flux less dissipative for
88 : * those variables.
89 : */
90 1 : class Hll final : public evolution::BoundaryCorrection {
91 : public:
92 0 : struct LargestOutgoingCharSpeed : db::SimpleTag {
93 0 : using type = Scalar<DataVector>;
94 : };
95 0 : struct LargestIngoingCharSpeed : db::SimpleTag {
96 0 : using type = Scalar<DataVector>;
97 : };
98 :
99 0 : struct MagneticFieldMagnitudeForHydro {
100 0 : static constexpr Options::String help = {
101 : "When the magnetic field is below this value we use the hydro "
102 : "characteristic speeds."};
103 0 : using type = double;
104 : };
105 0 : struct LightSpeedDensityCutoff {
106 0 : static constexpr Options::String help = {
107 : "When the density is below this value we just use the light speed for "
108 : "the characteristic speeds."};
109 0 : using type = double;
110 : };
111 0 : using options =
112 : tmpl::list<MagneticFieldMagnitudeForHydro, LightSpeedDensityCutoff>;
113 0 : static constexpr Options::String help = {
114 : "Computes the HLL boundary correction term for the GRMHD system."};
115 :
116 0 : Hll() = default;
117 0 : Hll(const Hll&) = default;
118 0 : Hll& operator=(const Hll&) = default;
119 0 : Hll(Hll&&) = default;
120 0 : Hll& operator=(Hll&&) = default;
121 0 : ~Hll() override = default;
122 :
123 0 : Hll(double magnetic_field_magnitude_for_hydro,
124 : double light_speed_density_cutoff);
125 :
126 : /// \cond
127 : explicit Hll(CkMigrateMessage* /*unused*/);
128 : using PUP::able::register_constructor;
129 : WRAPPED_PUPable_decl_template(Hll); // NOLINT
130 : /// \endcond
131 0 : void pup(PUP::er& p) override; // NOLINT
132 :
133 0 : std::unique_ptr<BoundaryCorrection> get_clone() const override;
134 :
135 0 : using dg_package_field_tags =
136 : tmpl::list<Tags::TildeD, Tags::TildeYe, Tags::TildeTau,
137 : Tags::TildeS<Frame::Inertial>, Tags::TildeB<Frame::Inertial>,
138 : Tags::TildePhi, ::Tags::NormalDotFlux<Tags::TildeD>,
139 : ::Tags::NormalDotFlux<Tags::TildeYe>,
140 : ::Tags::NormalDotFlux<Tags::TildeTau>,
141 : ::Tags::NormalDotFlux<Tags::TildeS<Frame::Inertial>>,
142 : ::Tags::NormalDotFlux<Tags::TildeB<Frame::Inertial>>,
143 : ::Tags::NormalDotFlux<Tags::TildePhi>,
144 : LargestOutgoingCharSpeed, LargestIngoingCharSpeed>;
145 0 : using dg_package_data_temporary_tags = tmpl::list<
146 : gr::Tags::Lapse<DataVector>, gr::Tags::Shift<DataVector, 3>,
147 : hydro::Tags::SpatialVelocityOneForm<DataVector, 3, Frame::Inertial>>;
148 0 : using dg_package_data_primitive_tags =
149 : tmpl::list<hydro::Tags::RestMassDensity<DataVector>,
150 : hydro::Tags::ElectronFraction<DataVector>,
151 : hydro::Tags::Temperature<DataVector>,
152 : hydro::Tags::SpatialVelocity<DataVector, 3>>;
153 0 : using dg_package_data_volume_tags =
154 : tmpl::list<hydro::Tags::GrmhdEquationOfState>;
155 0 : using dg_boundary_terms_volume_tags = tmpl::list<>;
156 :
157 0 : double dg_package_data(
158 : gsl::not_null<Scalar<DataVector>*> packaged_tilde_d,
159 : gsl::not_null<Scalar<DataVector>*> packaged_tilde_ye,
160 : gsl::not_null<Scalar<DataVector>*> packaged_tilde_tau,
161 : gsl::not_null<tnsr::i<DataVector, 3, Frame::Inertial>*> packaged_tilde_s,
162 : gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*> packaged_tilde_b,
163 : gsl::not_null<Scalar<DataVector>*> packaged_tilde_phi,
164 : gsl::not_null<Scalar<DataVector>*> packaged_normal_dot_flux_tilde_d,
165 : gsl::not_null<Scalar<DataVector>*> packaged_normal_dot_flux_tilde_ye,
166 : gsl::not_null<Scalar<DataVector>*> packaged_normal_dot_flux_tilde_tau,
167 : gsl::not_null<tnsr::i<DataVector, 3, Frame::Inertial>*>
168 : packaged_normal_dot_flux_tilde_s,
169 : gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*>
170 : packaged_normal_dot_flux_tilde_b,
171 : gsl::not_null<Scalar<DataVector>*> packaged_normal_dot_flux_tilde_phi,
172 : gsl::not_null<Scalar<DataVector>*> packaged_largest_outgoing_char_speed,
173 : gsl::not_null<Scalar<DataVector>*> packaged_largest_ingoing_char_speed,
174 :
175 : const Scalar<DataVector>& tilde_d, const Scalar<DataVector>& tilde_ye,
176 : const Scalar<DataVector>& tilde_tau,
177 : const tnsr::i<DataVector, 3, Frame::Inertial>& tilde_s,
178 : const tnsr::I<DataVector, 3, Frame::Inertial>& tilde_b,
179 : const Scalar<DataVector>& tilde_phi,
180 :
181 : const tnsr::I<DataVector, 3, Frame::Inertial>& flux_tilde_d,
182 : const tnsr::I<DataVector, 3, Frame::Inertial>& flux_tilde_ye,
183 : const tnsr::I<DataVector, 3, Frame::Inertial>& flux_tilde_tau,
184 : const tnsr::Ij<DataVector, 3, Frame::Inertial>& flux_tilde_s,
185 : const tnsr::IJ<DataVector, 3, Frame::Inertial>& flux_tilde_b,
186 : const tnsr::I<DataVector, 3, Frame::Inertial>& flux_tilde_phi,
187 :
188 : const Scalar<DataVector>& lapse,
189 : const tnsr::I<DataVector, 3, Frame::Inertial>& shift,
190 : const tnsr::i<DataVector, 3, Frame::Inertial>& spatial_velocity_one_form,
191 :
192 : const Scalar<DataVector>& rest_mass_density,
193 : const Scalar<DataVector>& electron_fraction,
194 : const Scalar<DataVector>& temperature,
195 : const tnsr::I<DataVector, 3, Frame::Inertial>& spatial_velocity,
196 :
197 : const tnsr::i<DataVector, 3, Frame::Inertial>& normal_covector,
198 : const tnsr::I<DataVector, 3, Frame::Inertial>& normal_vector,
199 : const std::optional<tnsr::I<DataVector, 3, Frame::Inertial>>&
200 : /*mesh_velocity*/,
201 : const std::optional<Scalar<DataVector>>& normal_dot_mesh_velocity,
202 : const EquationsOfState::EquationOfState<true, 3>& equation_of_state)
203 : const;
204 :
205 0 : static void dg_boundary_terms(
206 : gsl::not_null<Scalar<DataVector>*> boundary_correction_tilde_d,
207 : gsl::not_null<Scalar<DataVector>*> boundary_correction_tilde_ye,
208 : gsl::not_null<Scalar<DataVector>*> boundary_correction_tilde_tau,
209 : gsl::not_null<tnsr::i<DataVector, 3, Frame::Inertial>*>
210 : boundary_correction_tilde_s,
211 : gsl::not_null<tnsr::I<DataVector, 3, Frame::Inertial>*>
212 : boundary_correction_tilde_b,
213 : gsl::not_null<Scalar<DataVector>*> boundary_correction_tilde_phi,
214 : const Scalar<DataVector>& tilde_d_int,
215 : const Scalar<DataVector>& tilde_ye_int,
216 : const Scalar<DataVector>& tilde_tau_int,
217 : const tnsr::i<DataVector, 3, Frame::Inertial>& tilde_s_int,
218 : const tnsr::I<DataVector, 3, Frame::Inertial>& tilde_b_int,
219 : const Scalar<DataVector>& tilde_phi_int,
220 : const Scalar<DataVector>& normal_dot_flux_tilde_d_int,
221 : const Scalar<DataVector>& normal_dot_flux_tilde_ye_int,
222 : const Scalar<DataVector>& normal_dot_flux_tilde_tau_int,
223 : const tnsr::i<DataVector, 3, Frame::Inertial>&
224 : normal_dot_flux_tilde_s_int,
225 : const tnsr::I<DataVector, 3, Frame::Inertial>&
226 : normal_dot_flux_tilde_b_int,
227 : const Scalar<DataVector>& normal_dot_flux_tilde_phi_int,
228 : const Scalar<DataVector>& largest_outgoing_char_speed_int,
229 : const Scalar<DataVector>& largest_ingoing_char_speed_int,
230 : const Scalar<DataVector>& tilde_d_ext,
231 : const Scalar<DataVector>& tilde_ye_ext,
232 : const Scalar<DataVector>& tilde_tau_ext,
233 : const tnsr::i<DataVector, 3, Frame::Inertial>& tilde_s_ext,
234 : const tnsr::I<DataVector, 3, Frame::Inertial>& tilde_b_ext,
235 : const Scalar<DataVector>& tilde_phi_ext,
236 : const Scalar<DataVector>& normal_dot_flux_tilde_d_ext,
237 : const Scalar<DataVector>& normal_dot_flux_tilde_ye_ext,
238 : const Scalar<DataVector>& normal_dot_flux_tilde_tau_ext,
239 : const tnsr::i<DataVector, 3, Frame::Inertial>&
240 : normal_dot_flux_tilde_s_ext,
241 : const tnsr::I<DataVector, 3, Frame::Inertial>&
242 : normal_dot_flux_tilde_b_ext,
243 : const Scalar<DataVector>& normal_dot_flux_tilde_phi_ext,
244 : const Scalar<DataVector>& largest_outgoing_char_speed_ext,
245 : const Scalar<DataVector>& largest_ingoing_char_speed_ext,
246 : dg::Formulation dg_formulation);
247 :
248 : private:
249 0 : friend bool operator==(const Hll& lhs, const Hll& rhs);
250 :
251 0 : double magnetic_field_magnitude_for_hydro_{
252 : std::numeric_limits<double>::signaling_NaN()};
253 0 : double light_speed_density_cutoff_{
254 : std::numeric_limits<double>::signaling_NaN()};
255 : };
256 0 : bool operator!=(const Hll& lhs, const Hll& rhs);
257 : } // namespace grmhd::ValenciaDivClean::BoundaryCorrections
|
