Adapt (P)INSFV to functor props

It's perhaps worth mentioning by now that the current functor
implementation forces a paradigm in which 'interpolated' face
quantities are formed by composing atomic interpolated face quantities,
e.g. interp(u*v) becomes interp(u)*interp(v). I imagine we will
have some debate about whether this is satisfactory. My current
data point is that this change does not impact the order of
approximation error has demonstrated by continued passing of ins
and pins mms tests

framework/include/fvbcs/FVMatAdvectionOutflowBC.h

@@ -24,8 +24,4 @@ class FVMatAdvectionOutflowBC : public FVFluxBC
 
   /// The advected quantity on the elem
   const FunctorInterface<ADReal> & _adv_quant;
-
-  /// The interfacial velocity. We cache this in the residual computation in case a derived class
-  /// might want to use it
-  ADRealVectorValue _v;
 };

framework/include/fvkernels/FVMatAdvection.h

@@ -30,8 +30,4 @@ class FVMatAdvection : public FVFluxKernel
 
   /// The limiting method to use for the advected quantity
   std::unique_ptr<Moose::FV::Limiter> _limiter;
-
-  /// The interfacial velocity. We cache this in the residual computation in case a derived class
-  /// might want to use it
-  ADRealVectorValue _v;
 };

framework/src/fvbcs/FVMatAdvectionOutflowBC.C

@@ -47,8 +47,8 @@ FVMatAdvectionOutflowBC::FVMatAdvectionOutflowBC(const InputParameters & params)
 ADReal
 FVMatAdvectionOutflowBC::computeQpResidual()
 {
-  _v = _vel(std::make_tuple(_face_info, nullptr, true));
+  const auto v = _vel(std::make_tuple(_face_info, nullptr, true));
   const auto adv_quant_boundary =
-      _adv_quant(std::make_tuple(_face_info, nullptr, _v * _face_info->normal() > 0));
-  return _normal * _v * adv_quant_boundary;
+      _adv_quant(std::make_tuple(_face_info, nullptr, v * _face_info->normal() > 0));
+  return _normal * v * adv_quant_boundary;
 }

framework/src/fvkernels/FVMatAdvection.C

@@ -58,9 +58,10 @@ FVMatAdvection::computeQpResidual()
 {
   using namespace Moose::FV;
 
-  _v = _vel(std::make_tuple(_face_info, _cd_limiter.get(), /*this doesn't matter for cd*/ true));
+  const auto v =
+      _vel(std::make_tuple(_face_info, _cd_limiter.get(), /*this doesn't matter for cd*/ true));
   const auto adv_quant_interface =
-      _adv_quant(std::make_tuple(_face_info, _limiter.get(), _v * _face_info->normal() > 0));
+      _adv_quant(std::make_tuple(_face_info, _limiter.get(), v * _face_info->normal() > 0));
 
-  return _normal * _v * adv_quant_interface;
+  return _normal * v * adv_quant_interface;
 }

modules/navier_stokes/include/fvkernels/INSFVMomentumAdvection.h

@@ -36,10 +36,7 @@ class INSFVMomentumAdvection : public FVMatAdvection
   /**
    * interpolation overload for the velocity
    */
-  virtual void interpolate(Moose::FV::InterpMethod m,
-                           ADRealVectorValue & interp_v,
-                           const ADRealVectorValue & elem_v,
-                           const ADRealVectorValue & neighbor_v);
+  virtual void interpolate(Moose::FV::InterpMethod m, ADRealVectorValue & interp_v);
 
   virtual ADReal computeQpResidual() override;
 

modules/navier_stokes/include/fvkernels/PINSFVMomentumAdvection.h

@@ -25,20 +25,15 @@ class PINSFVMomentumAdvection : public INSFVMomentumAdvection
   /**
    * interpolation overload for the velocity
    */
