Maximize code-reuse for mass and momentum between INSFV and PINSFV

- Have PINSFVMassAdvection inherit from INSFVMassAdvection
- Add epsilon virtual method to INSFVMomentumAdvection

modules/navier_stokes/include/fvkernels/INSFVMomentumAdvection.h

@@ -24,19 +24,28 @@ class INSFVMomentumAdvection : public INSFVAdvectionKernel, public INSFVMomentum
   INSFVMomentumAdvection(const InputParameters & params);
   void gatherRCData(const Elem &) override final {}
   void gatherRCData(const FaceInfo & fi) override final;
+  void initialSetup() override;
 
 protected:
   virtual ADReal computeQpResidual() override;
 
+  /**
+   * A virtual method that allows us to reuse all the code from free-flow for porous
+   */
+  virtual const Moose::FunctorBase<ADReal> & epsilon() const { return _unity_functor; }
+
   /// Density
   const Moose::Functor<ADReal> & _rho;
 
   /// Our local momentum functor
-  const PiecewiseByBlockLambdaFunctor<ADReal> _rho_u;
+  std::unique_ptr<PiecewiseByBlockLambdaFunctor<ADReal>> _rho_u;
 
   /// The a coefficient for the element
   ADReal _ae = 0;
 
   /// The a coefficient for the neighbor
   ADReal _an = 0;
+
+  /// A unity functor used in the epsilon virtual method
+  const Moose::ConstantFunctor<ADReal> _unity_functor{1};
 };

modules/navier_stokes/include/fvkernels/PINSFVMassAdvection.h

@@ -9,20 +9,14 @@
 
 #pragma once
 
-#include "INSFVAdvectionKernel.h"
+#include "INSFVMassAdvection.h"
 
 /**
  * A flux kernel transporting mass in porous media across cell faces
  */
-class PINSFVMassAdvection : public INSFVAdvectionKernel
+class PINSFVMassAdvection : public INSFVMassAdvection
 {
 public:
   static InputParameters validParams();
   PINSFVMassAdvection(const InputParameters & params);
-
-protected:
-  ADReal computeQpResidual() override;
-
-  /// Density
-  const Moose::Functor<ADReal> & _rho;
 };

modules/navier_stokes/include/fvkernels/PINSFVMomentumAdvection.h

@@ -21,7 +21,7 @@ class PINSFVMomentumAdvection : public INSFVMomentumAdvection
   PINSFVMomentumAdvection(const InputParameters & params);
 
 protected:
-  virtual ADReal computeQpResidual() override;
+  const Moose::FunctorBase<ADReal> & epsilon() const override { return _eps; }
 
   /// porosity functor
   const Moose::Functor<ADReal> & _eps;

modules/navier_stokes/src/fvkernels/INSFVMomentumAdvection.C

@@ -27,13 +27,7 @@ INSFVMomentumAdvection::validParams()
 INSFVMomentumAdvection::INSFVMomentumAdvection(const InputParameters & params)
   : INSFVAdvectionKernel(params),
     INSFVMomentumResidualObject(*this),
-    _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())
+    _rho(getFunctor<ADReal>(NS::density))
 {
 #ifndef MOOSE_GLOBAL_AD_INDEXING
   mooseError("INSFV is not supported by local AD indexing. In order to use INSFV, please run the "
@@ -42,15 +36,31 @@ INSFVMomentumAdvection::INSFVMomentumAdvection(const InputParameters & params)
 #endif
 }
 
+void
+INSFVMomentumAdvection::initialSetup()
+{
+  INSFVAdvectionKernel::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);
+      },
+      std::set<ExecFlagType>({EXEC_ALWAYS}),
+      _mesh,
+      this->blockIDs());
+}
+
 ADReal
 INSFVMomentumAdvection::computeQpResidual()
 {
   using namespace Moose::FV;
 
   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,
+      interpCoeffsAndAdvected(*_rho_u,
                               makeFace(*_face_info,
                                        limiterType(_advected_interp_method),
                                        MetaPhysicL::raw_value(v) * _normal > 0,
@@ -62,11 +72,13 @@ INSFVMomentumAdvection::computeQpResidual()
   const auto neighbor_face = neighborFromFace();
 
   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;
+  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 * interp_coeffs.first;
+  _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 * interp_coeffs.second;
+  _an = -_normal * v * rho_neighbor / eps_neighbor * interp_coeffs.second;
 
   return _ae * var_elem - _an * var_neighbor;
 }

modules/navier_stokes/src/fvkernels/PINSFVMassAdvection.C

@@ -16,37 +16,17 @@ registerMooseObject("NavierStokesApp", PINSFVMassAdvection);
 InputParameters
 PINSFVMassAdvection::validParams()
 {
-  auto params = INSFVAdvectionKernel::validParams();
-  params.addRequiredParam<MooseFunctorName>(NS::density, "Density functor");
+  auto params = INSFVMassAdvection::validParams();
   params.addClassDescription("Object for advecting mass in porous media mass equation");
   return params;
 }
 
 PINSFVMassAdvection::PINSFVMassAdvection(const InputParameters & params)
-  : INSFVAdvectionKernel(params), _rho(getFunctor<ADReal>(NS::density))
+  : INSFVMassAdvection(params)
 {
 #ifndef MOOSE_GLOBAL_AD_INDEXING
   mooseError("PINSFV is not supported by local AD indexing. In order to use PINSFV, please run the "
              "configure script in the root MOOSE directory with the configure option "
              "'--with-ad-indexing-type=global'");
 #endif
 }
-
-ADReal
-PINSFVMassAdvection::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);
-  return _normal * v * rho_interface;
-}

modules/navier_stokes/src/fvkernels/PINSFVMomentumAdvection.C

@@ -29,21 +29,3 @@ PINSFVMomentumAdvection::PINSFVMomentumAdvection(const InputParameters & params)
   : INSFVMomentumAdvection(params), _eps(getFunctor<ADReal>(NS::porosity))
 {
 }
-
-ADReal
-PINSFVMomentumAdvection::computeQpResidual()
-{
-  const auto elem_face = elemFromFace();
-  const auto neighbor_face = neighborFromFace();
-
-  // Superficial velocity interpolation
-  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 / _eps(elem_face);
-  // Minus sign because we apply a minus sign to the residual in computeResidual
-  _an = -_normal * v * _rho(neighbor_face) * interp_coeffs.second / _eps(neighbor_face);
-
-  return _ae * _var(elem_face) - _an * _var(neighbor_face);
-}