https://mooseframework.inl.gov
Classes | Enumerations | Functions | Variables
Moose::FV Namespace Reference

Classes

class  CentralDifferenceLimiter
 Implements a limiter which reproduces a central-differencing scheme, defined by $(r_f) = 1$. More...
 
class  Limiter
 Base class for defining slope limiters for finite volume or potentially reconstructed Discontinuous-Galerkin applications. More...
 
struct  LimiterValueType
 
struct  LimiterValueType< ADReal >
 
struct  LimiterValueType< Real >
 
struct  LimiterValueType< T, typename std::enable_if< HasMemberType_value_type< T >::value >::type >
 
class  MinModLimiter
 The Min-Mod limiter function $(r_f)$ is defined as: More...
 
class  QUICKLimiter
 The QUICK (Quadratic Upstream Interpolation for Convective Kinematics) limiter function is derived from the following equations: More...
 
class  SOULimiter
 The SOU limiter is used for reproducing the second-order-upwind scheme. More...
 
class  UpwindLimiter
 Implements a limiter which reproduces the upwind scheme, defined by $(r_f) = 0$. More...
 
class  VanLeerLimiter
 The Van Leer limiter limiter function $(r_f)$ is defined as: More...
 
class  VenkatakrishnanLimiter
 The Venkatakrishnan limiter is derived from the following equations: More...
 

Enumerations

enum  LimiterType : int {
  LimiterType::VanLeer = 0, LimiterType::Upwind, LimiterType::CentralDifference, LimiterType::MinMod,
  LimiterType::SOU, LimiterType::QUICK, LimiterType::Venkatakrishnan
}
 