-  void interpolate(Moose::FV::InterpMethod m,
-                   ADRealVectorValue & interp_v,
-                   const ADRealVectorValue & elem_v,
-                   const ADRealVectorValue & neighbor_v) override;
+  void interpolate(Moose::FV::InterpMethod m, ADRealVectorValue & interp_v) override;
 
   virtual ADReal computeQpResidual() override;
   VectorValue<ADReal> coeffCalculator(const Elem & elem) const override;
 
   /// porosity variable to compute gradients
   const MooseVariableFV<Real> * const _eps_var;
   /// porosity in the current element
-  const VariableValue & _eps;
-  /// porosity in the neighbor element
-  const VariableValue & _eps_neighbor;
+  const FunctorInterface<ADReal> & _eps;
   /// Whether the porosity field is smooth or has discontinuities
   const bool _smooth_porosity;
 };

modules/navier_stokes/include/fvkernels/PINSFVMomentumDiffusion.h

@@ -10,6 +10,7 @@
 #pragma once
 
 #include "FVFluxKernel.h"
+#include "CentralDifferenceLimiter.h"
 
 /**
  * A flux kernel for diffusion of momentum in porous media across cell faces
@@ -24,26 +25,24 @@ class PINSFVMomentumDiffusion : public FVFluxKernel
   ADReal computeQpResidual() override;
 
   /// the current element viscosity
-  const ADMaterialProperty<Real> & _mu_elem;
-  /// the neighbor element viscosity
-  const ADMaterialProperty<Real> & _mu_neighbor;
+  const FunctorMaterialProperty<ADReal> & _mu;
 
   /// the porosity
-  const VariableValue & _eps;
-  /// the neighbor element porosity
-  const VariableValue & _eps_neighbor;
+  const FunctorInterface<ADReal> & _eps;
 
   // Parameters for the gradient diffusion term
   /// Which momentum component this kernel applies to
   const int _index;
 
   /// Velocity as material properties
-  const ADMaterialProperty<RealVectorValue> * _vel_elem;
-  const ADMaterialProperty<RealVectorValue> * _vel_neighbor;
+  const FunctorMaterialProperty<ADRealVectorValue> * const _vel;
 
   /// the porosity as a variable to be able to compute a face gradient
   const MooseVariableFVReal * const _eps_var;
 
   /// Whether to add the porosity gradient term, only for continuous porosity
   const bool _smooth_porosity;
+
+  /// The object used to perform average/central-difference interpolations
+  CentralDifferenceLimiter _cd_limiter;
 };

modules/navier_stokes/include/materials/INSFVMaterial.h

@@ -19,41 +19,39 @@ class INSFVMaterial : public Material
   INSFVMaterial(const InputParameters & parameters);
 
 protected:
-  void computeQpProperties() override;
-
   /// x-component velocity
-  const ADVariableValue & _u_vel;
+  const MooseVariableFVReal & _u_vel;
 
   /// y-component velocity
-  const ADVariableValue & _v_vel;
+  const MooseVariableFVReal * const _v_vel;
 
   /// z-component velocity
-  const ADVariableValue & _w_vel;
+  const MooseVariableFVReal * const _w_vel;
 
   /// pressure variable
-  const ADVariableValue & _p_var;
+  const MooseVariableFVReal & _p_var;
 
   /// The velocity as a vector
-  ADMaterialProperty<RealVectorValue> & _velocity;
+  FunctorMaterialProperty<ADRealVectorValue> & _velocity;
 
   /// The density times the x-velocity
-  ADMaterialProperty<Real> & _rho_u;
+  FunctorMaterialProperty<ADReal> & _rho_u;
 
   /// The density times the y-velocity
-  ADMaterialProperty<Real> & _rho_v;
+  FunctorMaterialProperty<ADReal> & _rho_v;
 
   /// The density times the z-velocity
-  ADMaterialProperty<Real> & _rho_w;
+  FunctorMaterialProperty<ADReal> & _rho_w;
 
   /// The pressure material property
-  ADMaterialProperty<Real> & _p;
+  FunctorMaterialProperty<ADReal> & _p;
 
   /// density
-  const Real & _rho;
+  const FunctorMaterialProperty<ADReal> & _rho;
 
   const bool _has_temperature;
 
-  const ADVariableValue * const _temperature;
-  const ADMaterialProperty<Real> * const _cp;
-  ADMaterialProperty<Real> * const _rho_cp_temp;
+  const MooseVariableFVReal * const _temperature;
+  const FunctorMaterialProperty<ADReal> * const _cp;
+  FunctorMaterialProperty<ADReal> * const _rho_cp_temp;
 };

modules/navier_stokes/src/fvbcs/INSFVMomentumAdvectionOutflowBC.C

@@ -67,14 +67,9 @@ INSFVMomentumAdvectionOutflowBC::computeQpResidual()
   if (_w_var)
     v(2) = _w_var->getBoundaryFaceValue(*_face_info);
 
-  ADReal adv_quant_boundary;
-  interpolate(_advected_interp_method,
-              adv_quant_boundary,
-              _adv_quant_elem[_qp],
-              _adv_quant_neighbor[_qp],
-              v,
-              *_face_info,
-              true);
+  const auto adv_quant_boundary =
+      _adv_quant(std::make_tuple(_face_info, nullptr, v * _face_info->normal() > 0));
+
   mooseAssert(_normal * v >= 0,
               "This boundary condition is for outflow but the flow is in the opposite direction of "
               "the boundary normal");

modules/navier_stokes/src/fvbcs/PINSFVMomentumAdvectionOutflowBC.C

@@ -70,14 +70,9 @@ PINSFVMomentumAdvectionOutflowBC::computeQpResidual()
   if (_w_var)
     v(2) = _w_var->getBoundaryFaceValue(*_face_info);
 
-  ADReal adv_quant_boundary;
-  interpolate(_advected_interp_method,
-              adv_quant_boundary,
-              _adv_quant_elem[_qp] / _eps[_qp],
-              _adv_quant_neighbor[_qp] / _eps_neighbor[_qp],
-              v,
-              *_face_info,
-              true);
+  const auto adv_quant_boundary =
+      _adv_quant(std::make_tuple(_face_info, nullptr, v * _face_info->normal() > 0));
+
   mooseAssert(_normal * v >= 0,
               "This boundary condition is for outflow but the flow is in the opposite direction of "
               "the boundary normal");

modules/navier_stokes/src/fvkernels/INSFVMassAdvection.C

@@ -34,27 +34,9 @@ ADReal
 INSFVMassAdvection::computeQpResidual()
 {
   ADRealVectorValue v;
-  ADReal face_rho;
-  Point normal;
 
-  this->interpolate(_velocity_interp_method, v, _vel_elem[_qp], _vel_neighbor[_qp]);
-  if (onBoundary(*_face_info))
-  {
-    const auto ft = _face_info->faceType(_var.name());
-    const bool out_of_elem = (ft == FaceInfo::VarFaceNeighbors::ELEM);
-    normal = out_of_elem ? _face_info->normal() : Point(-_face_info->normal());
-    face_rho = out_of_elem ? _rho(&_face_info->elem()) : _rho(_face_info->neighborPtr());
-  }
-  else
-  {
-    normal = _face_info->normal();
-    Moose::FV::interpolate(_advected_interp_method,
-                           face_rho,
-                           _rho(&_face_info->elem()),
-                           _rho(_face_info->neighborPtr()),
-                           v,
-                           *_face_info,
-                           true);
-  }
-  return normal * v * face_rho;
+  this->interpolate(_velocity_interp_method, v);
+  const auto rho_interface =
+      _rho(std::make_tuple(_face_info, _limiter.get(), v * _face_info->normal() > 0));
+  return _normal * v * rho_interface;
 }

modules/navier_stokes/src/fvkernels/INSFVMomentumAdvection.C

@@ -301,6 +301,8 @@ INSFVMomentumAdvection::coeffCalculator(const Elem & elem) const
                             const bool elem_has_info) {
     mooseAssert(fi, "We need a non-null FaceInfo");
     mooseAssert(&elem == &functor_elem, "Elems don't match");
+    mooseAssert((&elem == &fi->elem()) || (&elem == fi->neighborPtr()),
+                "Surely the element has to match one of the face information's elements right?");
 
     const Point normal = elem_has_info ? fi->normal() : Point(-fi->normal());
     const Point & rc_centroid = elem_has_info ? fi->elemCentroid() : fi->neighborCentroid();
@@ -314,8 +316,8 @@ INSFVMomentumAdvection::coeffCalculator(const Elem & elem) const
           "Let's make sure our normal is what we think it is");
 #endif
 
-    const auto elem_mu = _mu(&elem);
-    const auto elem_rho = _rho(&elem);
+    const auto face_mu = _mu(std::make_tuple(fi, _cd_limiter.get(), true));
+    const auto face_rho = _rho(std::make_tuple(fi, _cd_limiter.get(), true));
 
     // Unless specified otherwise, "elem" here refers to the element we're computing the
     // Rhie-Chow coefficient for. "neighbor" is the element across the current FaceInfo (fi)
@@ -331,7 +333,7 @@ INSFVMomentumAdvection::coeffCalculator(const Elem & elem) const
         {
           // Need to account for viscous shear stress from wall
           for (const auto i : make_range(_dim))
-            coeff(i) += elem_mu * surface_vector.norm() /
+            coeff(i) += face_mu * surface_vector.norm() /
                         std::abs((fi->faceCentroid() - rc_centroid) * normal) *
                         (1 - normal(i) * normal(i));
 
@@ -352,9 +354,7 @@ INSFVMomentumAdvection::coeffCalculator(const Elem & elem) const
           if (_w_var)
             face_velocity(2) = _w_var->getBoundaryFaceValue(*fi);
 
-          const auto advection_coeffs =
-              Moose::FV::interpCoeffs(_advected_interp_method, *fi, elem_has_info, face_velocity);
-          ADReal temp_coeff = elem_rho * face_velocity * surface_vector * advection_coeffs.first;
+          ADReal temp_coeff = face_rho * face_velocity * surface_vector;
 
           if (_fully_developed_flow_boundaries.find(bc_id) ==
               _fully_developed_flow_boundaries.end())
@@ -366,7 +366,7 @@ INSFVMomentumAdvection::coeffCalculator(const Elem & elem) const
             // equation 8.89. So relative to the internal face viscous term, we have substituted
             // eqn. 8.82 for 8.78
             temp_coeff +=
-                elem_mu * surface_vector.norm() / (fi->faceCentroid() - rc_centroid).norm();
+                face_mu * surface_vector.norm() / (fi->faceCentroid() - rc_centroid).norm();
 
           // For flow boundaries, the coefficient addition is the same for every velocity component
           for (const auto i : make_range(_dim))
@@ -379,7 +379,7 @@ INSFVMomentumAdvection::coeffCalculator(const Elem & elem) const
         {
           // Moukalled eqns. 15.154 - 15.156
           for (const auto i : make_range(_dim))
-            coeff(i) += 2. * elem_mu * surface_vector.norm() /
+            coeff(i) += 2. * face_mu * surface_vector.norm() /
                         std::abs((fi->faceCentroid() - rc_centroid) * normal) * normal(i) *
                         normal(i);
 
@@ -396,14 +396,8 @@ INSFVMomentumAdvection::coeffCalculator(const Elem & elem) const
 
     // Else we are on an internal face
 
-    const auto neighbor_mu = _mu(neighbor);
-    ADReal face_mu;
-    Moose::FV::interpolate(
-        Moose::FV::InterpMethod::Average, face_mu, elem_mu, neighbor_mu, *fi, elem_has_info);
-    const auto neighbor_rho = _rho(neighbor);
-    ADReal face_rho;
-    Moose::FV::interpolate(
-        Moose::FV::InterpMethod::Average, face_rho, elem_rho, neighbor_rho, *fi, elem_has_info);
+    mooseAssert((neighbor == &fi->elem()) || (neighbor == fi->neighborPtr()),
+                "Surely the neighbor has to match one of the face information's elements, right?");
 
     ADRealVectorValue neighbor_velocity(_u_var->getNeighborValue(neighbor, *fi, elem_velocity(0)));
     if (_v_var)
@@ -419,11 +413,19 @@ INSFVMomentumAdvection::coeffCalculator(const Elem & elem) const
                            *fi,
                            elem_has_info);
 
-    // we are only interested in the interpolation coefficient for the Rhie-Chow element,
-    // so we just use the 'first' member of the returned pair
+    // At the point that we support more advection interpolation strategies than upwind and average
+    // (central-differencing), we will need to pass phi_downwind, phi_upwind, and grad_phi_upwind
+    // arguments to interpCoeffs. However, because upwind and average interpolations are independent
+    // of the values of those arguments we can just pass in a dummy value for the phi values and a
+    // nullptr for the gradient
+    static const ADReal dummy = 0;
+    const bool fi_elem_is_upwind = interp_v * fi->normal() > 0;
     const auto advection_coeffs =
-        Moose::FV::interpCoeffs(_advected_interp_method, *fi, elem_has_info, interp_v);
-    ADReal temp_coeff = face_rho * interp_v * surface_vector * advection_coeffs.first;
+        Moose::FV::interpCoeffs(*_limiter, dummy, dummy, nullptr, *fi, fi_elem_is_upwind);
+    const bool functor_elem_is_upwind = (fi_elem_is_upwind == (&elem == &fi->elem()));
+    const auto & advection_coeff =
+        functor_elem_is_upwind ? advection_coeffs.first : advection_coeffs.second;
+    ADReal temp_coeff = face_rho * interp_v * surface_vector * advection_coeff;
 
     // Now add the viscous flux. Note that this includes only the orthogonal component! See
     // Moukalled equations 8.80, 8.78, and the orthogonal correction approach equation for
@@ -442,10 +444,7 @@ INSFVMomentumAdvection::coeffCalculator(const Elem & elem) const
 }
 
 void
-INSFVMomentumAdvection::interpolate(Moose::FV::InterpMethod m,
-                                    ADRealVectorValue & v,
-                                    const ADRealVectorValue & elem_v,
-                                    const ADRealVectorValue & neighbor_v)
+INSFVMomentumAdvection::interpolate(Moose::FV::InterpMethod m, ADRealVectorValue & v)
 {
   const Elem * const elem = &_face_info->elem();
   const Elem * const neighbor = _face_info->neighborPtr();
@@ -474,8 +473,7 @@ INSFVMomentumAdvection::interpolate(Moose::FV::InterpMethod m,
     return;
   }
 
-  Moose::FV::interpolate(
-      Moose::FV::InterpMethod::Average, v, elem_v, neighbor_v, *_face_info, true);
+  v = _vel(std::make_tuple(_face_info, _cd_limiter.get(), true));
 
   if (m == Moose::FV::InterpMethod::Average)
     return;
@@ -544,10 +542,7 @@ INSFVMomentumAdvection::coeffCalculator(const Elem &) const
 }
 
 void
-INSFVMomentumAdvection::interpolate(Moose::FV::InterpMethod,
-                                    ADRealVectorValue &,
-                                    const ADRealVectorValue &,
-                                    const ADRealVectorValue &)
+INSFVMomentumAdvection::interpolate(Moose::FV::InterpMethod, ADRealVectorValue &)
 {
   mooseError("INSFVMomentumAdvection only works with global AD indexing");
 }
@@ -557,16 +552,11 @@ ADReal
 INSFVMomentumAdvection::computeQpResidual()
 {
   ADRealVectorValue v;
-  ADReal adv_quant_interface;
-
-  this->interpolate(_velocity_interp_method, v, _vel_elem[_qp], _vel_neighbor[_qp]);
-  Moose::FV::interpolate(_advected_interp_method,
-                         adv_quant_interface,
-                         _adv_quant_elem[_qp],
-                         _adv_quant_neighbor[_qp],
-                         v,
-                         *_face_info,
-                         true);
+
+  this->interpolate(_velocity_interp_method, v);
+  const auto adv_quant_interface =
+      _adv_quant(std::make_tuple(_face_info, _limiter.get(), v * _face_info->normal() > 0));
+
   return _normal * v * adv_quant_interface;
 }