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 "MooseError.h"
13 : #include "FaceInfo.h"
14 : #include "Limiter.h"
15 : #include "MathUtils.h"
16 : #include "MooseFunctor.h"
17 : #include "metaphysicl/raw_type.h"
18 : #include "libmesh/compare_types.h"
19 : #include "libmesh/elem.h"
20 : #include <cmath>
21 : #include <limits>
22 : #include <tuple>
23 :
24 : template <typename>
25 : class MooseVariableFV;
26 : class FaceArgInterface;
27 :
28 : namespace Moose
29 : {
30 : namespace FV
31 : {
32 : /// This codifies a set of available ways to interpolate with elem+neighbor
33 : /// solution information to calculate values (e.g. solution, material
34 : /// properties, etc.) at the face (centroid). These methods are used in the
35 : /// class's interpolate functions. Some interpolation methods are only meant
36 : /// to be used with advective terms (e.g. upwind), others are more
37 : /// generically applicable.
38 : enum class InterpMethod
39 : {
40 : /// gc*elem+(1-gc)*neighbor
41 : Average,
42 : /// 1/(gc/elem+(1-gc)/neighbor)
43 : HarmonicAverage,
44 : /// (gc*elem+(1-gc)*neighbor)+gradient*(rf-rf')
45 : SkewCorrectedAverage,
46 : /// weighted
47 : Upwind,
48 : // Rhie-Chow
49 : RhieChow,
50 : VanLeer,
51 : MinMod,
52 : SOU,
53 : QUICK,
54 : Venkatakrishnan
55 : };
56 :
57 : enum class LinearFVComputationMode
58 : {
59 : RHS,
60 : Matrix,
61 : FullSystem
62 : };
63 :
64 : /**
65 : * Returns an enum with all the currently supported interpolation methods and the current default
66 : * for FV: first-order upwind
67 : * @return MooseEnum with all the face interpolation methods supported
68 : */
69 : MooseEnum interpolationMethods();
70 :
71 : /**
72 : * @return An input parameters object that contains the \p advected_interp_method parameter, e.g.
73 : * the interpolation method to use for an advected quantity
74 : */
75 : InputParameters advectedInterpolationParameter();
76 :
77 : /*
78 : * Converts from the interpolation method to the interpolation enum.
79 : * This routine is here in lieu of using a MooseEnum for InterpMethod
80 : * @param interp_method the name of the interpolation method
81 : * @return the interpolation method
82 : */
83 : InterpMethod selectInterpolationMethod(const std::string & interp_method);
84 :
85 : /**
86 : * Sets one interpolation method
87 : * @param obj The \p MooseObject with input parameters to query
88 : * @param interp_method The interpolation method we will set
89 : * @param param_name The name of the parameter setting this interpolation method
90 : * @return Whether the interpolation method has indicated that we will need more than the
91 : * default level of ghosting
92 : */
93 : bool setInterpolationMethod(const MooseObject & obj,
94 : Moose::FV::InterpMethod & interp_method,
95 : const std::string & param_name);
96 :
97 : /**
98 : * Produce the interpolation coefficients in the equation:
99 : *
100 : * \phi_f = c_1 * \phi_{F1} + c_2 * \phi_{F2}
101 : *
102 : * A couple of examples: if we are doing an average interpolation with
103 : * an orthogonal regular grid, then the pair will be (0.5, 0.5). If we are doing an
104 : * upwind interpolation with the velocity facing outward from the F1 element,
105 : * then the pair will be (1.0, 0.0).
106 : *
107 : * The template is needed in case the face flux carries derivatives in an AD setting.
108 : *
109 : * @param m The interpolation method
110 : * @param fi The face information
111 : * @param one_is_elem Whether fi.elem() == F1
112 : * @param face_flux The advecting face flux. Not relevant for an Average interpolation
113 : * @return a pair where the first Real is c_1 and the second Real is c_2
114 : */
115 : template <typename T = Real>
116 : std::pair<Real, Real>
117 169887091 : interpCoeffs(const InterpMethod m,
118 : const FaceInfo & fi,
119 : const bool one_is_elem,
120 : const T & face_flux = 0.0)
121 : {
122 169887091 : switch (m)
123 : {
124 169887091 : case InterpMethod::Average:
125 : case InterpMethod::SkewCorrectedAverage:
126 : {
127 169887091 : if (one_is_elem)
128 169887091 : return std::make_pair(fi.gC(), 1. - fi.gC());
129 : else
130 0 : return std::make_pair(1. - fi.gC(), fi.gC());
131 : }
132 :
133 0 : case InterpMethod::Upwind:
134 : {
135 0 : if ((face_flux > 0) == one_is_elem)
136 0 : return std::make_pair(1., 0.);
137 : else
138 0 : return std::make_pair(0., 1.);
139 : }
140 :
141 0 : default:
142 0 : mooseError("Unrecognized interpolation method");
143 : }
144 : }
145 :
146 : /**
147 : * A simple linear interpolation of values between cell centers to a cell face. The \p one_is_elem
148 : * parameter indicates whether value1 corresponds to the FaceInfo elem value; else it corresponds to
149 : * the FaceInfo neighbor value
150 : */
151 : template <typename T, typename T2>
152 : typename libMesh::CompareTypes<T, T2>::supertype
153 161326240 : linearInterpolation(const T & value1,
154 : const T2 & value2,
155 : const FaceInfo & fi,
156 : const bool one_is_elem,
157 : const InterpMethod interp_method = InterpMethod::Average)
158 : {
159 : mooseAssert(interp_method == InterpMethod::Average ||
160 : interp_method == InterpMethod::SkewCorrectedAverage,
161 : "The selected interpolation function only works with average or skewness-corrected "
162 : "average options!");
163 161326240 : const auto coeffs = interpCoeffs(interp_method, fi, one_is_elem);
164 161326240 : return coeffs.first * value1 + coeffs.second * value2;
165 : }
166 :
167 : /**
168 : * Computes the harmonic mean (1/(gc/value1+(1-gc)/value2)) of Reals, RealVectorValues and
169 : * RealTensorValues while accounting for the possibility that one or both of them are AD.
170 : * For tensors, we use a component-wise mean instead of the matrix-inverse based option.
171 : * @param value1 Reference to value1 in the (1/(gc/value1+(1-gc)/value2)) expression
172 : * @param value2 Reference to value2 in the (1/(gc/value1+(1-gc)/value2)) expression
173 : * @param fi Reference to the FaceInfo of the face on which the interpolation is requested
174 : * @param one_is_elem Boolean indicating if the interpolation weight on FaceInfo belongs to the
175 : * element corresponding to value1
176 : * @return The interpolated value
177 : */
178 : template <typename T1, typename T2>
179 : typename libMesh::CompareTypes<T1, T2>::supertype
180 8060920 : harmonicInterpolation(const T1 & value1,
181 : const T2 & value2,
182 : const FaceInfo & fi,
183 : const bool one_is_elem)
184 : {
185 : // We check if the base values of the given template types match, if not we throw a compile-time
186 : // error
187 : static_assert(std::is_same<typename MetaPhysicL::RawType<T1>::value_type,
188 : typename MetaPhysicL::RawType<T2>::value_type>::value,
189 : "The input values for harmonic interpolation need to have the same raw-value!");
190 :
191 : // Fetch the interpolation coefficients, we use the exact same coefficients as for a simple
192 : // weighted average
193 8060920 : const auto coeffs = interpCoeffs(InterpMethod::Average, fi, one_is_elem);
194 :
195 : // We check if the types are fit to compute the harmonic mean of. This is done compile-time
196 : // using constexpr. We start with Real/ADReal which is straightforward if the input values are
197 : // positive.
198 : if constexpr (libMesh::TensorTools::TensorTraits<T1>::rank == 0)
199 : {
200 : // The harmonic mean of mixed positive and negative numbers (and 0 as well) is not well-defined
201 : // so we assert that the input values shall be positive.
202 : #ifndef NDEBUG
203 : if (value1 * value2 <= 0)
204 : mooseWarning("Input values (" + Moose::stringify(MetaPhysicL::raw_value(value1)) + " & " +
205 : Moose::stringify(MetaPhysicL::raw_value(value2)) +
206 : ") must be of the same sign for harmonic interpolation");
207 : #endif
208 16044872 : return 1.0 / (coeffs.first / value1 + coeffs.second / value2);
209 : }
210 : // For vectors (ADRealVectorValue, VectorValue), we take the component-wise harmonic mean
211 : else if constexpr (libMesh::TensorTools::TensorTraits<T1>::rank == 1)
212 : {
213 38400 : typename libMesh::CompareTypes<T1, T2>::supertype result;
214 153600 : for (const auto i : make_range(Moose::dim))
215 : {
216 : #ifndef NDEBUG
217 : if (value1(i) * value2(i) <= 0)
218 : mooseWarning("Component " + std::to_string(i) + " of input values (" +
219 : Moose::stringify(MetaPhysicL::raw_value(value1(i))) + " & " +
220 : Moose::stringify(MetaPhysicL::raw_value(value2(i))) +
221 : ") must be of the same sign for harmonic interpolation");
222 : #endif
223 115200 : result(i) = 1.0 / (coeffs.first / value1(i) + coeffs.second / value2(i));
224 : }
225 76800 : return result;
226 38400 : }
227 : // For tensors (ADRealTensorValue, TensorValue), similarly to the vectors,
228 : // we take the component-wise harmonic mean instead of the matrix-inverse approach
229 : else if constexpr (libMesh::TensorTools::TensorTraits<T1>::rank == 2)
230 : {
231 : typename libMesh::CompareTypes<T1, T2>::supertype result;
232 : for (const auto i : make_range(Moose::dim))
233 : for (const auto j : make_range(Moose::dim))
234 : {
235 : #ifndef NDEBUG
236 : if (value1(i, j) * value2(i, j) <= 0)
237 : mooseWarning("Component (" + std::to_string(i) + "," + std::to_string(j) +
238 : ") of input values (" +
239 : Moose::stringify(MetaPhysicL::raw_value(value1(i, j))) + " & " +
240 : Moose::stringify(MetaPhysicL::raw_value(value2(i, j))) +
241 : ") must be of the same sign for harmonic interpolation");
242 : #endif
243 : result(i, j) = 1.0 / (coeffs.first / value1(i, j) + coeffs.second / value2(i, j));
244 : }
245 : return result;
246 : }
247 : // We ran out of options, harmonic mean is not supported for other types at the moment
248 : else
249 : // This line is supposed to throw an error when the user tries to compile this function with
250 : // types that are not supported. This is the reason we needed the always_false function. Hope as
251 : // C++ gets nicer, we can do this in a nicer way.
252 : static_assert(Moose::always_false<T1>,
253 : "Harmonic interpolation is not implemented for the used type!");
254 : }
255 :
256 : /**
257 : * Linear interpolation with skewness correction using the face gradient.
258 : * See more info in Moukalled Chapter 9. The correction involves a first order
259 : * Taylor expansion around the intersection of the cell face and the line
260 : * connecting the two cell centers.
261 : */
262 : template <typename T, typename T2, typename T3>
263 : typename libMesh::CompareTypes<T, T2>::supertype
264 0 : skewCorrectedLinearInterpolation(const T & value1,
265 : const T2 & value2,
266 : const T3 & face_gradient,
267 : const FaceInfo & fi,
268 : const bool one_is_elem)
269 : {
270 0 : const auto coeffs = interpCoeffs(InterpMethod::SkewCorrectedAverage, fi, one_is_elem);
271 :
272 0 : auto value = (coeffs.first * value1 + coeffs.second * value2) +
273 0 : face_gradient * fi.skewnessCorrectionVector();
274 0 : return value;
275 : }
276 :
277 : /// Provides interpolation of face values for non-advection-specific purposes (although it can/will
278 : /// still be used by advective kernels sometimes). The interpolated value is stored in result.
279 : /// This should be called when a face value needs to be computed from two neighboring
280 : /// cells/elements. value1 and value2 represent the cell property/values from which to compute the
281 : /// face value. The \p one_is_elem parameter indicates whether value1 corresponds to the FaceInfo
282 : /// elem value; else it corresponds to the FaceInfo neighbor value
283 : template <typename T, typename T2, typename T3>
284 : void
285 8672922 : interpolate(InterpMethod m,
286 : T & result,
287 : const T2 & value1,
288 : const T3 & value2,
289 : const FaceInfo & fi,
290 : const bool one_is_elem)
291 : {
292 8672922 : switch (m)
293 : {
294 612002 : case InterpMethod::Average:
295 : case InterpMethod::SkewCorrectedAverage:
296 612002 : result = linearInterpolation(value1, value2, fi, one_is_elem);
297 612002 : break;
298 8060920 : case InterpMethod::HarmonicAverage:
299 8060920 : result = harmonicInterpolation(value1, value2, fi, one_is_elem);
300 8060920 : break;
301 0 : default:
302 0 : mooseError("unsupported interpolation method for interpolate() function");
303 : }
304 8672922 : }
305 :
306 : /**
307 : * perform a possibly skew-corrected linear interpolation by evaluating the supplied functor with
308 : * the provided functor face argument
309 : */
310 : template <typename T, FunctorEvaluationKind FEK = FunctorEvaluationKind::Value>
311 : T
312 63898131 : linearInterpolation(const FunctorBase<T> & functor, const FaceArg & face, const StateArg & time)
313 : {
314 : static_assert((FEK == FunctorEvaluationKind::Value) || (FEK == FunctorEvaluationKind::Dot),
315 : "Only Value and Dot evaluations currently supported");
316 : mooseAssert(face.limiter_type == LimiterType::CentralDifference,
317 : "this interpolation method is meant for linear interpolations");
318 :
319 : mooseAssert(face.fi,
320 : "We must have a non-null face_info in order to prepare our ElemFromFace tuples");
321 :
322 63898131 : constexpr FunctorEvaluationKind GFEK = FunctorGradientEvaluationKind<FEK>::value;
323 :
324 63898131 : const auto elem_arg = face.makeElem();
325 63898131 : const auto neighbor_arg = face.makeNeighbor();
326 :
327 63898131 : const auto elem_eval = functor.template genericEvaluate<FEK>(elem_arg, time);
328 63898131 : const auto neighbor_eval = functor.template genericEvaluate<FEK>(neighbor_arg, time);
329 :
330 63898131 : if (face.correct_skewness)
331 : {
332 : // This condition ensures that the recursive algorithm (face_center->
333 : // face_gradient -> cell_gradient -> face_center -> ...) terminates after
334 : // one loop. It is hardcoded to one loop at this point since it yields
335 : // 2nd order accuracy on skewed meshes with the minimum additional effort.
336 0 : FaceArg new_face(face);
337 0 : new_face.correct_skewness = false;
338 0 : const auto surface_gradient = functor.template genericEvaluate<GFEK>(new_face, time);
339 :
340 0 : return skewCorrectedLinearInterpolation(
341 0 : elem_eval, neighbor_eval, surface_gradient, *face.fi, true);
342 0 : }
343 : else
344 63898131 : return linearInterpolation(elem_eval, neighbor_eval, *face.fi, true);
345 63896871 : }
346 :
347 : /**
348 : * Computes the product of the advected and the advector based on the given interpolation method
349 : */
350 : template <typename T1,
351 : typename T2,
352 : typename T3,
353 : template <typename> class Vector1,
354 : template <typename> class Vector2>
355 : void
356 : interpolate(InterpMethod m,
357 : Vector1<T1> & result,
358 : const T2 & fi_elem_advected,
359 : const T2 & fi_neighbor_advected,
360 : const Vector2<T3> & fi_elem_advector,
361 : const Vector2<T3> & fi_neighbor_advector,
362 : const FaceInfo & fi)
363 : {
364 : switch (m)
365 : {
366 : case InterpMethod::Average:
367 : result = linearInterpolation(fi_elem_advected * fi_elem_advector,
368 : fi_neighbor_advected * fi_neighbor_advector,
369 : fi,
370 : true);
371 : break;
372 : case InterpMethod::Upwind:
373 : {
374 : const auto face_advector = linearInterpolation(MetaPhysicL::raw_value(fi_elem_advector),
375 : MetaPhysicL::raw_value(fi_neighbor_advector),
376 : fi,
377 : true);
378 : Real elem_coeff, neighbor_coeff;
379 : if (face_advector * fi.normal() > 0)
380 : elem_coeff = 1, neighbor_coeff = 0;
381 : else
382 : elem_coeff = 0, neighbor_coeff = 1;
383 :
384 : result = elem_coeff * fi_elem_advected * fi_elem_advector +
385 : neighbor_coeff * fi_neighbor_advected * fi_neighbor_advector;
386 : break;
387 : }
388 : default:
389 : mooseError("unsupported interpolation method for FVFaceInterface::interpolate");
390 : }
391 : }
392 :
393 : /// Provides interpolation of face values for advective flux kernels. This should be called by
394 : /// advective kernels when a face value is needed from two neighboring cells/elements. The
395 : /// interpolated value is stored in result. value1 and value2 represent the two neighboring advected
396 : /// cell property/values. advector represents the vector quantity at the face that is doing the
397 : /// advecting (e.g. the flow velocity at the face); this value often will have been computed using a
398 : /// call to the non-advective interpolate function. The \p one_is_elem parameter indicates whether
399 : /// value1 corresponds to the FaceInfo elem value; else it corresponds to the FaceInfo neighbor
400 : /// value
401 : template <typename T, typename T2, typename T3, typename Vector>
402 : void
403 : interpolate(InterpMethod m,
404 : T & result,
405 : const T2 & value1,
406 : const T3 & value2,
407 : const Vector & advector,
408 : const FaceInfo & fi,
409 : const bool one_is_elem)
410 : {
411 : const auto coeffs = interpCoeffs(m, fi, one_is_elem, advector * fi.normal());
412 : result = coeffs.first * value1 + coeffs.second * value2;
413 : }
414 :
415 : /// Calculates and returns "grad_u dot normal" on the face to be used for
416 : /// diffusive terms. If using any cross-diffusion corrections, etc. all
417 : /// those calculations should be handled appropriately by this function.
418 : ADReal gradUDotNormal(const FaceInfo & face_info,
419 : const MooseVariableFV<Real> & fv_var,
420 : const Moose::StateArg & time,
421 : bool correct_skewness = false);
422 :
423 : /**
424 : * From Moukalled 12.30
425 : *
426 : * r_f = (phiC - phiU) / (phiD - phiC)
427 : *
428 : * However, this formula is only clear on grids where upwind locations can be readily determined,
429 : * which is not the case for unstructured meshes. So we leverage a virtual upwind location and
430 : * Moukalled 12.65
431 : *
432 : * phiD - phiU = 2 * grad(phi)_C * dCD ->
433 : * phiU = phiD - 2 * grad(phi)_C * dCD
434 : *
435 : * Combining the two equations and performing some algebraic manipulation yields this equation for
436 : * r_f:
437 : *
438 : * r_f = 2 * grad(phi)_C * dCD / (phiD - phiC) - 1
439 : *
440 : * This equation is clearly asymmetric considering the face between C and D because of the
441 : * subscript on grad(phi). Hence this method can be thought of as constructing an r associated with
442 : * the C side of the face.
443 : *
444 : * This is a branchless version, avoiding if-else statements to enable vectorization.
445 : */
446 : template <typename Scalar, typename Vector>
447 : Scalar
448 21498062 : rF(const Scalar & phiC, const Scalar & phiD, const Vector & gradC, const RealVectorValue & dCD)
449 : {
450 21498062 : const Scalar denom = phiD - phiC;
451 21498062 : const Scalar grad_dot = gradC * dCD;
452 :
453 : // This is an arbitrary number here, when we start seeing convergence issues we
454 : // can tune this but so far this has shown okay results.
455 21498062 : constexpr Real eps = 1e-10;
456 : const Real denom_sign =
457 21498062 : std::copysign(1.0,
458 21498062 : MetaPhysicL::raw_value(denom) +
459 21498062 : std::numeric_limits<Real>::min() * MetaPhysicL::raw_value(grad_dot));
460 21498062 : const Scalar safe = denom + denom_sign * eps;
461 21746852 : return 2.0 * grad_dot / safe - 1.0;
462 248790 : }
463 :
464 : /**
465 : * Produce the interpolation coefficients in the equation:
466 : *
467 : * \phi_f = c_upwind * \phi_{upwind} + c_downwind * \phi_{downwind}
468 : *
469 : * A couple of examples: if we are doing an average interpolation with
470 : * an orthogonal regular grid, then the pair will be (0.5, 0.5). If we are doing an
471 : * upwind interpolation then the pair will be (1.0, 0.0).
472 : *
473 : * @return a pair where the first Real is c_upwind and the second Real is c_downwind
474 : */
475 : template <typename T>
476 : std::pair<T, T>
477 3048046 : interpCoeffs(const Limiter<T> & limiter,
478 : const T & phi_upwind,
479 : const T & phi_downwind,
480 : const libMesh::VectorValue<T> * const grad_phi_upwind,
481 : const libMesh::VectorValue<T> * const grad_phi_face,
482 : const Real & max_value,
483 : const Real & min_value,
484 : const FaceInfo & fi,
485 : const bool fi_elem_is_upwind)
486 : {
487 : // Using beta, w_f, g nomenclature from Greenshields
488 3048046 : const auto beta = limiter(phi_upwind,
489 : phi_downwind,
490 : grad_phi_upwind,
491 : grad_phi_face,
492 3048046 : fi_elem_is_upwind ? fi.dCN() : Point(-fi.dCN()),
493 : max_value,
494 : min_value,
495 : &fi,
496 : fi_elem_is_upwind);
497 :
498 3048046 : const auto w_f = fi_elem_is_upwind ? fi.gC() : (1. - fi.gC());
499 :
500 3048046 : const auto g = beta * (1. - w_f);
501 :
502 5663558 : return std::make_pair(1. - g, g);
503 2615512 : }
504 :
505 : /**
506 : * Interpolates with a limiter
507 : */
508 : template <typename Scalar,
509 : typename Vector,
510 : typename Enable = typename std::enable_if<libMesh::ScalarTraits<Scalar>::value>::type>
511 : Scalar
512 2026407 : interpolate(const Limiter<Scalar> & limiter,
513 : const Scalar & phi_upwind,
514 : const Scalar & phi_downwind,
515 : const Vector * const grad_phi_upwind,
516 : const FaceInfo & fi,
517 : const bool fi_elem_is_upwind)
518 : {
519 4052790 : auto pr =
520 24 : interpCoeffs(limiter,
521 : phi_upwind,
522 : phi_downwind,
523 : grad_phi_upwind,
524 : /*grad_phi_face*/ static_cast<const libMesh::VectorValue<Scalar> *>(nullptr),
525 2026407 : /* max_value */ std::numeric_limits<Real>::max(),
526 2026407 : /* min_value */ std::numeric_limits<Real>::min(),
527 : fi,
528 : fi_elem_is_upwind);
529 4052802 : return pr.first * phi_upwind + pr.second * phi_downwind;
530 2026395 : }
531 :
532 : /**
533 : * Vector overload
534 : */
535 : template <typename Limiter, typename T, typename Tensor>
536 : libMesh::VectorValue<T>
537 4 : interpolate(const Limiter & limiter,
538 : const TypeVector<T> & phi_upwind,
539 : const TypeVector<T> & phi_downwind,
540 : const Tensor * const grad_phi_upwind,
541 : const FaceInfo & fi,
542 : const bool fi_elem_is_upwind)
543 : {
544 : mooseAssert(limiter.constant() || grad_phi_upwind,
545 : "Non-null gradient only supported for constant limiters.");
546 :
547 4 : const libMesh::VectorValue<T> * const gradient_example = nullptr;
548 4 : libMesh::VectorValue<T> ret;
549 16 : for (const auto i : make_range(unsigned(LIBMESH_DIM)))
550 : {
551 12 : if (grad_phi_upwind)
552 : {
553 12 : const libMesh::VectorValue<T> gradient = grad_phi_upwind->row(i);
554 12 : ret(i) =
555 12 : interpolate(limiter, phi_upwind(i), phi_downwind(i), &gradient, fi, fi_elem_is_upwind);
556 : }
557 : else
558 0 : ret(i) = interpolate(
559 : limiter, phi_upwind(i), phi_downwind(i), gradient_example, fi, fi_elem_is_upwind);
560 : }
561 :
562 4 : return ret;
563 : }
564 :
565 : /**
566 : * This function performs a full limited interpolation of a variable, taking into account
567 : * the values and gradients at both upwind and downwind locations, as well as geometric
568 : * information and limits. It applies the specified limiter to ensure the interpolation
569 : * adheres to the constraints imposed by the limiter.
570 : *
571 : * @tparam T The data type of the scalar values and the return type.
572 : * @param limiter The limiter object used to compute the flux limiting ratio.
573 : * @param phi_upwind The field value at the upwind location.
574 : * @param phi_downwind The field value at the downwind location.
575 : * @param grad_phi_upwind Pointer to the gradient at the upwind location.
576 : * @param grad_phi_face Pointer to the gradient at the face location.
577 : * @param fi FaceInfo object containing geometric details such as face centroid and cell centroids.
578 : * @param fi_elem_is_upwind Boolean indicating if the face info element is upwind.
579 : * @param max_value The maximum allowable value.
580 : * @param min_value The minimum allowable value.
581 : * @return The computed limited interpolation value.
582 : */
583 : template <typename T>
584 : T
585 4717997 : fullLimitedInterpolation(const Limiter<T> & limiter,
586 : const T & phi_upwind,
587 : const T & phi_downwind,
588 : const VectorValue<T> * const grad_phi_upwind,
589 : const VectorValue<T> * const grad_phi_face,
590 : const Real & max_value,
591 : const Real & min_value,
592 : const FaceArg & face)
593 : {
594 : // Compute the direction vector based on whether the current element is upwind
595 4717997 : const auto dCD = face.elem_is_upwind ? face.fi->dCN() : Point(-face.fi->dCN());
596 :
597 : // Compute the flux limiting ratio using the specified limiter
598 4717997 : const auto beta = limiter(phi_upwind,
599 : phi_downwind,
600 : grad_phi_upwind,
601 : grad_phi_face,
602 : dCD,
603 : max_value,
604 : min_value,
605 4717997 : face.fi,
606 4717997 : face.elem_is_upwind);
607 :
608 : // Determine the face centroid and the appropriate cell centroid
609 4717997 : const auto & face_centroid = face.fi->faceCentroid();
610 4717997 : const auto & cell_centroid =
611 4717997 : face.elem_is_upwind ? face.fi->elemCentroid() : face.fi->neighborCentroid();
612 :
613 : // Compute the delta value at the face
614 4717997 : T delta_face;
615 4717997 : if (face.limiter_type == LimiterType::Venkatakrishnan || face.limiter_type == LimiterType::SOU)
616 2636459 : delta_face = (*grad_phi_upwind) * (face_centroid - cell_centroid);
617 : else
618 2081538 : delta_face = (*grad_phi_face) * (face_centroid - cell_centroid);
619 :
620 : // Return the interpolated value
621 9435994 : return phi_upwind + beta * delta_face;
622 4717997 : }
623 :
624 : /**
625 : * This function calculates the minimum and maximum values within a two-cell stencil. The stencil
626 : * includes the immediate neighboring elements of the face's associated element and the neighboring
627 : * elements of those neighbors. It evaluates the values using a provided functor and accounts for
628 : * the valid (non-null) neighbors.
629 : *
630 : * @tparam T The data type for the values being computed. This is typically a scalar type.
631 : * @tparam FEK Enumeration type FunctorEvaluationKind with a default value of
632 : * FunctorEvaluationKind::Value. This determines the kind of evaluation the functor will perform.
633 : * @tparam Enable A type trait used for SFINAE (Substitution Failure Is Not An Error) to ensure that
634 : * T is a valid scalar type as determined by ScalarTraits<T>.
635 : *
636 : * @param functor An object of a functor class derived from FunctorBase<T>. This object provides the
637 : * genericEvaluate method to compute the value at a given element and time.
638 : * @param face An argument representing the face in the computational domain. The face provides
639 : * access to neighboring elements via neighbor_ptr_range().
640 : * @param time An argument representing the state or time at which the evaluation is performed.
641 : *
642 : * @return std::pair<T, T> A pair containing the minimum and maximum values computed across the
643 : * two-cell stencil. The first element is the maximum value, and the second element is the minimum
644 : * value.
645 : */
646 : template <typename T,
647 : FunctorEvaluationKind FEK = FunctorEvaluationKind::Value,
648 : typename Enable = typename std::enable_if<ScalarTraits<T>::value>::type>
649 : std::pair<Real, Real>
650 2636459 : computeMinMaxValue(const FunctorBase<T> & functor, const FaceArg & face, const StateArg & time)
651 : {
652 : // Initialize max_value to 0 and min_value to a large value
653 2636459 : Real max_value(std::numeric_limits<Real>::min()), min_value(std::numeric_limits<Real>::max());
654 :
655 : // Iterate over the direct neighbors of the element associated with the face
656 13163857 : for (const auto neighbor : (*face.fi).elem().neighbor_ptr_range())
657 : {
658 : // If not a valid neighbor, skip to the next one
659 10527398 : if (neighbor == nullptr)
660 743991 : continue;
661 :
662 : // Evaluate the functor at the neighbor and update max_value and min_value
663 9783407 : const ElemArg local_cell_arg = {neighbor, false};
664 9783407 : const auto value =
665 9783407 : MetaPhysicL::raw_value(functor.template genericEvaluate<FEK>(local_cell_arg, time));
666 9783407 : max_value = std::max(max_value, value);
667 9783407 : min_value = std::min(min_value, value);
668 : }
669 :
670 : // Iterate over the neighbors of the neighbor
671 13163857 : for (const auto neighbor : (*face.fi).neighbor().neighbor_ptr_range())
672 : {
673 : // If not a valid neighbor, skip to the next one
674 10527398 : if (neighbor == nullptr)
675 737183 : continue;
676 :
677 : // Evaluate the functor at the neighbor and update max_value and min_value
678 9790215 : const ElemArg local_cell_arg = {neighbor, false};
679 9790215 : const auto value =
680 9790215 : MetaPhysicL::raw_value(functor.template genericEvaluate<FEK>(local_cell_arg, time));
681 9790215 : max_value = std::max(max_value, value);
682 9790215 : min_value = std::min(min_value, value);
683 : }
684 :
685 : // Return a pair containing the computed maximum and minimum values
686 5272918 : return std::make_pair(max_value, min_value);
687 : }
688 :
689 : /**
690 : * This function calculates the minimum and maximum values of a specified component within a
691 : * two-cell stencil. The stencil includes the immediate neighboring elements of the face's
692 : * associated element and the neighboring elements of those neighbors. It evaluates the values using
693 : * a provided functor and accounts for the valid (non-null) neighbors.
694 : *
695 : * @tparam T The data type for the values being computed. This is typically a scalar type.
696 : *
697 : * @param functor An object of a functor class derived from FunctorBase<VectorValue<T>>. This object
698 : * provides the operator() method to compute the value at a given element and time.
699 : * @param face An argument representing the face in the computational domain. The face provides
700 : * access to neighboring elements via neighbor_ptr_range().
701 : * @param time An argument representing the state or time at which the evaluation is performed.
702 : * @param component An unsigned integer representing the specific component of the vector to be
703 : * evaluated.
704 : *
705 : * @return std::pair<T, T> A pair containing the minimum and maximum values of the specified
706 : * component computed across the two-cell stencil. The first element is the maximum value, and the
707 : * second element is the minimum value.
708 : *
709 : * Usage:
710 : * This function is typically used in the finite volume methods for min-max computations over
711 : * stencils (neighborhoods). It helps compute the limiting for limited second order upwind at the
712 : * faces.
713 : */
714 : template <typename T>
715 : std::pair<Real, Real>
716 0 : computeMinMaxValue(const FunctorBase<VectorValue<T>> & functor,
717 : const FaceArg & face,
718 : const StateArg & time,
719 : const unsigned int & component)
720 : {
721 : // Initialize max_value to 0 and min_value to a large value
722 0 : Real max_value(std::numeric_limits<Real>::min()), min_value(std::numeric_limits<Real>::max());
723 :
724 : // Iterate over the direct neighbors of the element associated with the face
725 0 : for (const auto neighbor : (*face.fi).elem().neighbor_ptr_range())
726 : {
727 : // If not a valid neighbor, skip to the next one
728 0 : if (neighbor == nullptr)
729 0 : continue;
730 :
731 : // Evaluate the functor at the neighbor for the specified component and update max_value and
732 : // min_value
733 0 : const ElemArg local_cell_arg = {neighbor, false};
734 0 : const auto value = MetaPhysicL::raw_value(functor(local_cell_arg, time)(component));
735 0 : max_value = std::max(max_value, value);
736 0 : min_value = std::min(min_value, value);
737 : }
738 :
739 : // Iterate over the neighbors of the neighbor associated with the face
740 0 : for (const auto neighbor : (*face.fi).neighbor().neighbor_ptr_range())
741 : {
742 : // If not a valid neighbor, skip to the next one
743 0 : if (neighbor == nullptr)
744 0 : continue;
745 :
746 : // Evaluate the functor at the neighbor for the specified component and update max_value and
747 : // min_value
748 0 : const ElemArg local_cell_arg = {neighbor, false};
749 0 : const auto value = MetaPhysicL::raw_value(functor(local_cell_arg, time)(component));
750 0 : max_value = std::max(max_value, value);
751 0 : min_value = std::min(min_value, value);
752 : }
753 :
754 : // Return a pair containing the computed maximum and minimum values
755 0 : return std::make_pair(max_value, min_value);
756 : }
757 :
758 : /**
759 : * This function interpolates values using a specified limiter and face argument. It evaluates the
760 : * values at upwind and downwind locations and computes interpolation coefficients and advected
761 : * values.
762 : *
763 : * @tparam T The data type for the values being interpolated. This is typically a scalar type.
764 : * @tparam FEK Enumeration type FunctorEvaluationKind with a default value of
765 : * FunctorEvaluationKind::Value. This determines the kind of evaluation the functor will perform.
766 : * @tparam Enable A type trait used for SFINAE (Substitution Failure Is Not An Error) to ensure that
767 : * T is a valid scalar type as determined by ScalarTraits<T>.
768 : *
769 : * @param functor An object of a functor class derived from FunctorBase<T>. This object provides the
770 : * genericEvaluate method to compute the value at a given element and time.
771 : * @param face An argument representing the face in the computational domain. The face provides
772 : * access to neighboring elements and limiter type.
773 : * @param time An argument representing the state or time at which the evaluation is performed.
774 : *
775 : * @return std::pair<std::pair<T, T>, std::pair<T, T>> A pair of pairs.
776 : * - The first pair corresponds to the interpolation coefficients, with the
777 : * first value corresponding to the face information element and the second to the face information
778 : * neighbor. This pair should sum to unity.
779 : * - The second pair corresponds to the face information functor element
780 : * value and neighbor value.
781 : *
782 : * Usage:
783 : * This function is used for interpolating values at faces in a finite volume method, ensuring that
784 : * the interpolation adheres to the constraints imposed by the limiter.
785 : */
786 : template <typename T,
787 : FunctorEvaluationKind FEK = FunctorEvaluationKind::Value,
788 : typename Enable = typename std::enable_if<libMesh::ScalarTraits<T>::value>::type>
789 : std::pair<std::pair<T, T>, std::pair<T, T>>
790 370435 : interpCoeffsAndAdvected(const FunctorBase<T> & functor, const FaceArg & face, const StateArg & time)
791 : {
792 : // Ensure only supported FunctorEvaluationKinds are used
793 : static_assert((FEK == FunctorEvaluationKind::Value) || (FEK == FunctorEvaluationKind::Dot),
794 : "Only Value and Dot evaluations currently supported");
795 :
796 : // Determine the gradient evaluation kind
797 370435 : constexpr FunctorEvaluationKind GFEK = FunctorGradientEvaluationKind<FEK>::value;
798 : typedef typename FunctorBase<T>::GradientType GradientType;
799 370435 : static const GradientType zero(0);
800 :
801 : mooseAssert(face.fi, "this must be non-null");
802 :
803 : // Construct the limiter based on the face limiter type
804 370435 : const auto limiter = Limiter<typename LimiterValueType<T>::value_type>::build(face.limiter_type);
805 :
806 : // Determine upwind and downwind arguments based on the face element
807 370435 : const auto upwind_arg = face.elem_is_upwind ? face.makeElem() : face.makeNeighbor();
808 370435 : const auto downwind_arg = face.elem_is_upwind ? face.makeNeighbor() : face.makeElem();
809 :
810 : // Evaluate the functor at the upwind and downwind locations
811 370435 : auto phi_upwind = functor.template genericEvaluate<FEK>(upwind_arg, time);
812 370435 : auto phi_downwind = functor.template genericEvaluate<FEK>(downwind_arg, time);
813 :
814 : // Initialize the interpolation coefficients
815 370435 : std::pair<T, T> interp_coeffs;
816 :
817 : // Compute interpolation coefficients for Upwind or CentralDifference limiters
818 370435 : if (face.limiter_type == LimiterType::Upwind ||
819 0 : face.limiter_type == LimiterType::CentralDifference)
820 370435 : interp_coeffs = interpCoeffs(*limiter,
821 : phi_upwind,
822 : phi_downwind,
823 : &zero,
824 : &zero,
825 370435 : std::numeric_limits<Real>::max(),
826 740870 : std::numeric_limits<Real>::min(),
827 370435 : *face.fi,
828 370435 : face.elem_is_upwind);
829 : else
830 : {
831 : // Determine the time argument for the limiter
832 0 : auto * time_arg = face.state_limiter;
833 0 : if (!time_arg)
834 : {
835 0 : static Moose::StateArg temp_state(0, Moose::SolutionIterationType::Time);
836 0 : time_arg = &temp_state;
837 : }
838 :
839 0 : Real max_value(std::numeric_limits<Real>::min()), min_value(std::numeric_limits<Real>::max());
840 :
841 : // Compute min-max values for min-max limiters
842 0 : if (face.limiter_type == LimiterType::Venkatakrishnan || face.limiter_type == LimiterType::SOU)
843 0 : std::tie(max_value, min_value) = computeMinMaxValue(functor, face, *time_arg);
844 :
845 : // Evaluate the gradient of the functor at the upwind and downwind locations
846 0 : const auto grad_phi_upwind = functor.template genericEvaluate<GFEK>(upwind_arg, *time_arg);
847 0 : const auto grad_phi_face = functor.template genericEvaluate<GFEK>(face, *time_arg);
848 :
849 : // Compute the interpolation coefficients using the specified limiter
850 0 : interp_coeffs = interpCoeffs(*limiter,
851 : phi_upwind,
852 : phi_downwind,
853 : &grad_phi_upwind,
854 : &grad_phi_face,
855 : max_value,
856 : min_value,
857 0 : *face.fi,
858 0 : face.elem_is_upwind);
859 0 : }
860 :
861 : // Return the interpolation coefficients and advected values
862 370435 : if (face.elem_is_upwind)
863 366307 : return std::make_pair(std::move(interp_coeffs),
864 732614 : std::make_pair(std::move(phi_upwind), std::move(phi_downwind)));
865 : else
866 0 : return std::make_pair(
867 4128 : std::make_pair(std::move(interp_coeffs.second), std::move(interp_coeffs.first)),
868 8256 : std::make_pair(std::move(phi_downwind), std::move(phi_upwind)));
869 370435 : }
870 :
871 : /**
872 : * This function interpolates values at faces in a computational grid using a specified functor,
873 : * face argument, and evaluation kind. It handles different limiter types and performs
874 : * interpolation accordingly.
875 : *
876 : * @tparam T The data type for the values being interpolated. This is typically a scalar type.
877 : * @tparam FEK Enumeration type FunctorEvaluationKind with a default value of
878 : * FunctorEvaluationKind::Value. This determines the kind of evaluation the functor will perform.
879 : * @tparam Enable A type trait used for SFINAE (Substitution Failure Is Not An Error) to ensure that
880 : * T is a valid scalar type as determined by ScalarTraits<T>.
881 : *
882 : * @param functor An object of a functor class derived from FunctorBase<T>. This object provides the
883 : * genericEvaluate method to compute the value at a given element and time.
884 : * @param face An argument representing the face in the computational domain. The face provides
885 : * access to neighboring elements and limiter type.
886 : * @param time An argument representing the state or time at which the evaluation is performed.
887 : *
888 : * @return T The interpolated value at the face.
889 : *
890 : * Usage:
891 : * This function is used for interpolating values at faces in a finite volume method, ensuring that
892 : * the interpolation adheres to the constraints imposed by the limiter.
893 : */
894 : template <typename T,
895 : FunctorEvaluationKind FEK = FunctorEvaluationKind::Value,
896 : typename Enable = typename std::enable_if<libMesh::ScalarTraits<T>::value>::type>
897 : T
898 68976987 : interpolate(const FunctorBase<T> & functor, const FaceArg & face, const StateArg & time)
899 : {
900 : // Ensure only supported FunctorEvaluationKinds are used
901 : static_assert((FEK == FunctorEvaluationKind::Value) || (FEK == FunctorEvaluationKind::Dot),
902 : "Only Value and Dot evaluations currently supported");
903 :
904 : // Special handling for central differencing as it is the only interpolation method which
905 : // currently supports skew correction
906 68976987 : if (face.limiter_type == LimiterType::CentralDifference)
907 63898131 : return linearInterpolation<T, FEK>(functor, face, time);
908 :
909 5078856 : if (face.limiter_type == LimiterType::Upwind ||
910 4708421 : face.limiter_type == LimiterType::CentralDifference)
911 : {
912 : // Compute interpolation coefficients and advected values
913 370435 : const auto [interp_coeffs, advected] = interpCoeffsAndAdvected<T, FEK>(functor, face, time);
914 : // Return the interpolated value
915 370435 : return interp_coeffs.first * advected.first + interp_coeffs.second * advected.second;
916 370435 : }
917 : else
918 : {
919 : // Determine the gradient evaluation kind
920 4708421 : constexpr FunctorEvaluationKind GFEK = FunctorGradientEvaluationKind<FEK>::value;
921 : typedef typename FunctorBase<T>::GradientType GradientType;
922 4708421 : static const GradientType zero(0);
923 9416842 : const auto & limiter =
924 4708421 : Limiter<typename LimiterValueType<T>::value_type>::build(face.limiter_type);
925 :
926 : // Determine upwind and downwind arguments based on the face element
927 4708421 : const auto & upwind_arg = face.elem_is_upwind ? face.makeElem() : face.makeNeighbor();
928 4708421 : const auto & downwind_arg = face.elem_is_upwind ? face.makeNeighbor() : face.makeElem();
929 4708421 : const auto & phi_upwind = functor.template genericEvaluate<FEK>(upwind_arg, time);
930 4708421 : const auto & phi_downwind = functor.template genericEvaluate<FEK>(downwind_arg, time);
931 :
932 : // Determine the time argument for the limiter
933 4708421 : auto * time_arg = face.state_limiter;
934 4708421 : if (!time_arg)
935 : {
936 35001 : static Moose::StateArg temp_state(0, Moose::SolutionIterationType::Time);
937 35001 : time_arg = &temp_state;
938 : }
939 :
940 : // Initialize min-max values
941 4708421 : Real max_value(std::numeric_limits<Real>::min()), min_value(std::numeric_limits<Real>::max());
942 4708421 : if (face.limiter_type == LimiterType::Venkatakrishnan || face.limiter_type == LimiterType::SOU)
943 2636459 : std::tie(max_value, min_value) = computeMinMaxValue(functor, face, *time_arg);
944 :
945 : // Evaluate the gradient of the functor at the upwind and downwind locations
946 4708421 : const auto & grad_phi_upwind = functor.template genericEvaluate<GFEK>(upwind_arg, *time_arg);
947 4708421 : const auto & grad_phi_face = functor.template genericEvaluate<GFEK>(face, *time_arg);
948 :
949 : // Perform full limited interpolation and return the interpolated value
950 4708421 : return fullLimitedInterpolation(*limiter,
951 : phi_upwind,
952 : phi_downwind,
953 : &grad_phi_upwind,
954 : &grad_phi_face,
955 : max_value,
956 : min_value,
957 4708421 : face);
958 4708421 : }
959 : }
960 :
961 : /**
962 : * This function interpolates vector values at faces in a computational grid using a specified
963 : * functor, face argument, and limiter type. It handles different limiter types and performs
964 : * interpolation accordingly.
965 : *
966 : * @tparam T The data type for the vector values being interpolated. This is typically a scalar
967 : * type.
968 : *
969 : * @param functor An object of a functor class derived from FunctorBase<VectorValue<T>>. This object
970 : * provides the operator() method to compute the value at a given element and time.
971 : * @param face An argument representing the face in the computational domain. The face provides
972 : * access to neighboring elements and limiter type.
973 : * @param time An argument representing the state or time at which the evaluation is performed.
974 : *
975 : * @return VectorValue<T> The interpolated vector value at the face.
976 : *
977 : * Usage:
978 : * This function is used for interpolating vector values at faces in a finite volume method,
979 : * ensuring that the interpolation adheres to the constraints imposed by the limiter.
980 : */
981 : template <typename T>
982 : libMesh::VectorValue<T>
983 220236 : interpolate(const FunctorBase<libMesh::VectorValue<T>> & functor,
984 : const FaceArg & face,
985 : const StateArg & time)
986 : {
987 : // Define a zero gradient vector for initialization
988 220236 : static const libMesh::VectorValue<T> grad_zero(0);
989 :
990 : mooseAssert(face.fi, "this must be non-null");
991 :
992 : // Construct the limiter based on the face limiter type
993 220236 : const auto limiter = Limiter<typename LimiterValueType<T>::value_type>::build(face.limiter_type);
994 :
995 : // Determine upwind and downwind arguments based on the face element
996 220236 : const auto upwind_arg = face.elem_is_upwind ? face.makeElem() : face.makeNeighbor();
997 220236 : const auto downwind_arg = face.elem_is_upwind ? face.makeNeighbor() : face.makeElem();
998 220236 : auto phi_upwind = functor(upwind_arg, time);
999 220236 : auto phi_downwind = functor(downwind_arg, time);
1000 :
1001 : // Initialize the return vector value
1002 220236 : libMesh::VectorValue<T> ret;
1003 76086 : T coeff_upwind, coeff_downwind;
1004 :
1005 : // Compute interpolation coefficients and advected values for Upwind or CentralDifference limiters
1006 220236 : if (face.limiter_type == LimiterType::Upwind ||
1007 218868 : face.limiter_type == LimiterType::CentralDifference)
1008 : {
1009 868176 : for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
1010 : {
1011 651132 : const auto &component_upwind = phi_upwind(i), component_downwind = phi_downwind(i);
1012 651132 : std::tie(coeff_upwind, coeff_downwind) = interpCoeffs(*limiter,
1013 : component_upwind,
1014 : component_downwind,
1015 : &grad_zero,
1016 : &grad_zero,
1017 651132 : std::numeric_limits<Real>::max(),
1018 651132 : std::numeric_limits<Real>::min(),
1019 651132 : *face.fi,
1020 651132 : face.elem_is_upwind);
1021 651132 : ret(i) = coeff_upwind * component_upwind + coeff_downwind * component_downwind;
1022 : }
1023 217044 : }
1024 : else
1025 : {
1026 : // Determine the time argument for the limiter
1027 3192 : auto * time_arg = face.state_limiter;
1028 3192 : if (!time_arg)
1029 : {
1030 3192 : static Moose::StateArg temp_state(0, Moose::SolutionIterationType::Time);
1031 3192 : time_arg = &temp_state;
1032 : }
1033 :
1034 : // Evaluate the gradient of the functor at the upwind and downwind locations
1035 3192 : const auto & grad_phi_upwind = functor.gradient(upwind_arg, *time_arg);
1036 3192 : const auto & grad_phi_downwind = functor.gradient(downwind_arg, *time_arg);
1037 :
1038 3192 : const auto coeffs = interpCoeffs(InterpMethod::Average, *face.fi, face.elem_is_upwind);
1039 :
1040 : // Compute interpolation coefficients and advected values for each component
1041 12768 : for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
1042 : {
1043 9576 : const auto &component_upwind = phi_upwind(i), component_downwind = phi_downwind(i);
1044 9576 : const libMesh::VectorValue<T> &grad_upwind = grad_phi_upwind.row(i),
1045 9576 : grad_face = coeffs.first * grad_phi_upwind.row(i) +
1046 0 : coeffs.second * grad_phi_downwind.row(i);
1047 :
1048 : // Initialize min-max values
1049 9576 : Real max_value(std::numeric_limits<Real>::min()), min_value(std::numeric_limits<Real>::max());
1050 9576 : if (face.limiter_type == LimiterType::Venkatakrishnan ||
1051 9576 : face.limiter_type == LimiterType::SOU)
1052 0 : std::tie(max_value, min_value) = computeMinMaxValue(functor, face, *time_arg, i);
1053 :
1054 : // Perform full limited interpolation for the component
1055 9576 : ret(i) = fullLimitedInterpolation(*limiter,
1056 : component_upwind,
1057 : component_downwind,
1058 : &grad_upwind,
1059 : &grad_face,
1060 : max_value,
1061 : min_value,
1062 : face);
1063 : }
1064 3192 : }
1065 :
1066 : // Return the interpolated vector value
1067 440472 : return ret;
1068 220236 : }
1069 :
1070 : /**
1071 : * This function interpolates container values at faces in a computational grid using a specified
1072 : * functor, face argument, and limiter type. It handles different limiter types and performs
1073 : * interpolation accordingly.
1074 : *
1075 : * @tparam T The data type for the container values being interpolated. This is typically a
1076 : * container type such as a vector or array.
1077 : *
1078 : * @param functor An object of a functor class derived from FunctorBase<T>. This object provides the
1079 : * operator() method to compute the value at a given element and time.
1080 : * @param face An argument representing the face in the computational domain. The face provides
1081 : * access to neighboring elements and limiter type.
1082 : * @param time An argument representing the state or time at which the evaluation is performed.
1083 : *
1084 : * @return T The interpolated container value at the face.
1085 : *
1086 : * Usage:
1087 : * This function is used for interpolating container values at faces in a finite volume method,
1088 : * ensuring that the interpolation adheres to the constraints imposed by the limiter.
1089 : */
1090 : template <typename T>
1091 : T
1092 72 : containerInterpolate(const FunctorBase<T> & functor, const FaceArg & face, const StateArg & time)
1093 : {
1094 : typedef typename FunctorBase<T>::GradientType ContainerGradientType;
1095 : typedef typename ContainerGradientType::value_type GradientType;
1096 72 : const GradientType * const example_gradient = nullptr;
1097 :
1098 : mooseAssert(face.fi, "this must be non-null");
1099 72 : const auto limiter = Limiter<typename T::value_type>::build(face.limiter_type);
1100 :
1101 72 : const auto upwind_arg = face.elem_is_upwind ? face.makeElem() : face.makeNeighbor();
1102 72 : const auto downwind_arg = face.elem_is_upwind ? face.makeNeighbor() : face.makeElem();
1103 72 : const auto phi_upwind = functor(upwind_arg, time);
1104 72 : const auto phi_downwind = functor(downwind_arg, time);
1105 :
1106 : // initialize in order to get proper size
1107 72 : T ret = phi_upwind;
1108 : typename T::value_type coeff_upwind, coeff_downwind;
1109 :
1110 72 : if (face.limiter_type == LimiterType::Upwind ||
1111 48 : face.limiter_type == LimiterType::CentralDifference)
1112 : {
1113 96 : for (const auto i : index_range(ret))
1114 : {
1115 48 : const auto &component_upwind = phi_upwind[i], component_downwind = phi_downwind[i];
1116 48 : std::tie(coeff_upwind, coeff_downwind) = interpCoeffs(*limiter,
1117 : component_upwind,
1118 : component_downwind,
1119 : example_gradient,
1120 : example_gradient,
1121 48 : std::numeric_limits<Real>::max(),
1122 48 : std::numeric_limits<Real>::min(),
1123 48 : *face.fi,
1124 48 : face.elem_is_upwind);
1125 48 : ret[i] = coeff_upwind * component_upwind + coeff_downwind * component_downwind;
1126 : }
1127 48 : }
1128 : else
1129 : {
1130 24 : const auto grad_phi_upwind = functor.gradient(upwind_arg, time);
1131 48 : for (const auto i : index_range(ret))
1132 : {
1133 24 : const auto &component_upwind = phi_upwind[i], component_downwind = phi_downwind[i];
1134 24 : const auto & grad = grad_phi_upwind[i];
1135 24 : std::tie(coeff_upwind, coeff_downwind) = interpCoeffs(*limiter,
1136 : component_upwind,
1137 : component_downwind,
1138 : &grad,
1139 : example_gradient,
1140 24 : std::numeric_limits<Real>::max(),
1141 24 : std::numeric_limits<Real>::min(),
1142 24 : *face.fi,
1143 24 : face.elem_is_upwind);
1144 24 : ret[i] = coeff_upwind * component_upwind + coeff_downwind * component_downwind;
1145 : }
1146 16 : }
1147 :
1148 120 : return ret;
1149 72 : }
1150 :
1151 : template <typename T>
1152 : std::vector<T>
1153 48 : interpolate(const FunctorBase<std::vector<T>> & functor,
1154 : const FaceArg & face,
1155 : const StateArg & time)
1156 : {
1157 48 : return containerInterpolate(functor, face, time);
1158 : }
1159 :
1160 : template <typename T, std::size_t N>
1161 : std::array<T, N>
1162 24 : interpolate(const FunctorBase<std::array<T, N>> & functor,
1163 : const FaceArg & face,
1164 : const StateArg & time)
1165 : {
1166 24 : return containerInterpolate(functor, face, time);
1167 : }
1168 :
1169 : /**
1170 : * Return whether the supplied face is on a boundary of the \p object's execution
1171 : */
1172 : template <typename SubdomainRestrictable>
1173 : bool
1174 18633415 : onBoundary(const SubdomainRestrictable & obj, const FaceInfo & fi)
1175 : {
1176 35622113 : const bool internal = fi.neighborPtr() && obj.hasBlocks(fi.elemSubdomainID()) &&
1177 16988698 : obj.hasBlocks(fi.neighborSubdomainID());
1178 18633415 : return !internal;
1179 : }
1180 :
1181 : /**
1182 : * Determine whether the passed-in face is on the boundary of an object that lives on the provided
1183 : * subdomains. Note that if the subdomain set is empty we consider that to mean that the object has
1184 : * no block restriction and lives on the entire mesh
1185 : */
1186 : bool onBoundary(const std::set<SubdomainID> & subs, const FaceInfo & fi);
1187 : }
1188 : }
|
