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 <boost/functional/hash.hpp>
8 : #include <cstddef>
9 : #include <optional>
10 : #include <string>
11 : #include <unordered_map>
12 : #include <utility>
13 :
14 : #include "DataStructures/Tensor/TypeAliases.hpp"
15 : #include "Domain/SizeOfElement.hpp"
16 : #include "Domain/Tags.hpp"
17 : #include "Evolution/DiscontinuousGalerkin/Limiters/Tags.hpp"
18 : #include "Evolution/DiscontinuousGalerkin/Limiters/Weno.hpp"
19 : #include "Evolution/DiscontinuousGalerkin/Limiters/WenoType.hpp"
20 : #include "Evolution/DiscontinuousGalerkin/NormalVectorTags.hpp"
21 : #include "Evolution/Systems/NewtonianEuler/Limiters/VariablesToLimit.hpp"
22 : #include "Evolution/Systems/NewtonianEuler/Tags.hpp"
23 : #include "Options/Auto.hpp"
24 : #include "Options/String.hpp"
25 : #include "PointwiseFunctions/Hydro/TagsDeclarations.hpp"
26 : #include "Utilities/TMPL.hpp"
27 :
28 : /// \cond
29 : class DataVector;
30 : template <size_t VolumeDim>
31 : class Direction;
32 : template <size_t VolumeDim>
33 : class Element;
34 : template <size_t VolumeDim>
35 : class ElementId;
36 : template <size_t VolumeDim>
37 : class Mesh;
38 : template <size_t VolumeDim>
39 : class OrientationMap;
40 :
41 : namespace EquationsOfState {
42 : template <bool IsRelativistic, size_t ThermodynamicDim>
43 : class EquationOfState;
44 : } // namespace EquationsOfState
45 :
46 : namespace gsl {
47 : template <typename T>
48 : class not_null;
49 : } // namespace gsl
50 :
51 : namespace PUP {
52 : class er;
53 : } // namespace PUP
54 : /// \endcond
55 :
56 : namespace NewtonianEuler {
57 : namespace Limiters {
58 :
59 : /// \ingroup LimitersGroup
60 : /// \brief A compact-stencil WENO limiter for the NewtonianEuler system.
61 : ///
62 : /// Implements the simple WENO limiter of \cite Zhong2013 and the Hermite WENO
63 : /// (HWENO) limiter of \cite Zhu2016. See the documentation of the
64 : /// system-agnostic ::Limiters::Weno limiter for a general discussion of the
65 : /// algorithm and the various options that control the action of the limiter.
66 : //
67 : /// This implemention is specialized to the NewtonianEuler evolution system.
68 : /// By specializing the limiter to the system, we can add a few features that
69 : /// improve its robustness:
70 : /// - the troubled-cell indicator (TCI) can be specialized to the features of
71 : /// the evolution system.
72 : /// - the limiter can be applied to the system's characteristic variables. This
73 : /// is the recommendation of the reference, because it reduces spurious
74 : /// oscillations in the post-limiter solution.
75 : /// - after limiting, the solution can be processed to remove any remaining
76 : /// unphysical values like negative densities and pressures. We do this by
77 : /// scaling the solution around its mean (a "flattener" or "bounds-preserving"
78 : /// filter). Note: the flattener is applied to all elements, including those
79 : /// where the limiter did not act to reduce the solution's slopes.
80 : ///
81 : /// The matrix of TCI, variables to limit, post-processing, etc. choices can
82 : /// rapidly grow large. Here we reduce the possibilities by tying the TCI to the
83 : /// limiter in keeping with each limiter's main reference: HWENO uses the KXRCF
84 : /// TCI and simple WENO uses the TVB TCI. To fully explore the matrix of
85 : /// possibilities, the source code could be generalized --- however, experience
86 : /// suggests it is unlikely that there exists one combination that will perform
87 : /// remarkably better than the others.
88 : template <size_t VolumeDim>
89 1 : class Weno {
90 : public:
91 0 : using ConservativeVarsWeno = ::Limiters::Weno<
92 : VolumeDim, tmpl::list<NewtonianEuler::Tags::MassDensityCons,
93 : NewtonianEuler::Tags::MomentumDensity<VolumeDim>,
94 : NewtonianEuler::Tags::EnergyDensity>>;
95 :
96 0 : struct VariablesToLimit {
97 0 : using type = NewtonianEuler::Limiters::VariablesToLimit;
98 0 : static type suggested_value() {
99 : return NewtonianEuler::Limiters::VariablesToLimit::Characteristic;
100 : }
101 0 : static constexpr Options::String help = {
102 : "Variable representation on which to apply the limiter"};
103 : };
104 : // Future design improvement: attach the TvbConstant/KxrcfConstant to the
105 : // limiter type, so that it isn't necessary to specify both (but with one
106 : // required to be 'None') in each input file.
107 0 : struct TvbConstant {
108 0 : using type = Options::Auto<double, Options::AutoLabel::None>;
109 0 : static constexpr Options::String help = {
110 : "Constant in RHS of the TVB minmod TCI, used when Type = SimpleWeno"};
111 : };
112 0 : struct KxrcfConstant {
113 0 : using type = Options::Auto<double, Options::AutoLabel::None>;
114 0 : static constexpr Options::String help = {
115 : "Constant in RHS of KXRCF TCI, used when Type = Hweno"};
116 : };
117 0 : struct ApplyFlattener {
118 0 : using type = bool;
119 0 : static constexpr Options::String help = {
120 : "Flatten after limiting to restore pointwise positivity"};
121 : };
122 0 : using options =
123 : tmpl::list<typename ConservativeVarsWeno::Type, VariablesToLimit,
124 : typename ConservativeVarsWeno::NeighborWeight, TvbConstant,
125 : KxrcfConstant, ApplyFlattener,
126 : typename ConservativeVarsWeno::DisableForDebugging>;
127 0 : static constexpr Options::String help = {
128 : "A WENO limiter specialized to the NewtonianEuler system"};
129 0 : static std::string name() { return "NewtonianEulerWeno"; };
130 :
131 0 : Weno(::Limiters::WenoType weno_type,
132 : NewtonianEuler::Limiters::VariablesToLimit vars_to_limit,
133 : double neighbor_linear_weight, std::optional<double> tvb_constant,
134 : std::optional<double> kxrcf_constant, bool apply_flattener,
135 : bool disable_for_debugging = false,
136 : const Options::Context& context = {});
137 :
138 0 : Weno() = default;
139 0 : Weno(const Weno& /*rhs*/) = default;
140 0 : Weno& operator=(const Weno& /*rhs*/) = default;
141 0 : Weno(Weno&& /*rhs*/) = default;
142 0 : Weno& operator=(Weno&& /*rhs*/) = default;
143 0 : ~Weno() = default;
144 :
145 : // NOLINTNEXTLINE(google-runtime-references)
146 0 : void pup(PUP::er& p);
147 :
148 0 : using PackagedData = typename ConservativeVarsWeno::PackagedData;
149 0 : using package_argument_tags =
150 : typename ConservativeVarsWeno::package_argument_tags;
151 :
152 : /// \brief Package data for sending to neighbor elements
153 1 : void package_data(gsl::not_null<PackagedData*> packaged_data,
154 : const Scalar<DataVector>& mass_density_cons,
155 : const tnsr::I<DataVector, VolumeDim>& momentum_density,
156 : const Scalar<DataVector>& energy_density,
157 : const Mesh<VolumeDim>& mesh,
158 : const std::array<double, VolumeDim>& element_size,
159 : const OrientationMap<VolumeDim>& orientation_map) const;
160 :
161 0 : using limit_tags =
162 : tmpl::list<NewtonianEuler::Tags::MassDensityCons,
163 : NewtonianEuler::Tags::MomentumDensity<VolumeDim>,
164 : NewtonianEuler::Tags::EnergyDensity>;
165 0 : using limit_argument_tags = tmpl::list<
166 : domain::Tags::Mesh<VolumeDim>, domain::Tags::Element<VolumeDim>,
167 : domain::Tags::SizeOfElement<VolumeDim>,
168 : domain::Tags::DetInvJacobian<Frame::ElementLogical, Frame::Inertial>,
169 : evolution::dg::Tags::NormalCovectorAndMagnitude<VolumeDim>,
170 : ::hydro::Tags::EquationOfState<false, 2>>;
171 :
172 : /// \brief Limit the solution on the element
173 1 : bool operator()(
174 : gsl::not_null<Scalar<DataVector>*> mass_density_cons,
175 : gsl::not_null<tnsr::I<DataVector, VolumeDim>*> momentum_density,
176 : gsl::not_null<Scalar<DataVector>*> energy_density,
177 : const Mesh<VolumeDim>& mesh, const Element<VolumeDim>& element,
178 : const std::array<double, VolumeDim>& element_size,
179 : const Scalar<DataVector>& det_inv_logical_to_inertial_jacobian,
180 : const typename evolution::dg::Tags::NormalCovectorAndMagnitude<
181 : VolumeDim>::type& normals_and_magnitudes,
182 : const EquationsOfState::EquationOfState<false, 2>& equation_of_state,
183 : const std::unordered_map<DirectionalId<VolumeDim>, PackagedData,
184 : boost::hash<DirectionalId<VolumeDim>>>&
185 : neighbor_data) const;
186 :
187 : private:
188 : template <size_t LocalDim>
189 : // NOLINTNEXTLINE(readability-redundant-declaration) false positive
190 0 : friend bool operator==(const Weno<LocalDim>& lhs, const Weno<LocalDim>& rhs);
191 :
192 0 : ::Limiters::WenoType weno_type_;
193 0 : NewtonianEuler::Limiters::VariablesToLimit vars_to_limit_;
194 0 : double neighbor_linear_weight_;
195 0 : std::optional<double> tvb_constant_;
196 0 : std::optional<double> kxrcf_constant_;
197 0 : bool apply_flattener_;
198 0 : bool disable_for_debugging_;
199 : // Note: conservative_vars_weno_ is always used for calls to package_data, and
200 : // is also used when limiting cons vars with the simple WENO algorithm. So we
201 : // construct conservative_vars_weno_ with the correct TVB constant when it is
202 : // used for limiting (precisely, when weno_type_ == SimpleWeno), and with a
203 : // dummy TVB constant value of NaN otherwise (weno_type_ == Hweno). This lets
204 : // us construct the conservative_vars_weno_ variable and use it for delegating
205 : // the package_data work even when the specialized limiter has no TVB constant
206 : // and is possible because package_data doesn't depend on the TCI.
207 0 : ConservativeVarsWeno conservative_vars_weno_;
208 : };
209 :
210 : template <size_t VolumeDim>
211 0 : bool operator!=(const Weno<VolumeDim>& lhs, const Weno<VolumeDim>& rhs);
212 :
213 : } // namespace Limiters
214 : } // namespace NewtonianEuler
|
