Merge pull request #20504 from lindsayad/general-advection-schemes-20493

Implement generalized advection schemes for RC NS implementation

framework/include/base/MooseFunctor.h

@@ -22,6 +22,8 @@
 #include "libmesh/remote_elem.h"
 #include "libmesh/tensor_tools.h"
 
+#include "metaphysicl/ct_types.h"
+
 #include <unordered_map>
 #include <functional>
 
@@ -36,12 +38,11 @@ struct ElemArg
 {
   const libMesh::Elem * elem;
   bool correct_skewness;
-  bool apply_gradient_to_skewness;
 
   friend bool operator<(const ElemArg & l, const ElemArg & r)
   {
-    return std::make_tuple(l.elem, l.correct_skewness, l.apply_gradient_to_skewness) <
-           std::make_tuple(r.elem, r.correct_skewness, r.apply_gradient_to_skewness);
+    return std::make_tuple(l.elem, l.correct_skewness) <
+           std::make_tuple(r.elem, r.correct_skewness);
   }
 };
 
@@ -62,14 +63,9 @@ struct ElemFromFaceArg
   /// information object will be used to help construct a ghost value evaluation
   const FaceInfo * fi;
 
-  /// Whether to apply skew correction weights
+  /// Whether to perform skew correction
   bool correct_skewness;
 
-  /// Whether to apply the face gradient when computing a skew corrected face value. A true value
-  /// for this data member in conjunction with a false value for \p correct_skewness does not make
-  /// sense
-  bool apply_gradient_to_skewness;
-
   /// a subdomain ID. This is useful when the functor is a material property and the user wants
   /// to indicate which material property definition should be used to evaluate the functor. For
   /// instance if we are using a flux kernel that is not defined on one side of the face, the
@@ -79,14 +75,12 @@ struct ElemFromFaceArg
   /**
    * Make a \p ElemArg from our data
    */
-  ElemArg makeElem() const { return {elem, correct_skewness, apply_gradient_to_skewness}; }
+  ElemArg makeElem() const { return {elem, correct_skewness}; }
 
   friend bool operator<(const ElemFromFaceArg & l, const ElemFromFaceArg & r)
   {
-    return std::make_tuple(
-               l.elem, l.fi, l.correct_skewness, l.apply_gradient_to_skewness, l.sub_id) <
-           std::make_tuple(
-               r.elem, r.fi, r.correct_skewness, r.apply_gradient_to_skewness, r.sub_id);
+    return std::make_tuple(l.elem, l.fi, l.correct_skewness, l.sub_id) <
+           std::make_tuple(r.elem, r.fi, r.correct_skewness, r.sub_id);
   }
 };
 
@@ -105,14 +99,9 @@ struct FaceArg
   /// a boolean which states whether the face information element is upwind of the face
   bool elem_is_upwind;
 
-  /// Whether to apply skew correction weights
+  /// Whether to perform skew correction
   bool correct_skewness;
 
-  /// Whether to apply the face gradient when computing a skew corrected face value. A true value
-  /// for this data member in conjunction with a false value for \p correct_skewness does not make
-  /// sense
-  bool apply_gradient_to_skewness;
-
   ///@{
   /**
    * a pair of subdomain IDs. These do not always correspond to the face info element subdomain
@@ -128,30 +117,24 @@ struct FaceArg
   /**
    * Make a \p ElemArg from our data using the face information element
    */
-  ElemArg makeElem() const { return {&fi->elem(), correct_skewness, apply_gradient_to_skewness}; }
+  ElemArg makeElem() const { return {&fi->elem(), correct_skewness}; }
 
   /**
    * Make a \p ElemArg from our data using the face information neighbor
    */
-  ElemArg makeNeighbor() const
-  {
-    return {fi->neighborPtr(), correct_skewness, apply_gradient_to_skewness};
-  }
+  ElemArg makeNeighbor() const { return {fi->neighborPtr(), correct_skewness}; }
 
   /**
    * Make a \p ElemFromFaceArg from our data using the face information element
    */
-  ElemFromFaceArg elemFromFace() const
-  {
-    return {&fi->elem(), fi, correct_skewness, apply_gradient_to_skewness, elem_sub_id};
-  }
+  ElemFromFaceArg elemFromFace() const { return {&fi->elem(), fi, correct_skewness, elem_sub_id}; }
 
   /**
    * Make a \p ElemFromFaceArg from our data using the face information neighbor
    */
   ElemFromFaceArg neighborFromFace() const
   {
-    return {fi->neighborPtr(), fi, correct_skewness, apply_gradient_to_skewness, neighbor_sub_id};
+    return {fi->neighborPtr(), fi, correct_skewness, neighbor_sub_id};
   }
 
   friend bool operator<(const FaceArg & l, const FaceArg & r)
@@ -160,13 +143,11 @@ struct FaceArg
                            l.limiter_type,
                            l.elem_is_upwind,
                            l.correct_skewness,
