Moose::FV Namespace Reference

Classes
class	CentralDifferenceLimiter
	Implements a limiter which reproduces a central-differencing scheme, defined by $\beta(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 $\beta(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 $\beta(r_f) = 0$. More...
class	VanLeerLimiter
	The Van Leer limiter limiter function $\beta(r_f)$ is defined as: More...
class	VenkatakrishnanLimiter
	The Venkatakrishnan limiter is derived from the following equations: More...

Enumerations
enum class	GradientLimiterType : int { None = -1 , Venkatakrishnan = 0 }
	Cell-gradient limiter variants used for MUSCL-style reconstructions. More...
enum class	LimiterType : int {   VanLeer = 0 , Upwind , CentralDifference , MinMod ,   SOU , QUICK , Venkatakrishnan }
enum class	InterpMethod {   Average , HarmonicAverage , SkewCorrectedAverage , Upwind ,   RhieChow , VanLeer , MinMod , SOU ,   QUICK , Venkatakrishnan }
	This codifies a set of available ways to interpolate with elem+neighbor solution information to calculate values (e.g. More...
enum class	LinearFVComputationMode { RHS , Matrix , FullSystem }

Functions
LimiterType	limiterType (InterpMethod interp_method)
	Return the limiter type associated with the supplied interpolation method.
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.
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.
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.
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.
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.
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.
template<typename T = Real>
std::pair< Real, Real >	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:
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.
template<typename T1 , typename T2 >
libMesh::CompareTypes< T1, T2 >::supertype	harmonicInterpolation (const T1 &value1, const T2 &value2, const Real weight1, const Real weight2)
	Computes a weighted harmonic mean of Reals, RealVectorValues, and RealTensorValues while accounting for the possibility that one or both values are AD.
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)) using the geometric weighting of a face.
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.
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).
template<typename T , FunctorEvaluationKind FEK = FunctorEvaluationKind::Value>
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 provided functor face argument
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.
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.
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.
template<typename Scalar , typename Vector >
Scalar	rF (const Scalar &phiC, const Scalar &phiD, const Vector &gradC, const RealVectorValue &dCD)
	From Moukalled 12.30.
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:
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.
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.
template<typename T >
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.
template<typename T , FunctorEvaluationKind FEK = FunctorEvaluationKind::Value, typename Enable = typename std::enable_if<ScalarTraits<T>::value>::type>
std::pair< Real, Real >	computeMinMaxValue (const FunctorBase< T > &functor, const FaceArg &face, const StateArg &time)
	This function calculates the minimum and maximum values within a two-cell stencil.
template<typename T >
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 stencil.
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.
template<typename T , FunctorEvaluationKind FEK = FunctorEvaluationKind::Value, typename Enable = typename std::enable_if<libMesh::ScalarTraits<T>::value>::type>
T	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.
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.
template<typename T >
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.
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.
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.

Variables
const MooseEnum	moose_limiter_type

Enumeration Type Documentation

GradientLimiterType

enum class Moose::FV::GradientLimiterType : int

strong

Cell-gradient limiter variants used for MUSCL-style reconstructions.

This is intentionally separate from face/value limiter selections. Gradient limiting typically has fewer supported options and uses different data/loops.

Enumerator
None	No gradient limiting.
Venkatakrishnan	Venkatakrishnan limiter (smooth, multidimensional).

Definition at line 22 of file GradientLimiterType.h.

{
  None = -1,
  Venkatakrishnan = 0
};

InterpMethod

enum class Moose::FV::InterpMethod

strong

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 38 of file MathFVUtils.h.

{
  Average,
  HarmonicAverage,
  SkewCorrectedAverage,
  Upwind,
  // Rhie-Chow
  RhieChow,
  VanLeer,
  MinMod,
  SOU,
  QUICK,
  Venkatakrishnan
};

LimiterType

enum class Moose::FV::LimiterType : int

strong

Enumerator
VanLeer
Upwind
CentralDifference
MinMod
SOU
QUICK
Venkatakrishnan

Definition at line 26 of file Limiter.h.

{
  VanLeer = 0,
  Upwind,
  CentralDifference,
  MinMod,
  SOU,
  QUICK,
  Venkatakrishnan
};

LinearFVComputationMode

enum class Moose::FV::LinearFVComputationMode

strong

Enumerator
RHS
Matrix
FullSystem

Definition at line 57 of file MathFVUtils.h.

{
  RHS,
  Matrix,
  FullSystem
};

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.

{
  auto params = emptyInputParameters();
  params.addParam<MooseEnum>("advected_interp_method",
                             interpolationMethods(),
                             "The interpolation to use for the advected quantity. Options are "
                             "'upwind', 'average', 'sou' (for second-order upwind), 'min_mod', "
                             "'vanLeer', 'quick', 'venkatakrishnan', and "
                             "'skewness-corrected' with the default being 'upwind'.");
  return params;
}

Referenced by FVAdvection::validParams().

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

T	The data type for the values being computed. This is typically a scalar type.
FEK	Enumeration type FunctorEvaluationKind with a default value of FunctorEvaluationKind::Value. This determines the kind of evaluation the functor will perform.
Enable	A 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

functor	An 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.
face	An argument representing the face in the computational domain. The face provides access to neighboring elements via neighbor_ptr_range().
time	An 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 663 of file MathFVUtils.h.

{
  // Initialize max_value to 0 and min_value to a large value
  Real max_value(std::numeric_limits<Real>::min()), min_value(std::numeric_limits<Real>::max());
 
  // Iterate over the direct neighbors of the element associated with the face
  for (const auto neighbor : (*face.fi).elem().neighbor_ptr_range())
  {
    // If not a valid neighbor, skip to the next one
    if (neighbor == nullptr)
      continue;
 
    // Evaluate the functor at the neighbor and update max_value and min_value
    const ElemArg local_cell_arg = {neighbor, false};
    const auto value =
        MetaPhysicL::raw_value(functor.template genericEvaluate<FEK>(local_cell_arg, time));
    max_value = std::max(max_value, value);
    min_value = std::min(min_value, value);
  }
 
  // Iterate over the neighbors of the neighbor
  for (const auto neighbor : (*face.fi).neighbor().neighbor_ptr_range())
  {
    // If not a valid neighbor, skip to the next one
    if (neighbor == nullptr)
      continue;
 
    // Evaluate the functor at the neighbor and update max_value and min_value
    const ElemArg local_cell_arg = {neighbor, false};
    const auto value =
        MetaPhysicL::raw_value(functor.template genericEvaluate<FEK>(local_cell_arg, time));
    max_value = std::max(max_value, value);
    min_value = std::min(min_value, value);
  }
 
  // Return a pair containing the computed maximum and minimum values
  return std::make_pair(max_value, min_value);
}

Referenced by interpCoeffsAndAdvected(), interpolate(), and interpolate().

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

T	The data type for the values being computed. This is typically a scalar type.

Parameters

functor	An 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.
face	An argument representing the face in the computational domain. The face provides access to neighboring elements via neighbor_ptr_range().
time	An argument representing the state or time at which the evaluation is performed.
component	An 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 729 of file MathFVUtils.h.

{
  // Initialize max_value to 0 and min_value to a large value
  Real max_value(std::numeric_limits<Real>::min()), min_value(std::numeric_limits<Real>::max());
 
  // Iterate over the direct neighbors of the element associated with the face
  for (const auto neighbor : (*face.fi).elem().neighbor_ptr_range())
  {
    // If not a valid neighbor, skip to the next one
    if (neighbor == nullptr)
      continue;
 
    // Evaluate the functor at the neighbor for the specified component and update max_value and
    // min_value
    const ElemArg local_cell_arg = {neighbor, false};
    const auto value = MetaPhysicL::raw_value(functor(local_cell_arg, time)(component));
    max_value = std::max(max_value, value);
    min_value = std::min(min_value, value);
  }
 
  // Iterate over the neighbors of the neighbor associated with the face
  for (const auto neighbor : (*face.fi).neighbor().neighbor_ptr_range())
  {
    // If not a valid neighbor, skip to the next one
    if (neighbor == nullptr)
      continue;
 
    // Evaluate the functor at the neighbor for the specified component and update max_value and
    // min_value
    const ElemArg local_cell_arg = {neighbor, false};
    const auto value = MetaPhysicL::raw_value(functor(local_cell_arg, time)(component));
    max_value = std::max(max_value, value);
    min_value = std::min(min_value, value);
  }
 
  // Return a pair containing the computed maximum and minimum values
  return std::make_pair(max_value, min_value);
}

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

T	The data type for the container values being interpolated. This is typically a container type such as a vector or array.

Parameters

functor	An 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.
face	An argument representing the face in the computational domain. The face provides access to neighboring elements and limiter type.
time	An 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 1105 of file MathFVUtils.h.

