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 <cmath>
8 : #include <cstddef>
9 :
10 : #include "NumericalAlgorithms/FiniteDifference/Reconstruct.hpp"
11 : #include "NumericalAlgorithms/FiniteDifference/Unlimited.hpp"
12 : #include "Utilities/ConstantExpressions.hpp"
13 : #include "Utilities/ErrorHandling/Assert.hpp"
14 : #include "Utilities/ForceInline.hpp"
15 : #include "Utilities/Gsl.hpp"
16 : #include "Utilities/Math.hpp"
17 :
18 : /// \cond
19 : class DataVector;
20 : template <size_t>
21 : class Direction;
22 : template <size_t Dim, typename T>
23 : class DirectionMap;
24 : template <size_t Dim>
25 : class Index;
26 : /// \endcond
27 :
28 : namespace fd::reconstruction {
29 : namespace detail {
30 : struct MonotonicityPreserving5Reconstructor {
31 : SPECTRE_ALWAYS_INLINE static std::array<double, 2> pointwise(
32 : const double* const q, const int stride, const double alpha,
33 : const double epsilon) {
34 : using std::abs;
35 : using std::max;
36 : using std::min;
37 :
38 : // define minmod function for 2 and 4 args
39 : const auto minmod2 = [](const double x, const double y) {
40 : return 0.5 * (sgn(x) + sgn(y)) * min(abs(x), abs(y));
41 : };
42 : const auto minmod4 = [](const double w, const double x, const double y,
43 : const double z) {
44 : const double sign_w = sgn(w);
45 : return 0.125 * (sign_w + sgn(x)) *
46 : abs((sign_w + sgn(y)) * (sign_w + sgn(z))) *
47 : min(abs(w), min(abs(x), min(abs(y), abs(z))));
48 : };
49 :
50 : // first, compute unlimited fifth-order finite difference reconstruction
51 : // values
52 : auto result = UnlimitedReconstructor<4>::pointwise(q, stride);
53 :
54 : // compute q_{j+1/2}
55 : const double q_mp_plus =
56 : q[0] + minmod2(q[stride] - q[0], alpha * (q[0] - q[-stride]));
57 : // compute q_{j-1/2}
58 : const double q_mp_minus =
59 : q[0] + minmod2(q[-stride] - q[0], alpha * (q[0] - q[stride]));
60 :
61 : const bool limit_q_plus =
62 : ((result[1] - q[0]) * (result[1] - q_mp_plus) > epsilon);
63 : const bool limit_q_minus =
64 : ((result[0] - q[0]) * (result[0] - q_mp_minus) > epsilon);
65 :
66 : if (limit_q_plus or limit_q_minus) {
67 : const double dp = q[2 * stride] + q[0] - 2.0 * q[stride];
68 : const double dj = q[stride] + q[-stride] - 2.0 * q[0];
69 : const double dm = q[0] + q[-2 * stride] - 2.0 * q[-stride];
70 : const double dm4_plus = minmod4(4.0 * dj - dp, 4 * dp - dj, dj, dp);
71 : const double dm4_minus = minmod4(4.0 * dj - dm, 4 * dm - dj, dj, dm);
72 :
73 : if (limit_q_plus) {
74 : const double q_ul = q[0] + alpha * (q[0] - q[-stride]);
75 : const double q_md =
76 : 0.5 * (q[0] + q[stride] - dm4_plus); // inline q^{AV}
77 : const double q_lc =
78 : q[0] + 0.5 * (q[0] - q[-stride]) + 1.3333333333333333 * dm4_minus;
79 : const double q_min =
80 : max(min(q[0], min(q[stride], q_md)), min(q[0], min(q_ul, q_lc)));
81 : const double q_max =
82 : min(max(q[0], max(q[stride], q_md)), max(q[0], max(q_ul, q_lc)));
83 :
84 : result[1] = result[1] + minmod2(q_min - result[1], q_max - result[1]);
85 : }
86 :
87 : if (limit_q_minus) {
88 : const double q_ul = q[0] + alpha * (q[0] - q[stride]);
89 : const double q_md =
90 : 0.5 * (q[0] + q[-stride] - dm4_minus); // inline q^{AV}
91 : const double q_lc =
92 : q[0] + 0.5 * (q[0] - q[stride]) + 1.3333333333333333 * dm4_plus;
93 : const double q_min =
94 : max(min(q[0], min(q[-stride], q_md)), min(q[0], min(q_ul, q_lc)));
95 : const double q_max =
96 : min(max(q[0], max(q[-stride], q_md)), max(q[0], max(q_ul, q_lc)));
97 :
98 : result[0] = result[0] + minmod2(q_min - result[0], q_max - result[0]);
99 : }
100 : }
101 :
102 : return result;
103 : }
104 :
105 : SPECTRE_ALWAYS_INLINE static constexpr size_t stencil_width() { return 5; }
106 : };
107 : } // namespace detail
108 :
109 : /*!
110 : * \ingroup FiniteDifferenceGroup
111 : * \brief Performs the fifth order monotonicity-preserving (MP5) reconstruction
112 : * \cite Suresh1997.
113 : *
114 : * First, calculate the original interface value \f$q_{j+1/2}^\text{OR}\f$ with
115 : * the (unlimited) fifth order finite difference reconstruction
116 : *
117 : * \f{align}
118 : * q_{j+1/2}^\text{OR} = \frac{1}{128}( 3 q_{j-2} - 20 q_{j-1} + 90 q_{j}
119 : * + 60 q_{j+1} - 5 q_{j+2}) .
120 : * \f}
121 : *
122 : * and compute
123 : *
124 : * \f{align}
125 : * q^\text{MP} = q_j + \text{minmod}(q_{j+1} - q_j, \alpha(q_j - q_{j-1}))
126 : * \f}
127 : *
128 : * for a given value of \f$\alpha\f$, which is usually set to \f$\alpha=4\f$.
129 : *
130 : * If \f$ (q_{j+1/2}^\text{OR} - q_j)(q_{j+1/2}^\text{OR} - q^\text{MP}) \leq
131 : * \epsilon\f$ where \f$\epsilon\f$ is a small tolerance value, use $q_{1/2}
132 : * = q_{j+1/2}^\text{OR}$.
133 : *
134 : * \note A proper value of \f$\epsilon\f$ may depend on the scale of the
135 : * quantity \f$q\f$ to be reconstructed; reference \cite Suresh1997 suggests
136 : * 1e-10. For hydro simulations with atmosphere treatment, setting
137 : * \f$\epsilon=0.0\f$ would be safe.
138 : *
139 : * Otherwise, calculate
140 : *
141 : * \f{align}
142 : * d_{j+1} & = q_{j+2} - 2q_{j+1} + q_{j} \\
143 : * d_j & = q_{j+1} - 2q_j + q_{j-1} \\
144 : * d_{j-1} & = q_{j} - 2q_{j-1} + q_{j-2} ,
145 : * \f}
146 : *
147 : * \f{align}
148 : * d^\text{M4}_{j+1/2} =
149 : * \text{minmod}(4d_j - d_{j+1}, 4d_{j+1}-d_j, d_j, d_{j+1}) \\
150 : * d^\text{M4}_{j-1/2} =
151 : * \text{minmod}(4d_j - d_{j-1}, 4d_{j-1}-d_j, d_j, d_{j-1}),
152 : * \f}
153 : *
154 : * \f{align}
155 : * q^\text{UL} & = q_j + \alpha(q_j - q_{j-1}) \\
156 : * q^\text{AV} & = (q_j + q_{j+1}) / 2 \\
157 : * q^\text{MD} & = q^\text{AV} - d^\text{M4}_{j+1/2} / 2 \\
158 : * q^\text{LC} & = q_j + (q_j - q_{j-1}) / 2 + (4/3) d^\text{M4}_{j-1/2}
159 : * \f}
160 : *
161 : * and
162 : * \f{align}
163 : * q^\text{min} & = \text{max}[ \text{min}(q_j, q_{j+1}, q^\text{MD}),
164 : * \text{min}(q_j, q^\text{UL}, q^\text{LC}) ] \\
165 : * q^\text{max} & = \text{min}[ \text{max}(q_j, q_{j+1}, q^\text{MD}),
166 : * \text{max}(q_j, q^\text{UL}, q^\text{LC}) ].
167 : * \f}
168 : *
169 : * The reconstructed value is given as
170 : * \f{align}
171 : * q_{i+1/2} = \text{median}(q_{j+1/2}^\text{OR}, q^\text{min}, q^\text{max})
172 : * \f}
173 : * where the median function can be expressed as
174 : * \f{align}
175 : * \text{median}(x,y,z) = x + \text{minmod}(y-x, z-x).
176 : * \f}
177 : *
178 : */
179 : template <size_t Dim>
180 1 : void monotonicity_preserving_5(
181 : const gsl::not_null<std::array<gsl::span<double>, Dim>*>
182 : reconstructed_upper_side_of_face_vars,
183 : const gsl::not_null<std::array<gsl::span<double>, Dim>*>
184 : reconstructed_lower_side_of_face_vars,
185 : const gsl::span<const double>& volume_vars,
186 : const DirectionMap<Dim, gsl::span<const double>>& ghost_cell_vars,
187 : const Index<Dim>& volume_extents, const size_t number_of_variables,
188 : const double alpha, const double epsilon) {
189 : detail::reconstruct<detail::MonotonicityPreserving5Reconstructor>(
190 : reconstructed_upper_side_of_face_vars,
191 : reconstructed_lower_side_of_face_vars, volume_vars, ghost_cell_vars,
192 : volume_extents, number_of_variables, alpha, epsilon);
193 : }
194 : } // namespace fd::reconstruction
|