enum  InterpMethod {
  InterpMethod::Average, InterpMethod::HarmonicAverage, InterpMethod::SkewCorrectedAverage, InterpMethod::Upwind,
  InterpMethod::RhieChow, InterpMethod::VanLeer, InterpMethod::MinMod, InterpMethod::SOU,
  InterpMethod::QUICK, InterpMethod::Venkatakrishnan
}
 This codifies a set of available ways to interpolate with elem+neighbor solution information to calculate values (e.g. More...
 
enum  LinearFVComputationMode { LinearFVComputationMode::RHS, LinearFVComputationMode::Matrix, LinearFVComputationMode::FullSystem }
 

Functions

LimiterType limiterType (InterpMethod interp_method)
 Return the limiter type associated with the supplied interpolation method. More...
 
bool elemHasFaceInfo (const Elem &elem, const Elem *const neighbor)
 This function infers based on elements if the faceinfo between them belongs to the element or not. More...
 
template<typename ActionFunctor >
void loopOverElemFaceInfo (const Elem &elem, const MooseMesh &mesh, ActionFunctor &act, const Moose::CoordinateSystemType coord_type, const unsigned int rz_radial_coord=libMesh::invalid_uint)
 
template<typename FVVar >
std::tuple< const Elem *, const Elem *, bool > determineElemOneAndTwo (const FaceInfo &fi, const FVVar &var)
 This utility determines element one and element two given a FaceInfo fi and variable var. More...
 
template<typename T , typename Enable = typename std::enable_if<libMesh::ScalarTraits<T>::value>::type>
libMesh::VectorValue< T > greenGaussGradient (const ElemArg &elem_arg, const StateArg &state_arg, const FunctorBase< T > &functor, const bool two_term_boundary_expansion, const MooseMesh &mesh, const bool force_green_gauss=false)
 Compute a cell gradient using the method of Green-Gauss. More...
 
template<typename T , typename Enable = typename std::enable_if<libMesh::ScalarTraits<T>::value>::type>
libMesh::VectorValue< T > greenGaussGradient (const FaceArg &face_arg, const StateArg &state_arg, const FunctorBase< T > &functor, const bool two_term_boundary_expansion, const MooseMesh &mesh)
 Compute a face gradient from Green-Gauss cell gradients, with orthogonality correction On the boundaries, the boundary element value is used. More...
 
template<typename T >
TensorValue< T > greenGaussGradient (const ElemArg &elem_arg, const StateArg &state_arg, const Moose::FunctorBase< libMesh::VectorValue< T >> &functor, const bool two_term_boundary_expansion, const MooseMesh &mesh)
 
template<typename T >
TensorValue< T > greenGaussGradient (const FaceArg &face_arg, const StateArg &state_arg, const Moose::FunctorBase< libMesh::VectorValue< T >> &functor, const bool two_term_boundary_expansion, const MooseMesh &mesh)
 
template<typename T >
Moose::FunctorBase< std::vector< T > >::GradientType greenGaussGradient (const ElemArg &elem_arg, const StateArg &state_arg, const Moose::FunctorBase< std::vector< T >> &functor, const bool two_term_boundary_expansion, const MooseMesh &mesh)
 
template<typename T >
Moose::FunctorBase< std::vector< T > >::GradientType greenGaussGradient (const FaceArg &face_arg, const StateArg &state_arg, const Moose::FunctorBase< std::vector< T >> &functor, const bool two_term_boundary_expansion, const MooseMesh &mesh)
 
template<typename T , std::size_t N>
Moose::FunctorBase< std::array< T, N > >::GradientType greenGaussGradient (const ElemArg &elem_arg, const StateArg &state_arg, const Moose::FunctorBase< std::array< T, N >> &functor, const bool two_term_boundary_expansion, const MooseMesh &mesh)
 
template<typename T , std::size_t N>
Moose::FunctorBase< std::array< T, N > >::GradientType greenGaussGradient (const FaceArg &face_arg, const StateArg &state_arg, const Moose::FunctorBase< std::array< T, N >> &functor, const bool two_term_boundary_expansion, const MooseMesh &mesh)
 
MooseEnum interpolationMethods ()
 Returns an enum with all the currently supported interpolation methods and the current default for FV: first-order upwind. More...
 
InputParameters advectedInterpolationParameter ()
 
InterpMethod selectInterpolationMethod (const std::string &interp_method)
 
bool setInterpolationMethod (const MooseObject &obj, Moose::FV::InterpMethod &interp_method, const std::string &param_name)
 Sets one interpolation method. More...
 
template<typename T = Real>
std::pair< Real, RealinterpCoeffs (const InterpMethod m, const FaceInfo &fi, const bool one_is_elem, const T &face_flux=0.0)
 Produce the interpolation coefficients in the equation: More...
 
template<typename T , typename T2 >
libMesh::CompareTypes< T, T2 >::supertype linearInterpolation (const T &value1, const T2 &value2, const FaceInfo &fi, const bool one_is_elem, const InterpMethod interp_method=InterpMethod::Average)
 A simple linear interpolation of values between cell centers to a cell face. More...
 
template<typename T1 , typename T2 >
libMesh::CompareTypes< T1, T2 >::supertype harmonicInterpolation (const T1 &value1, const T2 &value2, const FaceInfo &fi, const bool one_is_elem)
 Computes the harmonic mean (1/(gc/value1+(1-gc)/value2)) of Reals, RealVectorValues and RealTensorValues while accounting for the possibility that one or both of them are AD. More...
 
template<typename T , typename T2 , typename T3 >
libMesh::CompareTypes< T, T2 >::supertype skewCorrectedLinearInterpolation (const T &value1, const T2 &value2, const T3 &face_gradient, const FaceInfo &fi, const bool one_is_elem)
 Linear interpolation with skewness correction using the face gradient. More...
 
template<typename T , typename T2 , typename T3 >
void interpolate (InterpMethod m, T &result, const T2 &value1, const T3 &value2, const FaceInfo &fi, const bool one_is_elem)
 Provides interpolation of face values for non-advection-specific purposes (although it can/will still be used by advective kernels sometimes). More...
 
template<typename T , FunctorEvaluationKind FEK = FunctorEvaluationKind::Value>
linearInterpolation (const FunctorBase< T > &functor, const FaceArg &face, const StateArg &time)
 perform a possibly skew-corrected linear interpolation by evaluating the supplied functor with the provided functor face argument More...
 
template<typename T1 , typename T2 , typename T3 , template< typename > class Vector1, template< typename > class Vector2>
void interpolate (InterpMethod m, Vector1< T1 > &result, const T2 &fi_elem_advected, const T2 &fi_neighbor_advected, const Vector2< T3 > &fi_elem_advector, const Vector2< T3 > &fi_neighbor_advector, const FaceInfo &fi)
 Computes the product of the advected and the advector based on the given interpolation method. More...
 
template<typename T , typename T2 , typename T3 , typename Vector >
void interpolate (InterpMethod m, T &result, const T2 &value1, const T3 &value2, const Vector &advector, const FaceInfo &fi, const bool one_is_elem)
 Provides interpolation of face values for advective flux kernels. More...
 
ADReal gradUDotNormal (const FaceInfo &face_info, const MooseVariableFV< Real > &fv_var, const Moose::StateArg &time, bool correct_skewness=false)
 Calculates and returns "grad_u dot normal" on the face to be used for diffusive terms. More...
 
template<typename Scalar , typename Vector >
Scalar rF (const Scalar &phiC, const Scalar &phiD, const Vector &gradC, const RealVectorValue &dCD)
 From Moukalled 12.30. More...
 
template<typename T >
std::pair< T, T > interpCoeffs (const Limiter< T > &limiter, const T &phi_upwind, const T &phi_downwind, const libMesh::VectorValue< T > *const grad_phi_upwind, const libMesh::VectorValue< T > *const grad_phi_face, const Real &max_value, const Real &min_value, const FaceInfo &fi, const bool fi_elem_is_upwind)
 Produce the interpolation coefficients in the equation: More...
 
template<typename Scalar , typename Vector , typename Enable = typename std::enable_if<libMesh::ScalarTraits<Scalar>::value>::type>
Scalar interpolate (const Limiter< Scalar > &limiter, const Scalar &phi_upwind, const Scalar &phi_downwind, const Vector *const grad_phi_upwind, const FaceInfo &fi, const bool fi_elem_is_upwind)
 Interpolates with a limiter. More...
 
template<typename Limiter , typename T , typename Tensor >
libMesh::VectorValue< T > interpolate (const Limiter &limiter, const TypeVector< T > &phi_upwind, const TypeVector< T > &phi_downwind, const Tensor *const grad_phi_upwind, const FaceInfo &fi, const bool fi_elem_is_upwind)
 Vector overload. More...
 
template<typename T >
fullLimitedInterpolation (const Limiter< T > &limiter, const T &phi_upwind, const T &phi_downwind, const VectorValue< T > *const grad_phi_upwind, const VectorValue< T > *const grad_phi_face, const Real &max_value, const Real &min_value, const FaceArg &face)
 This function performs a full limited interpolation of a variable, taking into account the values and gradients at both upwind and downwind locations, as well as geometric information and limits. More...
 
template<typename T , FunctorEvaluationKind FEK = FunctorEvaluationKind::Value, typename Enable = typename std::enable_if<ScalarTraits<T>::value>::type>
std::pair< Real, RealcomputeMinMaxValue (const FunctorBase< T > &functor, const FaceArg &face, const StateArg &time)
 This function calculates the minimum and maximum values within a two-cell stencil. More...
 
template<typename T >
std::pair< Real, RealcomputeMinMaxValue (const FunctorBase< VectorValue< T >> &functor, const FaceArg &face, const StateArg &time, const unsigned int &component)
 This function calculates the minimum and maximum values of a specified component within a two-cell stencil. More...
 
template<typename T , FunctorEvaluationKind FEK = FunctorEvaluationKind::Value, typename Enable = typename std::enable_if<libMesh::ScalarTraits<T>::value>::type>
std::pair< std::pair< T, T >, std::pair< T, T > > interpCoeffsAndAdvected (const FunctorBase< T > &functor, const FaceArg &face, const StateArg &time)
 This function interpolates values using a specified limiter and face argument. More...
 
template<typename T , FunctorEvaluationKind FEK = FunctorEvaluationKind::Value, typename Enable = typename std::enable_if<libMesh::ScalarTraits<T>::value>::type>
interpolate (const FunctorBase< T > &functor, const FaceArg &face, const StateArg &time)
 This function interpolates values at faces in a computational grid using a specified functor, face argument, and evaluation kind. More...
 
template<typename T >
libMesh::VectorValue< T > interpolate (const FunctorBase< libMesh::VectorValue< T >> &functor, const FaceArg &face, const StateArg &time)
 This function interpolates vector values at faces in a computational grid using a specified functor, face argument, and limiter type. More...
 
template<typename T >
containerInterpolate (const FunctorBase< T > &functor, const FaceArg &face, const StateArg &time)
 This function interpolates container values at faces in a computational grid using a specified functor, face argument, and limiter type. More...
 
template<typename T >
std::vector< T > interpolate (const FunctorBase< std::vector< T >> &functor, const FaceArg &face, const StateArg &time)
 
template<typename T , std::size_t N>
std::array< T, N > interpolate (const FunctorBase< std::array< T, N >> &functor, const FaceArg &face, const StateArg &time)
 
template<typename SubdomainRestrictable >
bool onBoundary (const SubdomainRestrictable &obj, const FaceInfo &fi)
 Return whether the supplied face is on a boundary of the object's execution. More...
 
bool onBoundary (const std::set< SubdomainID > &subs, const FaceInfo &fi)
 Determine whether the passed-in face is on the boundary of an object that lives on the provided subdomains. More...
 
const MooseEnum moose_limiter_type ("vanLeer=0 upwind=1 central_difference=2 min_mod=3 sou=4 quick=5 venkatakrishnan=6", "upwind")
 

Variables

const MooseEnum moose_limiter_type
 

Enumeration Type Documentation

◆ InterpMethod

This codifies a set of available ways to interpolate with elem+neighbor solution information to calculate values (e.g.

solution, material properties, etc.) at the face (centroid). These methods are used in the class's interpolate functions. Some interpolation methods are only meant to be used with advective terms (e.g. upwind), others are more generically applicable.

Enumerator
Average 

gc*elem+(1-gc)*neighbor

HarmonicAverage 

1/(gc/elem+(1-gc)/neighbor)

SkewCorrectedAverage 

(gc*elem+(1-gc)*neighbor)+gradient*(rf-rf')

Upwind 

weighted

RhieChow 
VanLeer 
MinMod 
SOU 
QUICK 
Venkatakrishnan 

Definition at line 35 of file MathFVUtils.h.

◆ LimiterType

enum Moose::FV::LimiterType : int
strong
Enumerator
VanLeer 
Upwind 
CentralDifference 
MinMod 
SOU 
QUICK 
Venkatakrishnan 

Definition at line 26 of file Limiter.h.

◆ LinearFVComputationMode

Function Documentation

◆ advectedInterpolationParameter()

InputParameters Moose::FV::advectedInterpolationParameter ( )
Returns
An input parameters object that contains the advected_interp_method parameter, e.g. the interpolation method to use for an advected quantity

Definition at line 68 of file MathFVUtils.C.

Referenced by FVAdvection::validParams(), and LinearFVAdvection::validParams().

69 {
70  auto params = emptyInputParameters();
71  params.addParam<MooseEnum>("advected_interp_method",
73  "The interpolation to use for the advected quantity. Options are "
74  "'upwind', 'average', 'sou' (for second-order upwind), 'min_mod', "
75  "'vanLeer', 'quick', 'venkatakrishnan', and "
76  "'skewness-corrected' with the default being 'upwind'.");
77  return params;
78 }
MooseEnum interpolationMethods()
Returns an enum with all the currently supported interpolation methods and the current default for FV...
Definition: MathFVUtils.C:61
InputParameters emptyInputParameters()
This is a "smart" enum class intended to replace many of the shortcomings in the C++ enum type It sho...
Definition: MooseEnum.h:33

◆ computeMinMaxValue() [1/2]

template<typename T , FunctorEvaluationKind FEK = FunctorEvaluationKind::Value, typename Enable = typename std::enable_if<ScalarTraits<T>::value>::type>
std::pair<Real, Real> Moose::FV::computeMinMaxValue ( const FunctorBase< T > &  functor,
const FaceArg face,
const StateArg time 
)

This function calculates the minimum and maximum values within a two-cell stencil.

The stencil includes the immediate neighboring elements of the face's associated element and the neighboring elements of those neighbors. It evaluates the values using a provided functor and accounts for the valid (non-null) neighbors.

Template Parameters
TThe data type for the values being computed. This is typically a scalar type.
FEKEnumeration type FunctorEvaluationKind with a default value of FunctorEvaluationKind::Value. This determines the kind of evaluation the functor will perform.
EnableA type trait used for SFINAE (Substitution Failure Is Not An Error) to ensure that T is a valid scalar type as determined by ScalarTraits<T>.
Parameters
functorAn object of a functor class derived from FunctorBase<T>. This object provides the genericEvaluate method to compute the value at a given element and time.
faceAn argument representing the face in the computational domain. The face provides access to neighboring elements via neighbor_ptr_range().
timeAn argument representing the state or time at which the evaluation is performed.
Returns
std::pair<T, T> A pair containing the minimum and maximum values computed across the two-cell stencil. The first element is the maximum value, and the second element is the minimum value.

Definition at line 645 of file MathFVUtils.h.

Referenced by interpCoeffsAndAdvected(), and interpolate().

646 {
647  // Initialize max_value to 0 and min_value to a large value
649 
650  // Iterate over the direct neighbors of the element associated with the face
651  for (const auto neighbor : (*face.fi).elem().neighbor_ptr_range())
652  {
653  // If not a valid neighbor, skip to the next one
654  if (neighbor == nullptr)
655  continue;
656 
657  // Evaluate the functor at the neighbor and update max_value and min_value
658  const ElemArg local_cell_arg = {neighbor, false};
659  const auto value =
660  MetaPhysicL::raw_value(functor.template genericEvaluate<FEK>(local_cell_arg, time));
661  max_value = std::max(max_value, value);
662  min_value = std::min(min_value, value);
663  }
664 
665  // Iterate over the neighbors of the neighbor
666  for (const auto neighbor : (*face.fi).neighbor().neighbor_ptr_range())
667  {
668  // If not a valid neighbor, skip to the next one
669  if (neighbor == nullptr)
670  continue;
671 
672  // Evaluate the functor at the neighbor and update max_value and min_value
673  const ElemArg local_cell_arg = {neighbor, false};
674  const auto value =
675  MetaPhysicL::raw_value(functor.template genericEvaluate<FEK>(local_cell_arg, time));
676  max_value = std::max(max_value, value);
677  min_value = std::min(min_value, value);
678  }
679 
680  // Return a pair containing the computed maximum and minimum values
681  return std::make_pair(max_value, min_value);
682 }
auto raw_value(const Eigen::Map< T > &in)
Definition: EigenADReal.h:73
auto max(const L &left, const R &right)
Real value(unsigned n, unsigned alpha, unsigned beta, Real x)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
auto min(const L &left, const R &right)

◆ computeMinMaxValue() [2/2]

template<typename T >
std::pair<Real, Real> Moose::FV::computeMinMaxValue ( const FunctorBase< VectorValue< T >> &  functor,
const FaceArg face,
const StateArg time,
const unsigned int component 
)

This function calculates the minimum and maximum values of a specified component within a two-cell stencil.

The stencil includes the immediate neighboring elements of the face's associated element and the neighboring elements of those neighbors. It evaluates the values using a provided functor and accounts for the valid (non-null) neighbors.

Template Parameters
TThe data type for the values being computed. This is typically a scalar type.
Parameters
functorAn object of a functor class derived from FunctorBase<VectorValue<T>>. This object provides the operator() method to compute the value at a given element and time.
faceAn argument representing the face in the computational domain. The face provides access to neighboring elements via neighbor_ptr_range().
timeAn argument representing the state or time at which the evaluation is performed.
componentAn unsigned integer representing the specific component of the vector to be evaluated.
Returns
std::pair<T, T> A pair containing the minimum and maximum values of the specified component computed across the two-cell stencil. The first element is the maximum value, and the second element is the minimum value.

Usage: This function is typically used in the finite volume methods for min-max computations over stencils (neighborhoods). It helps compute the limiting for limited second order upwind at the faces.

Definition at line 711 of file MathFVUtils.h.

715 {
716  // Initialize max_value to 0 and min_value to a large value
718 
719  // Iterate over the direct neighbors of the element associated with the face
720  for (const auto neighbor : (*face.fi).elem().neighbor_ptr_range())
721  {
722  // If not a valid neighbor, skip to the next one
723  if (neighbor == nullptr)
724  continue;
725 
726  // Evaluate the functor at the neighbor for the specified component and update max_value and
727  // min_value
728  const ElemArg local_cell_arg = {neighbor, false};
729  const auto value = MetaPhysicL::raw_value(functor(local_cell_arg, time)(component));
730  max_value = std::max(max_value, value);
731  min_value = std::min(min_value, value);
732  }
733 
734  // Iterate over the neighbors of the neighbor associated with the face
735  for (const auto neighbor : (*face.fi).neighbor().neighbor_ptr_range())
736  {
737  // If not a valid neighbor, skip to the next one
738  if (neighbor == nullptr)
739  continue;
740 
741  // Evaluate the functor at the neighbor for the specified component and update max_value and
742  // min_value
743  const ElemArg local_cell_arg = {neighbor, false};
744  const auto value = MetaPhysicL::raw_value(functor(local_cell_arg, time)(component));
745  max_value = std::max(max_value, value);
746  min_value = std::min(min_value, value);
747  }
748 
749  // Return a pair containing the computed maximum and minimum values
750  return std::make_pair(max_value, min_value);
751 }
auto raw_value(const Eigen::Map< T > &in)
Definition: EigenADReal.h:73
auto max(const L &left, const R &right)
Real value(unsigned n, unsigned alpha, unsigned beta, Real x)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
auto min(const L &left, const R &right)

◆ containerInterpolate()

template<typename T >
T Moose::FV::containerInterpolate ( const FunctorBase< T > &  functor,
const FaceArg face,
const StateArg time 
)

This function interpolates container values at faces in a computational grid using a specified functor, face argument, and limiter type.

It handles different limiter types and performs interpolation accordingly.

Template Parameters
TThe data type for the container values being interpolated. This is typically a container type such as a vector or array.
Parameters
functorAn object of a functor class derived from FunctorBase<T>. This object provides the operator() method to compute the value at a given element and time.
faceAn argument representing the face in the computational domain. The face provides access to neighboring elements and limiter type.
timeAn argument representing the state or time at which the evaluation is performed.
Returns
T The interpolated container value at the face.

Usage: This function is used for interpolating container values at faces in a finite volume method, ensuring that the interpolation adheres to the constraints imposed by the limiter.

Definition at line 1087 of file MathFVUtils.h.

Referenced by interpolate().

1088 {
1089  typedef typename FunctorBase<T>::GradientType ContainerGradientType;
1090  typedef typename ContainerGradientType::value_type GradientType;
1091  const GradientType * const example_gradient = nullptr;
1092 
1093  mooseAssert(face.fi, "this must be non-null");
1094  const auto limiter = Limiter<typename T::value_type>::build(face.limiter_type);
1095 
1096  const auto upwind_arg = face.elem_is_upwind ? face.makeElem() : face.makeNeighbor();
1097  const auto downwind_arg = face.elem_is_upwind ? face.makeNeighbor() : face.makeElem();
1098  const auto phi_upwind = functor(upwind_arg, time);
1099  const auto phi_downwind = functor(downwind_arg, time);
1100 
1101  // initialize in order to get proper size
1102  T ret = phi_upwind;
1103  typename T::value_type coeff_upwind, coeff_downwind;
1104 
1105  if (face.limiter_type == LimiterType::Upwind ||
1106  face.limiter_type == LimiterType::CentralDifference)
1107  {
1108  for (const auto i : index_range(ret))
1109  {
1110  const auto &component_upwind = phi_upwind[i], component_downwind = phi_downwind[i];
1111  std::tie(coeff_upwind, coeff_downwind) = interpCoeffs(*limiter,
1112  component_upwind,
1113  component_downwind,
1114  example_gradient,
1115  example_gradient,
1118  *face.fi,
1119  face.elem_is_upwind);
1120  ret[i] = coeff_upwind * component_upwind + coeff_downwind * component_downwind;
1121  }
1122  }
1123  else
1124  {
1125  const auto grad_phi_upwind = functor.gradient(upwind_arg, time);
1126  for (const auto i : index_range(ret))
1127  {
1128  const auto &component_upwind = phi_upwind[i], component_downwind = phi_downwind[i];
1129  const auto & grad = grad_phi_upwind[i];
1130  std::tie(coeff_upwind, coeff_downwind) = interpCoeffs(*limiter,
1131  component_upwind,
1132  component_downwind,
1133  &grad,
1134  example_gradient,
1137  *face.fi,
1138  face.elem_is_upwind);
1139  ret[i] = coeff_upwind * component_upwind + coeff_downwind * component_downwind;
1140  }
1141  }
1142 
1143  return ret;
1144 }
auto max(const L &left, const R &right)
std::pair< T, T > interpCoeffs(const Limiter< T > &limiter, const T &phi_upwind, const T &phi_downwind, const libMesh::VectorValue< T > *const grad_phi_upwind, const libMesh::VectorValue< T > *const grad_phi_face, const Real &max_value, const Real &min_value, const FaceInfo &fi, const bool fi_elem_is_upwind)
Produce the interpolation coefficients in the equation:
Definition: MathFVUtils.h:472
auto min(const L &left, const R &right)
auto index_range(const T &sizable)

◆ determineElemOneAndTwo()

template<typename FVVar >
std::tuple<const Elem *, const Elem *, bool> Moose::FV::determineElemOneAndTwo ( const FaceInfo fi,
const FVVar &  var 
)

This utility determines element one and element two given a FaceInfo fi and variable var.

You may ask what in the world "element one" and "element two" means, and that would be a very good question. What it means is: a variable will always have degrees of freedom on element one. A variable may or may not have degrees of freedom on element two. So we are introducing a second terminology here. FaceInfo geometric objects have element-neighbor pairs. These element-neighbor pairs are purely geometric and have no relation to the algebraic system of variables. The elem1-elem2 notation introduced here is based on dof/algebraic information and may very well be different from variable to variable, e.g. elem1 may correspond to the FaceInfo elem for one variable (and correspondingly elem2 will be the FaceInfo neighbor), but elem1 may correspond to the FaceInfo neighbor for another variable (and correspondingly for that variable elem2 will be the FaceInfo elem).

Returns
A tuple, where the first item is elem1, the second item is elem2, and the third item is a boolean indicating whether elem1 corresponds to the FaceInfo elem

Definition at line 130 of file FVUtils.h.

Referenced by MooseVariableFV< Real >::getExtrapolatedBoundaryFaceValue().

131 {
132  auto ft = fi.faceType(std::make_pair(var.number(), var.sys().number()));
133  mooseAssert(ft == FaceInfo::VarFaceNeighbors::BOTH
134  ? var.hasBlocks(fi.elem().subdomain_id()) && fi.neighborPtr() &&
135  var.hasBlocks(fi.neighborPtr()->subdomain_id())
136  : true,
137  "Finite volume variable " << var.name()
138  << " does not exist on both sides of the face despite "
139  "what the FaceInfo is telling us.");
140  mooseAssert(ft == FaceInfo::VarFaceNeighbors::ELEM
141  ? var.hasBlocks(fi.elem().subdomain_id()) &&
142  (!fi.neighborPtr() || !var.hasBlocks(fi.neighborPtr()->subdomain_id()))
143  : true,
144  "Finite volume variable " << var.name()
145  << " does not exist on or only on the elem side of the "
146  "face despite what the FaceInfo is telling us.");
147  mooseAssert(ft == FaceInfo::VarFaceNeighbors::NEIGHBOR
148  ? fi.neighborPtr() && var.hasBlocks(fi.neighborPtr()->subdomain_id()) &&
149  !var.hasBlocks(fi.elem().subdomain_id())
150  : true,
151  "Finite volume variable " << var.name()
152  << " does not exist on or only on the neighbor side of the "
153  "face despite what the FaceInfo is telling us.");
154 
155  const bool one_is_elem =
157  const Elem * const elem_one = one_is_elem ? &fi.elem() : fi.neighborPtr();
158  mooseAssert(elem_one, "This elem should be non-null!");
159  const Elem * const elem_two = one_is_elem ? fi.neighborPtr() : &fi.elem();
160 
161  return std::make_tuple(elem_one, elem_two, one_is_elem);
162 }
const Elem & elem() const
Definition: FaceInfo.h:81
const Elem * neighborPtr() const
Definition: FaceInfo.h:84
VarFaceNeighbors faceType(const std::pair< unsigned int, unsigned int > &var_sys) const
Returns which side(s) the given variable-system number pair is defined on for this face...
Definition: FaceInfo.h:225

◆ elemHasFaceInfo()

bool Moose::FV::elemHasFaceInfo ( const Elem &  elem,
const Elem *const  neighbor 
)

This function infers based on elements if the faceinfo between them belongs to the element or not.

Parameters
elemReference to an element
neighborPointer to the neighbor of the element
Returns
If the element (first argument) is the owner of the faceinfo between the two elements

Definition at line 21 of file FVUtils.C.

Referenced by MooseMesh::buildFiniteVolumeInfo(), and loopOverElemFaceInfo().

22 {
23  // The face info belongs to elem:
24  // * at all mesh boundaries (i.e. where there is no neighbor)
25  // * if the element faces a neighbor which is on a lower refinement level
26  // * if the element is active and it has a lower ID than its neighbor
27  if (!neighbor)
28  return true;
29  else if (elem.level() != neighbor->level())
30  return neighbor->level() < elem.level();
31  else if (!neighbor->active())
32  return false;
33  else
34  return elem.id() < neighbor->id();
35 }

◆ fullLimitedInterpolation()

template<typename T >
T Moose::FV::fullLimitedInterpolation ( const Limiter< T > &  limiter,
const T &  phi_upwind,
const T &  phi_downwind,
const VectorValue< T > *const  grad_phi_upwind,
const VectorValue< T > *const  grad_phi_face,
const Real max_value,
const Real min_value,
const FaceArg face 
)

This function performs a full limited interpolation of a variable, taking into account the values and gradients at both upwind and downwind locations, as well as geometric information and limits.

It applies the specified limiter to ensure the interpolation adheres to the constraints imposed by the limiter.

Template Parameters
TThe data type of the scalar values and the return type.
Parameters
limiterThe limiter object used to compute the flux limiting ratio.
phi_upwindThe field value at the upwind location.
phi_downwindThe field value at the downwind location.
grad_phi_upwindPointer to the gradient at the upwind location.
grad_phi_facePointer to the gradient at the face location.
fiFaceInfo object containing geometric details such as face centroid and cell centroids.
fi_elem_is_upwindBoolean indicating if the face info element is upwind.
max_valueThe maximum allowable value.
min_valueThe minimum allowable value.
Returns
The computed limited interpolation value.

Definition at line 580 of file MathFVUtils.h.

Referenced by interpolate().

588 {
589  // Compute the direction vector based on whether the current element is upwind
590  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  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  face.fi,
601  face.elem_is_upwind);
602 
603  // Determine the face centroid and the appropriate cell centroid
604  const auto & face_centroid = face.fi->faceCentroid();
605  const auto & cell_centroid =
606  face.elem_is_upwind ? face.fi->elemCentroid() : face.fi->neighborCentroid();
607 
608  // Compute the delta value at the face
609  T delta_face;
610  if (face.limiter_type == LimiterType::Venkatakrishnan || face.limiter_type == LimiterType::SOU)
611  delta_face = (*grad_phi_upwind) * (face_centroid - cell_centroid);
612  else
613  delta_face = (*grad_phi_face) * (face_centroid - cell_centroid);
614 
615  // Return the interpolated value
616  return phi_upwind + beta * delta_face;
617 }

◆ gradUDotNormal()

ADReal Moose::FV::gradUDotNormal ( const FaceInfo face_info,
const MooseVariableFV< Real > &  fv_var,
const Moose::StateArg time,
bool  correct_skewness = false 
)

Calculates and returns "grad_u dot normal" on the face to be used for diffusive terms.

If using any cross-diffusion corrections, etc. all those calculations should be handled appropriately by this function.

Definition at line 18 of file MathFVUtils.C.

Referenced by FVFluxKernel::gradUDotNormal().

23 {
24  return fv_var.adGradSln(face_info, time, correct_skewness) * face_info.normal();
25 }
const Point & normal() const
Returns the unit normal vector for the face oriented outward from the face&#39;s elem element...
Definition: FaceInfo.h:68
const ADTemplateVariableGradient< OutputType > & adGradSln() const override
AD grad solution getter.

◆ greenGaussGradient() [1/8]

template<typename T , typename Enable = typename std::enable_if<libMesh::ScalarTraits<T>::value>::type>
libMesh::VectorValue<T> Moose::FV::greenGaussGradient ( const ElemArg elem_arg,
const StateArg state_arg,
const FunctorBase< T > &  functor,
const bool  two_term_boundary_expansion,
const MooseMesh mesh,
const bool  force_green_gauss = false 
)

Compute a cell gradient using the method of Green-Gauss.

Parameters
elem_argAn element argument specifying the current element and whether to perform skew corrections
state_argA state argument that indicates what temporal / solution iteration data to use when evaluating the provided functor
functorThe functor that will provide information such as cell and face value evaluations necessary to construct the cell gradient
two_term_boundary_expansionWhether to perform a two-term expansion to compute extrapolated boundary face values. If this is true, then an implicit system has to be solved. If false, then the cell center value will be used as the extrapolated boundary face value
meshThe mesh on which we are computing the gradient
Returns
The computed cell gradient

Definition at line 41 of file GreenGaussGradient.h.

Referenced by MooseVariableFV< Real >::adGradSln(), PiecewiseByBlockLambdaFunctor< T >::evaluateGradient(), and greenGaussGradient().

47 {
48  mooseAssert(elem_arg.elem, "This should be non-null");
49  const auto coord_type = mesh.getCoordSystem(elem_arg.elem->subdomain_id());
50  const auto rz_radial_coord = mesh.getAxisymmetricRadialCoord();
51 
52  const T elem_value = functor(elem_arg, state_arg);
53 
54  // We'll count the number of extrapolated boundary faces (ebfs)
55  unsigned int num_ebfs = 0;
56 
57  // Gradient to be returned
58  VectorValue<T> grad;
59 
60  if (!elem_arg.correct_skewness || force_green_gauss) // Do Green-Gauss
61  {
62  try
63  {
64  bool volume_set = false;
65  Real volume = 0;
66 
67  // If we are performing a two term Taylor expansion for extrapolated boundary faces (faces on
68  // boundaries that do not have associated Dirichlet conditions), then the element gradient
69  // depends on the boundary face value and the boundary face value depends on the element
70  // gradient, so we have a system of equations to solve. Here is the system:
71  //
72  // \nabla \phi_C - \frac{1}{V} \sum_{ebf} \phi_{ebf} \vec{S_f} =
73  // \frac{1}{V} \sum_{of} \phi_{of} \vec{S_f} eqn. 1
74  //
75  // \phi_{ebf} - \vec{d_{Cf}} \cdot \nabla \phi_C = \phi_C eqn. 2
76  //
77  // where $C$ refers to the cell centroid, $ebf$ refers to an extrapolated boundary face, $of$
78  // refers to "other faces", e.g. non-ebf faces, and $f$ is a general face. $d_{Cf}$ is the
79  // vector drawn from the element centroid to the face centroid, and $\vec{S_f}$ is the surface
80  // vector, e.g. the face area times the outward facing normal
81 
82  // ebf eqns: element gradient coefficients, e.g. eqn. 2, LHS term 2 coefficient
83  std::vector<libMesh::VectorValue<Real>> ebf_grad_coeffs;
84  // ebf eqns: rhs b values. These will actually correspond to the elem_value so we can use a
85  // pointer and avoid copying. This is the RHS of eqn. 2
86  std::vector<const T *> ebf_b;
87 
88  // elem grad eqns: ebf coefficients, e.g. eqn. 1, LHS term 2 coefficients
89  std::vector<libMesh::VectorValue<Real>> grad_ebf_coeffs;
90  // elem grad eqns: rhs b value, e.g. eqn. 1 RHS
91  libMesh::VectorValue<T> grad_b = 0;
92 
93  auto action_functor = [&volume_set,
94  &volume,
95  &elem_value,
96  &elem_arg,
97  &num_ebfs,
98  &ebf_grad_coeffs,
99  &ebf_b,
100  &grad_ebf_coeffs,
101  &grad_b,
102  &state_arg,
103  &functor,
104  two_term_boundary_expansion,
105  coord_type,
106  rz_radial_coord](const Elem & libmesh_dbg_var(functor_elem),
107  const Elem *,
108  const FaceInfo * const fi,
109  const Point & surface_vector,
110  Real coord,
111  const bool elem_has_info)
112  {
113  mooseAssert(fi, "We need a FaceInfo for this action_functor");
114  mooseAssert(
115  elem_arg.elem == &functor_elem,
116  "Just a sanity check that the element being passed in is the one we passed out.");
117 
118  if (functor.isExtrapolatedBoundaryFace(*fi, elem_arg.elem, state_arg))
119  {
120  if (two_term_boundary_expansion)
121  {
122  num_ebfs += 1;
123 
124  // eqn. 2
125  ebf_grad_coeffs.push_back(-1. * (elem_has_info
126  ? (fi->faceCentroid() - fi->elemCentroid())
127  : (fi->faceCentroid() - fi->neighborCentroid())));
128  ebf_b.push_back(&elem_value);
129 
130  // eqn. 1
131  grad_ebf_coeffs.push_back(-surface_vector);
132  }
133  else
134  // We are doing a one-term expansion for the extrapolated boundary faces, in which case
135  // we have no eqn. 2 and we have no second term in the LHS of eqn. 1. Instead we apply
136  // the element centroid value as the face value (one-term expansion) in the RHS of eqn.
137  // 1
138  grad_b += surface_vector * elem_value;
139  }
140  else
141  grad_b +=
142  surface_vector * functor(Moose::FaceArg{fi,
144  true,
145  elem_arg.correct_skewness,
146  elem_arg.elem,
147  nullptr},
148  state_arg);
149 
150  if (!volume_set)
151  {
152  // We use the FaceInfo volumes because those values have been pre-computed and cached.
153  // An explicit call to elem->volume() here would incur unnecessary expense
154  if (elem_has_info)
155  {
157  fi->elemCentroid(), coord, coord_type, rz_radial_coord);
158  volume = fi->elemVolume() * coord;
159  }
160  else
161  {
163  fi->neighborCentroid(), coord, coord_type, rz_radial_coord);
164  volume = fi->neighborVolume() * coord;
165  }
166 
167  volume_set = true;
168  }
169  };
170 
172  *elem_arg.elem, mesh, action_functor, coord_type, rz_radial_coord);
173 
174  mooseAssert(volume_set && volume > 0, "We should have set the volume");
175  grad_b /= volume;
176 
177  if (coord_type == Moose::CoordinateSystemType::COORD_RZ)
178  {
179  mooseAssert(rz_radial_coord != libMesh::invalid_uint, "rz_radial_coord must be set");
180  grad_b(rz_radial_coord) -= elem_value / elem_arg.elem->vertex_average()(rz_radial_coord);
181  }
182 
183  mooseAssert(
185  "We have not yet implemented the correct translation from gradient to divergence for "
186  "spherical coordinates yet.");
187 
188  // test for simple case
189  if (num_ebfs == 0)
190  grad = grad_b;
191  else
192  {
193  // We have to solve a system
194  const unsigned int sys_dim = Moose::dim + num_ebfs;
195  DenseVector<T> x(sys_dim), b(sys_dim);
196  DenseMatrix<T> A(sys_dim, sys_dim);
197 
198  // Let's make i refer to Moose::dim indices, and j refer to num_ebfs indices
199 
200  // eqn. 1
201  for (const auto i : make_range(Moose::dim))
202  {
203  // LHS term 1 coeffs
204  A(i, i) = 1;
205 
206  // LHS term 2 coeffs
207  for (const auto j : make_range(num_ebfs))
208  A(i, Moose::dim + j) = grad_ebf_coeffs[j](i) / volume;
209 
210  // RHS
211  b(i) = grad_b(i);
212  }
213 
214  // eqn. 2
215  for (const auto j : make_range(num_ebfs))
216  {
217  // LHS term 1 coeffs
218  A(Moose::dim + j, Moose::dim + j) = 1;
219 
220  // LHS term 2 coeffs
221  for (const auto i : make_range(unsigned(Moose::dim)))
222  A(Moose::dim + j, i) = ebf_grad_coeffs[j](i);
223 
224  // RHS
225  b(Moose::dim + j) = *ebf_b[j];
226  }
227 
228  A.lu_solve(b, x);
229  // libMesh is generous about what it considers nonsingular. Let's check a little more
230  // strictly
231  if (MooseUtils::absoluteFuzzyEqual(MetaPhysicL::raw_value(A(sys_dim - 1, sys_dim - 1)), 0))
232  throw libMesh::LogicError("Matrix A is singular!");
233  for (const auto i : make_range(Moose::dim))
234  grad(i) = x(i);
235  }
236  }
237  catch (libMesh::LogicError & e)
238  {
239  // Retry without two-term
240  if (!two_term_boundary_expansion)
241  mooseError(
242  "I believe we should only get singular systems when two-term boundary expansion is "
243  "being used. The error thrown during the computation of the gradient: ",
244  e.what());
245 
246  grad = greenGaussGradient(elem_arg, state_arg, functor, false, mesh, true);
247  }
248  }
249  else // Do Least-Squares
250  {
251  try
252  {
253 
254  // The least squares method aims to find the gradient at the element centroid by minimizing
255  // the difference between the estimated and actual value differences across neighboring cells.
256  // The least squares formulation is:
257  //
258  // Minimize J = \sum_{n} [ w_n ((\nabla \phi_C \cdot \delta \vec{x}_n) - \delta \phi_n) ]^2
259  //
260  // where:
261  // - \(\nabla \phi_C\) is the gradient at the element centroid C (unknown).
262  // - \(\delta \vec{x}_n = \vec{x}_n - \vec{x}_C\) is the vector from the element centroid to
263  // neighbor \(n\).
264  // - \(\delta \phi_n = \phi_n - \phi_C\) is the difference in the scalar value between
265  // neighbor \(n\) and element \(C\).
266  // - w_n = 1/||\vec{x}_n|| is a vector of weights that is used to handle large aspect ratio
267  // cells
268  //
269  // To handle extrapolated boundary faces (faces on boundaries without Dirichlet conditions),
270  // we include additional unknowns and equations in the least squares system.
271  // For ebfs, we may not know the value \(\phi_n\), so we include it as an unknown.
272  // This results in an augmented system that simultaneously solves for the gradient and the ebf
273  // values.
274 
275  // Get estimated number of faces in a cell
276  const auto ptr_range = elem_arg.elem->neighbor_ptr_range();
277  const size_t num_faces = std::distance(ptr_range.begin(), ptr_range.end());
278 
279  // Lists to store position differences and value differences for least squares computation
280  std::vector<Point> delta_x_list; // List of position differences between neighbor centroids
281  // and element centroid
282  delta_x_list.reserve(num_faces);
283 
284  std::vector<T>
285  delta_phi_list; // List of value differences between neighbor values and element value
286  delta_phi_list.reserve(num_faces);
287 
288  // Variables to handle extrapolated boundary faces (ebfs) in the least squares method
289  std::vector<Point> delta_x_ebf_list; // Position differences for ebfs
290  delta_phi_list.reserve(num_faces);
291 
292  std::vector<VectorValue<Real>>
293  ebf_grad_coeffs; // Coefficients for the gradient in ebf equations
294  delta_phi_list.reserve(num_faces);
295 
296  std::vector<const T *> ebf_b; // RHS values for ebf equations (pointers to avoid copying)
297  delta_phi_list.reserve(num_faces);
298 
299  unsigned int num_ebfs = 0; // Number of extrapolated boundary faces
300 
301  // Action functor to collect data from each face of the element
302  auto action_functor = [&elem_value,
303  &elem_arg,
304  &num_ebfs,
305  &ebf_grad_coeffs,
306  &ebf_b,
307  &delta_x_list,
308  &delta_phi_list,
309  &delta_x_ebf_list,
310  &state_arg,
311  &functor,
312  two_term_boundary_expansion](const Elem & libmesh_dbg_var(loc_elem),
313  const Elem * loc_neigh,
314  const FaceInfo * const fi,
315  const Point & /*surface_vector*/,
316  Real /*coord*/,
317  const bool elem_has_info)
318  {
319  mooseAssert(fi, "We need a FaceInfo for this action_functor");
320  mooseAssert(
321  elem_arg.elem == &loc_elem,
322  "Just a sanity check that the element being passed in is the one we passed out.");
323 
324  if (functor.isExtrapolatedBoundaryFace(*fi, elem_arg.elem, state_arg))
325  {
326  // Extrapolated Boundary Face (ebf)
327  const Point delta_x = elem_has_info ? (fi->faceCentroid() - fi->elemCentroid())
328  : (fi->faceCentroid() - fi->neighborCentroid());
329  delta_x_list.push_back(delta_x);
330 
331  T delta_phi;
332 
333  if (two_term_boundary_expansion)
334  {
335  // Two-term expansion: include ebf values as unknowns in the augmented system
336  num_ebfs += 1;
337 
338  // Coefficient for the gradient in the ebf equation
339  ebf_grad_coeffs.push_back(-delta_x);
340  // RHS value for the ebf equation
341  ebf_b.push_back(&elem_value);
342 
343  delta_phi = -elem_value;
344  delta_x_ebf_list.push_back(delta_x);
345  }
346  else
347  {
348  // One-term expansion: assume zero difference across the boundary (\delta \phi = 0)
349  delta_phi = T(0);
350  }
351  delta_phi_list.push_back(delta_phi);
352  }
353  else
354  {
355  // Internal Face or Boundary Face with Dirichlet condition
356  Point delta_x;
357  T neighbor_value;
358  if (functor.isInternalFace(*fi))
359  {
360  // Internal face with a neighboring element
361  delta_x = loc_neigh->vertex_average() - elem_arg.elem->vertex_average();
362 
363  const ElemArg neighbor_arg = {loc_neigh, elem_arg.correct_skewness};
364  neighbor_value = functor(neighbor_arg, state_arg);
365  }
366  else
367  {
368  // Boundary face with Dirichlet condition
369  delta_x = elem_has_info ? (fi->faceCentroid() - fi->elemCentroid())
370  : (fi->faceCentroid() - fi->neighborCentroid());
371  neighbor_value = functor(Moose::FaceArg{fi,
373  true,
374  elem_arg.correct_skewness,
375  elem_arg.elem,
376  nullptr},
377  state_arg);
378  }
379 
380  delta_x_list.push_back(delta_x);
381  delta_phi_list.push_back(neighbor_value - elem_value);
382  }
383  };
384 
385  // Loop over element faces and apply the action_functor
387  *elem_arg.elem, mesh, action_functor, coord_type, rz_radial_coord);
388 
389  // Compute Least Squares gradient
390  const unsigned int n = delta_x_list.size();
391  mooseAssert(n >= dim, "Not enough neighbors to perform least squares");
392 
393  DenseMatrix<T> ATA(dim, dim);
394  DenseVector<T> ATb(dim);
395 
396  // Assemble the normal equations for the least squares method
397  // ATA = \sum_n (\delta \vec{x}_n^T * \delta \vec{x}_n)
398  // ATb = \sum_n (\delta \vec{x}_n^T * \delta \phi_n)
399  ATA.zero();
400  ATb.zero();
401 
402  for (unsigned int i = 0; i < n; ++i)
403  {
404  const Point & delta_x = delta_x_list[i];
405  const T & delta_phi = delta_phi_list[i];
406 
407  for (unsigned int d1 = 0; d1 < dim; ++d1)
408  {
409  const Real dx_d1 = delta_x(d1);
410  const Real dx_norm = delta_x.norm();
411  const Real dx_d1_norm = dx_d1 / dx_norm;
412 
413  for (unsigned int d2 = 0; d2 < dim; ++d2)
414  {
415  const Real dx_d2_norm = delta_x(d2) / dx_norm;
416  ATA(d1, d2) += dx_d1_norm * dx_d2_norm;
417  }
418 
419  ATb(d1) += dx_d1_norm * delta_phi / dx_norm;
420  }
421  }
422 
423  // Add regularization to ATA to prevent singularity
424  T epsilon = 1e-12; // Small regularization parameter
425  for (unsigned int d = 0; d < dim; ++d)
426  {
427  ATA(d, d) += epsilon;
428  }
429 
430  if (num_ebfs == 0)
431  {
432  // Solve and store the least squares gradient
433  DenseVector<T> grad_ls(dim);
434  ATA.lu_solve(ATb, grad_ls);
435 
436  // Store the least squares gradient
437  for (unsigned int d = 0; d < dim; ++d)
438  {
439  grad(d) = grad_ls(d);
440  }
441  }
442  else
443  {
444  // We have to solve a system of equations
445  const unsigned int sys_dim = Moose::dim + num_ebfs;
446  DenseVector<T> x(sys_dim), b(sys_dim);
447  DenseMatrix<T> A(sys_dim, sys_dim);
448 
449  // Let's make i refer to Moose::dim indices, and j refer to num_ebfs indices
450 
451  // eqn. 1: Element gradient equations
452  for (unsigned int d1 = 0; d1 < dim; ++d1)
453  {
454  // LHS term 1 coefficients (gradient components)
455  for (unsigned int d2 = 0; d2 < dim; ++d2)
456  A(d1, d2) = ATA(d1, d2);
457 
458  // RHS
459  b(d1) = ATb(d1);
460  }
461 
462  // LHS term 2 coefficients (ebf contributions)
463  for (const auto i : make_range(num_ebfs))
464  {
465  const Point & delta_x = delta_x_ebf_list[i];
466  for (unsigned int d1 = 0; d1 < dim; ++d1)
467  {
468  const Real dx_d1 = delta_x(d1);
469  A(d1, Moose::dim + i) -= dx_d1 / delta_x.norm() / delta_x.norm();
470  }
471  }
472 
473  // eqn. 2: Extrapolated boundary face equations
474  for (const auto j : make_range(num_ebfs))
475  {
476  // LHS term 1 coefficients (ebf values)
477  A(Moose::dim + j, Moose::dim + j) = 1;
478 
479  // LHS term 2 coefficients (gradient components)
480  for (const auto i : make_range(unsigned(Moose::dim)))
481  A(Moose::dim + j, i) = ebf_grad_coeffs[j](i);
482 
483  // RHS
484  b(Moose::dim + j) = *ebf_b[j];
485  }
486 
487  // Solve the system
488  A.lu_solve(b, x);
489 
490  // Check for singularity
491  if (MooseUtils::absoluteFuzzyEqual(MetaPhysicL::raw_value(A(sys_dim - 1, sys_dim - 1)), 0))
492  throw libMesh::LogicError("Matrix A is singular!");
493 
494  // Extract the gradient components
495  for (const auto i : make_range(Moose::dim))
496  grad(i) = x(i);
497  }
498 
499  mooseAssert(
501  "We have not yet implemented the correct translation from gradient to divergence for "
502  "spherical coordinates yet.");
503  }
504  catch (libMesh::LogicError & e)
505  {
506  // Log warning and default to simple green Gauss
507  mooseWarning(
508  "Singular matrix encountered in least squares gradient computation. Falling back "
509  "to Green-Gauss gradient.");
510 
511  grad = greenGaussGradient(
512  elem_arg, state_arg, functor, false, mesh, /* force_green_gauss */ true);
513  }
514  }
515 
516  return grad;
517 
518  // Notes to future developer:
519  // Note 1:
520  // For the least squares gradient, the LS matrix could be precomputed and stored for every cell
521  // I tried doing this on October 2024, but the element lookup for these matrices is too slow
522  // and seems better to compute weights on the fly.
523  // Consider building a map from elem_id to these matrices and speed up lookup with octree
524  // Note 2:
525  // In the future one would like to have a hybrid gradient scheme, where:
526  // \nabla \phi = \beta (\nabla \phi)_{LS} + (1 - \beta) (\nabla \phi)_{GG}
527  // Then optize \beta based on mesh heuristics
528 }
void loopOverElemFaceInfo(const Elem &elem, const MooseMesh &mesh, ActionFunctor &act, const Moose::CoordinateSystemType coord_type, const unsigned int rz_radial_coord=libMesh::invalid_uint)
Definition: FVUtils.h:42
bool absoluteFuzzyEqual(const T &var1, const T2 &var2, const T3 &tol=libMesh::TOLERANCE *libMesh::TOLERANCE)
Function to check whether two variables are equal within an absolute tolerance.
Definition: MooseUtils.h:380
const unsigned int invalid_uint
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:302
void coordTransformFactor(const P &point, C &factor, const Moose::CoordinateSystemType coord_type, const unsigned int rz_radial_coord=libMesh::invalid_uint)
compute a coordinate transformation factor
void mooseWarning(Args &&... args)
Emit a warning message with the given stringified, concatenated args.
Definition: MooseError.h:336
MeshBase & mesh
auto raw_value(const Eigen::Map< T > &in)
Definition: EigenADReal.h:73
static constexpr std::size_t dim
This is the dimension of all vector and tensor datastructures used in MOOSE.
Definition: Moose.h:153
This data structure is used to store geometric and variable related metadata about each cell face in ...
Definition: FaceInfo.h:36
Moose::FunctorBase< std::array< T, N > >::GradientType greenGaussGradient(const FaceArg &face_arg, const StateArg &state_arg, const Moose::FunctorBase< std::array< T, N >> &functor, const bool two_term_boundary_expansion, const MooseMesh &mesh)
A structure defining a "face" evaluation calling argument for Moose functors.
Real volume(const MeshBase &mesh, unsigned int dim=libMesh::invalid_uint)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
IntRange< T > make_range(T beg, T end)