{
  typedef typename FunctorBase<T>::GradientType ContainerGradientType;
  typedef typename ContainerGradientType::value_type GradientType;
  const GradientType * const example_gradient = nullptr;
 
  mooseAssert(face.fi, "this must be non-null");
  const auto limiter = Limiter<typename T::value_type>::build(face.limiter_type);
 
  const auto upwind_arg = face.elem_is_upwind ? face.makeElem() : face.makeNeighbor();
  const auto downwind_arg = face.elem_is_upwind ? face.makeNeighbor() : face.makeElem();
  const auto phi_upwind = functor(upwind_arg, time);
  const auto phi_downwind = functor(downwind_arg, time);
 
  // initialize in order to get proper size
  T ret = phi_upwind;
  typename T::value_type coeff_upwind, coeff_downwind;
 
  if (face.limiter_type == LimiterType::Upwind ||
      face.limiter_type == LimiterType::CentralDifference)
  {
    for (const auto i : index_range(ret))
    {
      const auto &component_upwind = phi_upwind[i], component_downwind = phi_downwind[i];
      std::tie(coeff_upwind, coeff_downwind) = interpCoeffs(*limiter,
                                                            component_upwind,
                                                            component_downwind,
                                                            example_gradient,
                                                            example_gradient,
                                                            std::numeric_limits<Real>::max(),
                                                            std::numeric_limits<Real>::min(),
                                                            *face.fi,
                                                            face.elem_is_upwind);
      ret[i] = coeff_upwind * component_upwind + coeff_downwind * component_downwind;
    }
  }
  else
  {
    const auto grad_phi_upwind = functor.gradient(upwind_arg, time);
    for (const auto i : index_range(ret))
    {
      const auto &component_upwind = phi_upwind[i], component_downwind = phi_downwind[i];
      const auto & grad = grad_phi_upwind[i];
      std::tie(coeff_upwind, coeff_downwind) = interpCoeffs(*limiter,
                                                            component_upwind,
                                                            component_downwind,
                                                            &grad,
                                                            example_gradient,
                                                            std::numeric_limits<Real>::max(),
                                                            std::numeric_limits<Real>::min(),
                                                            *face.fi,
                                                            face.elem_is_upwind);
      ret[i] = coeff_upwind * component_upwind + coeff_downwind * component_downwind;
    }
  }
 
  return ret;
}

Referenced by interpolate(), and interpolate().

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.

{
  auto ft = fi.faceType(std::make_pair(var.number(), var.sys().number()));
  mooseAssert(ft == FaceInfo::VarFaceNeighbors::BOTH
                  ? var.hasBlocks(fi.elem().subdomain_id()) && fi.neighborPtr() &&
                        var.hasBlocks(fi.neighborPtr()->subdomain_id())
                  : true,
              "Finite volume variable " << var.name()
                                        << " does not exist on both sides of the face despite "
                                           "what the FaceInfo is telling us.");
  mooseAssert(ft == FaceInfo::VarFaceNeighbors::ELEM
                  ? var.hasBlocks(fi.elem().subdomain_id()) &&
                        (!fi.neighborPtr() || !var.hasBlocks(fi.neighborPtr()->subdomain_id()))
                  : true,
              "Finite volume variable " << var.name()
                                        << " does not exist on or only on the elem side of the "
                                           "face despite what the FaceInfo is telling us.");
  mooseAssert(ft == FaceInfo::VarFaceNeighbors::NEIGHBOR
                  ? fi.neighborPtr() && var.hasBlocks(fi.neighborPtr()->subdomain_id()) &&
                        !var.hasBlocks(fi.elem().subdomain_id())
                  : true,
              "Finite volume variable " << var.name()
                                        << " does not exist on or only on the neighbor side of the "
                                           "face despite what the FaceInfo is telling us.");
 
  const bool one_is_elem =
      ft == FaceInfo::VarFaceNeighbors::BOTH || ft == FaceInfo::VarFaceNeighbors::ELEM;
  const Elem * const elem_one = one_is_elem ? &fi.elem() : fi.neighborPtr();
  mooseAssert(elem_one, "This elem should be non-null!");
  const Elem * const elem_two = one_is_elem ? fi.neighborPtr() : &fi.elem();
 
  return std::make_tuple(elem_one, elem_two, one_is_elem);
}

Referenced by MooseVariableFV< RealEigenVector >::getDirichletBoundaryFaceValue().

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

elem	Reference to an element
neighbor	Pointer 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.

{
  // The face info belongs to elem:
  //  * at all mesh boundaries (i.e. where there is no neighbor)
  //  * if the element faces a neighbor which is on a lower refinement level
  //  * if the element is active and it has a lower ID than its neighbor
  if (!neighbor)
    return true;
  else if (elem.level() != neighbor->level())
    return neighbor->level() < elem.level();
  else if (!neighbor->active())
    return false;
  else
    return elem.id() < neighbor->id();
}

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

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

T	The data type of the scalar values and the return type.

Parameters

limiter	The limiter object used to compute the flux limiting ratio.
phi_upwind	The field value at the upwind location.
phi_downwind	The field value at the downwind location.
grad_phi_upwind	Pointer to the gradient at the upwind location.
grad_phi_face	Pointer to the gradient at the face location.
fi	FaceInfo object containing geometric details such as face centroid and cell centroids.
fi_elem_is_upwind	Boolean indicating if the face info element is upwind.
max_value	The maximum allowable value.
min_value	The minimum allowable value.

Returns

The computed limited interpolation value.

Definition at line 598 of file MathFVUtils.h.

{
  // Compute the direction vector based on whether the current element is upwind
  const auto dCD = face.elem_is_upwind ? face.fi->dCN() : Point(-face.fi->dCN());
 
  // Compute the flux limiting ratio using the specified limiter
  const auto beta = limiter(phi_upwind,
                            phi_downwind,
                            grad_phi_upwind,
                            grad_phi_face,
                            dCD,
                            max_value,
                            min_value,
                            face.fi,
                            face.elem_is_upwind);
 
  // Determine the face centroid and the appropriate cell centroid
  const auto & face_centroid = face.fi->faceCentroid();
  const auto & cell_centroid =
      face.elem_is_upwind ? face.fi->elemCentroid() : face.fi->neighborCentroid();
 
  // Compute the delta value at the face
  T delta_face;
  if (face.limiter_type == LimiterType::Venkatakrishnan || face.limiter_type == LimiterType::SOU)
    delta_face = (*grad_phi_upwind) * (face_centroid - cell_centroid);
  else
    delta_face = (*grad_phi_face) * (face_centroid - cell_centroid);
 
  // Return the interpolated value
  return phi_upwind + beta * delta_face;
}

Referenced by interpolate(), and interpolate().

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.

{
  return fv_var.adGradSln(face_info, time, correct_skewness) * face_info.normal();
}

Referenced by FVFluxKernel::gradUDotNormal().

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_arg	An element argument specifying the current element and whether to perform skew corrections
state_arg	A state argument that indicates what temporal / solution iteration data to use when evaluating the provided functor
functor	The functor that will provide information such as cell and face value evaluations necessary to construct the cell gradient
two_term_boundary_expansion	Whether 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
mesh	The mesh on which we are computing the gradient

Returns

The computed cell gradient

Definition at line 44 of file GreenGaussGradient.h.

