Use SingleSidedFaceArg in NS advection kernels

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

Other notes:
- 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

Refs joe61vette/ASP#8

framework/include/base/MooseFunctor.h

@@ -205,6 +205,25 @@ struct SingleSidedFaceArg
   /// 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 on 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, apply_gradient_to_skewness};
+  }
+
   /**
    * Make a \p ElemArg from our data using the face information element
    */

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 evluations 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/src/fvbcs/FVBoundaryCondition.C

@@ -92,14 +92,7 @@ FVBoundaryCondition::singleSidedFaceArg(const FaceInfo * fi,
   if (!fi)
     fi = _face_info;
   const bool use_elem = fi->faceType(_var.name()) == FaceInfo::VarFaceNeighbors::ELEM;
+  const auto sub_id = use_elem ? fi->elem().subdomain_id() : fi->neighborPtr()->subdomain_id();
 
-  if (use_elem)
-    return {fi, limiter_type, true, correct_skewness, correct_skewness, fi->elem().subdomain_id()};
-  else
-    return {fi,
-            limiter_type,
-            true,
-            correct_skewness,
-            correct_skewness,
-            fi->neighborPtr()->subdomain_id()};
+  return {fi, limiter_type, true, correct_skewness, correct_skewness, sub_id};
 }

framework/src/fvkernels/FVFluxKernel.C

@@ -338,3 +338,16 @@ FVFluxKernel::faceArgSubdomains(const FaceInfo * face_info) const
 
   return Moose::FV::faceArgSubdomains(*this, *face_info);
 }
+
+Moose::SingleSidedFaceArg
+FVFluxKernel::singleSidedFaceArg(const FaceInfo * fi,
+                                 const Moose::FV::LimiterType limiter_type,
+                                 const bool correct_skewness) const
+{
+  if (!fi)
+    fi = _face_info;
+  const bool use_elem = fi->faceType(_var.name()) == FaceInfo::VarFaceNeighbors::ELEM;
+  const auto sub_id = use_elem ? fi->elem().subdomain_id() : fi->neighborPtr()->subdomain_id();
+
+  return {fi, limiter_type, true, correct_skewness, correct_skewness, sub_id};
+}

modules/navier_stokes/include/fvkernels/INSFVAdvectionKernel.h

@@ -23,6 +23,18 @@ class INSFVAdvectionKernel : public FVFluxKernel, public INSFVBCInterface
   INSFVAdvectionKernel(const InputParameters & params);
   void initialSetup() override;
 
+  /**
+   * Sets the advection and velocity interpolation methods
+   * @param obj The \p MooseObject with input parameters to query
+   * @param advected_interp_method The advected interpolation method we will set
+   * @param velocity_interp_method The velocity interpolation method we will set
+   * @return Whether the interpolation methods have indicated that we will need more than our base
+   * level of ghosting
+   */
+  static bool setInterpolationMethods(const MooseObject & obj,
+                                      Moose::FV::InterpMethod & advected_interp_method,
+                                      Moose::FV::InterpMethod & velocity_interp_method);
+
 protected:
   bool skipForBoundary(const FaceInfo & fi) const override;
 

modules/navier_stokes/src/fvkernels/INSFVAdvectionKernel.C

@@ -36,6 +36,51 @@ INSFVAdvectionKernel::validParams()
   return params;
 }
 