◆ greenGaussGradient() [2/8]

template<typename T , typename Enable = typename std::enable_if<libMesh::ScalarTraits<T>::value>::type>
libMesh::VectorValue<T> Moose::FV::greenGaussGradient ( const FaceArg face_arg,
const StateArg state_arg,
const FunctorBase< T > &  functor,
const bool  two_term_boundary_expansion,
const MooseMesh mesh 
)

Compute a face gradient from Green-Gauss cell gradients, with orthogonality correction On the boundaries, the boundary element value is used.

Parameters
face_argA face argument specifying the current faceand whether to perform skew corrections
state_argA state argument that indicates what temporal / solution iteration data to use when evaluating the provided functor
functorThe functor that will provide information such as cell and face value evaluations necessary to construct the cell gradient
two_term_boundary_expansionWhether to perform a two-term expansion to compute extrapolated boundary face values. If this is true, then an implicit system has to be solved. If false, then the cell center value will be used as the extrapolated boundary face value
meshThe mesh on which we are computing the gradient
Returns
The computed cell gradient

Definition at line 548 of file GreenGaussGradient.h.

553 {
554  mooseAssert(face_arg.fi, "We should have a face info to compute a face gradient");
555  const auto & fi = *(face_arg.fi);
556  const auto & elem_arg = face_arg.makeElem();
557  const auto & neighbor_arg = face_arg.makeNeighbor();
558  const bool defined_on_elem = functor.hasBlocks(fi.elemSubdomainID());
559  const bool defined_on_neighbor = fi.neighborPtr() && functor.hasBlocks(fi.neighborSubdomainID());
560 
561  if (defined_on_elem && defined_on_neighbor)
562  {
563  const auto & value_elem = functor(elem_arg, state_arg);
564  const auto & value_neighbor = functor(neighbor_arg, state_arg);
565 
566  // This is the component of the gradient which is parallel to the line connecting
567  // the cell centers. Therefore, we can use our second order, central difference
568  // scheme to approximate it.
569  libMesh::VectorValue<T> face_gradient = (value_neighbor - value_elem) / fi.dCNMag() * fi.eCN();
570 
571  // We only need nonorthogonal correctors in 2+ dimensions
572  if (mesh.dimension() > 1)
573  {
574  // We are using an orthogonal approach for the non-orthogonal correction, for more information
575  // see Hrvoje Jasak's PhD Thesis (Imperial College, 1996)
576  libMesh::VectorValue<T> interpolated_gradient;
577 
578  // Compute the gradients in the two cells on both sides of the face
579  const auto & grad_elem =
580  greenGaussGradient(elem_arg, state_arg, functor, two_term_boundary_expansion, mesh);
581  const auto & grad_neighbor =
582  greenGaussGradient(neighbor_arg, state_arg, functor, two_term_boundary_expansion, mesh);
583 
585  interpolated_gradient,
586  grad_elem,
587  grad_neighbor,
588  fi,
589  true);
590 
591  face_gradient += interpolated_gradient - (interpolated_gradient * fi.eCN()) * fi.eCN();
592  }
593 
594  return face_gradient;
595  }
596  else if (defined_on_elem)
597  return functor.gradient(elem_arg, state_arg);
598  else if (defined_on_neighbor)
599  return functor.gradient(neighbor_arg, state_arg);
600  else
601  mooseError("The functor must be defined on one of the sides");
602 }
gc*elem+(1-gc)*neighbor
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:302
MeshBase & mesh
Moose::FunctorBase< std::array< T, N > >::GradientType greenGaussGradient(const FaceArg &face_arg, const StateArg &state_arg, const Moose::FunctorBase< std::array< T, N >> &functor, const bool two_term_boundary_expansion, const MooseMesh &mesh)
void interpolate(InterpMethod m, T &result, const T2 &value1, const T3 &value2, const FaceInfo &fi, const bool one_is_elem)
Provides interpolation of face values for non-advection-specific purposes (although it can/will still...
Definition: MathFVUtils.h:282

◆ greenGaussGradient() [3/8]

template<typename T >
TensorValue<T> Moose::FV::greenGaussGradient ( const ElemArg elem_arg,
const StateArg state_arg,
const Moose::FunctorBase< libMesh::VectorValue< T >> &  functor,
const bool  two_term_boundary_expansion,
const MooseMesh mesh 
)

Definition at line 606 of file GreenGaussGradient.h.

611 {
612  TensorValue<T> ret;
613  for (const auto i : make_range(Moose::dim))
614  {
615  VectorComponentFunctor<T> scalar_functor(functor, i);
616  const auto row_gradient =
617  greenGaussGradient(elem_arg, state_arg, scalar_functor, two_term_boundary_expansion, mesh);
618  for (const auto j : make_range(unsigned(Moose::dim)))
619  ret(i, j) = row_gradient(j);
620  }
621 
622  return ret;
623 }
static constexpr std::size_t dim
This is the dimension of all vector and tensor datastructures used in MOOSE.
Definition: Moose.h:153
Moose::FunctorBase< std::array< T, N > >::GradientType greenGaussGradient(const FaceArg &face_arg, const StateArg &state_arg, const Moose::FunctorBase< std::array< T, N >> &functor, const bool two_term_boundary_expansion, const MooseMesh &mesh)
IntRange< T > make_range(T beg, T end)

◆ greenGaussGradient() [4/8]

template<typename T >
TensorValue<T> Moose::FV::greenGaussGradient ( const FaceArg face_arg,
const StateArg state_arg,
const Moose::FunctorBase< libMesh::VectorValue< T >> &  functor,
const bool  two_term_boundary_expansion,
const MooseMesh mesh 
)

Definition at line 627 of file GreenGaussGradient.h.

632 {
633  TensorValue<T> ret;
634  for (const auto i : make_range(unsigned(Moose::dim)))
635  {
636  VectorComponentFunctor<T> scalar_functor(functor, i);
637  const auto row_gradient =
638  greenGaussGradient(face_arg, state_arg, scalar_functor, two_term_boundary_expansion, mesh);
639  for (const auto j : make_range(unsigned(Moose::dim)))
640  ret(i, j) = row_gradient(j);
641  }
642 
643  return ret;
644 }
static constexpr std::size_t dim
This is the dimension of all vector and tensor datastructures used in MOOSE.
Definition: Moose.h:153
Moose::FunctorBase< std::array< T, N > >::GradientType greenGaussGradient(const FaceArg &face_arg, const StateArg &state_arg, const Moose::FunctorBase< std::array< T, N >> &functor, const bool two_term_boundary_expansion, const MooseMesh &mesh)
IntRange< T > make_range(T beg, T end)

◆ greenGaussGradient() [5/8]

template<typename T >
Moose::FunctorBase<std::vector<T> >::GradientType Moose::FV::greenGaussGradient ( const ElemArg elem_arg,
const StateArg state_arg,
const Moose::FunctorBase< std::vector< T >> &  functor,
const bool  two_term_boundary_expansion,
const MooseMesh mesh 
)

Definition at line 648 of file GreenGaussGradient.h.

653 {
654  // Determine the size of the container
655  const auto vals = functor(elem_arg, state_arg);
657  GradientType ret(vals.size());
658  for (const auto i : index_range(ret))
659  {
660  // Note that this can be very inefficient. Within the scalar greenGaussGradient routine we're
661  // going to do value type evaluations of the array functor from scalar_functor and we will be
662  // discarding all the value type evaluations other than the one corresponding to i
663  ArrayComponentFunctor<T, FunctorBase<std::vector<T>>> scalar_functor(functor, i);
664  ret[i] =
665  greenGaussGradient(elem_arg, state_arg, scalar_functor, two_term_boundary_expansion, mesh);
666  }
667 
668  return ret;
669 }
Base class template for functor objects.
Definition: MooseFunctor.h:137
Moose::FunctorBase< std::array< T, N > >::GradientType greenGaussGradient(const FaceArg &face_arg, const StateArg &state_arg, const Moose::FunctorBase< std::array< T, N >> &functor, const bool two_term_boundary_expansion, const MooseMesh &mesh)
auto index_range(const T &sizable)

◆ greenGaussGradient() [6/8]

template<typename T >
Moose::FunctorBase<std::vector<T> >::GradientType Moose::FV::greenGaussGradient ( const FaceArg face_arg,
const StateArg state_arg,
const Moose::FunctorBase< std::vector< T >> &  functor,
const bool  two_term_boundary_expansion,
const MooseMesh mesh 
)

Definition at line 673 of file GreenGaussGradient.h.

678 {
679  // Determine the size of the container
680  const auto vals = functor(face_arg, state_arg);
682  GradientType ret(vals.size());
683  for (const auto i : index_range(ret))
684  {
685  // Note that this can be very inefficient. Within the scalar greenGaussGradient routine we're
686  // going to do value type evaluations of the array functor from scalar_functor and we will be
687  // discarding all the value type evaluations other than the one corresponding to i
688  ArrayComponentFunctor<T, FunctorBase<std::vector<T>>> scalar_functor(functor, i);
689  ret[i] =
690  greenGaussGradient(face_arg, state_arg, scalar_functor, two_term_boundary_expansion, mesh);
691  }
692 
693  return ret;
694 }
Base class template for functor objects.
Definition: MooseFunctor.h:137
Moose::FunctorBase< std::array< T, N > >::GradientType greenGaussGradient(const FaceArg &face_arg, const StateArg &state_arg, const Moose::FunctorBase< std::array< T, N >> &functor, const bool two_term_boundary_expansion, const MooseMesh &mesh)
auto index_range(const T &sizable)

◆ greenGaussGradient() [7/8]

template<typename T , std::size_t N>
Moose::FunctorBase<std::array<T, N> >::GradientType Moose::FV::greenGaussGradient ( const ElemArg elem_arg,
const StateArg state_arg,
const Moose::FunctorBase< std::array< T, N >> &  functor,
const bool  two_term_boundary_expansion,
const MooseMesh mesh 
)

Definition at line 698 of file GreenGaussGradient.h.

703 {
705  GradientType ret;
706  for (const auto i : make_range(N))
707  {
708  // Note that this can be very inefficient. Within the scalar greenGaussGradient routine we're
709  // going to do value type evaluations of the array functor from scalar_functor and we will be
710  // discarding all the value type evaluations other than the one corresponding to i
711  ArrayComponentFunctor<T, FunctorBase<std::array<T, N>>> scalar_functor(functor, i);
712  ret[i] =
713  greenGaussGradient(elem_arg, state_arg, scalar_functor, two_term_boundary_expansion, mesh);
714  }
715 
716  return ret;
717 }
Base class template for functor objects.
Definition: MooseFunctor.h:137
Moose::FunctorBase< std::array< T, N > >::GradientType greenGaussGradient(const FaceArg &face_arg, const StateArg &state_arg, const Moose::FunctorBase< std::array< T, N >> &functor, const bool two_term_boundary_expansion, const MooseMesh &mesh)
IntRange< T > make_range(T beg, T end)

◆ greenGaussGradient() [8/8]

template<typename T , std::size_t N>
Moose::FunctorBase<std::array<T, N> >::GradientType Moose::FV::greenGaussGradient ( const FaceArg face_arg,
const StateArg state_arg,
const Moose::FunctorBase< std::array< T, N >> &  functor,
const bool  two_term_boundary_expansion,
const MooseMesh mesh 
)

Definition at line 721 of file GreenGaussGradient.h.

726 {
728  GradientType ret;
729  for (const auto i : make_range(N))
730  {
731  // Note that this can be very inefficient. Within the scalar greenGaussGradient routine we're
732  // going to do value type evaluations of the array functor from scalar_functor and we will be
733  // discarding all the value type evaluations other than the one corresponding to i
734  ArrayComponentFunctor<T, FunctorBase<std::array<T, N>>> scalar_functor(functor, i);
735  ret[i] =
736  greenGaussGradient(face_arg, state_arg, scalar_functor, two_term_boundary_expansion, mesh);
737  }
738 
739  return ret;
740 }
Base class template for functor objects.
Definition: MooseFunctor.h:137
Moose::FunctorBase< std::array< T, N > >::GradientType greenGaussGradient(const FaceArg &face_arg, const StateArg &state_arg, const Moose::FunctorBase< std::array< T, N >> &functor, const bool two_term_boundary_expansion, const MooseMesh &mesh)
IntRange< T > make_range(T beg, T end)

◆ harmonicInterpolation()

template<typename T1 , typename T2 >
libMesh::CompareTypes<T1, T2>::supertype Moose::FV::harmonicInterpolation ( const T1 &  value1,
const T2 &  value2,
const FaceInfo fi,
const bool  one_is_elem 
)

Computes the harmonic mean (1/(gc/value1+(1-gc)/value2)) of Reals, RealVectorValues and RealTensorValues while accounting for the possibility that one or both of them are AD.

For tensors, we use a component-wise mean instead of the matrix-inverse based option.

Parameters
value1Reference to value1 in the (1/(gc/value1+(1-gc)/value2)) expression
value2Reference to value2 in the (1/(gc/value1+(1-gc)/value2)) expression
fiReference to the FaceInfo of the face on which the interpolation is requested
one_is_elemBoolean indicating if the interpolation weight on FaceInfo belongs to the element corresponding to value1
Returns
The interpolated value

Definition at line 177 of file MathFVUtils.h.

Referenced by interpolate().

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  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.
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)) + " & " +
203  ") must be of the same sign for harmonic interpolation");
204 #endif
205  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  {
211  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))) + " & " +
218  ") must be of the same sign for harmonic interpolation");
219 #endif
220  result(i) = 1.0 / (coeffs.first / value1(i) + coeffs.second / value2(i));
221  }
222  return result;
223  }
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  {
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))) + " & " +
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 }
void mooseWarning(Args &&... args)
Emit a warning message with the given stringified, concatenated args.
Definition: MooseError.h:336
auto raw_value(const Eigen::Map< T > &in)
Definition: EigenADReal.h:73
static constexpr std::size_t dim
This is the dimension of all vector and tensor datastructures used in MOOSE.
Definition: Moose.h:153
std::string stringify(const T &t)
conversion to string
Definition: Conversion.h:64
std::pair< T, T > interpCoeffs(const Limiter< T > &limiter, const T &phi_upwind, const T &phi_downwind, const libMesh::VectorValue< T > *const grad_phi_upwind, const libMesh::VectorValue< T > *const grad_phi_face, const Real &max_value, const Real &min_value, const FaceInfo &fi, const bool fi_elem_is_upwind)
Produce the interpolation coefficients in the equation:
Definition: MathFVUtils.h:472
IntRange< T > make_range(T beg, T end)

