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