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 :
11 : #include "DataStructures/DataBox/PrefixHelpers.hpp"
12 : #include "DataStructures/DataBox/Prefixes.hpp"
13 : #include "DataStructures/DataVector.hpp"
14 : #include "DataStructures/Tensor/TypeAliases.hpp"
15 : #include "Domain/Structure/DirectionalIdMap.hpp"
16 : #include "Domain/Tags.hpp"
17 : #include "Evolution/DgSubcell/GhostData.hpp"
18 : #include "Evolution/DgSubcell/Tags/GhostDataForReconstruction.hpp"
19 : #include "Evolution/DgSubcell/Tags/Mesh.hpp"
20 : #include "Evolution/Systems/NewtonianEuler/FiniteDifference/Reconstructor.hpp"
21 : #include "Evolution/Systems/NewtonianEuler/Tags.hpp"
22 : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
23 : #include "PointwiseFunctions/Hydro/Tags.hpp"
24 :
25 1 : namespace NewtonianEuler::fd {
26 : /*!
27 : * \brief Adaptive-order WENO reconstruction hybridizing orders 5 and 3. See
28 : * ::fd::reconstruction::aoweno_53() for details.
29 : */
30 : template <size_t Dim>
31 1 : class AoWeno53Prim : public Reconstructor<Dim> {
32 : private:
33 : // Conservative vars tags
34 0 : using MassDensityCons = NewtonianEuler::Tags::MassDensityCons;
35 0 : using EnergyDensity = NewtonianEuler::Tags::EnergyDensity;
36 0 : using MomentumDensity = NewtonianEuler::Tags::MomentumDensity<Dim>;
37 :
38 : // Primitive vars tags
39 0 : using MassDensity = hydro::Tags::RestMassDensity<DataVector>;
40 0 : using Velocity = hydro::Tags::SpatialVelocity<DataVector, Dim>;
41 0 : using SpecificInternalEnergy =
42 : hydro::Tags::SpecificInternalEnergy<DataVector>;
43 0 : using Pressure = hydro::Tags::Pressure<DataVector>;
44 :
45 0 : using prims_tags =
46 : tmpl::list<MassDensity, Velocity, SpecificInternalEnergy, Pressure>;
47 0 : using cons_tags = tmpl::list<MassDensityCons, MomentumDensity, EnergyDensity>;
48 0 : using flux_tags = db::wrap_tags_in<::Tags::Flux, cons_tags, tmpl::size_t<Dim>,
49 : Frame::Inertial>;
50 0 : using prim_tags_for_reconstruction =
51 : tmpl::list<MassDensity, Velocity, Pressure>;
52 :
53 : public:
54 0 : struct GammaHi {
55 0 : using type = double;
56 0 : static constexpr Options::String help = {
57 : "The linear weight for the 5th-order stencil."};
58 : };
59 0 : struct GammaLo {
60 0 : using type = double;
61 0 : static constexpr Options::String help = {
62 : "The linear weight for the central 3rd-order stencil."};
63 : };
64 0 : struct Epsilon {
65 0 : using type = double;
66 0 : static constexpr Options::String help = {
67 : "The parameter added to the oscillation indicators to avoid division "
68 : "by zero"};
69 : };
70 0 : struct NonlinearWeightExponent {
71 0 : using type = size_t;
72 0 : static constexpr Options::String help = {
73 : "The exponent q to which the oscillation indicators are raised"};
74 : };
75 :
76 0 : using options =
77 : tmpl::list<GammaHi, GammaLo, Epsilon, NonlinearWeightExponent>;
78 0 : static constexpr Options::String help{
79 : "Adaptive-order WENO reconstruction hybridizing orders 5 and 3 using "
80 : "primitive variables."};
81 :
82 0 : AoWeno53Prim() = default;
83 0 : AoWeno53Prim(AoWeno53Prim&&) = default;
84 0 : AoWeno53Prim& operator=(AoWeno53Prim&&) = default;
85 0 : AoWeno53Prim(const AoWeno53Prim&) = default;
86 0 : AoWeno53Prim& operator=(const AoWeno53Prim&) = default;
87 0 : ~AoWeno53Prim() override = default;
88 :
89 0 : AoWeno53Prim(double gamma_hi, double gamma_lo, double epsilon,
90 : size_t nonlinear_weight_exponent);
91 :
92 0 : explicit AoWeno53Prim(CkMigrateMessage* msg);
93 :
94 0 : WRAPPED_PUPable_decl_base_template(Reconstructor<Dim>, AoWeno53Prim);
95 :
96 0 : auto get_clone() const -> std::unique_ptr<Reconstructor<Dim>> override;
97 :
98 0 : void pup(PUP::er& p) override;
99 :
100 0 : size_t ghost_zone_size() const override { return 3; }
101 :
102 0 : using reconstruction_argument_tags =
103 : tmpl::list<::Tags::Variables<prims_tags>,
104 : hydro::Tags::EquationOfState<false, 2>,
105 : domain::Tags::Element<Dim>,
106 : evolution::dg::subcell::Tags::GhostDataForReconstruction<Dim>,
107 : evolution::dg::subcell::Tags::Mesh<Dim>>;
108 :
109 : template <typename TagsList>
110 0 : void reconstruct(
111 : gsl::not_null<std::array<Variables<TagsList>, Dim>*> vars_on_lower_face,
112 : gsl::not_null<std::array<Variables<TagsList>, Dim>*> vars_on_upper_face,
113 : const Variables<prims_tags>& volume_prims,
114 : const EquationsOfState::EquationOfState<false, 2>& eos,
115 : const Element<Dim>& element,
116 : const DirectionalIdMap<Dim, evolution::dg::subcell::GhostData>&
117 : ghost_data,
118 : const Mesh<Dim>& subcell_mesh) const;
119 :
120 : /// Called by an element doing DG when the neighbor is doing subcell.
121 : template <typename TagsList>
122 1 : void reconstruct_fd_neighbor(
123 : gsl::not_null<Variables<TagsList>*> vars_on_face,
124 : const Variables<prims_tags>& subcell_volume_prims,
125 : const EquationsOfState::EquationOfState<false, 2>& eos,
126 : const Element<Dim>& element,
127 : const DirectionalIdMap<Dim, evolution::dg::subcell::GhostData>&
128 : ghost_data,
129 : const Mesh<Dim>& subcell_mesh,
130 : const Direction<Dim> direction_to_reconstruct) const;
131 :
132 : private:
133 : template <size_t LocalDim>
134 : // NOLINTNEXTLINE(readability-redundant-declaration)
135 0 : friend bool operator==(const AoWeno53Prim<LocalDim>& lhs,
136 : const AoWeno53Prim<LocalDim>& rhs);
137 :
138 0 : double gamma_hi_ = std::numeric_limits<double>::signaling_NaN();
139 0 : double gamma_lo_ = std::numeric_limits<double>::signaling_NaN();
140 0 : double epsilon_ = std::numeric_limits<double>::signaling_NaN();
141 0 : size_t nonlinear_weight_exponent_ = 0;
142 :
143 0 : void (*reconstruct_)(gsl::not_null<std::array<gsl::span<double>, Dim>*>,
144 : gsl::not_null<std::array<gsl::span<double>, Dim>*>,
145 : const gsl::span<const double>&,
146 : const DirectionMap<Dim, gsl::span<const double>>&,
147 : const Index<Dim>&, size_t, double, double, double);
148 0 : void (*reconstruct_lower_neighbor_)(gsl::not_null<DataVector*>,
149 : const DataVector&, const DataVector&,
150 : const Index<Dim>&, const Index<Dim>&,
151 : const Direction<Dim>&, const double&,
152 : const double&, const double&);
153 0 : void (*reconstruct_upper_neighbor_)(gsl::not_null<DataVector*>,
154 : const DataVector&, const DataVector&,
155 : const Index<Dim>&, const Index<Dim>&,
156 : const Direction<Dim>&, const double&,
157 : const double&, const double&);
158 : };
159 :
160 : template <size_t Dim>
161 0 : bool operator!=(const AoWeno53Prim<Dim>& lhs, const AoWeno53Prim<Dim>& rhs) {
162 : return not(lhs == rhs);
163 : }
164 : } // namespace NewtonianEuler::fd
|
