Use SingleSidedFaceArg in NS advection kernels and always use skew we…

…ights

This forces atomic evaluation of components of aggregate functors (such as
functor material properties), which allows direct substitution of Dirichlet
information when component(s) are variables. This allows more accurate setting
of incoming fluxes

This dramatically improves our results on skewed meshes, as demonstrated
on one of Sinan's cask setups

This prevents zero mixing lengths and zero eddy viscosity.

Additionally: change default for fv = true variable two-term boundary
expansion to `true`. If users select finite volume variables through
the `fv = true` flag as opposed to specifying `type = MooseVariableFVReal`,
then we need to make sure they get the same default for two term boundary
expansion, which is `true`.

Except when doing vertex based interpolation. This necessitated adding
a special handler in `Moose::FV::interpolate(FunctorBase, FaceArg)` as
the `Moose::FV::linearInterpolation(FunctorBase, FaceArg)` function is
the only global function that currently handles skewness correction
correctly.

I decided to make this change for MooseVariableFV::evaluate(FaceArg)
after trigging a debug mode assertion when interpolating passive scalars
in the 2d-mixing-length test in INSFV. When I had changed the scalar
field advection kernel to functor evaluation based off Face arguments,
then I actually changed the evaluation behind-the-scenes from upwind
for this test to central-differencing. This is because all other face
evaluations in INSFV (I think) go through functor material properties
and end up using the global `Moose::FV::interpolate` functions which apply
the correct interpolation/limiting method supplies by the user. However,
prior to this commit variable evaluations did their own handling without
dispatching to the global `Moose::FV::interpolate` functions and always
applied some form of central-differencing on internal faces. That is
now changed so the mixing length test is back to upwinding its passive scalar
transport, which necessitated another re-golding

- We make sure to extend the change in advection kernel evaluation at boundaries
  to the VolumetricFlowRate postprocessor.
- Evaluation of incoming momentum flux is more accurate with these changes which
  puts more tax on the wall-function tests. This is because before when we upwinded
  the momentum at the boundaries we would actually apply positive y-velocities
  at the incoming face near the wall which led to a non-monotone y-velocity
  profile near the wall moving from the inlet down the channel. Now we
  actually apply the 0 y-velocity at the inlet face which makes for a more
  challenging solve due to the large gradient in y-velocity the solver wants
  to impose due to the large viscous forces. (This is my hypothesis anyway). To
  make things easier on the solver, I've reduced the Reynolds number. Note that
  if you psuedo-step to steady-state, convergence is excellent with the original
  Reynolds number
- Created a couple of global functions for help in creating interpolation parameters
  and for using those parameters to set advection and velocity interpolation object
  data members
- Limiter/interpolation coverage enhancements
- An additional change in this commit is making PINSFVEnergyAdvection
  a trivial derivative of INSFVEnergyAdvection such that we reuse all
  the code, e.g. computeQpResidual
- Test vector and array functor material properties
- Allow construction of VectorCompositeFunctor with 1-3 functors
- Convert mixing length tests into transient tests

Refs joe61vette/ASP#8

framework/include/base/MooseFunctor.h

@@ -38,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);
   }
 };
 
@@ -64,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
@@ -81,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);
   }
 };
 
@@ -107,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
@@ -130,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)
@@ -162,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);
   }
@@ -194,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);
   }
 };
 
@@ -268,6 +249,13 @@ class FunctorBase
   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
+  /// 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;

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

modules/navier_stokes/include/utils/ArrayComponentFunctor.h → framework/include/utils/ArrayComponentFunctor.h

@@ -24,7 +24,9 @@ class ArrayComponentFunctor : public Moose::FunctorBase<T>
   using typename Moose::FunctorBase<T>::DotType;
 
   ArrayComponentFunctor(const ArrayTypeFunctor & array, const unsigned int component)
-    : _array(array), _component(component)
+    : Moose::FunctorBase<T>(array.functorName() + "_" + std::to_string(component)),
+      _array(array),
+      _component(component)
   {
   }
 
@@ -73,4 +75,11 @@ class ArrayComponentFunctor : public Moose::FunctorBase<T>
   {
     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];
+  }
 };