{
  mooseAssert(elem_arg.elem, "This should be non-null");
  const auto coord_type = mesh.getCoordSystem(elem_arg.elem->subdomain_id());
  const auto rz_radial_coord = mesh.getAxisymmetricRadialCoord();
 
  const T elem_value = functor(elem_arg, state_arg);
 
  // We'll count the number of extrapolated boundary faces (ebfs)
  unsigned int num_ebfs = 0;
 
  // Gradient to be returned
  VectorValue<T> grad;
 
  if (!elem_arg.correct_skewness || force_green_gauss) // Do Green-Gauss
  {
    try
    {
      bool volume_set = false;
      Real volume = 0;
 
      // If we are performing a two term Taylor expansion for extrapolated boundary faces (faces on
      // boundaries that do not have associated Dirichlet conditions), then the element gradient
      // depends on the boundary face value and the boundary face value depends on the element
      // gradient, so we have a system of equations to solve. Here is the system:
      //
      // \nabla \phi_C - \frac{1}{V} \sum_{ebf} \phi_{ebf} \vec{S_f} =
      //   \frac{1}{V} \sum_{of} \phi_{of} \vec{S_f}                       eqn. 1
      //
      // \phi_{ebf} - \vec{d_{Cf}} \cdot \nabla \phi_C = \phi_C            eqn. 2
      //
      // where $C$ refers to the cell centroid, $ebf$ refers to an extrapolated boundary face, $of$
      // refers to "other faces", e.g. non-ebf faces, and $f$ is a general face. $d_{Cf}$ is the
      // vector drawn from the element centroid to the face centroid, and $\vec{S_f}$ is the surface
      // vector, e.g. the face area times the outward facing normal
 
      // ebf eqns: element gradient coefficients, e.g. eqn. 2, LHS term 2 coefficient
      std::vector<libMesh::VectorValue<Real>> ebf_grad_coeffs;
      // ebf eqns: rhs b values. These will actually correspond to the elem_value so we can use a
      // pointer and avoid copying. This is the RHS of eqn. 2
      std::vector<const T *> ebf_b;
 
      // elem grad eqns: ebf coefficients, e.g. eqn. 1, LHS term 2 coefficients
      std::vector<libMesh::VectorValue<Real>> grad_ebf_coeffs;
      // elem grad eqns: rhs b value, e.g. eqn. 1 RHS
      libMesh::VectorValue<T> grad_b = 0;
 
      auto action_functor = [&volume_set,
                             &volume,
                             &elem_value,
                             &elem_arg,
                             &num_ebfs,
                             &ebf_grad_coeffs,
                             &ebf_b,
                             &grad_ebf_coeffs,
                             &grad_b,
                             &state_arg,
                             &functor,
                             two_term_boundary_expansion,
                             coord_type,
                             rz_radial_coord](const Elem & libmesh_dbg_var(functor_elem),
                                              const Elem *,
                                              const FaceInfo * const fi,
                                              const Point & surface_vector,
                                              Real coord,
                                              const bool elem_has_info)
      {
        mooseAssert(fi, "We need a FaceInfo for this action_functor");
        mooseAssert(
            elem_arg.elem == &functor_elem,
            "Just a sanity check that the element being passed in is the one we passed out.");
 
        if (functor.isExtrapolatedBoundaryFace(*fi, elem_arg.elem, state_arg))
        {
          if (two_term_boundary_expansion)
          {
            num_ebfs += 1;
 
            // eqn. 2
            ebf_grad_coeffs.push_back(-1. * (elem_has_info
                                                 ? (fi->faceCentroid() - fi->elemCentroid())
                                                 : (fi->faceCentroid() - fi->neighborCentroid())));
            ebf_b.push_back(&elem_value);
 
            // eqn. 1
            grad_ebf_coeffs.push_back(-surface_vector);
          }
          else
            // We are doing a one-term expansion for the extrapolated boundary faces, in which case
            // we have no eqn. 2 and we have no second term in the LHS of eqn. 1. Instead we apply
            // the element centroid value as the face value (one-term expansion) in the RHS of eqn.
            // 1
            grad_b += surface_vector * elem_value;
        }
        else
          grad_b +=
              surface_vector * functor(Moose::FaceArg{fi,
                                                      Moose::FV::LimiterType::CentralDifference,
                                                      true,
                                                      elem_arg.correct_skewness,
                                                      elem_arg.elem,
                                                      nullptr},
                                       state_arg);
 
        if (!volume_set)
        {
          // We use the FaceInfo volumes because those values have been pre-computed and cached.
          // An explicit call to elem->volume() here would incur unnecessary expense
          if (elem_has_info)
          {
            MooseMeshUtils::coordTransformFactor(
                fi->elemCentroid(), coord, coord_type, rz_radial_coord);
            volume = fi->elemVolume() * coord;
          }
          else
          {
            MooseMeshUtils::coordTransformFactor(
                fi->neighborCentroid(), coord, coord_type, rz_radial_coord);
            volume = fi->neighborVolume() * coord;
          }
 
          volume_set = true;
        }
      };
 
      Moose::FV::loopOverElemFaceInfo(
          *elem_arg.elem, mesh, action_functor, coord_type, rz_radial_coord);
 
      mooseAssert(volume_set && volume > 0, "We should have set the volume");
      grad_b /= volume;
 
      if (coord_type == Moose::CoordinateSystemType::COORD_RZ)
      {
        mooseAssert(rz_radial_coord != libMesh::invalid_uint, "rz_radial_coord must be set");
        grad_b(rz_radial_coord) -= elem_value / elem_arg.elem->vertex_average()(rz_radial_coord);
      }
 
      mooseAssert(
          coord_type != Moose::CoordinateSystemType::COORD_RSPHERICAL,
          "We have not yet implemented the correct translation from gradient to divergence for "
          "spherical coordinates yet.");
 
      // test for simple case
      if (num_ebfs == 0)
        grad = grad_b;
      else
      {
        // We have to solve a system
        const unsigned int sys_dim = Moose::dim + num_ebfs;
        DenseVector<T> x(sys_dim), b(sys_dim);
        DenseMatrix<T> A(sys_dim, sys_dim);
 
        // Let's make i refer to Moose::dim indices, and j refer to num_ebfs indices
 
        // eqn. 1
        for (const auto i : make_range(Moose::dim))
        {
          // LHS term 1 coeffs
          A(i, i) = 1;
 
          // LHS term 2 coeffs
          for (const auto j : make_range(num_ebfs))
            A(i, Moose::dim + j) = grad_ebf_coeffs[j](i) / volume;
 
          // RHS
          b(i) = grad_b(i);
        }
 
        // eqn. 2
        for (const auto j : make_range(num_ebfs))
        {
          // LHS term 1 coeffs
          A(Moose::dim + j, Moose::dim + j) = 1;
 
          // LHS term 2 coeffs
          for (const auto i : make_range(unsigned(Moose::dim)))
            A(Moose::dim + j, i) = ebf_grad_coeffs[j](i);
 
          // RHS
          b(Moose::dim + j) = *ebf_b[j];
        }
 
        A.lu_solve(b, x);
        // libMesh is generous about what it considers nonsingular. Let's check a little more
        // strictly
        if (MooseUtils::absoluteFuzzyEqual(MetaPhysicL::raw_value(A(sys_dim - 1, sys_dim - 1)), 0))
          throw libMesh::LogicError("Matrix A is singular!");
        for (const auto i : make_range(Moose::dim))
          grad(i) = x(i);
      }
    }
    catch (std::exception & e)
    {
      // Don't ignore any *real* errors; we only handle matrix
      // singularity (LogicError in older libMesh, DegenerateMap in
      // newer) here
      if (!strstr(e.what(), "singular"))
        throw;
 
      // Retry without two-term
      if (!two_term_boundary_expansion)
        mooseError(
            "I believe we should only get singular systems when two-term boundary expansion is "
            "being used. The error thrown during the computation of the gradient: ",
            e.what());
 
      grad = greenGaussGradient(elem_arg, state_arg, functor, false, mesh, true);
    }
  }
  else // Do Least-Squares
  {
    try
    {
 
      // The least squares method aims to find the gradient at the element centroid by minimizing
      // the difference between the estimated and actual value differences across neighboring cells.
      // The least squares formulation is:
      //
      // Minimize J = \sum_{n} [ w_n ((\nabla \phi_C \cdot \delta \vec{x}_n) - \delta \phi_n) ]^2
      //
      // where:
      // - \‍(\nabla \phi_C\‍) is the gradient at the element centroid C (unknown).
      // - \‍(\delta \vec{x}_n = \vec{x}_n - \vec{x}_C\‍) is the vector from the element centroid to
      // neighbor \‍(n\‍).
      // - \‍(\delta \phi_n = \phi_n - \phi_C\‍) is the difference in the scalar value between
      // neighbor \‍(n\‍) and element \‍(C\‍).
      // - w_n = 1/||\vec{x}_n|| is a vector of weights that is used to handle large aspect ratio
      // cells
      //
      // To handle extrapolated boundary faces (faces on boundaries without Dirichlet conditions),
      // we include additional unknowns and equations in the least squares system.
      // For ebfs, we may not know the value \‍(\phi_n\‍), so we include it as an unknown.
      // This results in an augmented system that simultaneously solves for the gradient and the ebf
      // values.
 
      // Get estimated number of faces in a cell
      const auto ptr_range = elem_arg.elem->neighbor_ptr_range();
      const size_t num_faces = std::distance(ptr_range.begin(), ptr_range.end());
 
      // Lists to store position differences and value differences for least squares computation
      std::vector<Point> delta_x_list; // List of position differences between neighbor centroids
                                       // and element centroid
      delta_x_list.reserve(num_faces);
 
      std::vector<T>
          delta_phi_list; // List of value differences between neighbor values and element value
      delta_phi_list.reserve(num_faces);
 
      // Variables to handle extrapolated boundary faces (ebfs) in the least squares method
      std::vector<Point> delta_x_ebf_list; // Position differences for ebfs
      delta_phi_list.reserve(num_faces);
 
      std::vector<VectorValue<Real>>
          ebf_grad_coeffs; // Coefficients for the gradient in ebf equations
      delta_phi_list.reserve(num_faces);
 
      std::vector<const T *> ebf_b; // RHS values for ebf equations (pointers to avoid copying)
      delta_phi_list.reserve(num_faces);
 
      unsigned int num_ebfs = 0; // Number of extrapolated boundary faces
 
      // Action functor to collect data from each face of the element
      auto action_functor = [&elem_value,
                             &elem_arg,
                             &num_ebfs,
                             &ebf_grad_coeffs,
                             &ebf_b,
                             &delta_x_list,
                             &delta_phi_list,
                             &delta_x_ebf_list,
                             &state_arg,
                             &functor,
                             two_term_boundary_expansion](const Elem & libmesh_dbg_var(loc_elem),
                                                          const Elem * loc_neigh,
                                                          const FaceInfo * const fi,
                                                          const Point & /*surface_vector*/,
                                                          Real /*coord*/,
                                                          const bool elem_has_info)
      {
        mooseAssert(fi, "We need a FaceInfo for this action_functor");
        mooseAssert(
            elem_arg.elem == &loc_elem,
            "Just a sanity check that the element being passed in is the one we passed out.");
 
        if (functor.isExtrapolatedBoundaryFace(*fi, elem_arg.elem, state_arg))
        {
          // Extrapolated Boundary Face (ebf)
          const Point delta_x = elem_has_info ? (fi->faceCentroid() - fi->elemCentroid())
                                              : (fi->faceCentroid() - fi->neighborCentroid());
          delta_x_list.push_back(delta_x);
 
          T delta_phi;
 
          if (two_term_boundary_expansion)
          {
            // Two-term expansion: include ebf values as unknowns in the augmented system
            num_ebfs += 1;
 
            // Coefficient for the gradient in the ebf equation
            ebf_grad_coeffs.push_back(-delta_x);
            // RHS value for the ebf equation
            ebf_b.push_back(&elem_value);
 
            delta_phi = -elem_value;
            delta_x_ebf_list.push_back(delta_x);
          }
          else
          {
            // One-term expansion: assume zero difference across the boundary (\delta \phi = 0)
            delta_phi = T(0);
          }
          delta_phi_list.push_back(delta_phi);
        }
        else
        {
          // Internal Face or Boundary Face with Dirichlet condition
          Point delta_x;
          T neighbor_value;
          if (functor.isInternalFace(*fi))
          {
            // Internal face with a neighboring element
            delta_x = loc_neigh->vertex_average() - elem_arg.elem->vertex_average();
 
            const ElemArg neighbor_arg = {loc_neigh, elem_arg.correct_skewness};
            neighbor_value = functor(neighbor_arg, state_arg);
          }
          else
          {
            // Boundary face with Dirichlet condition
            delta_x = elem_has_info ? (fi->faceCentroid() - fi->elemCentroid())
                                    : (fi->faceCentroid() - fi->neighborCentroid());
            neighbor_value = functor(Moose::FaceArg{fi,
                                                    Moose::FV::LimiterType::CentralDifference,
                                                    true,
                                                    elem_arg.correct_skewness,
                                                    elem_arg.elem,
                                                    nullptr},
                                     state_arg);
          }
 
          delta_x_list.push_back(delta_x);
          delta_phi_list.push_back(neighbor_value - elem_value);
        }
      };
 
      // Loop over element faces and apply the action_functor
      Moose::FV::loopOverElemFaceInfo(
          *elem_arg.elem, mesh, action_functor, coord_type, rz_radial_coord);
 
      // Compute Least Squares gradient
      const unsigned int n = delta_x_list.size();
      mooseAssert(n >= dim, "Not enough neighbors to perform least squares");
 
      DenseMatrix<T> ATA(dim, dim);
      DenseVector<T> ATb(dim);
 
      // Assemble the normal equations for the least squares method
      // ATA = \sum_n (\delta \vec{x}_n^T * \delta \vec{x}_n)
      // ATb = \sum_n (\delta \vec{x}_n^T * \delta \phi_n)
      ATA.zero();
      ATb.zero();
 
      for (unsigned int i = 0; i < n; ++i)
      {
        const Point & delta_x = delta_x_list[i];
        const T & delta_phi = delta_phi_list[i];
 
        for (unsigned int d1 = 0; d1 < dim; ++d1)
        {
          const Real dx_d1 = delta_x(d1);
          const Real dx_norm = delta_x.norm();
          const Real dx_d1_norm = dx_d1 / dx_norm;
 
          for (unsigned int d2 = 0; d2 < dim; ++d2)
          {
            const Real dx_d2_norm = delta_x(d2) / dx_norm;
            ATA(d1, d2) += dx_d1_norm * dx_d2_norm;
          }
 
          ATb(d1) += dx_d1_norm * delta_phi / dx_norm;
        }
      }
 
      // Add regularization to ATA to prevent singularity
      T epsilon = 1e-12; // Small regularization parameter
      for (unsigned int d = 0; d < dim; ++d)
      {
        ATA(d, d) += epsilon;
      }
 
      if (num_ebfs == 0)
      {
        // Solve and store the least squares gradient
        DenseVector<T> grad_ls(dim);
        ATA.lu_solve(ATb, grad_ls);
 
        // Store the least squares gradient
        for (unsigned int d = 0; d < dim; ++d)
        {
          grad(d) = grad_ls(d);
        }
      }
      else
      {
        // We have to solve a system of equations
        const unsigned int sys_dim = Moose::dim + num_ebfs;
        DenseVector<T> x(sys_dim), b(sys_dim);
        DenseMatrix<T> A(sys_dim, sys_dim);
 
        // Let's make i refer to Moose::dim indices, and j refer to num_ebfs indices
 
        // eqn. 1: Element gradient equations
        for (unsigned int d1 = 0; d1 < dim; ++d1)
        {
          // LHS term 1 coefficients (gradient components)
          for (unsigned int d2 = 0; d2 < dim; ++d2)
            A(d1, d2) = ATA(d1, d2);
 
          // RHS
          b(d1) = ATb(d1);
        }
 
        // LHS term 2 coefficients (ebf contributions)
        for (const auto i : make_range(num_ebfs))
        {
          const Point & delta_x = delta_x_ebf_list[i];
          for (unsigned int d1 = 0; d1 < dim; ++d1)
          {
            const Real dx_d1 = delta_x(d1);
            A(d1, Moose::dim + i) -= dx_d1 / delta_x.norm() / delta_x.norm();
          }
        }
 
        // eqn. 2: Extrapolated boundary face equations
        for (const auto j : make_range(num_ebfs))
        {
          // LHS term 1 coefficients (ebf values)
          A(Moose::dim + j, Moose::dim + j) = 1;
 
          // LHS term 2 coefficients (gradient components)
          for (const auto i : make_range(unsigned(Moose::dim)))
            A(Moose::dim + j, i) = ebf_grad_coeffs[j](i);
 
          // RHS
          b(Moose::dim + j) = *ebf_b[j];
        }
 
        // Solve the system
        A.lu_solve(b, x);
 
        // Check for singularity
        if (MooseUtils::absoluteFuzzyEqual(MetaPhysicL::raw_value(A(sys_dim - 1, sys_dim - 1)), 0))
          throw libMesh::LogicError("Matrix A is singular!");
 
        // Extract the gradient components
        for (const auto i : make_range(Moose::dim))
          grad(i) = x(i);
      }
 
      mooseAssert(
          coord_type != Moose::CoordinateSystemType::COORD_RSPHERICAL,
          "We have not yet implemented the correct translation from gradient to divergence for "
          "spherical coordinates yet.");
    }
    catch (std::exception & e)
    {
      // Don't ignore any *real* errors; we only handle matrix
      // singularity (LogicError in older libMesh, DegenerateMap in
      // newer) here
      if (!strstr(e.what(), "singular"))
        throw;
 
      // Log warning and default to simple green Gauss
      mooseWarning(
          "Singular matrix encountered in least squares gradient computation. Falling back "
          "to Green-Gauss gradient.");
 
      grad = greenGaussGradient(
          elem_arg, state_arg, functor, false, mesh, /* force_green_gauss */ true);
    }
  }
 
  return grad;
 
  // Notes to future developer:
  // Note 1:
  // For the least squares gradient, the LS matrix could be precomputed and stored for every cell
  // I tried doing this on October 2024, but the element lookup for these matrices is too slow
  // and seems better to compute weights on the fly.
  // Consider building a map from elem_id to these matrices and speed up lookup with octree
  // Note 2:
  // In the future one would like to have a hybrid gradient scheme, where:
  // \nabla \phi = \beta (\nabla \phi)_{LS} + (1 - \beta) (\nabla \phi)_{GG}
  // Then optize \beta based on mesh heuristics
}

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

