Merge pull request #26774 from tanoret/fv_dt_faces

extending MooseVariableFV to face and qp for transient two-phase flow…

framework/include/base/MooseFunctor.h

@@ -29,6 +29,92 @@
 
 namespace Moose
 {
+/**
+ * An enumeration of possible functor evaluation kinds. The available options are value, gradient,
+ * time derivative (dot), and gradient of time derivative (gradDot)
+ */
+enum class FunctorEvaluationKind
+{
+  Value,
+  Gradient,
+  Dot,
+  GradDot
+};
+
+/**
+ * A structure that defines the return type of a functor based on the type of the functor and the
+ * requested evaluation kind, e.g. value, gradient, time derivative, or gradient of time derivative
+ */
+template <typename, FunctorEvaluationKind>
+struct FunctorReturnType;
+
+/**
+ * The return type for a value evaluation is just the type of the functor
+ */
+template <typename T>
+struct FunctorReturnType<T, FunctorEvaluationKind::Value>
+{
+  typedef T type;
+};
+
+/**
+ * The return type of a gradient evaluation is the rank increment of a value return type. So if the
+ * value type is Real, then a gradient will be a VectorValue<Real>. This also allows for containers
+ * of mathematical types. So if a value type is std::vector<Real>, then the gradient type will be
+ * std::vector<VectorValue<Real>>
+ */
+template <typename T>
+struct FunctorReturnType<T, FunctorEvaluationKind::Gradient>
+{
+  typedef typename MetaPhysicL::ReplaceAlgebraicType<
+      T,
+      typename TensorTools::IncrementRank<typename MetaPhysicL::ValueType<T>::type>::type>::type
+      type;
+};
+
+/**
+ * The return type of a time derivative evaluation is the same as the value type
+ */
+template <typename T>
+struct FunctorReturnType<T, FunctorEvaluationKind::Dot>
+{
+  typedef T type;
+};
+
+/**
+ * The return type of a gradient of time derivative evaluation is the same as the gradient type
+ */
+template <typename T>
+struct FunctorReturnType<T, FunctorEvaluationKind::GradDot>
+{
+  typedef typename FunctorReturnType<T, FunctorEvaluationKind::Gradient>::type type;
+};
+
+/**
+ * This structure takes an evaluation kind as a template argument and defines a constant expression
+ * indicating the associated gradient kind
+ */
+template <FunctorEvaluationKind>
+struct FunctorGradientEvaluationKind;
+
+/**
+ * The gradient kind associated with a value is simply the gradient
+ */
+template <>
+struct FunctorGradientEvaluationKind<FunctorEvaluationKind::Value>
+{
+  static constexpr FunctorEvaluationKind value = FunctorEvaluationKind::Gradient;
+};
+
+/**
+ * The gradient kind associated with a time derivative is the gradient of the time derivative
+ */
+template <>
+struct FunctorGradientEvaluationKind<FunctorEvaluationKind::Dot>
+{
+  static constexpr FunctorEvaluationKind value = FunctorEvaluationKind::GradDot;
+};
+
 /**
  * Abstract base class that can be used to hold collections of functors
  */
@@ -52,7 +138,6 @@ class FunctorBase : public FunctorAbstract
 {
 public:
   using FunctorType = FunctorBase<T>;
-  using FunctorReturnType = T;
   using ValueType = T;
   /// This rigmarole makes it so that a user can create functors that return containers (std::vector,
   /// std::array). This logic will make it such that if a user requests a functor type T that is a
@@ -61,9 +146,7 @@ class FunctorBase : public FunctorAbstract
   /// std::vector<Real>, then GradientType will be std::vector<VectorValue<Real>>. As another
   /// example: T = std::array<VectorValue<Real>, 1> -> GradientType = std::array<TensorValue<Real>,
   /// 1>
-  using GradientType = typename MetaPhysicL::ReplaceAlgebraicType<
-      T,
-      typename TensorTools::IncrementRank<typename MetaPhysicL::ValueType<T>::type>::type>::type;
+  using GradientType = typename FunctorReturnType<T, FunctorEvaluationKind::Gradient>::type;
   using DotType = ValueType;
 
   virtual ~FunctorBase() = default;
@@ -73,6 +156,14 @@ class FunctorBase : public FunctorAbstract
   {
   }
 
+  /**
+   * Perform a generic evaluation based on the supplied template argument \p FET and supplied
+   * spatial and temporal arguments
+   */
+  template <FunctorEvaluationKind FET, typename Space, typename State>
+  typename FunctorReturnType<T, FET>::type genericEvaluate(const Space & r,
+                                                           const State & state) const;
+
   /// Return the functor name
   const MooseFunctorName & functorName() const { return _functor_name; }
 
@@ -870,6 +961,21 @@ FunctorBase<T>::hasFaceSide(const FaceInfo & fi, const bool fi_elem_side) const
     return fi.neighborPtr() && hasBlocks(fi.neighbor().subdomain_id());
 }
 
+template <typename T>
+template <FunctorEvaluationKind FET, typename Space, typename State>
+typename FunctorReturnType<T, FET>::type
+FunctorBase<T>::genericEvaluate(const Space & r, const State & state) const
+{
+  if constexpr (FET == FunctorEvaluationKind::Value)
+    return (*this)(r, state);
+  else if constexpr (FET == FunctorEvaluationKind::Gradient)
+    return gradient(r, state);
+  else if constexpr (FET == FunctorEvaluationKind::Dot)
+    return dot(r, state);
+  else
+    return gradDot(r, state);
+}
+
 /**
  * A non-templated base class for functors that allow an owner object to hold
  * different class template instantiations of \p Functor in a single container
@@ -1139,7 +1245,6 @@ class ConstantFunctor final : public FunctorBase<T>
 {
 public:
   using typename FunctorBase<T>::FunctorType;
-  using typename FunctorBase<T>::FunctorReturnType;
   using typename FunctorBase<T>::ValueType;
   using typename FunctorBase<T>::GradientType;
   using typename FunctorBase<T>::DotType;
@@ -1202,7 +1307,6 @@ class NullFunctor final : public FunctorBase<T>
 {
 public:
   using typename FunctorBase<T>::FunctorType;
-  using typename FunctorBase<T>::FunctorReturnType;
   using typename FunctorBase<T>::ValueType;
   using typename FunctorBase<T>::GradientType;
   using typename FunctorBase<T>::DotType;

framework/include/utils/MathFVUtils.h

@@ -294,19 +294,26 @@ interpolate(InterpMethod m,
  * perform a possibly skew-corrected linear interpolation by evaluating the supplied functor with
  * the provided functor face argument
  */
-template <typename T>
+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->
@@ -315,14 +322,13 @@ linearInterpolation(const FunctorBase<T> & functor, const FaceArg & face, const
     // 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.gradient(new_face, time);
+    const auto surface_gradient = functor.template genericEvaluate<GFEK>(new_face, time);
 
     return skewCorrectedLinearInterpolation(
-        functor(elem_arg, time), functor(neighbor_arg, time), surface_gradient, *face.fi, true);
+        elem_eval, neighbor_eval, surface_gradient, *face.fi, true);
   }
   else
