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 <cstddef>
8 : #include <limits>
9 : #include <memory>
10 : #include <utility>
11 :
12 : #include "DataStructures/VariablesTag.hpp"
13 : #include "Domain/Structure/DirectionalIdMap.hpp"
14 : #include "Domain/Tags.hpp"
15 : #include "Evolution/DgSubcell/Tags/GhostDataForReconstruction.hpp"
16 : #include "Evolution/DgSubcell/Tags/Mesh.hpp"
17 : #include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/ReconstructWork.hpp"
18 : #include "Evolution/Systems/GrMhd/ValenciaDivClean/FiniteDifference/Reconstructor.hpp"
19 : #include "Evolution/VariableFixing/FixToAtmosphere.hpp"
20 : #include "Evolution/VariableFixing/Tags.hpp"
21 : #include "NumericalAlgorithms/FiniteDifference/FallbackReconstructorType.hpp"
22 : #include "Options/Auto.hpp"
23 : #include "Options/Context.hpp"
24 : #include "Options/String.hpp"
25 : #include "PointwiseFunctions/Hydro/Tags.hpp"
26 : #include "Utilities/TMPL.hpp"
27 :
28 : /// \cond
29 : class DataVector;
30 : template <size_t Dim>
31 : class Direction;
32 : template <size_t Dim>
33 : class Element;
34 : template <size_t Dim>
35 : class ElementId;
36 : namespace EquationsOfState {
37 : template <bool IsRelativistic, size_t ThermodynamicDim>
38 : class EquationOfState;
39 : } // namespace EquationsOfState
40 : template <size_t Dim>
41 : class Mesh;
42 : namespace gsl {
43 : template <typename T>
44 : class not_null;
45 : } // namespace gsl
46 : namespace PUP {
47 : class er;
48 : } // namespace PUP
49 : template <typename TagsList>
50 : class Variables;
51 : namespace evolution::dg::subcell {
52 : class GhostData;
53 : } // namespace evolution::dg::subcell
54 : /// \endcond
55 :
56 : namespace grmhd::ValenciaDivClean::fd {
57 : /*!
58 : * \brief Positivity-preserving adaptive order reconstruction. See
59 : * ::fd::reconstruction::positivity_preserving_adaptive_order() for details.
60 : *
61 : * The rest mass density, electron fraction, and the pressure are kept positive.
62 : */
63 1 : class PositivityPreservingAdaptiveOrderPrim : public Reconstructor {
64 : private:
65 0 : using prims_to_reconstruct_tags =
66 : tmpl::list<hydro::Tags::RestMassDensity<DataVector>,
67 : hydro::Tags::ElectronFraction<DataVector>,
68 : hydro::Tags::Temperature<DataVector>,
69 : hydro::Tags::LorentzFactorTimesSpatialVelocity<DataVector, 3>,
70 : hydro::Tags::MagneticField<DataVector, 3>,
71 : hydro::Tags::DivergenceCleaningField<DataVector>>;
72 :
73 0 : using positivity_preserving_tags =
74 : tmpl::list<hydro::Tags::RestMassDensity<DataVector>,
75 : hydro::Tags::ElectronFraction<DataVector>,
76 : hydro::Tags::Temperature<DataVector>>;
77 0 : using non_positive_tags =
78 : tmpl::list<hydro::Tags::LorentzFactorTimesSpatialVelocity<DataVector, 3>,
79 : hydro::Tags::MagneticField<DataVector, 3>,
80 : hydro::Tags::DivergenceCleaningField<DataVector>>;
81 :
82 0 : using FallbackReconstructorType =
83 : ::fd::reconstruction::FallbackReconstructorType;
84 :
85 : public:
86 0 : static constexpr size_t dim = 3;
87 :
88 0 : struct Alpha5 {
89 0 : using type = double;
90 0 : static constexpr Options::String help = {
91 : "The alpha parameter in the Persson convergence measurement. 4 is the "
92 : "right value, but anything in the range of 3-5 is 'reasonable'. "
93 : "Smaller values allow for more oscillations."};
94 : };
95 0 : struct Alpha7 {
96 0 : using type = Options::Auto<double, Options::AutoLabel::None>;
97 0 : static constexpr Options::String help = {
98 : "The alpha parameter in the Persson convergence measurement. 4 is the "
99 : "right value, but anything in the range of 3-5 is 'reasonable'. "
100 : "Smaller values allow for more oscillations. If specified to None, "
101 : "then 7th-order reconstruction is not used."};
102 : };
103 0 : struct Alpha9 {
104 0 : using type = Options::Auto<double, Options::AutoLabel::None>;
105 0 : static constexpr Options::String help = {
106 : "The alpha parameter in the Persson convergence measurement. 4 is the "
107 : "right value, but anything in the range of 3-5 is 'reasonable'. "
108 : "Smaller values allow for more oscillations. If specified to None, "
109 : "then 9th-order reconstruction is not used."};
110 : };
111 0 : struct LowOrderReconstructor {
112 0 : using type = FallbackReconstructorType;
113 0 : static constexpr Options::String help = {
114 : "The 2nd/3rd-order reconstruction scheme to use if unlimited 5th-order "
115 : "isn't okay."};
116 : };
117 0 : struct ReconstructRhoTimesTemperature {
118 0 : using type = bool;
119 0 : static constexpr Options::String help = {
120 : "If 'true' then we reconstruct the rho*T, if 'false' we reconstruct "
121 : "T."};
122 : };
123 :
124 0 : using options = tmpl::list<Alpha5, Alpha7, Alpha9, LowOrderReconstructor,
125 : ReconstructRhoTimesTemperature>;
126 :
127 0 : static constexpr Options::String help{
128 : "Positivity-preserving adaptive-order reconstruction."};
129 :
130 0 : PositivityPreservingAdaptiveOrderPrim() = default;
131 0 : PositivityPreservingAdaptiveOrderPrim(
132 : PositivityPreservingAdaptiveOrderPrim&&) = default;
133 0 : PositivityPreservingAdaptiveOrderPrim& operator=(
134 : PositivityPreservingAdaptiveOrderPrim&&) = default;
135 0 : PositivityPreservingAdaptiveOrderPrim(
136 : const PositivityPreservingAdaptiveOrderPrim&) = default;
137 0 : PositivityPreservingAdaptiveOrderPrim& operator=(
138 : const PositivityPreservingAdaptiveOrderPrim&) = default;
139 0 : ~PositivityPreservingAdaptiveOrderPrim() override = default;
140 :
141 0 : PositivityPreservingAdaptiveOrderPrim(
142 : double alpha_5, std::optional<double> alpha_7,
143 : std::optional<double> alpha_9,
144 : FallbackReconstructorType low_order_reconstructor,
145 : bool reconstruct_rho_times_temperature,
146 : const Options::Context& context = {});
147 :
148 0 : explicit PositivityPreservingAdaptiveOrderPrim(CkMigrateMessage* msg);
149 :
150 0 : WRAPPED_PUPable_decl_base_template(Reconstructor,
151 : PositivityPreservingAdaptiveOrderPrim);
152 :
153 0 : auto get_clone() const -> std::unique_ptr<Reconstructor> override;
154 :
155 0 : static constexpr bool use_adaptive_order = true;
156 0 : bool supports_adaptive_order() const override { return use_adaptive_order; }
157 :
158 0 : void pup(PUP::er& p) override;
159 :
160 0 : size_t ghost_zone_size() const override {
161 : return eight_to_the_alpha_9_.has_value()
162 : ? 5
163 : : (six_to_the_alpha_7_.has_value() ? 4 : 3);
164 : }
165 :
166 0 : using reconstruction_argument_tags =
167 : tmpl::list<::Tags::Variables<hydro::grmhd_tags<DataVector>>,
168 : hydro::Tags::GrmhdEquationOfState, domain::Tags::Element<dim>,
169 : evolution::dg::subcell::Tags::GhostDataForReconstruction<dim>,
170 : evolution::dg::subcell::Tags::Mesh<dim>,
171 : ::Tags::VariableFixer<VariableFixing::FixToAtmosphere<dim>>>;
172 :
173 : template <size_t ThermodynamicDim>
174 0 : void reconstruct(
175 : gsl::not_null<std::array<Variables<tags_list_for_reconstruct>, dim>*>
176 : vars_on_lower_face,
177 : gsl::not_null<std::array<Variables<tags_list_for_reconstruct>, dim>*>
178 : vars_on_upper_face,
179 : gsl::not_null<std::optional<std::array<gsl::span<std::uint8_t>, dim>>*>
180 : reconstruction_order,
181 : const Variables<hydro::grmhd_tags<DataVector>>& volume_prims,
182 : const EquationsOfState::EquationOfState<true, ThermodynamicDim>& eos,
183 : const Element<dim>& element,
184 : const DirectionalIdMap<dim, evolution::dg::subcell::GhostData>&
185 : ghost_data,
186 : const Mesh<dim>& subcell_mesh,
187 : const VariableFixing::FixToAtmosphere<dim>& fix_to_atmosphere) const;
188 :
189 : template <size_t ThermodynamicDim>
190 0 : void reconstruct_fd_neighbor(
191 : gsl::not_null<Variables<tags_list_for_reconstruct>*> vars_on_face,
192 : const Variables<hydro::grmhd_tags<DataVector>>& subcell_volume_prims,
193 : const EquationsOfState::EquationOfState<true, ThermodynamicDim>& eos,
194 : const Element<dim>& element,
195 : const DirectionalIdMap<dim, evolution::dg::subcell::GhostData>&
196 : ghost_data,
197 : const Mesh<dim>& subcell_mesh,
198 : const VariableFixing::FixToAtmosphere<dim>& fix_to_atmosphere,
199 : Direction<dim> direction_to_reconstruct) const;
200 :
201 0 : bool reconstruct_rho_times_temperature() const override;
202 :
203 : private:
204 : // NOLINTNEXTLINE(readability-redundant-declaration)
205 0 : friend bool operator==(const PositivityPreservingAdaptiveOrderPrim& lhs,
206 : const PositivityPreservingAdaptiveOrderPrim& rhs);
207 0 : friend bool operator!=(const PositivityPreservingAdaptiveOrderPrim& lhs,
208 : const PositivityPreservingAdaptiveOrderPrim& rhs);
209 0 : void set_function_pointers();
210 :
211 0 : double four_to_the_alpha_5_ = std::numeric_limits<double>::signaling_NaN();
212 0 : std::optional<double> six_to_the_alpha_7_{};
213 0 : std::optional<double> eight_to_the_alpha_9_{};
214 0 : FallbackReconstructorType low_order_reconstructor_ =
215 : FallbackReconstructorType::None;
216 0 : bool reconstruct_rho_times_temperature_{false};
217 :
218 0 : using PointerReconsOrder = void (*)(
219 : gsl::not_null<std::array<gsl::span<double>, dim>*>,
220 : gsl::not_null<std::array<gsl::span<double>, dim>*>,
221 : gsl::not_null<std::optional<std::array<gsl::span<std::uint8_t>, dim>>*>,
222 : const gsl::span<const double>&,
223 : const DirectionMap<dim, gsl::span<const double>>&, const Index<dim>&,
224 : size_t, double, double, double);
225 0 : using PointerRecons =
226 : void (*)(gsl::not_null<std::array<gsl::span<double>, dim>*>,
227 : gsl::not_null<std::array<gsl::span<double>, dim>*>,
228 : const gsl::span<const double>&,
229 : const DirectionMap<dim, gsl::span<const double>>&,
230 : const Index<dim>&, size_t, double, double, double);
231 0 : PointerRecons reconstruct_ = nullptr;
232 0 : PointerReconsOrder pp_reconstruct_ = nullptr;
233 :
234 0 : using PointerNeighbor = void (*)(gsl::not_null<DataVector*>,
235 : const DataVector&, const DataVector&,
236 : const Index<dim>&, const Index<dim>&,
237 : const Direction<dim>&, const double&,
238 : const double&, const double&);
239 0 : PointerNeighbor reconstruct_lower_neighbor_ = nullptr;
240 0 : PointerNeighbor reconstruct_upper_neighbor_ = nullptr;
241 0 : PointerNeighbor pp_reconstruct_lower_neighbor_ = nullptr;
242 0 : PointerNeighbor pp_reconstruct_upper_neighbor_ = nullptr;
243 : };
244 :
245 : } // namespace grmhd::ValenciaDivClean::fd
|