greenGaussGradient() [2/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 621 of file GreenGaussGradient.h.

{
  TensorValue<T> ret;
  for (const auto i : make_range(Moose::dim))
  {
    VectorComponentFunctor<T> scalar_functor(functor, i);
    const auto row_gradient =
        greenGaussGradient(elem_arg, state_arg, scalar_functor, two_term_boundary_expansion, mesh);
    for (const auto j : make_range(unsigned(Moose::dim)))
      ret(i, j) = row_gradient(j);
  }
 
  return ret;
}

greenGaussGradient() [3/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 713 of file GreenGaussGradient.h.

{
  typedef typename Moose::FunctorBase<std::array<T, N>>::GradientType GradientType;
  GradientType ret;
  for (const auto i : make_range(N))
  {
    // Note that this can be very inefficient. Within the scalar greenGaussGradient routine we're
    // going to do value type evaluations of the array functor from scalar_functor and we will be
    // discarding all the value type evaluations other than the one corresponding to i
    ArrayComponentFunctor<T, FunctorBase<std::array<T, N>>> scalar_functor(functor, i);
    ret[i] =
        greenGaussGradient(elem_arg, state_arg, scalar_functor, two_term_boundary_expansion, mesh);
  }
 
  return ret;
}

greenGaussGradient() [4/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 663 of file GreenGaussGradient.h.

{
  // Determine the size of the container
  const auto vals = functor(elem_arg, state_arg);
  typedef typename Moose::FunctorBase<std::vector<T>>::GradientType GradientType;
  GradientType ret(vals.size());
  for (const auto i : index_range(ret))
  {
    // Note that this can be very inefficient. Within the scalar greenGaussGradient routine we're
    // going to do value type evaluations of the array functor from scalar_functor and we will be
    // discarding all the value type evaluations other than the one corresponding to i
    ArrayComponentFunctor<T, FunctorBase<std::vector<T>>> scalar_functor(functor, i);
    ret[i] =
        greenGaussGradient(elem_arg, state_arg, scalar_functor, two_term_boundary_expansion, mesh);
  }
 
  return ret;
}

greenGaussGradient() [5/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_arg	A face argument specifying the current faceand whether to perform skew corrections
state_arg	A state argument that indicates what temporal / solution iteration data to use when evaluating the provided functor
functor	The functor that will provide information such as cell and face value evaluations necessary to construct the cell gradient
two_term_boundary_expansion	Whether 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
mesh	The mesh on which we are computing the gradient

Returns

The computed cell gradient

Definition at line 563 of file GreenGaussGradient.h.

{
  mooseAssert(face_arg.fi, "We should have a face info to compute a face gradient");
  const auto & fi = *(face_arg.fi);
  const auto & elem_arg = face_arg.makeElem();
  const auto & neighbor_arg = face_arg.makeNeighbor();
  const bool defined_on_elem = functor.hasBlocks(fi.elemSubdomainID());
  const bool defined_on_neighbor = fi.neighborPtr() && functor.hasBlocks(fi.neighborSubdomainID());
 
  if (defined_on_elem && defined_on_neighbor)
  {
    const auto & value_elem = functor(elem_arg, state_arg);
    const auto & value_neighbor = functor(neighbor_arg, state_arg);
 
    // This is the component of the gradient which is parallel to the line connecting
    // the cell centers. Therefore, we can use our second order, central difference
    // scheme to approximate it.
    libMesh::VectorValue<T> face_gradient = (value_neighbor - value_elem) / fi.dCNMag() * fi.eCN();
 
    // We only need nonorthogonal correctors in 2+ dimensions
    if (mesh.dimension() > 1)
    {
      // We are using an orthogonal approach for the non-orthogonal correction, for more information
      // see Hrvoje Jasak's PhD Thesis (Imperial College, 1996)
      libMesh::VectorValue<T> interpolated_gradient;
 
      // Compute the gradients in the two cells on both sides of the face
      const auto & grad_elem =
          greenGaussGradient(elem_arg, state_arg, functor, two_term_boundary_expansion, mesh);
      const auto & grad_neighbor =
          greenGaussGradient(neighbor_arg, state_arg, functor, two_term_boundary_expansion, mesh);
 
      Moose::FV::interpolate(Moose::FV::InterpMethod::Average,
                             interpolated_gradient,
                             grad_elem,
                             grad_neighbor,
                             fi,
                             true);
 
      face_gradient += interpolated_gradient - (interpolated_gradient * fi.eCN()) * fi.eCN();
    }
 
    return face_gradient;
  }
  else if (defined_on_elem)
    return functor.gradient(elem_arg, state_arg);
  else if (defined_on_neighbor)
    return functor.gradient(neighbor_arg, state_arg);
  else
    mooseError("The functor must be defined on one of the sides");
}

greenGaussGradient() [6/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 642 of file GreenGaussGradient.h.

{
  TensorValue<T> ret;
  for (const auto i : make_range(unsigned(Moose::dim)))
  {
    VectorComponentFunctor<T> scalar_functor(functor, i);
    const auto row_gradient =
        greenGaussGradient(face_arg, state_arg, scalar_functor, two_term_boundary_expansion, mesh);
    for (const auto j : make_range(unsigned(Moose::dim)))
      ret(i, j) = row_gradient(j);
  }
 
  return ret;
}

greenGaussGradient() [7/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 736 of file GreenGaussGradient.h.

{
  typedef typename Moose::FunctorBase<std::array<T, N>>::GradientType GradientType;
  GradientType ret;
  for (const auto i : make_range(N))
  {
    // Note that this can be very inefficient. Within the scalar greenGaussGradient routine we're
    // going to do value type evaluations of the array functor from scalar_functor and we will be
    // discarding all the value type evaluations other than the one corresponding to i
    ArrayComponentFunctor<T, FunctorBase<std::array<T, N>>> scalar_functor(functor, i);
    ret[i] =
        greenGaussGradient(face_arg, state_arg, scalar_functor, two_term_boundary_expansion, mesh);
  }
 
  return ret;
}

greenGaussGradient() [8/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 688 of file GreenGaussGradient.h.

{
  // Determine the size of the container
  const auto vals = functor(face_arg, state_arg);
  typedef typename Moose::FunctorBase<std::vector<T>>::GradientType GradientType;
  GradientType ret(vals.size());
  for (const auto i : index_range(ret))
  {
    // Note that this can be very inefficient. Within the scalar greenGaussGradient routine we're
    // going to do value type evaluations of the array functor from scalar_functor and we will be
    // discarding all the value type evaluations other than the one corresponding to i
    ArrayComponentFunctor<T, FunctorBase<std::vector<T>>> scalar_functor(functor, i);
    ret[i] =
        greenGaussGradient(face_arg, state_arg, scalar_functor, two_term_boundary_expansion, mesh);
  }
 
  return ret;
}

harmonicInterpolation() [1/2]

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)) using the geometric weighting of a face.

Parameters

value1	Reference to value1 in the (1/(gc/value1+(1-gc)/value2)) expression
value2	Reference to value2 in the (1/(gc/value1+(1-gc)/value2)) expression
fi	Reference to the FaceInfo of the face on which the interpolation is requested
one_is_elem	Boolean indicating if the interpolation weight on FaceInfo belongs to the element corresponding to value1

Returns

The interpolated value

Definition at line 260 of file MathFVUtils.h.

{
  const auto coeffs = interpCoeffs(InterpMethod::Average, fi, one_is_elem);
  return harmonicInterpolation(value1, value2, coeffs.first, coeffs.second);
}

harmonicInterpolation() [2/2]

template<typename T1 , typename T2 >

libMesh::CompareTypes< T1, T2 >::supertype Moose::FV::harmonicInterpolation	(	const T1 &	value1,
		const T2 &	value2,
		const Real	weight1,
		const Real	weight2
	)

Computes a weighted harmonic mean of Reals, RealVectorValues, and RealTensorValues while accounting for the possibility that one or both values are AD.

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

Parameters

value1	Reference to value1 in the 1/(weight1/value1 + weight2/value2) expression
value2	Reference to value2 in the 1/(weight1/value1 + weight2/value2) expression
weight1	Weight corresponding to value1
weight2	Weight corresponding to value2

Returns

The interpolated value

Definition at line 179 of file MathFVUtils.h.

{
  // We check if the base values of the given template types match, if not we throw a compile-time
  // error
  static_assert(std::is_same<typename MetaPhysicL::RawType<T1>::value_type,
                             typename MetaPhysicL::RawType<T2>::value_type>::value,
                "The input values for harmonic interpolation need to have the same raw-value!");
 
  // We check if the types are fit to compute the harmonic mean of. This is done compile-time
  // using constexpr. We start with Real/ADReal which is straightforward if the input values are
  // positive.
  if constexpr (libMesh::TensorTools::TensorTraits<T1>::rank == 0)
  {
    // The harmonic mean of mixed positive and negative numbers (and 0 as well) is not well-defined
    // so we assert that the input values shall be positive.
#ifndef NDEBUG
    if (value1 * value2 <= 0)
      mooseWarning("Input values (" + Moose::stringify(MetaPhysicL::raw_value(value1)) + " & " +
                   Moose::stringify(MetaPhysicL::raw_value(value2)) +
                   ") must be of the same sign for harmonic interpolation");
#endif
    return 1.0 / (weight1 / value1 + weight2 / value2);
  }
  // For vectors (ADRealVectorValue, VectorValue), we take the component-wise harmonic mean
  else if constexpr (libMesh::TensorTools::TensorTraits<T1>::rank == 1)
  {
    typename libMesh::CompareTypes<T1, T2>::supertype result;
    for (const auto i : make_range(Moose::dim))
    {
#ifndef NDEBUG
      if (value1(i) * value2(i) <= 0)
        mooseWarning("Component " + std::to_string(i) + " of input values (" +
                     Moose::stringify(MetaPhysicL::raw_value(value1(i))) + " & " +
                     Moose::stringify(MetaPhysicL::raw_value(value2(i))) +
                     ") must be of the same sign for harmonic interpolation");
#endif
      result(i) = 1.0 / (weight1 / value1(i) + weight2 / value2(i));
    }
    return result;
  }
  // For tensors (ADRealTensorValue, TensorValue), similarly to the vectors,
  // we take the component-wise harmonic mean instead of the matrix-inverse approach
  else if constexpr (libMesh::TensorTools::TensorTraits<T1>::rank == 2)
  {
    typename libMesh::CompareTypes<T1, T2>::supertype result;
    for (const auto i : make_range(Moose::dim))
      for (const auto j : make_range(Moose::dim))
      {
#ifndef NDEBUG
        if (value1(i, j) * value2(i, j) <= 0)
          mooseWarning("Component (" + std::to_string(i) + "," + std::to_string(j) +
                       ") of input values (" +
                       Moose::stringify(MetaPhysicL::raw_value(value1(i, j))) + " & " +
                       Moose::stringify(MetaPhysicL::raw_value(value2(i, j))) +
                       ") must be of the same sign for harmonic interpolation");
#endif
        result(i, j) = 1.0 / (weight1 / value1(i, j) + weight2 / value2(i, j));
      }
    return result;
  }
  // We ran out of options, harmonic mean is not supported for other types at the moment
  else
    // This line is supposed to throw an error when the user tries to compile this function with
    // types that are not supported. This is the reason we needed the always_false function. Hope as
    // C++ gets nicer, we can do this in a nicer way.
    static_assert(Moose::always_false<T1>,
                  "Harmonic interpolation is not implemented for the used type!");
}

Referenced by harmonicInterpolation(), and interpolate().

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:

\phi_f = c_1 * \phi_{F1} + c_2 * \phi_{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

m	The interpolation method
fi	The face information
one_is_elem	Whether fi.elem() == F1
face_flux	The 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 117 of file MathFVUtils.h.

{
  switch (m)
  {
    case InterpMethod::Average:
    case InterpMethod::SkewCorrectedAverage:
    {
      if (one_is_elem)
        return std::make_pair(fi.gC(), 1. - fi.gC());
      else
        return std::make_pair(1. - fi.gC(), fi.gC());
    }
 
    case InterpMethod::Upwind:
    {
      if ((face_flux > 0) == one_is_elem)
        return std::make_pair(1., 0.);
      else
        return std::make_pair(0., 1.);
    }
 
    default:
      mooseError("Unrecognized interpolation method");
  }
}

Referenced by containerInterpolate(), harmonicInterpolation(), interpCoeffsAndAdvected(), interpolate(), interpolate(), interpolate(), linearInterpolation(), and skewCorrectedLinearInterpolation().

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:

\phi_f = c_upwind * \phi_{upwind} + c_downwind * \phi_{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 490 of file MathFVUtils.h.

{
  // Using beta, w_f, g nomenclature from Greenshields
  const auto beta = limiter(phi_upwind,
                            phi_downwind,
                            grad_phi_upwind,
                            grad_phi_face,
                            fi_elem_is_upwind ? fi.dCN() : Point(-fi.dCN()),
                            max_value,
                            min_value,
                            &fi,
                            fi_elem_is_upwind);
 
  const auto w_f = fi_elem_is_upwind ? fi.gC() : (1. - fi.gC());
 
  const auto g = beta * (1. - w_f);
 
  return std::make_pair(1. - g, g);
}

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

T	The data type for the values being interpolated. This is typically a scalar type.
FEK	Enumeration type FunctorEvaluationKind with a default value of FunctorEvaluationKind::Value. This determines the kind of evaluation the functor will perform.
Enable	A 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

functor	An 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.
face	An argument representing the face in the computational domain. The face provides access to neighboring elements and limiter type.
time	An 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 803 of file MathFVUtils.h.

{
  // Ensure only supported FunctorEvaluationKinds are used
  static_assert((FEK == FunctorEvaluationKind::Value) || (FEK == FunctorEvaluationKind::Dot),
                "Only Value and Dot evaluations currently supported");
 
  // Determine the gradient evaluation kind
  constexpr FunctorEvaluationKind GFEK = FunctorGradientEvaluationKind<FEK>::value;
  typedef typename FunctorBase<T>::GradientType GradientType;
  static const GradientType zero(0);
 
  mooseAssert(face.fi, "this must be non-null");
 
  // Construct the limiter based on the face limiter type
  const auto limiter = Limiter<typename LimiterValueType<T>::value_type>::build(face.limiter_type);
 
  // Determine upwind and downwind arguments based on the face element
  const auto upwind_arg = face.elem_is_upwind ? face.makeElem() : face.makeNeighbor();
  const auto downwind_arg = face.elem_is_upwind ? face.makeNeighbor() : face.makeElem();
 
  // Evaluate the functor at the upwind and downwind locations
  auto phi_upwind = functor.template genericEvaluate<FEK>(upwind_arg, time);
  auto phi_downwind = functor.template genericEvaluate<FEK>(downwind_arg, time);
 
  // Initialize the interpolation coefficients
  std::pair<T, T> interp_coeffs;
 
  // Compute interpolation coefficients for Upwind or CentralDifference limiters
  if (face.limiter_type == LimiterType::Upwind ||
      face.limiter_type == LimiterType::CentralDifference)
    interp_coeffs = interpCoeffs(*limiter,
                                 phi_upwind,
                                 phi_downwind,
                                 &zero,
                                 &zero,
                                 std::numeric_limits<Real>::max(),
                                 std::numeric_limits<Real>::min(),
                                 *face.fi,
                                 face.elem_is_upwind);
  else
  {
    // Determine the time argument for the limiter
    auto * time_arg = face.state_limiter;
    if (!time_arg)
    {
      static Moose::StateArg temp_state(0, Moose::SolutionIterationType::Time);
      time_arg = &temp_state;
    }
 
    Real max_value(std::numeric_limits<Real>::min()), min_value(std::numeric_limits<Real>::max());
 
    // Compute min-max values for min-max limiters
    if (face.limiter_type == LimiterType::Venkatakrishnan || face.limiter_type == LimiterType::SOU)
      std::tie(max_value, min_value) = computeMinMaxValue(functor, face, *time_arg);
 
    // Evaluate the gradient of the functor at the upwind and downwind locations
    const auto grad_phi_upwind = functor.template genericEvaluate<GFEK>(upwind_arg, *time_arg);
    const auto grad_phi_face = functor.template genericEvaluate<GFEK>(face, *time_arg);
 
    // Compute the interpolation coefficients using the specified limiter
    interp_coeffs = interpCoeffs(*limiter,
                                 phi_upwind,
                                 phi_downwind,
                                 &grad_phi_upwind,
                                 &grad_phi_face,
                                 max_value,
                                 min_value,
                                 *face.fi,
                                 face.elem_is_upwind);
  }
 
  // Return the interpolation coefficients and advected values
  if (face.elem_is_upwind)
    return std::make_pair(std::move(interp_coeffs),
                          std::make_pair(std::move(phi_upwind), std::move(phi_downwind)));
  else
    return std::make_pair(
        std::make_pair(std::move(interp_coeffs.second), std::move(interp_coeffs.first)),
        std::make_pair(std::move(phi_downwind), std::move(phi_upwind)));
}

interpolate() [1/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

T	The data type for the vector values being interpolated. This is typically a scalar type.

Parameters

functor	An 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.
face	An argument representing the face in the computational domain. The face provides access to neighboring elements and limiter type.
time	An 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 996 of file MathFVUtils.h.

{
  // Define a zero gradient vector for initialization
  static const libMesh::VectorValue<T> grad_zero(0);
 
  mooseAssert(face.fi, "this must be non-null");
 
  // Construct the limiter based on the face limiter type
  const auto limiter = Limiter<typename LimiterValueType<T>::value_type>::build(face.limiter_type);
 
  // Determine upwind and downwind arguments based on the face element
  const auto upwind_arg = face.elem_is_upwind ? face.makeElem() : face.makeNeighbor();
  const auto downwind_arg = face.elem_is_upwind ? face.makeNeighbor() : face.makeElem();
  auto phi_upwind = functor(upwind_arg, time);
  auto phi_downwind = functor(downwind_arg, time);
 
  // Initialize the return vector value
  libMesh::VectorValue<T> ret;
  T coeff_upwind, coeff_downwind;
 
  // Compute interpolation coefficients and advected values for Upwind or CentralDifference limiters
  if (face.limiter_type == LimiterType::Upwind ||
      face.limiter_type == LimiterType::CentralDifference)
  {
    for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
    {
      const auto &component_upwind = phi_upwind(i), component_downwind = phi_downwind(i);
      std::tie(coeff_upwind, coeff_downwind) = interpCoeffs(*limiter,
                                                            component_upwind,
                                                            component_downwind,
                                                            &grad_zero,
                                                            &grad_zero,
                                                            std::numeric_limits<Real>::max(),
                                                            std::numeric_limits<Real>::min(),
                                                            *face.fi,
                                                            face.elem_is_upwind);
      ret(i) = coeff_upwind * component_upwind + coeff_downwind * component_downwind;
    }
  }
  else
  {
    // Determine the time argument for the limiter
    auto * time_arg = face.state_limiter;
    if (!time_arg)
    {
      static Moose::StateArg temp_state(0, Moose::SolutionIterationType::Time);
      time_arg = &temp_state;
    }
 
    // Evaluate the gradient of the functor at the upwind and downwind locations
    const auto & grad_phi_upwind = functor.gradient(upwind_arg, *time_arg);
    const auto & grad_phi_downwind = functor.gradient(downwind_arg, *time_arg);
 
    const auto coeffs = interpCoeffs(InterpMethod::Average, *face.fi, face.elem_is_upwind);
 
    // Compute interpolation coefficients and advected values for each component
    for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
    {
      const auto &component_upwind = phi_upwind(i), component_downwind = phi_downwind(i);
      const libMesh::VectorValue<T> &grad_upwind = grad_phi_upwind.row(i),
                                    grad_face = coeffs.first * grad_phi_upwind.row(i) +
                                                coeffs.second * grad_phi_downwind.row(i);
 
      // Initialize min-max values
      Real max_value(std::numeric_limits<Real>::min()), min_value(std::numeric_limits<Real>::max());
      if (face.limiter_type == LimiterType::Venkatakrishnan ||
          face.limiter_type == LimiterType::SOU)
        std::tie(max_value, min_value) = computeMinMaxValue(functor, face, *time_arg, i);
 
      // Perform full limited interpolation for the component
      ret(i) = fullLimitedInterpolation(*limiter,
                                        component_upwind,
                                        component_downwind,
                                        &grad_upwind,
                                        &grad_face,
                                        max_value,
                                        min_value,
                                        face);
    }
  }
 
  // Return the interpolated vector value
  return ret;
}

interpolate() [2/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 1175 of file MathFVUtils.h.

{
  return containerInterpolate(functor, face, time);
}

interpolate() [3/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 1166 of file MathFVUtils.h.

{
  return containerInterpolate(functor, face, time);
}

interpolate() [4/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

T	The data type for the values being interpolated. This is typically a scalar type.
FEK	Enumeration type FunctorEvaluationKind with a default value of FunctorEvaluationKind::Value. This determines the kind of evaluation the functor will perform.
Enable	A 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

functor	An 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.
face	An argument representing the face in the computational domain. The face provides access to neighboring elements and limiter type.
time	An 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 911 of file MathFVUtils.h.

{
  // Ensure only supported FunctorEvaluationKinds are used
  static_assert((FEK == FunctorEvaluationKind::Value) || (FEK == FunctorEvaluationKind::Dot),
                "Only Value and Dot evaluations currently supported");
 
  // Special handling for central differencing as it is the only interpolation method which
  // currently supports skew correction
  if (face.limiter_type == LimiterType::CentralDifference)
    return linearInterpolation<T, FEK>(functor, face, time);
 
  if (face.limiter_type == LimiterType::Upwind ||
      face.limiter_type == LimiterType::CentralDifference)
  {
    // Compute interpolation coefficients and advected values
    const auto [interp_coeffs, advected] = interpCoeffsAndAdvected<T, FEK>(functor, face, time);
    // Return the interpolated value
    return interp_coeffs.first * advected.first + interp_coeffs.second * advected.second;
  }
  else
  {
    // Determine the gradient evaluation kind
    constexpr FunctorEvaluationKind GFEK = FunctorGradientEvaluationKind<FEK>::value;
    typedef typename FunctorBase<T>::GradientType GradientType;
    static const GradientType zero(0);
    const auto & limiter =
        Limiter<typename LimiterValueType<T>::value_type>::build(face.limiter_type);
 
    // Determine upwind and downwind arguments based on the face element
    const auto & upwind_arg = face.elem_is_upwind ? face.makeElem() : face.makeNeighbor();
    const auto & downwind_arg = face.elem_is_upwind ? face.makeNeighbor() : face.makeElem();
    const auto & phi_upwind = functor.template genericEvaluate<FEK>(upwind_arg, time);
    const auto & phi_downwind = functor.template genericEvaluate<FEK>(downwind_arg, time);
 
    // Determine the time argument for the limiter
    auto * time_arg = face.state_limiter;
    if (!time_arg)
    {
      static Moose::StateArg temp_state(0, Moose::SolutionIterationType::Time);
      time_arg = &temp_state;
    }
 
    // Initialize min-max values
    Real max_value(std::numeric_limits<Real>::min()), min_value(std::numeric_limits<Real>::max());
    if (face.limiter_type == LimiterType::Venkatakrishnan || face.limiter_type == LimiterType::SOU)
      std::tie(max_value, min_value) = computeMinMaxValue(functor, face, *time_arg);
 
    // Evaluate the gradient of the functor at the upwind and downwind locations
    const auto & grad_phi_upwind = functor.template genericEvaluate<GFEK>(upwind_arg, *time_arg);
    const auto & grad_phi_face = functor.template genericEvaluate<GFEK>(face, *time_arg);
 
    // Perform full limited interpolation and return the interpolated value
    return fullLimitedInterpolation(*limiter,
                                    phi_upwind,
                                    phi_downwind,
                                    &grad_phi_upwind,
                                    &grad_phi_face,
                                    max_value,
                                    min_value,
                                    face);
  }
}

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 550 of file MathFVUtils.h.

{
  mooseAssert(limiter.constant() || grad_phi_upwind,
              "Non-null gradient only supported for constant limiters.");
 
  const libMesh::VectorValue<T> * const gradient_example = nullptr;
  libMesh::VectorValue<T> ret;
  for (const auto i : make_range(unsigned(LIBMESH_DIM)))
  {
    if (grad_phi_upwind)
    {
      const libMesh::VectorValue<T> gradient = grad_phi_upwind->row(i);
      ret(i) =
          interpolate(limiter, phi_upwind(i), phi_downwind(i), &gradient, fi, fi_elem_is_upwind);
    }
    else
      ret(i) = interpolate(
          limiter, phi_upwind(i), phi_downwind(i), gradient_example, fi, fi_elem_is_upwind);
  }
 
  return ret;
}

interpolate() [6/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 525 of file MathFVUtils.h.

{
  auto pr =
      interpCoeffs(limiter,
                   phi_upwind,
                   phi_downwind,
                   grad_phi_upwind,
                   /*grad_phi_face*/ static_cast<const libMesh::VectorValue<Scalar> *>(nullptr),
                   /* max_value */ std::numeric_limits<Real>::max(),
                   /* min_value */ std::numeric_limits<Real>::min(),
                   fi,
                   fi_elem_is_upwind);
  return pr.first * phi_upwind + pr.second * phi_downwind;
}

interpolate() [7/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 298 of file MathFVUtils.h.

{
  switch (m)
  {
    case InterpMethod::Average:
    case InterpMethod::SkewCorrectedAverage:
      result = linearInterpolation(value1, value2, fi, one_is_elem);
      break;
    case InterpMethod::HarmonicAverage:
      result = harmonicInterpolation(value1, value2, fi, one_is_elem);
      break;
    default:
      mooseError("unsupported interpolation method for interpolate() function");
  }
}

Referenced by FVOrthogonalDiffusion::computeQpResidual(), MooseVariableFV< OutputType >::evaluate(), greenGaussGradient(), interpolate(), and FVInterfaceKernel::interpolateValue().

interpolate() [8/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 416 of file MathFVUtils.h.

{
  const auto coeffs = interpCoeffs(m, fi, one_is_elem, advector * fi.normal());
  result = coeffs.first * value1 + coeffs.second * value2;
}

interpolate() [9/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 369 of file MathFVUtils.h.

{
  switch (m)
  {
    case InterpMethod::Average:
      result = linearInterpolation(fi_elem_advected * fi_elem_advector,
                                   fi_neighbor_advected * fi_neighbor_advector,
                                   fi,
                                   true);
      break;
    case InterpMethod::Upwind:
    {
      const auto face_advector = linearInterpolation(MetaPhysicL::raw_value(fi_elem_advector),
                                                     MetaPhysicL::raw_value(fi_neighbor_advector),
                                                     fi,
                                                     true);
      Real elem_coeff, neighbor_coeff;
      if (face_advector * fi.normal() > 0)
        elem_coeff = 1, neighbor_coeff = 0;
      else
        elem_coeff = 0, neighbor_coeff = 1;
 
      result = elem_coeff * fi_elem_advected * fi_elem_advector +
               neighbor_coeff * fi_neighbor_advected * fi_neighbor_advector;
      break;
    }
    default:
      mooseError("unsupported interpolation method for FVFaceInterface::interpolate");
  }
}

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.

{
  return MooseEnum("average upwind sou min_mod vanLeer quick venkatakrishnan skewness-corrected",
                   "upwind");
}

Referenced by advectedInterpolationParameter().

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.

{
  switch (interp_method)
  {
    case InterpMethod::Average:
    case InterpMethod::SkewCorrectedAverage:
      return LimiterType::CentralDifference;
 
    case InterpMethod::Upwind:
      return LimiterType::Upwind;
 
    case InterpMethod::VanLeer:
      return LimiterType::VanLeer;
 
    case InterpMethod::MinMod:
      return LimiterType::MinMod;
 
    case InterpMethod::SOU:
      return LimiterType::SOU;
 
    case InterpMethod::QUICK:
      return LimiterType::QUICK;
 
    case InterpMethod::Venkatakrishnan:
      return LimiterType::Venkatakrishnan;
 
    default:
      mooseError("Unrecognized interpolation method type.");
  }
}

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

linearInterpolation() [1/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 325 of file MathFVUtils.h.

{
  static_assert((FEK == FunctorEvaluationKind::Value) || (FEK == FunctorEvaluationKind::Dot),
                "Only Value and Dot evaluations currently supported");
  mooseAssert(face.limiter_type == LimiterType::CentralDifference,
              "this interpolation method is meant for linear interpolations");
 
  mooseAssert(face.fi,
              "We must have a non-null face_info in order to prepare our ElemFromFace tuples");
 
  constexpr FunctorEvaluationKind GFEK = FunctorGradientEvaluationKind<FEK>::value;
 
  const auto elem_arg = face.makeElem();
  const auto neighbor_arg = face.makeNeighbor();
 
  const auto elem_eval = functor.template genericEvaluate<FEK>(elem_arg, time);
  const auto neighbor_eval = functor.template genericEvaluate<FEK>(neighbor_arg, time);
 
  if (face.correct_skewness)
  {
    // This condition ensures that the recursive algorithm (face_center->
    // face_gradient -> cell_gradient -> face_center -> ...) terminates after
    // one loop. It is hardcoded to one loop at this point since it yields
    // 2nd order accuracy on skewed meshes with the minimum additional effort.
    FaceArg new_face(face);
    new_face.correct_skewness = false;
    const auto surface_gradient = functor.template genericEvaluate<GFEK>(new_face, time);
 
    return skewCorrectedLinearInterpolation(
        elem_eval, neighbor_eval, surface_gradient, *face.fi, true);
  }
  else
    return linearInterpolation(elem_eval, neighbor_eval, *face.fi, true);
}

linearInterpolation() [2/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 153 of file MathFVUtils.h.

{
  mooseAssert(interp_method == InterpMethod::Average ||
                  interp_method == InterpMethod::SkewCorrectedAverage,
              "The selected interpolation function only works with average or skewness-corrected "
              "average options!");
  const auto coeffs = interpCoeffs(interp_method, fi, one_is_elem);
  return coeffs.first * value1 + coeffs.second * value2;
}

Referenced by LinearFVGradientReader::gradient(), FVGreenGaussGradient::internalFaceValues(), interpolate(), interpolate(), linearInterpolation(), and MooseVariableFV< OutputType >::uncorrectedAdGradSln().

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.

{
  mooseAssert(elem.active(), "We should never call this method with an inactive element");
 
  for (const auto side : elem.side_index_range())
  {
    const Elem * const candidate_neighbor = elem.neighbor_ptr(side);
 
    bool elem_has_info = elemHasFaceInfo(elem, candidate_neighbor);
 
    std::set<const Elem *> neighbors;
 
    const bool inactive_neighbor_detected =
        candidate_neighbor ? !candidate_neighbor->active() : false;
 
    // See MooseMesh::buildFaceInfo for corresponding checks/additions of FaceInfo
    if (inactive_neighbor_detected)
    {
      // We must be next to an element that has been refined
      mooseAssert(candidate_neighbor->has_children(), "We should have children");
 
      const auto candidate_neighbor_side = candidate_neighbor->which_neighbor_am_i(&elem);
 
      for (const auto child_num : make_range(candidate_neighbor->n_children()))
        if (candidate_neighbor->is_child_on_side(child_num, candidate_neighbor_side))
        {
          const Elem * const child = candidate_neighbor->child_ptr(child_num);
          mooseAssert(child->level() - elem.level() == 1, "The math doesn't work out here.");
          mooseAssert(child->has_neighbor(&elem), "Elem should be a neighbor of this child.");
          mooseAssert(child->active(),
                      "We shouldn't have greater than a face mismatch level of one");
          neighbors.insert(child);
        }
    }
    else
      neighbors.insert(candidate_neighbor);
 
    for (const Elem * const neighbor : neighbors)
    {
      const FaceInfo * const fi =
          elem_has_info ? mesh.faceInfo(&elem, side)
                        : mesh.faceInfo(neighbor, neighbor->which_neighbor_am_i(&elem));
 
      mooseAssert(fi, "We should have found a FaceInfo");
      mooseAssert(elem_has_info ? &elem == &fi->elem() : &elem == fi->neighborPtr(),
                  "Doesn't seem like we understand how this FaceInfo thing is working");
      if (neighbor)
      {
        mooseAssert(neighbor != libMesh::remote_elem,
                    "Remote element detected. This indicates that you have insufficient geometric "
                    "ghosting. Please contact your application developers.");
        mooseAssert(elem_has_info ? neighbor == fi->neighborPtr() : neighbor == &fi->elem(),
                    "Doesn't seem like we understand how this FaceInfo thing is working");
      }
 
      const Point elem_normal = elem_has_info ? fi->normal() : Point(-fi->normal());
 
      Real coord;
      MooseMeshUtils::coordTransformFactor(fi->faceCentroid(), coord, coord_type, rz_radial_coord);
 
      const Point surface_vector = elem_normal * fi->faceArea() * coord;
 
      act(elem, neighbor, fi, surface_vector, coord, elem_has_info);
    }
  }
}

Referenced by greenGaussGradient().

onBoundary() [1/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.

{
  if (!fi.neighborPtr())
    // We're on the exterior boundary
    return true;
 
  if (subs.empty())
    // The face is internal and our functor lives on all subdomains
    return false;
 
  const auto sub_count =
      subs.count(fi.elem().subdomain_id()) + subs.count(fi.neighbor().subdomain_id());
 
  switch (sub_count)
  {
    case 0:
      mooseError("We should not be calling isExtrapolatedBoundaryFace on a functor that doesn't "
                 "live on either of the face information's neighboring elements");
 
    case 1:
      // We only live on one of the subs
      return true;
 
    case 2:
      // We live on both of the subs
      return false;
 
    default:
      mooseError("There should be no other sub_count options");
  }
}

onBoundary() [2/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 1187 of file MathFVUtils.h.

{
  const bool internal = fi.neighborPtr() && obj.hasBlocks(fi.elemSubdomainID()) &&
                        obj.hasBlocks(fi.neighborSubdomainID());
  return !internal;
}

Referenced by FVFluxKernel::onBoundary().

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.

This is a branchless version, avoiding if-else statements to enable vectorization.

Definition at line 461 of file MathFVUtils.h.

{
  const Scalar denom = phiD - phiC;
  const Scalar grad_dot = gradC * dCD;
 
  // This is an arbitrary number here, when we start seeing convergence issues we
  // can tune this but so far this has shown okay results.
  constexpr Real eps = 1e-10;
  const Real denom_sign =
      std::copysign(1.0,
                    MetaPhysicL::raw_value(denom) +
                        std::numeric_limits<Real>::min() * MetaPhysicL::raw_value(grad_dot));
  const Scalar safe = denom + denom_sign * eps;
  return 2.0 * grad_dot / safe - 1.0;
}

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

selectInterpolationMethod()

InterpMethod Moose::FV::selectInterpolationMethod

(

const std::string &

interp_method

)

Definition at line 81 of file MathFVUtils.C.

{
  if (interp_method == "average")
    return InterpMethod::Average;
  else if (interp_method == "harmonic")
    return InterpMethod::HarmonicAverage;
  else if (interp_method == "skewness-corrected")
    return InterpMethod::SkewCorrectedAverage;
  else if (interp_method == "upwind")
    return InterpMethod::Upwind;
  else if (interp_method == "rc")
    return InterpMethod::RhieChow;
  else if (interp_method == "vanLeer")
    return InterpMethod::VanLeer;
  else if (interp_method == "min_mod")
    return InterpMethod::MinMod;
  else if (interp_method == "sou")
    return InterpMethod::SOU;
  else if (interp_method == "quick")
    return InterpMethod::QUICK;
  else if (interp_method == "venkatakrishnan")
    return InterpMethod::Venkatakrishnan;
  else
    mooseError("Interpolation method ",
               interp_method,
               " is not currently an option in Moose::FV::selectInterpolationMethod");
}

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

setInterpolationMethod()

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

Sets one interpolation method.

Parameters

obj	The MooseObject with input parameters to query
interp_method	The interpolation method we will set
param_name	The 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.

{
  bool need_more_ghosting = false;
 
  const auto & interp_method_in = obj.getParam<MooseEnum>(param_name);
  interp_method = selectInterpolationMethod(interp_method_in);
 
  if (interp_method == InterpMethod::SOU || interp_method == InterpMethod::MinMod ||
      interp_method == InterpMethod::VanLeer || interp_method == InterpMethod::QUICK ||
      interp_method == InterpMethod::Venkatakrishnan)
    need_more_ghosting = true;
 
  return need_more_ghosting;
}

Referenced by FVAdvection::FVAdvection().

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 277 of file MathFVUtils.h.

{
  const auto coeffs = interpCoeffs(InterpMethod::SkewCorrectedAverage, fi, one_is_elem);
 
  auto value = (coeffs.first * value1 + coeffs.second * value2) +
               face_gradient * fi.skewnessCorrectionVector();
  return value;
}

Referenced by linearInterpolation().

Variable Documentation

moose_limiter_type

const MooseEnum Moose::FV::moose_limiter_type

extern