◆ interpCoeffs() [1/2]

template<typename T = Real>
std::pair<Real, Real> Moose::FV::interpCoeffs ( const InterpMethod  m,
const FaceInfo fi,
const bool  one_is_elem,
const T &  face_flux = 0.0 
)

Produce the interpolation coefficients in the equation:

= c_1 * {F1} + c_2 * {F2}

A couple of examples: if we are doing an average interpolation with an orthogonal regular grid, then the pair will be (0.5, 0.5). If we are doing an upwind interpolation with the velocity facing outward from the F1 element, then the pair will be (1.0, 0.0).

The template is needed in case the face flux carries derivatives in an AD setting.

Parameters
mThe interpolation method
fiThe face information
one_is_elemWhether fi.elem() == F1
face_fluxThe advecting face flux. Not relevant for an Average interpolation
Returns
a pair where the first Real is c_1 and the second Real is c_2

Definition at line 114 of file MathFVUtils.h.

Referenced by LinearFVAdvection::computeElemMatrixContribution(), LinearFVAnisotropicDiffusion::computeFluxRHSContribution(), LinearFVDiffusion::computeFluxRHSContribution(), LinearFVAdvection::computeNeighborMatrixContribution(), containerInterpolate(), harmonicInterpolation(), interpCoeffsAndAdvected(), interpolate(), linearInterpolation(), and skewCorrectedLinearInterpolation().

