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