-                           l.apply_gradient_to_skewness,
                            l.elem_sub_id,
                            l.neighbor_sub_id) < std::make_tuple(r.fi,
                                                                 r.limiter_type,
                                                                 r.elem_is_upwind,
                                                                 r.correct_skewness,
-                                                                r.apply_gradient_to_skewness,
                                                                 r.elem_sub_id,
                                                                 l.neighbor_sub_id);
   }
@@ -192,43 +173,45 @@ struct SingleSidedFaceArg
   /// a boolean which states whether the face information element is upwind of the face
   bool elem_is_upwind;
 
-  /// Whether to apply skew correction weights
+  /// Whether to perform skew correction
   bool correct_skewness;
 
-  /// Whether to apply the face gradient when computing a skew corrected face value. A true value
-  /// for this data member in conjunction with a false value for \p correct_skewness does not make
-  /// sense
-  bool apply_gradient_to_skewness;
-
   /// The subdomain ID which denotes the side of the face information we are evaluating on
   SubdomainID sub_id;
 
+  /**
+   * Make an element argument corresponding to the subdomain ID on which we are defined
+   */
+  ElemArg makeSidedElem() const
+  {
+    mooseAssert(fi, "The face must be non-null");
+#ifndef NDEBUG
+    const unsigned short matches = (sub_id == fi->elem().subdomain_id()) +
+                                   (fi->neighborPtr() && (sub_id == fi->neighbor().subdomain_id()));
+    mooseAssert(matches == 1,
+                "We should only have one match. If we have more or less, then this is not an "
+                "appropriate use of SingleSidedFaceArg");
+#endif
+
+    const Elem * const ret_elem =
+        sub_id == fi->elem().subdomain_id() ? &fi->elem() : fi->neighborPtr();
+    return {ret_elem, correct_skewness};
+  }
+
   /**
    * Make a \p ElemArg from our data using the face information element
    */
-  ElemArg makeElem() const { return {&fi->elem(), correct_skewness, apply_gradient_to_skewness}; }
+  ElemArg makeElem() const { return {&fi->elem(), correct_skewness}; }
 
   /**
    * Make a \p ElemArg from our data using the face information neighbor
    */
-  ElemArg makeNeighbor() const
-  {
-    return {fi->neighborPtr(), correct_skewness, apply_gradient_to_skewness};
-  }
+  ElemArg makeNeighbor() const { return {fi->neighborPtr(), correct_skewness}; }
 
   friend bool operator<(const SingleSidedFaceArg & l, const SingleSidedFaceArg & r)
   {
-    return std::make_tuple(l.fi,
-                           l.limiter_type,
-                           l.elem_is_upwind,
-                           l.correct_skewness,
-                           l.apply_gradient_to_skewness,
-                           l.sub_id) < std::make_tuple(r.fi,
-                                                       r.limiter_type,
-                                                       r.elem_is_upwind,
-                                                       r.correct_skewness,
-                                                       r.apply_gradient_to_skewness,
-                                                       r.sub_id);
+    return std::make_tuple(l.fi, l.limiter_type, l.elem_is_upwind, l.correct_skewness, l.sub_id) <
+           std::make_tuple(r.fi, r.limiter_type, r.elem_is_upwind, r.correct_skewness, r.sub_id);
   }
 };
 
@@ -266,7 +249,16 @@ class FunctorBase
   using FunctorType = FunctorBase<T>;
   using FunctorReturnType = T;
   using ValueType = T;
-  using GradientType = typename libMesh::TensorTools::IncrementRank<T>::type;
+  /// 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
+  /// container of algebraic types, for example Reals, then the GradientType will be a container of
+  /// the gradients of those algebraic types, in this example VectorValue<Reals>. So if T is
+  /// 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 DotType = ValueType;
 
   virtual ~FunctorBase() = default;
@@ -328,9 +320,15 @@ class FunctorBase
   void setCacheClearanceSchedule(const std::set<ExecFlagType> & clearance_schedule);
 
   /**
-   * Returns whether this face is an extrapolated boundary face for this functor
+   * Returns a pair where the first member is whether this face is an extrapolated boundary face for
+   * this functor. The second member is the element on which this functor is defined if this is an
+   * extrapolated boundary face (if it is not an extrapolated boundary face, then we just return the
+   * face information \p &elem())
    */
-  virtual bool isExtrapolatedBoundaryFace(const FaceInfo &) const { mooseError("not implemented"); }
+  virtual std::pair<bool, const Elem *> isExtrapolatedBoundaryFace(const FaceInfo &) const
+  {
+    mooseError("not implemented");
+  }
 
   /**
    * Returns true if this functor is a constant

framework/include/fvkernels/FVFluxKernel.h

@@ -110,6 +110,20 @@ class FVFluxKernel : public FVKernel,
    */
   std::pair<SubdomainID, SubdomainID> faceArgSubdomains(const FaceInfo * face_info = nullptr) const;
 