+bool
+INSFVAdvectionKernel::setInterpolationMethods(const MooseObject & obj,
+                                              Moose::FV::InterpMethod & advected_interp_method,
+                                              Moose::FV::InterpMethod & velocity_interp_method)
+{
+  using namespace Moose::FV;
+
+  bool need_more_ghosting = false;
+
+  const auto & advected_interp_method_in = obj.getParam<MooseEnum>("advected_interp_method");
+  if (advected_interp_method_in == "average")
+    advected_interp_method = InterpMethod::Average;
+  else if (advected_interp_method_in == "skewness-corrected")
+    advected_interp_method = Moose::FV::InterpMethod::SkewCorrectedAverage;
+  else if (advected_interp_method_in == "upwind")
+    advected_interp_method = InterpMethod::Upwind;
+  else
+  {
+    if (advected_interp_method_in == "sou")
+      advected_interp_method = InterpMethod::SOU;
+    else if (advected_interp_method_in == "min_mod")
+      advected_interp_method = InterpMethod::MinMod;
+    else if (advected_interp_method_in == "vanLeer")
+      advected_interp_method = InterpMethod::VanLeer;
+    else if (advected_interp_method_in == "quick")
+      advected_interp_method = InterpMethod::QUICK;
+    else
+      obj.mooseError("Unrecognized interpolation type ",
+                     static_cast<std::string>(advected_interp_method_in));
+
+    need_more_ghosting = true;
+  }
+
+  const auto & velocity_interp_method_in = obj.getParam<MooseEnum>("velocity_interp_method");
+  if (velocity_interp_method_in == "average")
+    velocity_interp_method = InterpMethod::Average;
+  else if (velocity_interp_method_in == "rc")
+    velocity_interp_method = InterpMethod::RhieChow;
+  else
+    obj.mooseError("Unrecognized interpolation type ",
+                   static_cast<std::string>(velocity_interp_method_in));
+
+  return need_more_ghosting;
+}
+
 INSFVAdvectionKernel::INSFVAdvectionKernel(const InputParameters & params)
   : FVFluxKernel(params),
     _rc_vel_provider(getUserObject<INSFVRhieChowInterpolator>("rhie_chow_user_object"))
@@ -45,66 +90,33 @@ INSFVAdvectionKernel::INSFVAdvectionKernel(const InputParameters & params)
              "configure script in the root MOOSE directory with the configure option "
              "'--with-ad-indexing-type=global'");
 #endif
-  using namespace Moose::FV;
-
-  const auto & advected_interp_method = getParam<MooseEnum>("advected_interp_method");
-  if (advected_interp_method == "average")
-    _advected_interp_method = InterpMethod::Average;
-  else if (advected_interp_method == "skewness-corrected")
-    _advected_interp_method = Moose::FV::InterpMethod::SkewCorrectedAverage;
-  else if (advected_interp_method == "upwind")
-    _advected_interp_method = InterpMethod::Upwind;
-  else
+  const bool need_more_ghosting =
+      setInterpolationMethods(*this, _advected_interp_method, _velocity_interp_method);
+  if (need_more_ghosting && _tid == 0)
   {
-    if (advected_interp_method == "sou")
-      _advected_interp_method = InterpMethod::SOU;
-    else if (advected_interp_method == "min_mod")
-      _advected_interp_method = InterpMethod::MinMod;
-    else if (advected_interp_method == "vanLeer")
-      _advected_interp_method = InterpMethod::VanLeer;
-    else if (advected_interp_method == "quick")
-      _advected_interp_method = InterpMethod::QUICK;
-    else
-      mooseError("Unrecognized interpolation type ",
-                 static_cast<std::string>(advected_interp_method));
-
-    if (_tid == 0)
-    {
-      auto & factory = _app.getFactory();
-
-      auto rm_params = factory.getValidParams("ElementSideNeighborLayers");
-
-      rm_params.set<std::string>("for_whom") = name();
-      rm_params.set<MooseMesh *>("mesh") = &const_cast<MooseMesh &>(_mesh);
-      rm_params.set<Moose::RelationshipManagerType>("rm_type") =
-          Moose::RelationshipManagerType::GEOMETRIC | Moose::RelationshipManagerType::ALGEBRAIC |
-          Moose::RelationshipManagerType::COUPLING;
-      FVKernel::setRMParams(
-          _pars,
-          rm_params,
-          std::max((unsigned short)(3), _pars.get<unsigned short>("ghost_layers")));
-      mooseAssert(rm_params.areAllRequiredParamsValid(),
-                  "All relationship manager parameters should be valid.");
-
-      auto rm_obj = factory.create<RelationshipManager>(
-          "ElementSideNeighborLayers", name() + "_skew_correction", rm_params);
-
-      // Delete the resources created on behalf of the RM if it ends up not being added to the
-      // App.
-      if (!_app.addRelationshipManager(rm_obj))
-        factory.releaseSharedObjects(*rm_obj);
-    }
+    auto & factory = _app.getFactory();
+
+    auto rm_params = factory.getValidParams("ElementSideNeighborLayers");
+
+    rm_params.set<std::string>("for_whom") = name();
+    rm_params.set<MooseMesh *>("mesh") = &const_cast<MooseMesh &>(_mesh);
+    rm_params.set<Moose::RelationshipManagerType>("rm_type") =
+        Moose::RelationshipManagerType::GEOMETRIC | Moose::RelationshipManagerType::ALGEBRAIC |
+        Moose::RelationshipManagerType::COUPLING;
+    FVKernel::setRMParams(
+        _pars, rm_params, std::max((unsigned short)(3), _pars.get<unsigned short>("ghost_layers")));
+    mooseAssert(rm_params.areAllRequiredParamsValid(),
+                "All relationship manager parameters should be valid.");
+
+    auto rm_obj = factory.create<RelationshipManager>(
+        "ElementSideNeighborLayers", name() + "_skew_correction", rm_params);
+
+    // Delete the resources created on behalf of the RM if it ends up not being added to the
+    // App.
+    if (!_app.addRelationshipManager(rm_obj))
+      factory.releaseSharedObjects(*rm_obj);
   }
 
