MathFVUtils.h

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#pragma once

#include "MooseError.h"
#include "FaceInfo.h"
#include "Limiter.h"
#include "MathUtils.h"
#include "MooseFunctor.h"
#include "metaphysicl/raw_type.h"
#include "libmesh/compare_types.h"
#include "libmesh/elem.h"
#include <cmath>
#include <limits>
#include <tuple>

template <typename>
class MooseVariableFV;
class FaceArgInterface;

namespace Moose
{
namespace FV
{
enum class InterpMethod
{
  Average,
  HarmonicAverage,
  SkewCorrectedAverage,
  Upwind,
  // Rhie-Chow
  RhieChow,
  VanLeer,
  MinMod,
  SOU,
  QUICK,
  Venkatakrishnan
};

enum class LinearFVComputationMode
{
  RHS,
  Matrix,
  FullSystem
};

MooseEnum interpolationMethods();

InputParameters advectedInterpolationParameter();

/*
 * Converts from the interpolation method to the interpolation enum.
 * This routine is here in lieu of using a MooseEnum for InterpMethod
 * @param interp_method the name of the interpolation method
 * @return the interpolation method
 */
InterpMethod selectInterpolationMethod(const std::string & interp_method);

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

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)
{
  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");
  }
}

template <typename T, typename T2>
typename 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)
{
  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;
}

template <typename T1, typename T2>
typename libMesh::CompareTypes<T1, T2>::supertype
harmonicInterpolation(const T1 & value1,
                      const T2 & value2,
                      const FaceInfo & fi,
                      const bool one_is_elem)
{
  // 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!");

  // Fetch the interpolation coefficients, we use the exact same coefficients as for a simple
  // weighted average
  const auto coeffs = interpCoeffs(InterpMethod::Average, fi, one_is_elem);

  // 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 / (coeffs.first / value1 + coeffs.second / 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 / (coeffs.first / value1(i) + coeffs.second / 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 / (coeffs.first / value1(i, j) + coeffs.second / 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!");
}

template <typename T, typename T2, typename T3>
typename libMesh::CompareTypes<T, T2>::supertype
skewCorrectedLinearInterpolation(const T & value1,
                                 const T2 & value2,
                                 const T3 & face_gradient,
                                 const FaceInfo & fi,
                                 const bool one_is_elem)
{
  const auto coeffs = interpCoeffs(InterpMethod::SkewCorrectedAverage, fi, one_is_elem);

  auto value = (coeffs.first * value1 + coeffs.second * value2) +
               face_gradient * fi.skewnessCorrectionVector();
  return value;
}

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)
{
  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");
  }
}

template <typename T, FunctorEvaluationKind FEK = FunctorEvaluationKind::Value>
T
linearInterpolation(const FunctorBase<T> & functor, const FaceArg & face, const StateArg & time)
{
  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);
}

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)
{
  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");
  }
}

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)
{
  const auto coeffs = interpCoeffs(m, fi, one_is_elem, advector * fi.normal());
  result = coeffs.first * value1 + coeffs.second * value2;
}

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

template <typename Scalar, typename Vector>
Scalar
rF(const Scalar & phiC, const Scalar & phiD, const Vector & gradC, const RealVectorValue & dCD)
{
  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;
}

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)
{
  // 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);
}

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)
{
  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;
}

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)
{
  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;
}

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)
{
  // 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;
}

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)
{
  // 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);
}

template <typename T>
std::pair<Real, Real>
computeMinMaxValue(const FunctorBase<VectorValue<T>> & functor,
                   const FaceArg & face,
                   const StateArg & time,
                   const unsigned int & component)
{
  // 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);
}

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)
{
  // 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)));
}

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)
{
  // 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);
  }
}

template <typename T>
libMesh::VectorValue<T>
interpolate(const FunctorBase<libMesh::VectorValue<T>> & functor,
            const FaceArg & face,
            const StateArg & time)
{
  // 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;
}

template <typename T>
T
containerInterpolate(const FunctorBase<T> & functor, const FaceArg & face, const StateArg & time)
{
  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;
}

template <typename T>
std::vector<T>
interpolate(const FunctorBase<std::vector<T>> & functor,
            const FaceArg & face,
            const StateArg & time)
{
  return containerInterpolate(functor, face, 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)
{
  return containerInterpolate(functor, face, time);
}

template <typename SubdomainRestrictable>
bool
onBoundary(const SubdomainRestrictable & obj, const FaceInfo & fi)
{
  const bool internal = fi.neighborPtr() && obj.hasBlocks(fi.elemSubdomainID()) &&
                        obj.hasBlocks(fi.neighborSubdomainID());
  return !internal;
}

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