+  /**
+   * Determine the single sided face argument when evaluating a functor on a face.
+   * This is used to perform evaluations of material properties with the actual face values of
+   * their dependences, rather than interpolate the material property to the boundary.
+   * @param fi the FaceInfo for this face
+   * @param limiter_type the limiter type, to be specified if more than the default average
+   *        interpolation is required for the parameters of the functor
+   * @param correct_skewness whether to perform skew correction at the face
+   */
+  Moose::SingleSidedFaceArg singleSidedFaceArg(
+      const FaceInfo * fi = nullptr,
+      Moose::FV::LimiterType limiter_type = Moose::FV::LimiterType::CentralDifference,
+      bool correct_skewness = false) const;
+
   const bool _force_boundary_execution;
 
   std::unordered_set<BoundaryID> _boundaries_to_force;

framework/include/mesh/FaceInfo.h

@@ -165,9 +165,6 @@ class FaceInfo
   /// Return the geometric weighting factor
   Real gC() const { return _gc; }
 
-  /// Return the weighting factor for skewed element-pairs
-  Real gCSkewed() const { return _gc_skewed; }
-
   /**
    * @return the distance vector drawn from centroid C to F, or in terms of MOOSE implementation,
    * the distance vector obtained from subtracting the element centroid from the neighbor centroid
@@ -241,9 +238,6 @@ class FaceInfo
   const Point _neighbor_centroid;
   const Real _neighbor_volume;
 
-  /// Geometric weighting factor for face value interpolation
-  const Real _gc;
-
   /// the distance vector between neighbor and element centroids
   const RealVectorValue _d_cf;
 
@@ -256,9 +250,8 @@ class FaceInfo
   /// The vector to the intersection of d_{CF} and the face.
   Point _r_intersection;
 
-  /// Geometric weighting factor for face value interpolation in case of skewed
-  /// cell-connections
-  Real _gc_skewed;
+  /// Geometric weighting factor for face value interpolation
+  Real _gc;
 
   /// cached locations of variables in solution vectors
   /// TODO: make this more efficient by not using a map if possible

framework/include/utils/ArrayComponentFunctor.h

@@ -0,0 +1,85 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* 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 "MooseFunctor.h"
+
+/**
+ * This is essentially a forwarding functor that forwards the spatial and temporal evaluation
+ * arguments to the parent array functor and then returns the result indexed at a given component.
+ */
+template <typename T, typename ArrayTypeFunctor>
+class ArrayComponentFunctor : public Moose::FunctorBase<T>
+{
+public:
+  using typename Moose::FunctorBase<T>::ValueType;
+  using typename Moose::FunctorBase<T>::GradientType;
+  using typename Moose::FunctorBase<T>::DotType;
+
+  ArrayComponentFunctor(const ArrayTypeFunctor & array, const unsigned int component)
+    : Moose::FunctorBase<T>(array.functorName() + "_" + std::to_string(component)),
+      _array(array),
+      _component(component)
+  {
+  }
+
+  std::pair<bool, const Elem *> isExtrapolatedBoundaryFace(const FaceInfo & fi) const override
+  {
+    return _array.isExtrapolatedBoundaryFace(fi);
+  }
+
+private:
+  /// The parent array functor
+  const ArrayTypeFunctor & _array;
+
+  /// The component at which we'll index the parent array functor evaluation result
+  const unsigned int _component;
+
+  ValueType evaluate(const Moose::ElemArg & elem, const unsigned int state) const override final
+  {
+    return _array(elem, state)[_component];
+  }
+
+  ValueType evaluate(const Moose::ElemFromFaceArg & elem_from_face,
+                     const unsigned int state) const override final
+  {
+    return _array(elem_from_face, state)[_component];
+  }
+
+  ValueType evaluate(const Moose::FaceArg & face, const unsigned int state) const override final
+  {
+    return _array(face, state)[_component];
+  }
+
+  ValueType evaluate(const Moose::SingleSidedFaceArg & face,
+                     const unsigned int state) const override final
+  {
+    return _array(face, state)[_component];
+  }
+
+  ValueType evaluate(const Moose::ElemQpArg & elem_qp,
+                     const unsigned int state) const override final
+  {
+    return _array(elem_qp, state)[_component];
+  }
+
+  ValueType evaluate(const Moose::ElemSideQpArg & elem_side_qp,
+                     const unsigned int state) const override final
+  {
+    return _array(elem_side_qp, state)[_component];
+  }
+
+  using Moose::FunctorBase<T>::evaluateGradient;
+  GradientType evaluateGradient(const Moose::ElemArg & elem,
+                                const unsigned int state) const override final
+  {
+    return _array.gradient(elem, state)[_component];
+  }
+};