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
  */
@@ -530,6 +541,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

@@ -230,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 typename libMesh::TensorTools::IncrementRank<T>::type example_gradient(0);
-
-  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>

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/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;
 }

modules/navier_stokes/src/fvkernels/INSFVMomentumAdvection.C

@@ -27,7 +27,13 @@ INSFVMomentumAdvection::validParams()
 INSFVMomentumAdvection::INSFVMomentumAdvection(const InputParameters & params)
   : INSFVAdvectionKernel(params),
     INSFVMomentumResidualObject(*this),
-    _rho(getFunctor<ADReal>(NS::density))
+    _rho(getFunctor<ADReal>(NS::density)),
+    _rho_u(
+        "rho_u",
+        [this](const auto & r, const auto & t) -> ADReal { return _rho(r, t) * _var(r, t); },
+        std::set<ExecFlagType>({EXEC_ALWAYS}),
+        _mesh,
+        this->blockIDs())
 {
 #ifndef MOOSE_GLOBAL_AD_INDEXING
   mooseError("INSFV is not supported by local AD indexing. In order to use INSFV, please run the "
@@ -39,16 +45,27 @@ INSFVMomentumAdvection::INSFVMomentumAdvection(const InputParameters & params)
 ADReal
 INSFVMomentumAdvection::computeQpResidual()
 {
+  using namespace Moose::FV;
+
+  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()));
+
   const auto elem_face = elemFromFace();
   const auto neighbor_face = neighborFromFace();
 
-  const auto v = _rc_vel_provider.getVelocity(_velocity_interp_method, *_face_info, _tid);
-  const auto interp_coeffs = Moose::FV::interpCoeffs(_advected_interp_method, *_face_info, true, v);
-  _ae = _normal * v * _rho(elem_face) * interp_coeffs.first;
+  const auto rho_elem = _rho(elem_face), rho_neighbor = _rho(neighbor_face);
+  const auto var_elem = advected.first / rho_elem, var_neighbor = advected.second / rho_neighbor;
+
+  _ae = _normal * v * rho_elem * interp_coeffs.first;
   // Minus sign because we apply a minus sign to the residual in computeResidual
-  _an = -_normal * v * _rho(neighbor_face) * interp_coeffs.second;
+  _an = -_normal * v * rho_neighbor * interp_coeffs.second;
 
-  return _ae * _var(elem_face) - _an * _var(neighbor_face);
+  return _ae * var_elem - _an * var_neighbor;
 }
 
 void