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