Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #pragma once
11 :
12 : #include "MooseObject.h"
13 : #include "MooseEnum.h"
14 : #include "MooseTypes.h"
15 : #include "InputParameters.h"
16 : #include "THMEnums.h"
17 :
18 : #include "libmesh/elem.h"
19 : #include "libmesh/vector_value.h"
20 : #include "libmesh/point.h"
21 :
22 : /**
23 : * Interface class for 1-D slope reconstruction
24 : *
25 : * This class provides interfaces for generating slopes for an arbitrary set
26 : * of variables. A number of reconstruction options are provided, including a
27 : * number of TVD slope limiters.
28 : */
29 : template <bool is_ad>
30 : class SlopeReconstruction1DInterface
31 : {
32 : public:
33 : SlopeReconstruction1DInterface(const MooseObject * moose_object);
34 :
35 : /// Slope reconstruction type
36 : enum ESlopeReconstructionType
37 : {
38 : None, ///< No reconstruction; Godunov scheme
39 : Full, ///< Full reconstruction; no limitation
40 : Minmod, ///< Minmod slope limiter
41 : MC, ///< Monotonized Central-Difference slope limiter
42 : Superbee ///< Superbee slope limiter
43 : };
44 : /// Map of slope reconstruction type string to enum
45 : static const std::map<std::string, ESlopeReconstructionType> _slope_reconstruction_type_to_enum;
46 :
47 : /**
48 : * Gets a MooseEnum for slope reconstruction type
49 : *
50 : * @param[in] name default value
51 : * @returns MooseEnum for slope reconstruction type
52 : */
53 : static MooseEnum getSlopeReconstructionMooseEnum(const std::string & name = "");
54 :
55 : protected:
56 : /**
57 : * Gets the primitive solution vector and position of neighbor(s)
58 : *
59 : * @param[in] elem Element for which to get slopes
60 : * @param[in] W_neighbor Primitive solution vector for each neighbor
61 : * @param[in] x_neighbor Position for each neighbor
62 : */
63 : virtual void
64 : getNeighborPrimitiveVariables(const Elem * elem,
65 : std::vector<std::vector<GenericReal<is_ad>>> & W_neighbor,
66 : std::vector<Point> & x_neighbor) const;
67 :
68 : /**
69 : * Gets limited slopes for the primitive variables in the 1-D direction
70 : *
71 : * @param[in] elem Element for which to get slopes
72 : *
73 : * @returns Vector of slopes for the element in the 1-D direction
74 : */
75 : std::vector<GenericReal<is_ad>> getElementSlopes(const Elem * elem) const;
76 :
77 : /**
78 : * Gets limited slopes for the primitive variables in the 1-D direction for boundary element
79 : *
80 : * @param[in] W_elem Primitive solution for element
81 : * @param[in] x_elem Position for element
82 : * @param[in] dir Direction for element
83 : * @param[in] W_neighbor Primitive solution vector for each neighbor
84 : * @param[in] x_neighbor Position for each neighbor
85 : * @param[in] W_boundary Primitive solution vector for boundary
86 : *
87 : * @returns Vector of slopes for the element in the 1-D direction
88 : */
89 : std::vector<GenericReal<is_ad>>
90 : getBoundaryElementSlopes(const std::vector<GenericReal<is_ad>> & W_elem,
91 : const Point & x_elem,
92 : const RealVectorValue & dir,
93 : std::vector<std::vector<GenericReal<is_ad>>> W_neighbor,
94 : std::vector<Point> x_neighbor,
95 : const std::vector<GenericReal<is_ad>> & W_boundary) const;
96 :
97 : /**
98 : * Gets limited slopes for the primitive variables in the 1-D direction
99 : *
100 : * @param[in] W_elem Primitive solution for element
101 : * @param[in] x_elem Position for element
102 : * @param[in] dir Direction for element
103 : * @param[in] W_neighbor Primitive solution vector for each neighbor
104 : * @param[in] x_neighbor Position for each neighbor
105 : *
106 : * @returns Vector of slopes for the element in the 1-D direction
107 : */
108 : std::vector<GenericReal<is_ad>>
109 : getElementSlopes(const std::vector<GenericReal<is_ad>> & W_elem,
110 : const Point & x_elem,
111 : const RealVectorValue & dir,
112 : const std::vector<std::vector<GenericReal<is_ad>>> & W_neighbor,
113 : const std::vector<Point> & x_neighbor) const;
114 :
115 : /**
116 : * Computes the cell-average primitive variable values for an element
117 : *
118 : * @param[in] elem Element for which to get values
119 : *
120 : * @returns Vector of values on element
121 : */
122 : virtual std::vector<GenericReal<is_ad>>
123 : computeElementPrimitiveVariables(const Elem * elem) const = 0;
124 :
125 : /// MooseObject this interface is extending
126 : const MooseObject * const _moose_object;
127 :
128 : /// Slope reconstruction scheme
129 : const ESlopeReconstructionType _scheme;
130 :
131 : public:
132 : static InputParameters validParams();
133 :
134 : protected:
135 : /// Number of sides
136 : static const unsigned int _n_side;
137 : /// Number of elemental values in stencil for computing slopes
138 : static const unsigned int _n_sten;
139 : };
140 :
141 : namespace THM
142 : {
143 : template <>
144 : SlopeReconstruction1DInterface<true>::ESlopeReconstructionType
145 : stringToEnum<SlopeReconstruction1DInterface<true>::ESlopeReconstructionType>(const std::string & s);
146 :
147 : template <>
148 : SlopeReconstruction1DInterface<false>::ESlopeReconstructionType
149 : stringToEnum<SlopeReconstruction1DInterface<false>::ESlopeReconstructionType>(
150 : const std::string & s);
151 : }
152 :
153 : template <bool is_ad>
154 : const std::map<std::string,
155 : typename SlopeReconstruction1DInterface<is_ad>::ESlopeReconstructionType>
156 : SlopeReconstruction1DInterface<is_ad>::_slope_reconstruction_type_to_enum{
157 : {"NONE", None}, {"FULL", Full}, {"MINMOD", Minmod}, {"MC", MC}, {"SUPERBEE", Superbee}};
158 :
159 : template <bool is_ad>
160 : MooseEnum
161 : SlopeReconstruction1DInterface<is_ad>::getSlopeReconstructionMooseEnum(const std::string & name)
162 : {
163 : return THM::getMooseEnum<SlopeReconstruction1DInterface<is_ad>::ESlopeReconstructionType>(
164 28194 : name, _slope_reconstruction_type_to_enum);
165 : }
166 :
167 : template <bool is_ad>
168 : const unsigned int SlopeReconstruction1DInterface<is_ad>::_n_side = 2;
169 :
170 : template <bool is_ad>
171 : const unsigned int SlopeReconstruction1DInterface<is_ad>::_n_sten = 3;
172 :
173 : template <bool is_ad>
174 : InputParameters
175 19453 : SlopeReconstruction1DInterface<is_ad>::validParams()
176 : {
177 19453 : InputParameters params = emptyInputParameters();
178 :
179 77812 : params.addRequiredParam<MooseEnum>(
180 : "scheme",
181 : SlopeReconstruction1DInterface::getSlopeReconstructionMooseEnum("None"),
182 : "Slope reconstruction scheme");
183 :
184 19453 : return params;
185 0 : }
186 :
187 : template <bool is_ad>
188 15044 : SlopeReconstruction1DInterface<is_ad>::SlopeReconstruction1DInterface(
189 : const MooseObject * moose_object)
190 15044 : : _moose_object(moose_object),
191 15044 : _scheme(THM::stringToEnum<ESlopeReconstructionType>(
192 15044 : moose_object->parameters().get<MooseEnum>("scheme")))
193 : {
194 15044 : }
195 :
196 : template <bool is_ad>
197 : void
198 4709918 : SlopeReconstruction1DInterface<is_ad>::getNeighborPrimitiveVariables(
199 : const Elem * elem,
200 : std::vector<std::vector<GenericReal<is_ad>>> & W_neighbor,
201 : std::vector<Point> & x_neighbor) const
202 : {
203 4709918 : W_neighbor.clear();
204 4709918 : x_neighbor.clear();
205 14129754 : for (unsigned int i_side = 0; i_side < _n_side; i_side++)
206 : {
207 : auto neighbor = elem->neighbor_ptr(i_side);
208 9419836 : if (neighbor && (neighbor->processor_id() == _moose_object->processor_id()))
209 : {
210 8935490 : x_neighbor.push_back(neighbor->vertex_average());
211 17870980 : W_neighbor.push_back(computeElementPrimitiveVariables(neighbor));
212 : }
213 : }
214 4709918 : }
215 :
216 : template <bool is_ad>
217 : std::vector<GenericReal<is_ad>>
218 4653824 : SlopeReconstruction1DInterface<is_ad>::getElementSlopes(const Elem * elem) const
219 : {
220 : mooseAssert(elem, "The supplied element is a nullptr.");
221 :
222 4653824 : const auto W_elem = computeElementPrimitiveVariables(elem);
223 4653824 : const Point x_elem = elem->vertex_average();
224 4653824 : const RealVectorValue dir = (elem->node_ref(1) - elem->node_ref(0)).unit();
225 :
226 : std::vector<Point> x_neighbor;
227 : std::vector<std::vector<GenericReal<is_ad>>> W_neighbor;
228 4653824 : getNeighborPrimitiveVariables(elem, W_neighbor, x_neighbor);
229 :
230 9307648 : return getElementSlopes(W_elem, x_elem, dir, W_neighbor, x_neighbor);
231 4653824 : }
232 :
233 : template <bool is_ad>
234 : std::vector<GenericReal<is_ad>>
235 65038 : SlopeReconstruction1DInterface<is_ad>::getBoundaryElementSlopes(
236 : const std::vector<GenericReal<is_ad>> & W_elem,
237 : const Point & x_elem,
238 : const RealVectorValue & dir,
239 : std::vector<std::vector<GenericReal<is_ad>>> W_neighbor,
240 : std::vector<Point> x_neighbor,
241 : const std::vector<GenericReal<is_ad>> & W_boundary) const
242 : {
243 65038 : if (W_neighbor.size() == 1)
244 : {
245 64974 : W_neighbor.push_back(W_boundary);
246 :
247 : // The boundary point will be assumed to be the same distance away as neighbor
248 : const Point dx = x_elem - x_neighbor[0];
249 : const Point x_boundary = x_elem + dx;
250 64974 : x_neighbor.push_back(x_boundary);
251 : }
252 :
253 65038 : return getElementSlopes(W_elem, x_elem, dir, W_neighbor, x_neighbor);
254 : }
255 :
256 : template <bool is_ad>
257 : std::vector<GenericReal<is_ad>>
258 4718862 : SlopeReconstruction1DInterface<is_ad>::getElementSlopes(
259 : const std::vector<GenericReal<is_ad>> & W_elem,
260 : const Point & x_elem,
261 : const RealVectorValue & dir,
262 : const std::vector<std::vector<GenericReal<is_ad>>> & W_neighbor,
263 : const std::vector<Point> & x_neighbor) const
264 : {
265 : mooseAssert(x_neighbor.size() == W_neighbor.size(),
266 : "Neighbor positions size must equal neighbor solutions size.");
267 :
268 : using std::abs, std::min, std::max;
269 :
270 : // get the number of slopes to be stored
271 4718862 : const unsigned int n_slopes = W_elem.size();
272 :
273 : // compute one-sided slope(s)
274 : std::vector<std::vector<GenericReal<is_ad>>> slopes_one_sided;
275 13728246 : for (unsigned int i = 0; i < W_neighbor.size(); i++)
276 : {
277 9009384 : const Real dx = (x_elem - x_neighbor[i]) * dir;
278 :
279 9009384 : std::vector<GenericReal<is_ad>> slopes(n_slopes, 0.0);
280 36037536 : for (unsigned int m = 0; m < n_slopes; m++)
281 54056304 : slopes[m] = (W_elem[m] - W_neighbor[i][m]) / dx;
282 :
283 9009384 : slopes_one_sided.push_back(slopes);
284 : }
285 :
286 : // Fill in any missing one-sided slopes and compute central slope
287 4718862 : std::vector<GenericReal<is_ad>> slopes_central(n_slopes, 0.0);
288 4718862 : if (W_neighbor.size() == 2)
289 : {
290 4290586 : const Real dx = (x_neighbor[0] - x_neighbor[1]) * dir;
291 17162344 : for (unsigned int m = 0; m < n_slopes; m++)
292 25743516 : slopes_central[m] = (W_neighbor[0][m] - W_neighbor[1][m]) / dx;
293 : }
294 428276 : else if (W_neighbor.size() == 1)
295 : {
296 428212 : slopes_one_sided.push_back(slopes_one_sided[0]);
297 428212 : slopes_central = slopes_one_sided[0];
298 : }
299 : else // only one element; use zero slopes
300 : {
301 64 : slopes_one_sided.push_back(slopes_central);
302 64 : slopes_one_sided.push_back(slopes_central);
303 : }
304 :
305 : // vector for the (possibly limited) slopes
306 4718862 : std::vector<GenericReal<is_ad>> slopes_limited(n_slopes, 0.0);
307 :
308 : // limit the slopes
309 4718862 : switch (_scheme)
310 : {
311 : // first-order, zero slope
312 : case None:
313 : break;
314 :
315 : // full reconstruction; no limitation
316 1030 : case Full:
317 :
318 1030 : slopes_limited = slopes_central;
319 : break;
320 :
321 : // minmod limiter
322 : case Minmod:
323 :
324 17155736 : for (unsigned int m = 0; m < n_slopes; m++)
325 : {
326 25733604 : if ((slopes_one_sided[0][m] * slopes_one_sided[1][m]) > 0.0)
327 : {
328 10749859 : if (abs(slopes_one_sided[0][m]) < abs(slopes_one_sided[1][m]))
329 4366878 : slopes_limited[m] = slopes_one_sided[0][m];
330 : else
331 6382981 : slopes_limited[m] = slopes_one_sided[1][m];
332 : }
333 : }
334 : break;
335 :
336 : // MC (monotonized central-difference) limiter
337 : case MC:
338 :
339 780800 : for (unsigned int m = 0; m < n_slopes; m++)
340 : {
341 585600 : if (slopes_central[m] > 0.0 && slopes_one_sided[0][m] > 0.0 && slopes_one_sided[1][m] > 0.0)
342 3856 : slopes_limited[m] =
343 3856 : min(slopes_central[m], 2.0 * min(slopes_one_sided[0][m], slopes_one_sided[1][m]));
344 583672 : else if (slopes_central[m] < 0.0 && slopes_one_sided[0][m] < 0.0 &&
345 0 : slopes_one_sided[1][m] < 0.0)
346 3408 : slopes_limited[m] =
347 3408 : max(slopes_central[m], 2.0 * max(slopes_one_sided[0][m], slopes_one_sided[1][m]));
348 : }
349 : break;
350 :
351 : // superbee limiter
352 : case Superbee:
353 :
354 780800 : for (unsigned int m = 0; m < n_slopes; m++)
355 : {
356 585600 : GenericReal<is_ad> slope1 = 0.0;
357 585600 : GenericReal<is_ad> slope2 = 0.0;
358 :
359 : // calculate slope1 with minmod
360 585600 : if (slopes_one_sided[1][m] > 0.0 && slopes_one_sided[0][m] > 0.0)
361 2992 : slope1 = min(slopes_one_sided[1][m], 2.0 * slopes_one_sided[0][m]);
362 584104 : else if (slopes_one_sided[1][m] < 0.0 && slopes_one_sided[0][m] < 0.0)
363 3776 : slope1 = max(slopes_one_sided[1][m], 2.0 * slopes_one_sided[0][m]);
364 :
365 : // calculate slope2 with minmod
366 585600 : if (slopes_one_sided[1][m] > 0.0 && slopes_one_sided[0][m] > 0.0)
367 2992 : slope2 = min(2.0 * slopes_one_sided[1][m], slopes_one_sided[0][m]);
368 584104 : else if (slopes_one_sided[1][m] < 0.0 && slopes_one_sided[0][m] < 0.0)
369 3776 : slope2 = max(2.0 * slopes_one_sided[1][m], slopes_one_sided[0][m]);
370 :
371 : // calculate slope with maxmod
372 585600 : if (slope1 > 0.0 && slope2 > 0.0)
373 1496 : slopes_limited[m] = max(slope1, slope2);
374 584104 : else if (slope1 < 0.0 && slope2 < 0.0)
375 1888 : slopes_limited[m] = min(slope1, slope2);
376 : }
377 : break;
378 :
379 0 : default:
380 0 : mooseError("Unknown slope reconstruction scheme");
381 : break;
382 : }
383 4718862 : return slopes_limited;
384 4718862 : }
|