118 {
119  switch (m)
120  {
121  case InterpMethod::Average:
122  case InterpMethod::SkewCorrectedAverage:
123  {
124  if (one_is_elem)
125  return std::make_pair(fi.gC(), 1. - fi.gC());
126  else
127  return std::make_pair(1. - fi.gC(), fi.gC());
128  }
129 
130  case InterpMethod::Upwind:
131  {
132  if ((face_flux > 0) == one_is_elem)
133  return std::make_pair(1., 0.);
134  else
135  return std::make_pair(0., 1.);
136  }
137 
138  default:
139  mooseError("Unrecognized interpolation method");
140  }
141 }
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:302
Real gC() const
Return the geometric weighting factor.
Definition: FaceInfo.h:132

◆ interpCoeffs() [2/2]

template<typename T >
std::pair<T, T> Moose::FV::interpCoeffs ( const Limiter< T > &  limiter,
const T &  phi_upwind,
const T &  phi_downwind,
const libMesh::VectorValue< T > *const  grad_phi_upwind,
const libMesh::VectorValue< T > *const  grad_phi_face,
const Real max_value,
const Real min_value,
const FaceInfo fi,
const bool  fi_elem_is_upwind 
)

Produce the interpolation coefficients in the equation:

= c_upwind * {upwind} + c_downwind * {downwind}

