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 4710482 : 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 : W_neighbor.clear();
204 : x_neighbor.clear();
205 14131446 : for (unsigned int i_side = 0; i_side < _n_side; i_side++)
206 : {
207 : auto neighbor = elem->neighbor_ptr(i_side);
208 9420964 : if (neighbor && (neighbor->processor_id() == _moose_object->processor_id()))
209 : {
210 8936534 : x_neighbor.push_back(neighbor->vertex_average());
211 17873068 : W_neighbor.push_back(computeElementPrimitiveVariables(neighbor));
212 : }
213 : }
214 4710482 : }
215 :
216 : template <bool is_ad>
217 : std::vector<GenericReal<is_ad>>
218 4654356 : SlopeReconstruction1DInterface<is_ad>::getElementSlopes(const Elem * elem) const
219 : {
220 : mooseAssert(elem, "The supplied element is a nullptr.");
221 :
222 4654356 : const auto W_elem = computeElementPrimitiveVariables(elem);
223 4654356 : const Point x_elem = elem->vertex_average();
224 4654356 : 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 4654356 : getNeighborPrimitiveVariables(elem, W_neighbor, x_neighbor);
229 :
230 9308712 : return getElementSlopes(W_elem, x_elem, dir, W_neighbor, x_neighbor);
231 4654356 : }
232 :
233 : template <bool is_ad>
234 : std::vector<GenericReal<is_ad>>
235 65102 : 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 65102 : if (W_neighbor.size() == 1)
244 : {
245 65038 : 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 65038 : x_neighbor.push_back(x_boundary);
251 : }
252 :
253 65102 : 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 4719458 : 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 : // get the number of slopes to be stored
269 4719458 : const unsigned int n_slopes = W_elem.size();
270 :
271 : // compute one-sided slope(s)
272 : std::vector<std::vector<GenericReal<is_ad>>> slopes_one_sided;
273 13729982 : for (unsigned int i = 0; i < W_neighbor.size(); i++)
274 : {
275 9010524 : const Real dx = (x_elem - x_neighbor[i]) * dir;
276 :
277 9010524 : std::vector<GenericReal<is_ad>> slopes(n_slopes, 0.0);
278 36042096 : for (unsigned int m = 0; m < n_slopes; m++)
279 54063144 : slopes[m] = (W_elem[m] - W_neighbor[i][m]) / dx;
280 :
281 9010524 : slopes_one_sided.push_back(slopes);
282 : }
283 :
284 : // Fill in any missing one-sided slopes and compute central slope
285 4719458 : std::vector<GenericReal<is_ad>> slopes_central(n_slopes, 0.0);
286 4719458 : if (W_neighbor.size() == 2)
287 : {
288 4291130 : const Real dx = (x_neighbor[0] - x_neighbor[1]) * dir;
289 17164520 : for (unsigned int m = 0; m < n_slopes; m++)
290 25746780 : slopes_central[m] = (W_neighbor[0][m] - W_neighbor[1][m]) / dx;
291 : }
292 428328 : else if (W_neighbor.size() == 1)
293 : {
294 428264 : slopes_one_sided.push_back(slopes_one_sided[0]);
295 428264 : slopes_central = slopes_one_sided[0];
296 : }
297 : else // only one element; use zero slopes
298 : {
299 64 : slopes_one_sided.push_back(slopes_central);
300 64 : slopes_one_sided.push_back(slopes_central);
301 : }
302 :
303 : // vector for the (possibly limited) slopes
304 4719458 : std::vector<GenericReal<is_ad>> slopes_limited(n_slopes, 0.0);
305 :
306 : // limit the slopes
307 4719458 : switch (_scheme)
308 : {
309 : // first-order, zero slope
310 : case None:
311 : break;
312 :
313 : // full reconstruction; no limitation
314 1030 : case Full:
315 :
316 1030 : slopes_limited = slopes_central;
317 : break;
318 :
319 : // minmod limiter
320 : case Minmod:
321 :
322 17158120 : for (unsigned int m = 0; m < n_slopes; m++)
323 : {
324 25737180 : if ((slopes_one_sided[0][m] * slopes_one_sided[1][m]) > 0.0)
325 : {
326 10751639 : if (std::abs(slopes_one_sided[0][m]) < std::abs(slopes_one_sided[1][m]))
327 4367843 : slopes_limited[m] = slopes_one_sided[0][m];
328 : else
329 6383796 : slopes_limited[m] = slopes_one_sided[1][m];
330 : }
331 : }
332 : break;
333 :
334 : // MC (monotonized central-difference) limiter
335 : case MC:
336 :
337 780800 : for (unsigned int m = 0; m < n_slopes; m++)
338 : {
339 585600 : if (slopes_central[m] > 0.0 && slopes_one_sided[0][m] > 0.0 && slopes_one_sided[1][m] > 0.0)
340 3856 : slopes_limited[m] = std::min(
341 3856 : slopes_central[m], 2.0 * std::min(slopes_one_sided[0][m], slopes_one_sided[1][m]));
342 583672 : else if (slopes_central[m] < 0.0 && slopes_one_sided[0][m] < 0.0 &&
343 0 : slopes_one_sided[1][m] < 0.0)
344 3408 : slopes_limited[m] = std::max(
345 3408 : slopes_central[m], 2.0 * std::max(slopes_one_sided[0][m], slopes_one_sided[1][m]));
346 : }
347 : break;
348 :
349 : // superbee limiter
350 : case Superbee:
351 :
352 780800 : for (unsigned int m = 0; m < n_slopes; m++)
353 : {
354 585600 : GenericReal<is_ad> slope1 = 0.0;
355 585600 : GenericReal<is_ad> slope2 = 0.0;
356 :
357 : // calculate slope1 with minmod
358 585600 : if (slopes_one_sided[1][m] > 0.0 && slopes_one_sided[0][m] > 0.0)
359 2992 : slope1 = std::min(slopes_one_sided[1][m], 2.0 * slopes_one_sided[0][m]);
360 584104 : else if (slopes_one_sided[1][m] < 0.0 && slopes_one_sided[0][m] < 0.0)
361 3776 : slope1 = std::max(slopes_one_sided[1][m], 2.0 * slopes_one_sided[0][m]);
362 :
363 : // calculate slope2 with minmod
364 585600 : if (slopes_one_sided[1][m] > 0.0 && slopes_one_sided[0][m] > 0.0)
365 2992 : slope2 = std::min(2.0 * slopes_one_sided[1][m], slopes_one_sided[0][m]);
366 584104 : else if (slopes_one_sided[1][m] < 0.0 && slopes_one_sided[0][m] < 0.0)
367 3776 : slope2 = std::max(2.0 * slopes_one_sided[1][m], slopes_one_sided[0][m]);
368 :
369 : // calculate slope with maxmod
370 585600 : if (slope1 > 0.0 && slope2 > 0.0)
371 1496 : slopes_limited[m] = std::max(slope1, slope2);
372 584104 : else if (slope1 < 0.0 && slope2 < 0.0)
373 1888 : slopes_limited[m] = std::min(slope1, slope2);
374 : }
375 : break;
376 :
377 0 : default:
378 0 : mooseError("Unknown slope reconstruction scheme");
379 : break;
380 : }
381 4719458 : return slopes_limited;
382 4719458 : }
|