-  const auto & velocity_interp_method = getParam<MooseEnum>("velocity_interp_method");
-  if (velocity_interp_method == "average")
-    _velocity_interp_method = InterpMethod::Average;
-  else if (velocity_interp_method == "rc")
-    _velocity_interp_method = InterpMethod::RhieChow;
-  else
-    mooseError("Unrecognized interpolation type ",
-               static_cast<std::string>(velocity_interp_method));
-
   auto param_check = [&params, this](const auto & param_name)
   {
     if (params.isParamSetByUser(param_name))

modules/navier_stokes/src/fvkernels/INSFVMassAdvection.C

@@ -35,11 +35,15 @@ ADReal
 INSFVMassAdvection::computeQpResidual()
 {
   const auto v = _rc_vel_provider.getVelocity(_velocity_interp_method, *_face_info, _tid);
-  const auto rho_interface =
-      Moose::FV::interpolate(_rho,
-                             Moose::FV::makeFace(*_face_info,
-                                                 limiterType(_advected_interp_method),
-                                                 MetaPhysicL::raw_value(v) * _normal > 0,
-                                                 faceArgSubdomains()));
-  return _normal * v * rho_interface;
+  ADReal rho_face{};
+
+  if (onBoundary(*_face_info))
+    rho_face = _rho(singleSidedFaceArg());
+  else
+    rho_face = Moose::FV::interpolate(_rho,
+                                      Moose::FV::makeFace(*_face_info,
+                                                          limiterType(_advected_interp_method),
+                                                          MetaPhysicL::raw_value(v) * _normal > 0,
+                                                          faceArgSubdomains()));
+  return _normal * v * rho_face;
 }

modules/navier_stokes/src/fvkernels/INSFVMomentumAdvection.C

@@ -11,6 +11,7 @@
 #include "NS.h"
 #include "FVUtils.h"
 #include "INSFVRhieChowInterpolator.h"
+#include "SystemBase.h"
 
 registerMooseObject("NavierStokesApp", INSFVMomentumAdvection);
 
@@ -43,9 +44,8 @@ INSFVMomentumAdvection::initialSetup()
 
   _rho_u = std::make_unique<PiecewiseByBlockLambdaFunctor<ADReal>>(
       "rho_u",
-      [this](const auto & r, const auto & t) -> ADReal {
-        return _rho(r, t) * _var(r, t) / epsilon()(r, t);
-      },
+      [this](const auto & r, const auto & t) -> ADReal
+      { return _rho(r, t) * _var(r, t) / epsilon()(r, t); },
       std::set<ExecFlagType>({EXEC_ALWAYS}),
       _mesh,
       this->blockIDs());
@@ -59,28 +59,48 @@ INSFVMomentumAdvection::computeQpResidual()
   const bool correct_skewness = _advected_interp_method == InterpMethod::SkewCorrectedAverage;
   const auto v = _rc_vel_provider.getVelocity(_velocity_interp_method, *_face_info, _tid);
 
-  const auto [interp_coeffs, advected] =
-      interpCoeffsAndAdvected(*_rho_u,
-                              makeFace(*_face_info,
-                                       limiterType(_advected_interp_method),
-                                       MetaPhysicL::raw_value(v) * _normal > 0,
-                                       faceArgSubdomains(),
-                                       correct_skewness,
-                                       correct_skewness));
-
-  const auto elem_face = elemFromFace();
-  const auto neighbor_face = neighborFromFace();
-
-  const auto rho_elem = _rho(elem_face), rho_neighbor = _rho(neighbor_face);
-  const auto eps_elem = epsilon()(elem_face), eps_neighbor = epsilon()(neighbor_face);
-  const auto var_elem = advected.first / rho_elem * eps_elem,
-             var_neighbor = advected.second / rho_neighbor * eps_neighbor;
-
-  _ae = _normal * v * rho_elem / eps_elem * interp_coeffs.first;
-  // Minus sign because we apply a minus sign to the residual in computeResidual
-  _an = -_normal * v * rho_neighbor / eps_neighbor * interp_coeffs.second;
-
-  return _ae * var_elem - _an * var_neighbor;
+  if (onBoundary(*_face_info))
+  {
+    const auto ssf = singleSidedFaceArg();
+    const Elem * const sided_elem = ssf.makeSidedElem().elem;
+    const auto dof_number = sided_elem->dof_number(_sys.number(), _var.number(), 0);
+    const auto rhof = _rho(ssf);
+    const auto epsf = epsilon()(ssf);
+    const auto uf = _var(ssf);
+    const Real duf_du = uf.derivatives()[dof_number];
+    const auto coeff = _normal * v * rhof / epsf;
+    if (sided_elem == &_face_info->elem())
+      _ae = coeff * duf_du;
+    else
+      _an = -coeff * duf_du;
+
+    return coeff * uf;
+  }
+  else
+  {
+    const auto [interp_coeffs, advected] =
+        interpCoeffsAndAdvected(*_rho_u,
+                                makeFace(*_face_info,
+                                         limiterType(_advected_interp_method),
+                                         MetaPhysicL::raw_value(v) * _normal > 0,
+                                         faceArgSubdomains(),
+                                         correct_skewness,
+                                         correct_skewness));
+
+    const auto elem_face = elemFromFace();
+    const auto neighbor_face = neighborFromFace();
+
+    const auto rho_elem = _rho(elem_face), rho_neighbor = _rho(neighbor_face);
+    const auto eps_elem = epsilon()(elem_face), eps_neighbor = epsilon()(neighbor_face);
+    const auto var_elem = advected.first / rho_elem * eps_elem,
+               var_neighbor = advected.second / rho_neighbor * eps_neighbor;
+
+    _ae = _normal * v * rho_elem / eps_elem * interp_coeffs.first;
+    // Minus sign because we apply a minus sign to the residual in computeResidual
+    _an = -_normal * v * rho_neighbor / eps_neighbor * interp_coeffs.second;
+
+    return _ae * var_elem - _an * var_neighbor;
+  }
 }
 
 void