A couple of examples: if we are doing an average interpolation with an orthogonal regular grid, then the pair will be (0.5, 0.5). If we are doing an upwind interpolation then the pair will be (1.0, 0.0).

Returns
a pair where the first Real is c_upwind and the second Real is c_downwind

Definition at line 472 of file MathFVUtils.h.

481 {
482  // Using beta, w_f, g nomenclature from Greenshields
483  const auto beta = limiter(phi_upwind,
484  phi_downwind,
485  grad_phi_upwind,
486  grad_phi_face,
487  fi_elem_is_upwind ? fi.dCN() : Point(-fi.dCN()),
488  max_value,
489  min_value,
490  &fi,
491  fi_elem_is_upwind);
492 
493  const auto w_f = fi_elem_is_upwind ? fi.gC() : (1. - fi.gC());
494 
495  const auto g = beta * (1. - w_f);
496 
497  return std::make_pair(1. - g, g);
498 }
Real gC() const
Return the geometric weighting factor.
Definition: FaceInfo.h:132
const Point & dCN() const
Definition: FaceInfo.h:138

◆ interpCoeffsAndAdvected()

template<typename T , FunctorEvaluationKind FEK = FunctorEvaluationKind::Value, typename Enable = typename std::enable_if<libMesh::ScalarTraits<T>::value>::type>
std::pair<std::pair<T, T>, std::pair<T, T> > Moose::FV::interpCoeffsAndAdvected ( const FunctorBase< T > &  functor,
const FaceArg face,
const StateArg time 
)

This function interpolates values using a specified limiter and face argument.

It evaluates the values at upwind and downwind locations and computes interpolation coefficients and advected values.

Template Parameters
TThe data type for the values being interpolated. This is typically a scalar type.
FEKEnumeration type FunctorEvaluationKind with a default value of FunctorEvaluationKind::Value. This determines the kind of evaluation the functor will perform.
EnableA type trait used for SFINAE (Substitution Failure Is Not An Error) to ensure that T is a valid scalar type as determined by ScalarTraits<T>.
Parameters
functorAn object of a functor class derived from FunctorBase<T>. This object provides the genericEvaluate method to compute the value at a given element and time.
faceAn argument representing the face in the computational domain. The face provides access to neighboring elements and limiter type.
timeAn argument representing the state or time at which the evaluation is performed.
Returns
std::pair<std::pair<T, T>, std::pair<T, T>> A pair of pairs.
  • The first pair corresponds to the interpolation coefficients, with the first value corresponding to the face information element and the second to the face information neighbor. This pair should sum to unity.
  • The second pair corresponds to the face information functor element value and neighbor value.

Usage: This function is used for interpolating values at faces in a finite volume method, ensuring that the interpolation adheres to the constraints imposed by the limiter.

Definition at line 785 of file MathFVUtils.h.

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  constexpr FunctorEvaluationKind GFEK = FunctorGradientEvaluationKind<FEK>::value;
793  typedef typename FunctorBase<T>::GradientType GradientType;
794  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  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  const auto upwind_arg = face.elem_is_upwind ? face.makeElem() : face.makeNeighbor();
803  const auto downwind_arg = face.elem_is_upwind ? face.makeNeighbor() : face.makeElem();
804 
805  // Evaluate the functor at the upwind and downwind locations
806  auto phi_upwind = functor.template genericEvaluate<FEK>(upwind_arg, time);
807  auto phi_downwind = functor.template genericEvaluate<FEK>(downwind_arg, time);
808 
809  // Initialize the interpolation coefficients
810  std::pair<T, T> interp_coeffs;
811 
812  // Compute interpolation coefficients for Upwind or CentralDifference limiters
813  if (face.limiter_type == LimiterType::Upwind ||
814  face.limiter_type == LimiterType::CentralDifference)
815  interp_coeffs = interpCoeffs(*limiter,
816  phi_upwind,
817  phi_downwind,
818  &zero,
819  &zero,
822  *face.fi,
823  face.elem_is_upwind);
824  else
825  {
826  // Determine the time argument for the limiter
827  auto * time_arg = face.state_limiter;
828  if (!time_arg)
829  {
831  time_arg = &temp_state;
832  }
833 
835 
836  // Compute min-max values for min-max limiters
837  if (face.limiter_type == LimiterType::Venkatakrishnan || face.limiter_type == LimiterType::SOU)
838  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  const auto grad_phi_upwind = functor.template genericEvaluate<GFEK>(upwind_arg, *time_arg);
842  const auto grad_phi_face = functor.template genericEvaluate<GFEK>(face, *time_arg);
843 
844  // Compute the interpolation coefficients using the specified limiter
845  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  *face.fi,
853  face.elem_is_upwind);
854  }
855 
856  // Return the interpolation coefficients and advected values
857  if (face.elem_is_upwind)
858  return std::make_pair(std::move(interp_coeffs),
859  std::make_pair(std::move(phi_upwind), std::move(phi_downwind)));
860  else
861  return std::make_pair(
862  std::make_pair(std::move(interp_coeffs.second), std::move(interp_coeffs.first)),
863  std::make_pair(std::move(phi_downwind), std::move(phi_upwind)));
864 }
const Number zero
auto max(const L &left, const R &right)
FunctorEvaluationKind
An enumeration of possible functor evaluation kinds.
Definition: MooseFunctor.h:36
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
std::pair< T, T > interpCoeffs(const Limiter< T > &limiter, const T &phi_upwind, const T &phi_downwind, const libMesh::VectorValue< T > *const grad_phi_upwind, const libMesh::VectorValue< T > *const grad_phi_face, const Real &max_value, const Real &min_value, const FaceInfo &fi, const bool fi_elem_is_upwind)
Produce the interpolation coefficients in the equation:
Definition: MathFVUtils.h:472
State argument for evaluating functors.
std::pair< Real, Real > computeMinMaxValue(const FunctorBase< VectorValue< T >> &functor, const FaceArg &face, const StateArg &time, const unsigned int &component)
This function calculates the minimum and maximum values of a specified component within a two-cell st...
Definition: MathFVUtils.h:711
auto min(const L &left, const R &right)

◆ interpolate() [1/9]

template<typename T , typename T2 , typename T3 >
void Moose::FV::interpolate ( InterpMethod  m,
T &  result,
const T2 &  value1,
const T3 &  value2,
const FaceInfo fi,
const bool  one_is_elem 
)

Provides interpolation of face values for non-advection-specific purposes (although it can/will still be used by advective kernels sometimes).

The interpolated value is stored in result. This should be called when a face value needs to be computed from two neighboring cells/elements. value1 and value2 represent the cell property/values from which to compute the face value. The one_is_elem parameter indicates whether value1 corresponds to the FaceInfo elem value; else it corresponds to the FaceInfo neighbor value

Definition at line 282 of file MathFVUtils.h.

Referenced by InterfaceDiffusiveFluxIntegralTempl< is_ad >::computeQpIntegral(), FVDiffusionInterface::computeQpResidual(), FVOneVarDiffusionInterface::computeQpResidual(), FVOrthogonalDiffusion::computeQpResidual(), FVAnisotropicDiffusion::computeQpResidual(), FVDiffusion::computeQpResidual(), PiecewiseByBlockLambdaFunctor< T >::evaluate(), MooseLinearVariableFV< Real >::evaluate(), MooseVariableFV< Real >::evaluate(), greenGaussGradient(), and interpolate().

288 {
289  switch (m)
290  {
291  case InterpMethod::Average:
292  case InterpMethod::SkewCorrectedAverage:
293  result = linearInterpolation(value1, value2, fi, one_is_elem);
294  break;
295  case InterpMethod::HarmonicAverage:
296  result = harmonicInterpolation(value1, value2, fi, one_is_elem);
297  break;
298  default:
299  mooseError("unsupported interpolation method for interpolate() function");
300  }
301 }
libMesh::CompareTypes< T1, T2 >::supertype harmonicInterpolation(const T1 &value1, const T2 &value2, const FaceInfo &fi, const bool one_is_elem)
Computes the harmonic mean (1/(gc/value1+(1-gc)/value2)) of Reals, RealVectorValues and RealTensorVal...
Definition: MathFVUtils.h:177
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:302
T linearInterpolation(const FunctorBase< T > &functor, const FaceArg &face, const StateArg &time)
perform a possibly skew-corrected linear interpolation by evaluating the supplied functor with the pr...
Definition: MathFVUtils.h:309

◆ interpolate() [2/9]

template<typename T1 , typename T2 , typename T3 , template< typename > class Vector1, template< typename > class Vector2>
void Moose::FV::interpolate ( InterpMethod  m,
Vector1< T1 > &  result,
const T2 &  fi_elem_advected,
const T2 &  fi_neighbor_advected,
const Vector2< T3 > &  fi_elem_advector,
const Vector2< T3 > &  fi_neighbor_advector,
const FaceInfo fi 
)

Computes the product of the advected and the advector based on the given interpolation method.

Definition at line 355 of file MathFVUtils.h.

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 }
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:302
auto raw_value(const Eigen::Map< T > &in)
Definition: EigenADReal.h:73
T linearInterpolation(const FunctorBase< T > &functor, const FaceArg &face, const StateArg &time)
perform a possibly skew-corrected linear interpolation by evaluating the supplied functor with the pr...
Definition: MathFVUtils.h:309
const Point & normal() const
Returns the unit normal vector for the face oriented outward from the face&#39;s elem element...
Definition: FaceInfo.h:68
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

◆ interpolate() [3/9]

template<typename T , typename T2 , typename T3 , typename Vector >
void Moose::FV::interpolate ( InterpMethod  m,
T &  result,
const T2 &  value1,
const T3 &  value2,
const Vector &  advector,
const FaceInfo fi,
const bool  one_is_elem 
)

Provides interpolation of face values for advective flux kernels.

This should be called by advective kernels when a face value is needed from two neighboring cells/elements. The interpolated value is stored in result. value1 and value2 represent the two neighboring advected cell property/values. advector represents the vector quantity at the face that is doing the advecting (e.g. the flow velocity at the face); this value often will have been computed using a call to the non-advective interpolate function. The one_is_elem parameter indicates whether value1 corresponds to the FaceInfo elem value; else it corresponds to the FaceInfo neighbor value

Definition at line 402 of file MathFVUtils.h.

409 {
410  const auto coeffs = interpCoeffs(m, fi, one_is_elem, advector * fi.normal());
411  result = coeffs.first * value1 + coeffs.second * value2;
412 }
const Point & normal() const
Returns the unit normal vector for the face oriented outward from the face&#39;s elem element...
Definition: FaceInfo.h:68
std::pair< T, T > interpCoeffs(const Limiter< T > &limiter, const T &phi_upwind, const T &phi_downwind, const libMesh::VectorValue< T > *const grad_phi_upwind, const libMesh::VectorValue< T > *const grad_phi_face, const Real &max_value, const Real &min_value, const FaceInfo &fi, const bool fi_elem_is_upwind)
Produce the interpolation coefficients in the equation:
Definition: MathFVUtils.h:472

◆ interpolate() [4/9]

template<typename Scalar , typename Vector , typename Enable = typename std::enable_if<libMesh::ScalarTraits<Scalar>::value>::type>
Scalar Moose::FV::interpolate ( const Limiter< Scalar > &  limiter,
const Scalar &  phi_upwind,
const Scalar &  phi_downwind,
const Vector *const  grad_phi_upwind,
const FaceInfo fi,
const bool  fi_elem_is_upwind 
)

Interpolates with a limiter.

Definition at line 507 of file MathFVUtils.h.

513 {
514  auto pr =
515  interpCoeffs(limiter,
516  phi_upwind,
517  phi_downwind,
518  grad_phi_upwind,
519  /*grad_phi_face*/ static_cast<const libMesh::VectorValue<Scalar> *>(nullptr),
520  /* max_value */ std::numeric_limits<Real>::max(),
521  /* min_value */ std::numeric_limits<Real>::min(),
522  fi,
523  fi_elem_is_upwind);
524  return pr.first * phi_upwind + pr.second * phi_downwind;
525 }
auto max(const L &left, const R &right)
std::pair< T, T > interpCoeffs(const Limiter< T > &limiter, const T &phi_upwind, const T &phi_downwind, const libMesh::VectorValue< T > *const grad_phi_upwind, const libMesh::VectorValue< T > *const grad_phi_face, const Real &max_value, const Real &min_value, const FaceInfo &fi, const bool fi_elem_is_upwind)
Produce the interpolation coefficients in the equation:
Definition: MathFVUtils.h:472
auto min(const L &left, const R &right)

◆ interpolate() [5/9]

template<typename Limiter , typename T , typename Tensor >
libMesh::VectorValue<T> Moose::FV::interpolate ( const Limiter limiter,
const TypeVector< T > &  phi_upwind,
const TypeVector< T > &  phi_downwind,
const Tensor *const  grad_phi_upwind,
const FaceInfo fi,
const bool  fi_elem_is_upwind 
)

Vector overload.

Definition at line 532 of file MathFVUtils.h.

538 {
539  mooseAssert(limiter.constant() || grad_phi_upwind,
540  "Non-null gradient only supported for constant limiters.");
541 
542  const libMesh::VectorValue<T> * const gradient_example = nullptr;
544  for (const auto i : make_range(unsigned(LIBMESH_DIM)))
545  {
546  if (grad_phi_upwind)
547  {
548  const libMesh::VectorValue<T> gradient = grad_phi_upwind->row(i);
549  ret(i) =
550  interpolate(limiter, phi_upwind(i), phi_downwind(i), &gradient, fi, fi_elem_is_upwind);
551  }
552  else
553  ret(i) = interpolate(
554  limiter, phi_upwind(i), phi_downwind(i), gradient_example, fi, fi_elem_is_upwind);
555  }
556 
557  return ret;
558 }
std::array< T, N > interpolate(const FunctorBase< std::array< T, N >> &functor, const FaceArg &face, const StateArg &time)
Definition: MathFVUtils.h:1157
IntRange< T > make_range(T beg, T end)

