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