Implement generalized advection schemes for RC NS implementation

Closes idaholab#20493

framework/include/utils/MathFVUtils.h

@@ -90,6 +90,23 @@ faceArgSubdomains(const SubdomainRestrictable & obj, const FaceInfo & fi)
                         makeSidedFace(obj, fi.neighborPtr(), fi).sub_id);
 }
 
+inline FaceArg
+makeFace(const FaceInfo & fi,
+         const LimiterType limiter_type,
+         const bool elem_is_upwind,
+         const std::pair<SubdomainID, SubdomainID> & subs,
+         const bool skewness_correction = false,
+         const bool apply_gradient_correction = false)
+{
+  return {&fi,
+          limiter_type,
+          elem_is_upwind,
+          skewness_correction,
+          apply_gradient_correction,
+          subs.first,
+          subs.second};
+}
+
 /**
  * Make a functor face argument with a central differencing limiter, e.g. compose a face argument
  * that will tell functors to perform (possibly skew-corrected) linear interpolations from cell
@@ -111,13 +128,12 @@ makeCDFace(const FaceInfo & fi,
            const bool skewness_correction = false,
            const bool apply_gradient_correction = false)
 {
-  return {&fi,
-          Moose::FV::LimiterType::CentralDifference,
-          true,
-          skewness_correction,
-          apply_gradient_correction,
-          subs.first,
-          subs.second};
+  return makeFace(fi,
+                  LimiterType::CentralDifference,
+                  true,
+                  subs,
+                  skewness_correction,
+                  apply_gradient_correction);
 }
 
 /**
@@ -482,11 +498,6 @@ interpCoeffs(const Limiter<T> & limiter,
   return std::make_pair(1. - g, g);
 }
 
-template <typename T>
-struct GradientType
-{
-};
-
 /**
  * Interpolates with a limiter
  */
@@ -598,6 +609,105 @@ interpolate(const Limiter & limiter,
   return ret;
 }
 
+/**
+ * Interpolates with a limiter and a face argument
+ * @return a pair of pairs. The first pair corresponds to the interpolation coefficients with the
+ * first corresponding to the face information element and the second corresponding to the face
+ * 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>
+std::pair<std::pair<T, T>, std::pair<T, T>>
+interpCoeffsAndAdvected(const FunctorBase<T> & functor, const FaceArg & face)
+{
+  typedef typename FunctorBase<T>::GradientType GradientType;
+  static const GradientType zero(0);
+
+  mooseAssert(face.fi, "this must be non-null");
+  const auto limiter = Limiter<typename LimiterValueType<T>::value_type>::build(face.limiter_type);
+
+  const auto upwind_arg = face.elem_is_upwind ? face.elemFromFace() : face.neighborFromFace();
+  const auto downwind_arg = face.elem_is_upwind ? face.neighborFromFace() : face.elemFromFace();
+  auto phi_upwind = functor(upwind_arg);
+  auto phi_downwind = functor(downwind_arg);
+
+  std::pair<T, T> interp_coeffs;
+  if (face.limiter_type == LimiterType::Upwind ||
+      face.limiter_type == LimiterType::CentralDifference)
+    interp_coeffs =
+        interpCoeffs(*limiter, phi_upwind, phi_downwind, &zero, *face.fi, face.elem_is_upwind);
+  else
+  {
+    const auto grad_phi_upwind = functor.gradient(upwind_arg);
+    interp_coeffs = interpCoeffs(
+        *limiter, phi_upwind, phi_downwind, &grad_phi_upwind, *face.fi, face.elem_is_upwind);
+  }
+
+  if (face.elem_is_upwind)
+    return std::make_pair(std::move(interp_coeffs),
+                          std::make_pair(std::move(phi_upwind), std::move(phi_downwind)));
+  else
+    return std::make_pair(
+        std::make_pair(std::move(interp_coeffs.second), std::move(interp_coeffs.first)),
+        std::make_pair(std::move(phi_downwind), std::move(phi_upwind)));
+}
+
+template <typename T, typename Enable = typename std::enable_if<ScalarTraits<T>::value>::type>
+T
+interpolate(const FunctorBase<T> & functor, const FaceArg & face)
+{
+  const auto [interp_coeffs, advected] = interpCoeffsAndAdvected(functor, face);
+  return interp_coeffs.first * advected.first + interp_coeffs.second * advected.second;
+}
+
+template <typename T>
+VectorValue<T>
+interpolate(const FunctorBase<VectorValue<T>> & functor, const FaceArg & face)
+{
+  static const VectorValue<T> grad_zero(0);
+
+  mooseAssert(face.fi, "this must be non-null");
+  const auto limiter = Limiter<typename LimiterValueType<T>::value_type>::build(face.limiter_type);
+
+  const auto upwind_arg = face.elem_is_upwind ? face.elemFromFace() : face.neighborFromFace();
+  const auto downwind_arg = face.elem_is_upwind ? face.neighborFromFace() : face.elemFromFace();
+  auto phi_upwind = functor(upwind_arg);
+  auto phi_downwind = functor(downwind_arg);
+
+  VectorValue<T> ret;
+  T coeff_upwind, coeff_downwind;
+
+  if (face.limiter_type == LimiterType::Upwind ||
+      face.limiter_type == LimiterType::CentralDifference)
+  {
+    for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
+    {
+      const auto &component_upwind = phi_upwind(i), component_downwind = phi_downwind(i);
+      std::tie(coeff_upwind, coeff_downwind) = interpCoeffs(*limiter,
+                                                            component_upwind,
+                                                            component_downwind,
+                                                            &grad_zero,
+                                                            *face.fi,
+                                                            face.elem_is_upwind);
+      ret(i) = coeff_upwind * component_upwind + coeff_downwind * component_downwind;
+    }
+  }
+  else
+  {
+    const auto grad_phi_upwind = functor.gradient(upwind_arg);
+    for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
+    {
+      const auto &component_upwind = phi_upwind(i), component_downwind = phi_downwind(i);
+      const VectorValue<T> grad = grad_phi_upwind.row(i);
+      std::tie(coeff_upwind, coeff_downwind) = interpCoeffs(
+          *limiter, component_upwind, component_downwind, &grad, *face.fi, face.elem_is_upwind);
+      ret(i) = coeff_upwind * component_upwind + coeff_downwind * component_downwind;
+    }
+  }
+
+  return ret;
+}
+
 /**
  * Return whether the supplied face is on a boundary of the \p object's execution
  */