-    return linearInterpolation(
-        functor(elem_arg, time), functor(neighbor_arg, time), *face.fi, true);
+    return linearInterpolation(elem_eval, neighbor_eval, *face.fi, true);
 }
 
 /**
@@ -529,10 +535,16 @@ interpolate(const Limiter & limiter,
  * information neighbor. This pair should sum to unity. The second pair corresponds to the face
  * information functor element value and neighbor
  */
-template <typename T, typename Enable = typename std::enable_if<ScalarTraits<T>::value>::type>
+template <typename T,
+          FunctorEvaluationKind FEK = FunctorEvaluationKind::Value,
+          typename Enable = typename std::enable_if<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)
 {
+  static_assert((FEK == FunctorEvaluationKind::Value) || (FEK == FunctorEvaluationKind::Dot),
+                "Only Value and Dot evaluations currently supported");
+
+  constexpr FunctorEvaluationKind GFEK = FunctorGradientEvaluationKind<FEK>::value;
   typedef typename FunctorBase<T>::GradientType GradientType;
   static const GradientType zero(0);
 
@@ -541,8 +553,8 @@ interpCoeffsAndAdvected(const FunctorBase<T> & functor, const FaceArg & face, co
 
   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);
+  auto phi_upwind = functor.template genericEvaluate<FEK>(upwind_arg, time);
+  auto phi_downwind = functor.template genericEvaluate<FEK>(downwind_arg, time);
 
   std::pair<T, T> interp_coeffs;
   if (face.limiter_type == LimiterType::Upwind ||
@@ -551,7 +563,7 @@ interpCoeffsAndAdvected(const FunctorBase<T> & functor, const FaceArg & face, co
         interpCoeffs(*limiter, phi_upwind, phi_downwind, &zero, *face.fi, face.elem_is_upwind);
   else
   {
-    const auto grad_phi_upwind = functor.gradient(upwind_arg, time);
+    const auto grad_phi_upwind = functor.template genericEvaluate<GFEK>(upwind_arg, time);
     interp_coeffs = interpCoeffs(
         *limiter, phi_upwind, phi_downwind, &grad_phi_upwind, *face.fi, face.elem_is_upwind);
   }
@@ -565,16 +577,21 @@ interpCoeffsAndAdvected(const FunctorBase<T> & functor, const FaceArg & face, co
         std::make_pair(std::move(phi_downwind), std::move(phi_upwind)));
 }
 
-template <typename T, typename Enable = typename std::enable_if<ScalarTraits<T>::value>::type>
+template <typename T,
+          FunctorEvaluationKind FEK = FunctorEvaluationKind::Value,
+          typename Enable = typename std::enable_if<ScalarTraits<T>::value>::type>
 T
 interpolate(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");
+
   // 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(functor, face, time);
+    return linearInterpolation<T, FEK>(functor, face, time);
 
-  const auto [interp_coeffs, advected] = interpCoeffsAndAdvected(functor, face, time);
+  const auto [interp_coeffs, advected] = interpCoeffsAndAdvected<T, FEK>(functor, face, time);
   return interp_coeffs.first * advected.first + interp_coeffs.second * advected.second;
 }
 

framework/include/utils/PiecewiseByBlockLambdaFunctor.h

@@ -62,7 +62,6 @@ class PiecewiseByBlockLambdaFunctor : public Moose::FunctorBase<T>
   using typename Moose::FunctorBase<T>::ValueType;
   using typename Moose::FunctorBase<T>::DotType;
   using typename Moose::FunctorBase<T>::GradientType;
-  using typename Moose::FunctorBase<T>::FunctorReturnType;
 
 protected:
   using ElemFn = std::function<T(const Moose::ElemArg &, const Moose::StateArg &)>;

framework/include/variables/MooseVariableFE.h

@@ -81,7 +81,6 @@ class MooseVariableFE : public MooseVariableField<OutputType>
   using DoFValue = typename MooseVariableField<OutputType>::DoFValue;
 
   using FunctorArg = typename Moose::ADType<OutputType>::type;
-  using typename Moose::FunctorBase<FunctorArg>::FunctorReturnType;
   using typename Moose::FunctorBase<FunctorArg>::ValueType;
   using typename Moose::FunctorBase<FunctorArg>::GradientType;
   using typename Moose::FunctorBase<FunctorArg>::DotType;

framework/include/variables/MooseVariableFV.h

@@ -548,6 +548,8 @@ class MooseVariableFV : public MooseVariableField<OutputType>
   GradientType evaluateGradient(const ElemArg & elem_arg, const StateArg &) const override final;
   GradientType evaluateGradient(const FaceArg & face, const StateArg &) const override final;
   DotType evaluateDot(const ElemArg & elem, const StateArg &) const override final;
+  DotType evaluateDot(const FaceArg & face, const StateArg &) const override final;
+  DotType evaluateDot(const ElemQpArg & elem_qp, const StateArg &) const override final;
 
   /**
    * Setup the boundary to Dirichlet BC map

framework/src/variables/MooseVariableFV.C

@@ -851,6 +851,39 @@ MooseVariableFV<Real>::evaluateDot(const ElemArg & elem_arg, const StateArg & st
     return (*_sys.solutionUDot())(dof_index);
 }
 
+template <>
+ADReal
+MooseVariableFV<Real>::evaluateDot(const FaceArg & face, const StateArg & state) const
+{
+  const FaceInfo * const fi = face.fi;
+  mooseAssert(fi, "The face information must be non-null");
+  if (isDirichletBoundaryFace(*fi, face.face_side, state))
+    return ADReal(0.0); // No time derivative if boundary value is set
+  else if (isExtrapolatedBoundaryFace(*fi, face.face_side, state))
+  {
+    mooseAssert(face.face_side && this->hasBlocks(face.face_side->subdomain_id()),
+                "If we are an extrapolated boundary face, then our FunctorBase::checkFace method "
+                "should have assigned a non-null element that we are defined on");
+    const auto elem_arg = ElemArg({face.face_side, face.correct_skewness});
+    // For extrapolated boundary faces, note that we take the value of the time derivative at the
+    // cell in contact with the face
+    return evaluateDot(elem_arg, state);
+  }
+  else
+  {
+    mooseAssert(this->isInternalFace(*fi),
+                "We must be either Dirichlet, extrapolated, or internal");
+    return Moose::FV::interpolate<ADReal, FunctorEvaluationKind::Dot>(*this, face, state);
+  }
+}
+
+template <>
+ADReal
+MooseVariableFV<Real>::evaluateDot(const ElemQpArg & elem_qp, const StateArg & state) const
+{
+  return evaluateDot(ElemArg({elem_qp.elem, /*correct_skewness*/ false}), state);
+}
+
 template <typename OutputType>
 void
 MooseVariableFV<OutputType>::prepareAux()