◆ interpolate() [6/9]

template<typename T , FunctorEvaluationKind FEK = FunctorEvaluationKind::Value, typename Enable = typename std::enable_if<libMesh::ScalarTraits<T>::value>::type>
T Moose::FV::interpolate ( const FunctorBase< T > &  functor,
const FaceArg face,
const StateArg time 
)

This function interpolates values at faces in a computational grid using a specified functor, face argument, and evaluation kind.

It handles different limiter types and performs interpolation accordingly.

Template Parameters
TThe data type for the values being interpolated. This is typically a scalar type.
FEKEnumeration type FunctorEvaluationKind with a default value of FunctorEvaluationKind::Value. This determines the kind of evaluation the functor will perform.
EnableA type trait used for SFINAE (Substitution Failure Is Not An Error) to ensure that T is a valid scalar type as determined by ScalarTraits<T>.
Parameters
functorAn object of a functor class derived from FunctorBase<T>. This object provides the genericEvaluate method to compute the value at a given element and time.
faceAn argument representing the face in the computational domain. The face provides access to neighboring elements and limiter type.
timeAn argument representing the state or time at which the evaluation is performed.
Returns
T The interpolated value at the face.

Usage: This function is used for interpolating values at faces in a finite volume method, ensuring that the interpolation adheres to the constraints imposed by the limiter.

Definition at line 893 of file MathFVUtils.h.

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  if (face.limiter_type == LimiterType::CentralDifference)
902  return linearInterpolation<T, FEK>(functor, face, time);
903 
904  if (face.limiter_type == LimiterType::Upwind ||
905  face.limiter_type == LimiterType::CentralDifference)
906  {
907  // Compute interpolation coefficients and advected values
908  const auto [interp_coeffs, advected] = interpCoeffsAndAdvected<T, FEK>(functor, face, time);
909  // Return the interpolated value
910  return interp_coeffs.first * advected.first + interp_coeffs.second * advected.second;
911  }
912  else
913  {
914  // Determine the gradient evaluation kind
915  constexpr FunctorEvaluationKind GFEK = FunctorGradientEvaluationKind<FEK>::value;
916  typedef typename FunctorBase<T>::GradientType GradientType;
917  static const GradientType zero(0);
918  const auto & limiter =
919  Limiter<typename LimiterValueType<T>::value_type>::build(face.limiter_type);
920 
921  // Determine upwind and downwind arguments based on the face element
922  const auto & upwind_arg = face.elem_is_upwind ? face.makeElem() : face.makeNeighbor();
923  const auto & downwind_arg = face.elem_is_upwind ? face.makeNeighbor() : face.makeElem();
924  const auto & phi_upwind = functor.template genericEvaluate<FEK>(upwind_arg, time);
925  const auto & phi_downwind = functor.template genericEvaluate<FEK>(downwind_arg, time);
926 
927  // Determine the time argument for the limiter
928  auto * time_arg = face.state_limiter;
929  if (!time_arg)
930  {
932  time_arg = &temp_state;
933  }
934 
935  // Initialize min-max values
937  if (face.limiter_type == LimiterType::Venkatakrishnan || face.limiter_type == LimiterType::SOU)
938  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  const auto & grad_phi_upwind = functor.template genericEvaluate<GFEK>(upwind_arg, *time_arg);
942  const auto & grad_phi_face = functor.template genericEvaluate<GFEK>(face, *time_arg);
943 
944  // Perform full limited interpolation and return the interpolated value
945  return fullLimitedInterpolation(*limiter,
946  phi_upwind,
947  phi_downwind,
948  &grad_phi_upwind,
949  &grad_phi_face,
950  max_value,
951  min_value,
952  face);
953  }
954 }
const Number zero
auto max(const L &left, const R &right)
FunctorEvaluationKind
An enumeration of possible functor evaluation kinds.
Definition: MooseFunctor.h:36
T fullLimitedInterpolation(const Limiter< T > &limiter, const T &phi_upwind, const T &phi_downwind, const VectorValue< T > *const grad_phi_upwind, const VectorValue< T > *const grad_phi_face, const Real &max_value, const Real &min_value, const FaceArg &face)
This function performs a full limited interpolation of a variable, taking into account the values and...
Definition: MathFVUtils.h:580
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
State argument for evaluating functors.
std::pair< Real, Real > computeMinMaxValue(const FunctorBase< VectorValue< T >> &functor, const FaceArg &face, const StateArg &time, const unsigned int &component)
This function calculates the minimum and maximum values of a specified component within a two-cell st...
Definition: MathFVUtils.h:711
auto min(const L &left, const R &right)

◆ interpolate() [7/9]

template<typename T >
libMesh::VectorValue<T> Moose::FV::interpolate ( const FunctorBase< libMesh::VectorValue< T >> &  functor,
const FaceArg face,
const StateArg time 
)

This function interpolates vector values at faces in a computational grid using a specified functor, face argument, and limiter type.

It handles different limiter types and performs interpolation accordingly.

Template Parameters
TThe data type for the vector values being interpolated. This is typically a scalar type.
Parameters
functorAn object of a functor class derived from FunctorBase<VectorValue<T>>. This object provides the operator() method to compute the value at a given element and time.
faceAn argument representing the face in the computational domain. The face provides access to neighboring elements and limiter type.
timeAn argument representing the state or time at which the evaluation is performed.
Returns
VectorValue<T> The interpolated vector value at the face.

Usage: This function is used for interpolating vector values at faces in a finite volume method, ensuring that the interpolation adheres to the constraints imposed by the limiter.

Definition at line 978 of file MathFVUtils.h.

981 {
982  // Define a zero gradient vector for initialization
983  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  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  const auto upwind_arg = face.elem_is_upwind ? face.makeElem() : face.makeNeighbor();
992  const auto downwind_arg = face.elem_is_upwind ? face.makeNeighbor() : face.makeElem();
993  auto phi_upwind = functor(upwind_arg, time);
994  auto phi_downwind = functor(downwind_arg, time);
995 
996  // Initialize the return vector value
998  T coeff_upwind, coeff_downwind;
999 
1000  // Compute interpolation coefficients and advected values for Upwind or CentralDifference limiters
1001  if (face.limiter_type == LimiterType::Upwind ||
1002  face.limiter_type == LimiterType::CentralDifference)
1003  {
1004  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
1005  {
1006  const auto &component_upwind = phi_upwind(i), component_downwind = phi_downwind(i);
1007  std::tie(coeff_upwind, coeff_downwind) = interpCoeffs(*limiter,
1008  component_upwind,
1009  component_downwind,
1010  &grad_zero,
1011  &grad_zero,
1014  *face.fi,
1015  face.elem_is_upwind);
1016  ret(i) = coeff_upwind * component_upwind + coeff_downwind * component_downwind;
1017  }
1018  }
1019  else
1020  {
1021  // Determine the time argument for the limiter
1022  auto * time_arg = face.state_limiter;
1023  if (!time_arg)
1024  {
1026  time_arg = &temp_state;
1027  }
1028 
1029  // Evaluate the gradient of the functor at the upwind and downwind locations
1030  const auto & grad_phi_upwind = functor.gradient(upwind_arg, *time_arg);
1031  const auto & grad_phi_downwind = functor.gradient(downwind_arg, *time_arg);
1032 
1033  const auto coeffs = interpCoeffs(InterpMethod::Average, *face.fi, face.elem_is_upwind);
1034 
1035  // Compute interpolation coefficients and advected values for each component
1036  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
1037  {
1038  const auto &component_upwind = phi_upwind(i), component_downwind = phi_downwind(i);
1039  const libMesh::VectorValue<T> &grad_upwind = grad_phi_upwind.row(i),
1040  grad_face = coeffs.first * grad_phi_upwind.row(i) +
1041  coeffs.second * grad_phi_downwind.row(i);
1042 
1043  // Initialize min-max values
1045  if (face.limiter_type == LimiterType::Venkatakrishnan ||
1046  face.limiter_type == LimiterType::SOU)
1047  std::tie(max_value, min_value) = computeMinMaxValue(functor, face, *time_arg, i);
1048 
1049  // Perform full limited interpolation for the component
1050  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  }
1060 
1061  // Return the interpolated vector value
1062  return ret;
1063 }
auto max(const L &left, const R &right)
T fullLimitedInterpolation(const Limiter< T > &limiter, const T &phi_upwind, const T &phi_downwind, const VectorValue< T > *const grad_phi_upwind, const VectorValue< T > *const grad_phi_face, const Real &max_value, const Real &min_value, const FaceArg &face)
This function performs a full limited interpolation of a variable, taking into account the values and...
Definition: MathFVUtils.h:580
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
std::pair< T, T > interpCoeffs(const Limiter< T > &limiter, const T &phi_upwind, const T &phi_downwind, const libMesh::VectorValue< T > *const grad_phi_upwind, const libMesh::VectorValue< T > *const grad_phi_face, const Real &max_value, const Real &min_value, const FaceInfo &fi, const bool fi_elem_is_upwind)
Produce the interpolation coefficients in the equation:
Definition: MathFVUtils.h:472
State argument for evaluating functors.
std::pair< Real, Real > computeMinMaxValue(const FunctorBase< VectorValue< T >> &functor, const FaceArg &face, const StateArg &time, const unsigned int &component)
This function calculates the minimum and maximum values of a specified component within a two-cell st...
Definition: MathFVUtils.h:711
auto min(const L &left, const R &right)

◆ interpolate() [8/9]

template<typename T >
std::vector<T> Moose::FV::interpolate ( const FunctorBase< std::vector< T >> &  functor,
const FaceArg face,
const StateArg time 
)

Definition at line 1148 of file MathFVUtils.h.

1151 {
1152  return containerInterpolate(functor, face, time);
1153 }
T containerInterpolate(const FunctorBase< T > &functor, const FaceArg &face, const StateArg &time)
This function interpolates container values at faces in a computational grid using a specified functo...
Definition: MathFVUtils.h:1087

◆ interpolate() [9/9]

template<typename T , std::size_t N>
std::array<T, N> Moose::FV::interpolate ( const FunctorBase< std::array< T, N >> &  functor,
const FaceArg face,
const StateArg time 
)

Definition at line 1157 of file MathFVUtils.h.

1160 {
1161  return containerInterpolate(functor, face, time);
1162 }
T containerInterpolate(const FunctorBase< T > &functor, const FaceArg &face, const StateArg &time)
This function interpolates container values at faces in a computational grid using a specified functo...
Definition: MathFVUtils.h:1087

◆ interpolationMethods()

MooseEnum Moose::FV::interpolationMethods ( )

Returns an enum with all the currently supported interpolation methods and the current default for FV: first-order upwind.

Returns
MooseEnum with all the face interpolation methods supported

Definition at line 61 of file MathFVUtils.C.

Referenced by advectedInterpolationParameter().

62 {
63  return MooseEnum("average upwind sou min_mod vanLeer quick venkatakrishnan skewness-corrected",
64  "upwind");
65 }
This is a "smart" enum class intended to replace many of the shortcomings in the C++ enum type It sho...
Definition: MooseEnum.h:33

◆ limiterType()

LimiterType Moose::FV::limiterType ( InterpMethod  interp_method)

Return the limiter type associated with the supplied interpolation method.

Definition at line 63 of file Limiter.C.

Referenced by FVAdvection::computeQpResidual(), FVMatAdvection::computeQpResidual(), and FVDivergence::computeQpResidual().

64 {
65  switch (interp_method)
66  {
67  case InterpMethod::Average:
68  case InterpMethod::SkewCorrectedAverage:
69  return LimiterType::CentralDifference;
70 
71  case InterpMethod::Upwind:
72  return LimiterType::Upwind;
73 
74  case InterpMethod::VanLeer:
75  return LimiterType::VanLeer;
76 
77  case InterpMethod::MinMod:
78  return LimiterType::MinMod;
79 
80  case InterpMethod::SOU:
81  return LimiterType::SOU;
82 
83  case InterpMethod::QUICK:
84  return LimiterType::QUICK;
85 
86  case InterpMethod::Venkatakrishnan:
87  return LimiterType::Venkatakrishnan;
88 
89  default:
90  mooseError("Unrecognized interpolation method type.");
91  }
92 }
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:302

◆ linearInterpolation() [1/2]

template<typename T , typename T2 >
libMesh::CompareTypes<T, T2>::supertype Moose::FV::linearInterpolation ( const T &  value1,
const T2 &  value2,
const FaceInfo fi,
const bool  one_is_elem,
const InterpMethod  interp_method = InterpMethod::Average 
)

A simple linear interpolation of values between cell centers to a cell face.

The one_is_elem parameter indicates whether value1 corresponds to the FaceInfo elem value; else it corresponds to the FaceInfo neighbor value

Definition at line 150 of file MathFVUtils.h.

Referenced by MooseLinearVariableFV< Real >::gradSln(), interpolate(), linearInterpolation(), ComputeLinearFVGreenGaussGradientFaceThread::operator()(), and MooseVariableFV< Real >::uncorrectedAdGradSln().

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  const auto coeffs = interpCoeffs(interp_method, fi, one_is_elem);
161  return coeffs.first * value1 + coeffs.second * value2;
162 }
std::pair< T, T > interpCoeffs(const Limiter< T > &limiter, const T &phi_upwind, const T &phi_downwind, const libMesh::VectorValue< T > *const grad_phi_upwind, const libMesh::VectorValue< T > *const grad_phi_face, const Real &max_value, const Real &min_value, const FaceInfo &fi, const bool fi_elem_is_upwind)
Produce the interpolation coefficients in the equation:
Definition: MathFVUtils.h:472

◆ linearInterpolation() [2/2]

template<typename T , FunctorEvaluationKind FEK = FunctorEvaluationKind::Value>
T Moose::FV::linearInterpolation ( const FunctorBase< T > &  functor,
const FaceArg face,
const StateArg time 
)