framework/include/utils/PiecewiseByBlockLambdaFunctor.h

@@ -78,6 +78,8 @@ class PiecewiseByBlockLambdaFunctor : public Moose::FunctorBase<T>
   using Moose::FunctorBase<T>::evaluateGradient;
   GradientType evaluateGradient(const Moose::ElemArg & elem_arg, unsigned int) const override;
   GradientType evaluateGradient(const Moose::FaceArg & face_arg, unsigned int) const override;
+  GradientType evaluateGradient(const Moose::ElemFromFaceArg & elem_from_face_arg,
+                                unsigned int) const override;
 
 private:
   /**
@@ -228,36 +230,12 @@ PiecewiseByBlockLambdaFunctor<T>::evaluate(const Moose::SingleSidedFaceArg & fac
 
 template <typename T>
 typename PiecewiseByBlockLambdaFunctor<T>::ValueType
-PiecewiseByBlockLambdaFunctor<T>::evaluate(const Moose::FaceArg & face, unsigned int state) const
+PiecewiseByBlockLambdaFunctor<T>::evaluate(const Moose::FaceArg & face,
+                                           unsigned int libmesh_dbg_var(state)) const
 {
   using namespace Moose::FV;
-
-  mooseAssert(face.fi,
-              "We must have a non-null face_info in order to prepare our ElemFromFace tuples");
-  static const GradientType example_gradient;
-
-  switch (face.limiter_type)
-  {
-    case LimiterType::CentralDifference:
-    case LimiterType::Upwind:
-    {
-      const auto elem_from_face = face.elemFromFace();
-      const auto neighbor_from_face = face.neighborFromFace();
-      const auto & upwind_elem = face.elem_is_upwind ? elem_from_face : neighbor_from_face;
-      const auto & downwind_elem = face.elem_is_upwind ? neighbor_from_face : elem_from_face;
-      auto limiter_obj =
-          Limiter<typename LimiterValueType<T>::value_type>::build(face.limiter_type);
-      return interpolate(*limiter_obj,
-                         evaluate(upwind_elem, state),
-                         evaluate(downwind_elem, state),
-                         &example_gradient,
-                         *face.fi,
-                         face.elem_is_upwind);
-    }
-
-    default:
-      mooseError("Unsupported limiter type in PiecewiseByBlockLambdaFunctor::evaluate(FaceArg)");
-  }
+  mooseAssert(state == 0, "Only current time state supported.");
+  return interpolate(*this, face);
 }
 
 template <typename T>
@@ -301,3 +279,15 @@ PiecewiseByBlockLambdaFunctor<T>::evaluateGradient(const Moose::FaceArg & face_a
 {
   return Moose::FV::greenGaussGradient(face_arg, *this, true, _mesh);
 }
+
+template <typename T>
+typename PiecewiseByBlockLambdaFunctor<T>::GradientType
+PiecewiseByBlockLambdaFunctor<T>::evaluateGradient(
+    const Moose::ElemFromFaceArg & elem_from_face_arg, unsigned int) const
+{
+  const auto elem_arg = elem_from_face_arg.makeElem();
+  if (elem_arg.elem)
+    return Moose::FV::greenGaussGradient(elem_arg, *this, true, _mesh);
+  else
+    return GradientType(0);
+}

framework/src/limiters/Limiter.C

@@ -66,6 +66,18 @@ limiterType(const InterpMethod interp_method)
     case InterpMethod::Upwind:
       return LimiterType::Upwind;
 
+    case InterpMethod::VanLeer:
+      return LimiterType::VanLeer;
+
+    case InterpMethod::MinMod:
+      return LimiterType::MinMod;
+
+    case InterpMethod::SOU:
+      return LimiterType::SOU;
+
+    case InterpMethod::QUICK:
+      return LimiterType::QUICK;
+
     default:
       mooseError("Unrecognized interpolation method type.");
   }

modules/navier_stokes/include/fvkernels/INSFVMomentumAdvection.h

@@ -11,6 +11,7 @@
 
 #include "INSFVAdvectionKernel.h"
 #include "INSFVMomentumResidualObject.h"
+#include "PiecewiseByBlockLambdaFunctor.h"
 
 /**
  * An advection kernel that implements interpolation schemes specific to Navier-Stokes flow
@@ -30,6 +31,9 @@ class INSFVMomentumAdvection : public INSFVAdvectionKernel, public INSFVMomentum
   /// Density
   const Moose::Functor<ADReal> & _rho;
 
+  /// Our local momentum functor
+  const PiecewiseByBlockLambdaFunctor<ADReal> _rho_u;
+
   /// The a coefficient for the element
   ADReal _ae = 0;
 

modules/navier_stokes/src/fvkernels/INSFVAdvectionKernel.C

@@ -10,12 +10,14 @@
 #include "INSFVAdvectionKernel.h"
 #include "NS.h"
 #include "MooseVariableFV.h"
+#include "RelationshipManager.h"
 
 InputParameters
 INSFVAdvectionKernel::validParams()
 {
   InputParameters params = FVFluxKernel::validParams();
-  MooseEnum advected_interp_method("average upwind skewness-corrected", "upwind");
+  MooseEnum advected_interp_method("average upwind sou min_mod vanLeer quick skewness-corrected",
+                                   "upwind");
   params.addParam<MooseEnum>(
       "advected_interp_method",
       advected_interp_method,
@@ -53,8 +55,46 @@ INSFVAdvectionKernel::INSFVAdvectionKernel(const InputParameters & params)
   else if (advected_interp_method == "upwind")
     _advected_interp_method = InterpMethod::Upwind;
   else
-    mooseError("Unrecognized interpolation type ",
-               static_cast<std::string>(advected_interp_method));
+  {
+    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);
+    }
+  }
 
   const auto & velocity_interp_method = getParam<MooseEnum>("velocity_interp_method");
   if (velocity_interp_method == "average")

modules/navier_stokes/src/fvkernels/INSFVEnergyAdvection.C

@@ -35,18 +35,12 @@ INSFVEnergyAdvection::INSFVEnergyAdvection(const InputParameters & params)
 ADReal
 INSFVEnergyAdvection::computeQpResidual()
 {
-  ADReal adv_quant_interface;
-
-  const auto elem_face = elemFromFace();
-  const auto neighbor_face = neighborFromFace();
-
   const auto v = _rc_vel_provider.getVelocity(_velocity_interp_method, *_face_info, _tid);
-  Moose::FV::interpolate(_advected_interp_method,
-                         adv_quant_interface,
-                         _adv_quant(elem_face),
-                         _adv_quant(neighbor_face),
-                         v,
-                         *_face_info,
-                         true);
+  const auto adv_quant_interface =
+      Moose::FV::interpolate(_adv_quant,
+                             Moose::FV::makeFace(*_face_info,
+                                                 limiterType(_advected_interp_method),
+                                                 MetaPhysicL::raw_value(v) * _normal > 0,
+                                                 faceArgSubdomains()));
   return _normal * v * adv_quant_interface;
 }

modules/navier_stokes/src/fvkernels/INSFVMassAdvection.C

@@ -34,18 +34,12 @@ INSFVMassAdvection::INSFVMassAdvection(const InputParameters & params)
 ADReal
 INSFVMassAdvection::computeQpResidual()
 {
-  ADReal rho_interface;
-
-  const auto elem_face = elemFromFace();
-  const auto neighbor_face = neighborFromFace();
-
   const auto v = _rc_vel_provider.getVelocity(_velocity_interp_method, *_face_info, _tid);
-  Moose::FV::interpolate(_advected_interp_method,
-                         rho_interface,
-                         _rho(elem_face),
-                         _rho(neighbor_face),
-                         v,
-                         *_face_info,
-                         true);
+  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;
 }