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