perform a possibly skew-corrected linear interpolation by evaluating the supplied functor with the provided functor face argument

Definition at line 309 of file MathFVUtils.h.

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  constexpr FunctorEvaluationKind GFEK = FunctorGradientEvaluationKind<FEK>::value;
320 
321  const auto elem_arg = face.makeElem();
322  const auto neighbor_arg = face.makeNeighbor();
323 
324  const auto elem_eval = functor.template genericEvaluate<FEK>(elem_arg, time);
325  const auto neighbor_eval = functor.template genericEvaluate<FEK>(neighbor_arg, time);
326 
327  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  FaceArg new_face(face);
334  new_face.correct_skewness = false;
335  const auto surface_gradient = functor.template genericEvaluate<GFEK>(new_face, time);
336 
338  elem_eval, neighbor_eval, surface_gradient, *face.fi, true);
339  }
340  else
341  return linearInterpolation(elem_eval, neighbor_eval, *face.fi, true);
342 }
FunctorEvaluationKind
An enumeration of possible functor evaluation kinds.
Definition: MooseFunctor.h:36
libMesh::CompareTypes< T, T2 >::supertype skewCorrectedLinearInterpolation(const T &value1, const T2 &value2, const T3 &face_gradient, const FaceInfo &fi, const bool one_is_elem)
Linear interpolation with skewness correction using the face gradient.
Definition: MathFVUtils.h:261
T linearInterpolation(const FunctorBase< T > &functor, const FaceArg &face, const StateArg &time)
perform a possibly skew-corrected linear interpolation by evaluating the supplied functor with the pr...
Definition: MathFVUtils.h:309

◆ loopOverElemFaceInfo()

template<typename ActionFunctor >
void Moose::FV::loopOverElemFaceInfo ( const Elem &  elem,
const MooseMesh mesh,
ActionFunctor &  act,
const Moose::CoordinateSystemType  coord_type,
const unsigned int  rz_radial_coord = libMesh::invalid_uint 
)

Definition at line 42 of file FVUtils.h.

Referenced by greenGaussGradient().

47 {
48  mooseAssert(elem.active(), "We should never call this method with an inactive element");
49 
50  for (const auto side : elem.side_index_range())
51  {
52  const Elem * const candidate_neighbor = elem.neighbor_ptr(side);
53 
54  bool elem_has_info = elemHasFaceInfo(elem, candidate_neighbor);
55 
56  std::set<const Elem *> neighbors;
57 
58  const bool inactive_neighbor_detected =
59  candidate_neighbor ? !candidate_neighbor->active() : false;
60 
61  // See MooseMesh::buildFaceInfo for corresponding checks/additions of FaceInfo
62  if (inactive_neighbor_detected)
63  {
64  // We must be next to an element that has been refined
65  mooseAssert(candidate_neighbor->has_children(), "We should have children");
66 
67  const auto candidate_neighbor_side = candidate_neighbor->which_neighbor_am_i(&elem);
68 
69  for (const auto child_num : make_range(candidate_neighbor->n_children()))
70  if (candidate_neighbor->is_child_on_side(child_num, candidate_neighbor_side))
71  {
72  const Elem * const child = candidate_neighbor->child_ptr(child_num);
73  mooseAssert(child->level() - elem.level() == 1, "The math doesn't work out here.");
74  mooseAssert(child->has_neighbor(&elem), "Elem should be a neighbor of this child.");
75  mooseAssert(child->active(),
76  "We shouldn't have greater than a face mismatch level of one");
77  neighbors.insert(child);
78  }
79  }
80  else
81  neighbors.insert(candidate_neighbor);
82 
83  for (const Elem * const neighbor : neighbors)
84  {
85  const FaceInfo * const fi =
86  elem_has_info ? mesh.faceInfo(&elem, side)
87  : mesh.faceInfo(neighbor, neighbor->which_neighbor_am_i(&elem));
88 
89  mooseAssert(fi, "We should have found a FaceInfo");
90  mooseAssert(elem_has_info ? &elem == &fi->elem() : &elem == fi->neighborPtr(),
91  "Doesn't seem like we understand how this FaceInfo thing is working");
92  if (neighbor)
93  {
94  mooseAssert(neighbor != libMesh::remote_elem,
95  "Remote element detected. This indicates that you have insufficient geometric "
96  "ghosting. Please contact your application developers.");
97  mooseAssert(elem_has_info ? neighbor == fi->neighborPtr() : neighbor == &fi->elem(),
98  "Doesn't seem like we understand how this FaceInfo thing is working");
99  }
100 
101  const Point elem_normal = elem_has_info ? fi->normal() : Point(-fi->normal());
102 
103  Real coord;
104  MooseMeshUtils::coordTransformFactor(fi->faceCentroid(), coord, coord_type, rz_radial_coord);
105 
106  const Point surface_vector = elem_normal * fi->faceArea() * coord;
107 
108  act(elem, neighbor, fi, surface_vector, coord, elem_has_info);
109  }
110  }
111 }
bool elemHasFaceInfo(const Elem &elem, const Elem *const neighbor)
This function infers based on elements if the faceinfo between them belongs to the element or not...
Definition: FVUtils.C:21
void coordTransformFactor(const P &point, C &factor, const Moose::CoordinateSystemType coord_type, const unsigned int rz_radial_coord=libMesh::invalid_uint)
compute a coordinate transformation factor
const Elem & elem() const
Definition: FaceInfo.h:81
const Point & faceCentroid() const
Returns the coordinates of the face centroid.
Definition: FaceInfo.h:71
MeshBase & mesh
Real faceArea() const
Returns the face area of face id.
Definition: FaceInfo.h:60
This data structure is used to store geometric and variable related metadata about each cell face in ...
Definition: FaceInfo.h:36
const Elem * neighborPtr() const
Definition: FaceInfo.h:84
const Point & normal() const
Returns the unit normal vector for the face oriented outward from the face&#39;s elem element...
Definition: FaceInfo.h:68
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
IntRange< T > make_range(T beg, T end)
const RemoteElem * remote_elem

◆ moose_limiter_type()

const MooseEnum Moose::FV::moose_limiter_type ( )

◆ onBoundary() [1/2]

template<typename SubdomainRestrictable >
bool Moose::FV::onBoundary ( const SubdomainRestrictable &  obj,
const FaceInfo fi 
)

Return whether the supplied face is on a boundary of the object's execution.

Definition at line 1169 of file MathFVUtils.h.

Referenced by FVFluxKernel::onBoundary(), ThreadedElementLoopBase< ConstElemPointerRange >::operator()(), and ThreadedFaceLoop< RangeType >::operator()().

1170 {
1171  const bool internal = fi.neighborPtr() && obj.hasBlocks(fi.elemSubdomainID()) &&
1172  obj.hasBlocks(fi.neighborSubdomainID());
1173  return !internal;
1174 }
SubdomainID elemSubdomainID() const
Definition: FaceInfo.h:102
const Elem * neighborPtr() const
Definition: FaceInfo.h:84
SubdomainID neighborSubdomainID() const
Definition: FaceInfo.h:252

◆ onBoundary() [2/2]

bool Moose::FV::onBoundary ( const std::set< SubdomainID > &  subs,
const FaceInfo fi 
)

Determine whether the passed-in face is on the boundary of an object that lives on the provided subdomains.

Note that if the subdomain set is empty we consider that to mean that the object has no block restriction and lives on the entire mesh

Definition at line 28 of file MathFVUtils.C.

29 {
30  if (!fi.neighborPtr())
31  // We're on the exterior boundary
32  return true;
33 
34  if (subs.empty())
35  // The face is internal and our functor lives on all subdomains
36  return false;
37 
38  const auto sub_count =
39  subs.count(fi.elem().subdomain_id()) + subs.count(fi.neighbor().subdomain_id());
40 
41  switch (sub_count)
42  {
43  case 0:
44  mooseError("We should not be calling isExtrapolatedBoundaryFace on a functor that doesn't "
45  "live on either of the face information's neighboring elements");
46 
47  case 1:
48  // We only live on one of the subs
49  return true;
50 
51  case 2:
52  // We live on both of the subs
53  return false;
54 
55  default:
56  mooseError("There should be no other sub_count options");
57  }
58 }
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:302
const Elem & elem() const
Definition: FaceInfo.h:81
const Elem * neighborPtr() const
Definition: FaceInfo.h:84
const Elem & neighbor() const
Definition: FaceInfo.h:216

◆ rF()

template<typename Scalar , typename Vector >
Scalar Moose::FV::rF ( const Scalar &  phiC,
const Scalar &  phiD,
const Vector &  gradC,
const RealVectorValue dCD 
)

From Moukalled 12.30.

r_f = (phiC - phiU) / (phiD - phiC)

However, this formula is only clear on grids where upwind locations can be readily determined, which is not the case for unstructured meshes. So we leverage a virtual upwind location and Moukalled 12.65

phiD - phiU = 2 * grad(phi)_C * dCD -> phiU = phiD - 2 * grad(phi)_C * dCD

Combining the two equations and performing some algebraic manipulation yields this equation for r_f:

r_f = 2 * grad(phi)_C * dCD / (phiD - phiC) - 1

This equation is clearly asymmetric considering the face between C and D because of the subscript on grad(phi). Hence this method can be thought of as constructing an r associated with the C side of the face

Definition at line 445 of file MathFVUtils.h.

Referenced by Moose::FV::MinModLimiter< T >::limit(), Moose::FV::VanLeerLimiter< T >::limit(), and Moose::FV::QUICKLimiter< T >::limit().

446 {
447  static const auto zero_vec = RealVectorValue(0);
448  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  return 1e6 * MathUtils::sign(gradC * dCD) + zero_vec * gradC;
455 
456  return 2. * gradC * dCD / (phiD - phiC) - 1.;
457 }
T sign(T x)
Definition: MathUtils.h:84

◆ selectInterpolationMethod()

InterpMethod Moose::FV::selectInterpolationMethod ( const std::string &  interp_method)

Definition at line 81 of file MathFVUtils.C.

Referenced by setInterpolationMethod(), FVRelationshipManagerInterface::setRMParamsAdvection(), and FVRelationshipManagerInterface::setRMParamsDiffusion().

82 {
83  if (interp_method == "average")
84  return InterpMethod::Average;
85  else if (interp_method == "harmonic")
86  return InterpMethod::HarmonicAverage;
87  else if (interp_method == "skewness-corrected")
88  return InterpMethod::SkewCorrectedAverage;
89  else if (interp_method == "upwind")
90  return InterpMethod::Upwind;
91  else if (interp_method == "rc")
92  return InterpMethod::RhieChow;
93  else if (interp_method == "vanLeer")
94  return InterpMethod::VanLeer;
95  else if (interp_method == "min_mod")
96  return InterpMethod::MinMod;
97  else if (interp_method == "sou")
98  return InterpMethod::SOU;
99  else if (interp_method == "quick")
100  return InterpMethod::QUICK;
101  else if (interp_method == "venkatakrishnan")
102  return InterpMethod::Venkatakrishnan;
103  else
104  mooseError("Interpolation method ",
105  interp_method,
106  " is not currently an option in Moose::FV::selectInterpolationMethod");
107 }
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:302

◆ setInterpolationMethod()

bool Moose::FV::setInterpolationMethod ( const MooseObject obj,
Moose::FV::InterpMethod interp_method,
const std::string &  param_name 
)

Sets one interpolation method.

Parameters
objThe MooseObject with input parameters to query
interp_methodThe interpolation method we will set
param_nameThe name of the parameter setting this interpolation method
Returns
Whether the interpolation method has indicated that we will need more than the default level of ghosting

Definition at line 110 of file MathFVUtils.C.

Referenced by FVAdvection::FVAdvection(), and LinearFVAdvection::LinearFVAdvection().

113 {
114  bool need_more_ghosting = false;
115 
116  const auto & interp_method_in = obj.getParam<MooseEnum>(param_name);
117  interp_method = selectInterpolationMethod(interp_method_in);
118 
119  if (interp_method == InterpMethod::SOU || interp_method == InterpMethod::MinMod ||
120  interp_method == InterpMethod::VanLeer || interp_method == InterpMethod::QUICK ||
121  interp_method == InterpMethod::Venkatakrishnan)
122  need_more_ghosting = true;
123 
124  return need_more_ghosting;
125 }
const T & getParam(const std::string &name) const
Retrieve a parameter for the object.
This is a "smart" enum class intended to replace many of the shortcomings in the C++ enum type It sho...
Definition: MooseEnum.h:33
InterpMethod selectInterpolationMethod(const std::string &interp_method)
Definition: MathFVUtils.C:81

◆ skewCorrectedLinearInterpolation()

template<typename T , typename T2 , typename T3 >
libMesh::CompareTypes<T, T2>::supertype Moose::FV::skewCorrectedLinearInterpolation ( const T &  value1,
const T2 &  value2,
const T3 &  face_gradient,
const FaceInfo fi,
const bool  one_is_elem 
)

Linear interpolation with skewness correction using the face gradient.

See more info in Moukalled Chapter 9. The correction involves a first order Taylor expansion around the intersection of the cell face and the line connecting the two cell centers.

Definition at line 261 of file MathFVUtils.h.

Referenced by linearInterpolation().

266 {
267  const auto coeffs = interpCoeffs(InterpMethod::SkewCorrectedAverage, fi, one_is_elem);
268 
269  auto value = (coeffs.first * value1 + coeffs.second * value2) +
270  face_gradient * fi.skewnessCorrectionVector();
271  return value;
272 }
Point skewnessCorrectionVector() const
Returns the skewness-correction vector (vector between the approximate and real face centroids)...
Definition: FaceInfo.C:80
Real value(unsigned n, unsigned alpha, unsigned beta, Real x)
std::pair< T, T > interpCoeffs(const Limiter< T > &limiter, const T &phi_upwind, const T &phi_downwind, const libMesh::VectorValue< T > *const grad_phi_upwind, const libMesh::VectorValue< T > *const grad_phi_face, const Real &max_value, const Real &min_value, const FaceInfo &fi, const bool fi_elem_is_upwind)
Produce the interpolation coefficients in the equation:
Definition: MathFVUtils.h:472

Variable Documentation

◆ moose_limiter_type

const MooseEnum Moose::FV::moose_limiter_type