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