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/Systems/RadiationTransport/M1Grey/BoundaryCorrections/BoundaryCorrection.hpp"
12 : #include "Evolution/Systems/RadiationTransport/M1Grey/Tags.hpp"
13 : #include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
14 : #include "Options/String.hpp"
15 : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.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 RadiationTransport::M1Grey::BoundaryCorrections {
33 :
34 : namespace Rusanov_detail {
35 : void dg_package_data_impl(
36 : gsl::not_null<Scalar<DataVector>*> packaged_tilde_e,
37 : gsl::not_null<tnsr::i<DataVector, 3, Frame::Inertial>*> packaged_tilde_s,
38 : gsl::not_null<Scalar<DataVector>*> packaged_normal_dot_flux_tilde_e,
39 : gsl::not_null<tnsr::i<DataVector, 3, Frame::Inertial>*>
40 : packaged_normal_dot_flux_tilde_s,
41 : const Scalar<DataVector>& tilde_e,
42 : const tnsr::i<DataVector, 3, Frame::Inertial>& tilde_s,
43 : const tnsr::I<DataVector, 3, Frame::Inertial>& flux_tilde_e,
44 : const tnsr::Ij<DataVector, 3, Frame::Inertial>& flux_tilde_s,
45 : const tnsr::i<DataVector, 3, Frame::Inertial>& normal_covector,
46 : const tnsr::I<DataVector, 3, Frame::Inertial>& normal_vector,
47 : const std::optional<tnsr::I<DataVector, 3, Frame::Inertial>>& mesh_velocity,
48 : const std::optional<Scalar<DataVector>>& normal_dot_mesh_velocity);
49 :
50 : void dg_boundary_terms_impl(
51 : gsl::not_null<Scalar<DataVector>*> boundary_correction_tilde_e,
52 : gsl::not_null<tnsr::i<DataVector, 3, Frame::Inertial>*>
53 : boundary_correction_tilde_s,
54 : const Scalar<DataVector>& tilde_e_int,
55 : const tnsr::i<DataVector, 3, Frame::Inertial>& tilde_s_int,
56 : const Scalar<DataVector>& normal_dot_flux_tilde_e_int,
57 : const tnsr::i<DataVector, 3, Frame::Inertial>& normal_dot_flux_tilde_s_int,
58 : const Scalar<DataVector>& tilde_e_ext,
59 : const tnsr::i<DataVector, 3, Frame::Inertial>& tilde_s_ext,
60 : const Scalar<DataVector>& normal_dot_flux_tilde_e_ext,
61 : const tnsr::i<DataVector, 3, Frame::Inertial>& normal_dot_flux_tilde_s_ext,
62 : dg::Formulation dg_formulation);
63 : } // namespace Rusanov_detail
64 :
65 : /// \cond
66 : template <typename NeutrinoSpeciesList>
67 : class Rusanov;
68 : /// \endcond
69 :
70 : /*!
71 : * \brief A Rusanov/local Lax-Friedrichs Riemann solver
72 : *
73 : * Let \f$U\f$ be the state vector of evolved variables, \f$F^i\f$ the
74 : * corresponding fluxes, and \f$n_i\f$ be the outward directed unit normal to
75 : * the interface. Denoting \f$F := n_i F^i\f$, the %Rusanov boundary correction
76 : * is
77 : *
78 : * \f{align*}
79 : * G_\text{Rusanov} = \frac{F_\text{int} - F_\text{ext}}{2} -
80 : * \frac{\text{max}\left(\{|\lambda_\text{int}|\},
81 : * \{|\lambda_\text{ext}|\}\right)}{2} \left(U_\text{ext} - U_\text{int}\right),
82 : * \f}
83 : *
84 : * where "int" and "ext" stand for interior and exterior, and
85 : * \f$\{|\lambda|\}\f$ is the set of characteristic/signal speeds. The minus
86 : * sign in front of the \f$F_{\text{ext}}\f$ is necessary because the outward
87 : * directed normal of the neighboring element has the opposite sign, i.e.
88 : * \f$n_i^{\text{ext}}=-n_i^{\text{int}}\f$.
89 : *
90 : * For radiation, the max characteristic/signal speed is 1 (the speed of light).
91 : * Note that the characteristic speeds of this system are *not* yet fully
92 : * implemented, see the warning in the documentation for the characteristics.
93 : *
94 : * \note In the strong form the `dg_boundary_terms` function returns
95 : * \f$G - F_\text{int}\f$
96 : */
97 : template <typename... NeutrinoSpecies>
98 1 : class Rusanov<tmpl::list<NeutrinoSpecies...>> final
99 : : public BoundaryCorrection<tmpl::list<NeutrinoSpecies...>> {
100 : public:
101 0 : using options = tmpl::list<>;
102 0 : static constexpr Options::String help = {
103 : "Computes the Rusanov or local Lax-Friedrichs boundary correction term "
104 : "for the M1Grey radiation transport system."};
105 :
106 0 : Rusanov() = default;
107 0 : Rusanov(const Rusanov&) = default;
108 0 : Rusanov& operator=(const Rusanov&) = default;
109 0 : Rusanov(Rusanov&&) = default;
110 0 : Rusanov& operator=(Rusanov&&) = default;
111 0 : ~Rusanov() override = default;
112 :
113 : /// \cond
114 : explicit Rusanov(CkMigrateMessage* msg)
115 : : BoundaryCorrection<tmpl::list<NeutrinoSpecies...>>(msg) {}
116 : using PUP::able::register_constructor;
117 : WRAPPED_PUPable_decl_template(Rusanov); // NOLINT
118 : /// \endcond
119 0 : void pup(PUP::er& p) override { // NOLINT
120 : BoundaryCorrection<tmpl::list<NeutrinoSpecies...>>::pup(p);
121 : }
122 :
123 : std::unique_ptr<BoundaryCorrection<tmpl::list<NeutrinoSpecies...>>>
124 0 : get_clone() const override {
125 : return std::make_unique<Rusanov>(*this);
126 : }
127 :
128 0 : using dg_package_field_tags = tmpl::list<
129 : Tags::TildeE<Frame::Inertial, NeutrinoSpecies>...,
130 : Tags::TildeS<Frame::Inertial, NeutrinoSpecies>...,
131 : ::Tags::NormalDotFlux<Tags::TildeE<Frame::Inertial, NeutrinoSpecies>>...,
132 : ::Tags::NormalDotFlux<Tags::TildeS<Frame::Inertial, NeutrinoSpecies>>...>;
133 0 : using dg_package_data_temporary_tags = tmpl::list<>;
134 0 : using dg_package_data_primitive_tags = tmpl::list<>;
135 0 : using dg_package_data_volume_tags = tmpl::list<>;
136 0 : using dg_boundary_terms_volume_tags = tmpl::list<>;
137 :
138 0 : double dg_package_data(
139 : const gsl::not_null<typename Tags::TildeE<
140 : Frame::Inertial, NeutrinoSpecies>::type*>... packaged_tilde_e,
141 : const gsl::not_null<typename Tags::TildeS<
142 : Frame::Inertial, NeutrinoSpecies>::type*>... packaged_tilde_s,
143 : const gsl::not_null<typename Tags::TildeE<
144 : Frame::Inertial,
145 : NeutrinoSpecies>::type*>... packaged_normal_dot_flux_tilde_e,
146 : const gsl::not_null<typename Tags::TildeS<
147 : Frame::Inertial,
148 : NeutrinoSpecies>::type*>... packaged_normal_dot_flux_tilde_s,
149 :
150 : const typename Tags::TildeE<Frame::Inertial,
151 : NeutrinoSpecies>::type&... tilde_e,
152 : const typename Tags::TildeS<Frame::Inertial,
153 : NeutrinoSpecies>::type&... tilde_s,
154 : const typename ::Tags::Flux<
155 : Tags::TildeE<Frame::Inertial, NeutrinoSpecies>, tmpl::size_t<3>,
156 : Frame::Inertial>::type&... flux_tilde_e,
157 : const typename ::Tags::Flux<
158 : Tags::TildeS<Frame::Inertial, NeutrinoSpecies>, tmpl::size_t<3>,
159 : Frame::Inertial>::type&... flux_tilde_s,
160 :
161 : const tnsr::i<DataVector, 3, Frame::Inertial>& normal_covector,
162 : const tnsr::I<DataVector, 3, Frame::Inertial>& normal_vector,
163 : const std::optional<tnsr::I<DataVector, 3, Frame::Inertial>>&
164 : mesh_velocity,
165 : const std::optional<Scalar<DataVector>>& normal_dot_mesh_velocity) const {
166 : EXPAND_PACK_LEFT_TO_RIGHT(Rusanov_detail::dg_package_data_impl(
167 : packaged_tilde_e, packaged_tilde_s, packaged_normal_dot_flux_tilde_e,
168 : packaged_normal_dot_flux_tilde_s, tilde_e, tilde_s, flux_tilde_e,
169 : flux_tilde_s, normal_covector, normal_vector, mesh_velocity,
170 : normal_dot_mesh_velocity));
171 : // max speed = speed of light
172 : return 1.0;
173 : }
174 :
175 0 : void dg_boundary_terms(
176 : const gsl::not_null<typename Tags::TildeE<
177 : Frame::Inertial,
178 : NeutrinoSpecies>::type*>... boundary_correction_tilde_e,
179 : const gsl::not_null<typename Tags::TildeS<
180 : Frame::Inertial,
181 : NeutrinoSpecies>::type*>... boundary_correction_tilde_s,
182 : const typename Tags::TildeE<Frame::Inertial,
183 : NeutrinoSpecies>::type&... tilde_e_int,
184 : const typename Tags::TildeS<Frame::Inertial,
185 : NeutrinoSpecies>::type&... tilde_s_int,
186 : const typename Tags::TildeE<Frame::Inertial, NeutrinoSpecies>::
187 : type&... normal_dot_flux_tilde_e_int,
188 : const typename Tags::TildeS<Frame::Inertial, NeutrinoSpecies>::
189 : type&... normal_dot_flux_tilde_s_int,
190 : const typename Tags::TildeE<Frame::Inertial,
191 : NeutrinoSpecies>::type&... tilde_e_ext,
192 : const typename Tags::TildeS<Frame::Inertial,
193 : NeutrinoSpecies>::type&... tilde_s_ext,
194 : const typename Tags::TildeE<Frame::Inertial, NeutrinoSpecies>::
195 : type&... normal_dot_flux_tilde_e_ext,
196 : const typename Tags::TildeS<Frame::Inertial, NeutrinoSpecies>::
197 : type&... normal_dot_flux_tilde_s_ext,
198 : dg::Formulation dg_formulation) const {
199 : EXPAND_PACK_LEFT_TO_RIGHT(Rusanov_detail::dg_boundary_terms_impl(
200 : boundary_correction_tilde_e, boundary_correction_tilde_s, tilde_e_int,
201 : tilde_s_int, normal_dot_flux_tilde_e_int, normal_dot_flux_tilde_s_int,
202 : tilde_e_ext, tilde_s_ext, normal_dot_flux_tilde_e_ext,
203 : normal_dot_flux_tilde_s_ext, dg_formulation));
204 : }
205 : };
206 :
207 : /// \cond
208 : template <typename... NeutrinoSpecies>
209 : // NOLINTNEXTLINE
210 : PUP::able::PUP_ID Rusanov<tmpl::list<NeutrinoSpecies...>>::my_PUP_ID = 0;
211 : /// \endcond
212 :
213 : } // namespace RadiationTransport::M1Grey::BoundaryCorrections
|