[WIP] Add porous medium treatment to linear FV segregated solvers.

modules/navier_stokes/doc/content/design/porous_linear_simple.md

@@ -0,0 +1,51 @@
+# Porous SIMPLE Formulation
+
+`PorousRhieChowMassFlux` and `LinearPNSFVMomentumPressure` extend the
+segregated linear finite-volume solver to the porous media equations
+
+$$
+\nabla \cdot (\epsilon \rho \mathbf{u} \mathbf{u})
+ - \nabla \cdot (\epsilon \mu \nabla \mathbf{u})
+ = - \epsilon \nabla p + \mathbf{f},
+\qquad
+\nabla \cdot (\epsilon \rho \mathbf{u}) = 0,
+$$
+
+where $\mathbf{u}$ is the interstitial velocity and $\epsilon$ is the porosity.
+The modifications relative to the non-porous SIMPLE implementation are:
+
+1.  The Rhie–Chow interpolation stores $\epsilon \rho$ in both the $H/A$ and
+    $A^{-1}$ functors, which makes the advective fluxes and the pressure
+    correction consistent with $\nabla \cdot (\epsilon \rho \mathbf{u}) = 0$.
+2.  The pressure gradient contribution in the momentum predictor is multiplied
+    by the local porosity.
+3.  Optional Bernoulli jumps and user-defined losses can be injected at block
+    interfaces to represent empirical head changes.
+
+Users should provide a viscosity functor that already contains the porosity
+factor (e.g. through an `ADParsedFunctorMaterial` or a parsed function)
+whenever the diffusion term is required.
+
+## Porosity Interpolation Strategy
+
+Face porosities are evaluated with harmonic weighting by default to maintain
+continuity of the superficial flux when the interstitial velocity is allowed
+to jump between blocks.  Arithmetic and geometric options are available
+through the `porosity_face_interpolation` parameter for cases where a different
+scheme is preferable.
+
+## Interface Pressure Jumps
+
+When a face is tagged (either automatically by detecting $\epsilon$ jumps
+between neighboring blocks or manually through `porosity_jump_sidesets`) the
+user object computes a Bernoulli correction
+
+$$ \Delta p = \frac{1}{2}\rho_d v_d^2 - \frac{1}{2}\rho_u v_u^2, $$
+
+where the subscripts denote upwind and downwind states along the face normal.
+Form factors supplied in `interface_form_factors` add irreversible losses based
+on the downwind superficial velocity for contractions and the upwind velocity
+for expansions, while `interface_pressure_jumps` allows the user to prescribe
+additional constant drops.  The resulting pressure jump is subtracted from the
+downwind pressure unknown before the face gradient is evaluated, which embeds
+the loss in both the predictor and corrector steps without special kernels.

modules/navier_stokes/doc/content/source/linearfvkernels/LinearPNSFVMomentumPressure.md

@@ -0,0 +1,29 @@
+!description /LinearFVKernels/LinearPNSFVMomentumPressure
+
+## Description
+
+`LinearPNSFVMomentumPressure` multiplies the pressure gradient source that is
+added to the segregated momentum predictor by the local porosity.  This kernel
+should be used together with `PorousRhieChowMassFlux` and the porous mass
+balance kernels when the unknown velocity is the interstitial velocity,
+yielding the porous momentum balance
+
+$$
+\nabla \cdot \left( \epsilon \mu \nabla \mathbf{u} \right) =
+ - \epsilon \nabla p + \cdots
+$$
+
+Only the right-hand side contribution is affected, so the kernel can be
+combined with the existing linear FV momentum flux kernels by supplying a
+porosity-weighted viscosity functor.
+
+## Example Input File
+
+!listing modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/porous_linear_basic.i
+  block=LinearFVKernels/momentum_pressure
+  caption=Porosity-weighted pressure gradient term.
+
+## Parameters
+
+!parameters /LinearFVKernels/LinearPNSFVMomentumPressure
+

modules/navier_stokes/doc/content/source/userobjects/PorousRhieChowMassFlux.md

@@ -0,0 +1,42 @@
+!description /UserObjects/PorousRhieChowMassFlux
+
+## Description
+
+`PorousRhieChowMassFlux` extends the linear finite-volume Rhie–Chow mass flux
+interpolator to porous-medium applications.  The user object computes face
+fluxes according to
+
+$$
+\dot{m}_f = -\left(\frac{H}{A}\right)_f + (A^{-1})_f \nabla p_f
+$$
+
+where the density terms in both $H/A$ and $A^{-1}$ include the local porosity,
+so that the advected mass flux is $\epsilon \rho \mathbf{u} \cdot \mathbf{n}$.
+It also exposes face-averaged porosity values in the `Ainv` functor that
+appear in the pressure diffusion equation.  Harmonic porosity weighting is
+used by default to maintain continuity of superficial fluxes when the
+interstitial velocity is discontinuous.
+
+In addition to the porosity weighting, `PorousRhieChowMassFlux` can inject
+manual or automatically-detected pressure jumps at block interfaces.  When
+porosity jumps across an internal face (or when the face belongs to one of
+`porosity_jump_sidesets`), a Bernoulli correction
+
+$$ \Delta p = \frac{1}{2}\rho_d v_d^2 - \frac{1}{2}\rho_u v_u^2 $$
+
+is applied to the downstream momentum equation.  Optional form factors and
+user-supplied pressure jumps may be added through the `interface_*`
+parameters to mimic irreversible losses and empirical head terms.  The manual
+faces are supplied by creating a sideset on the block interface (for example
+with `SideSetsBetweenSubdomainsGenerator`).
+
+## Example Input File
+
+!listing modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/porous_linear_basic.i
+  block=UserObjects
+  caption=Porous Rhie–Chow interpolator setup in a SIMPLE input.
+
+## Parameters
+
+!parameters /UserObjects/PorousRhieChowMassFlux
+

modules/navier_stokes/include/linearfvkernels/LinearPNSFVMomentumPressure.h

@@ -0,0 +1,31 @@
+//* This file is part of the MOOSE framework
+//* https://mooseframework.inl.gov
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "LinearFVMomentumPressure.h"
+
+#include "MooseFunctor.h"
+
+/**
+ * Pressure gradient kernel for porous SIMPLE momentum equations that multiplies
+ * the contribution by porosity.
+ */
+class LinearPNSFVMomentumPressure : public LinearFVMomentumPressure
+{
+public:
+  static InputParameters validParams();
+  LinearPNSFVMomentumPressure(const InputParameters & params);
+
+  virtual Real computeRightHandSideContribution() override;
+
+protected:
+  /// Porosity functor
+  const Moose::Functor<Real> & _eps;
+};

modules/navier_stokes/include/userobjects/PorousRhieChowMassFlux.h

@@ -0,0 +1,96 @@
+//* This file is part of the MOOSE framework
+//* https://mooseframework.inl.gov
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "RhieChowMassFlux.h"
+#include "MooseFunctor.h"
+
+#include <unordered_map>
+
+#include "MooseEnum.h"
+
+/**
+ * Rhie-Chow face flux provider that treats porous media by incorporating porosity
+ * into the face mass fluxes and enabling manual or automatically detected pressure
+ * jumps at interfaces.
+ */
+class PorousRhieChowMassFlux : public RhieChowMassFlux
+{
+public:
+  static InputParameters validParams();
+  PorousRhieChowMassFlux(const InputParameters & params);
+
+  virtual void initFaceMassFlux() override;
+  virtual void populateCouplingFunctors(
+      const std::vector<std::unique_ptr<NumericVector<Number>>> & raw_hbya,
+      const std::vector<std::unique_ptr<NumericVector<Number>>> & raw_Ainv) override;
+  virtual void computeFaceMassFlux() override;
+  virtual Real getVolumetricFaceFlux(const FaceInfo & fi) const override;
+
+protected:
+  /// Simple container for user-supplied pressure jump parameters
+  struct JumpOptions
+  {
+    Real constant = 0.0;
+    Real form_factor = 0.0;
+    Real bernoulli = 0.0;
+  };
+
+  /// Returns true if the face is tagged manually as an interface
+  bool isManualInterfaceFace(const FaceInfo & fi) const;
+
+  /// Returns true if either manual tagging or porosity jump detection marks the face
+  bool treatAsInterface(const FaceInfo & fi, const Moose::StateArg & time) const;
+
+  /**
+   * Compute the pressure jump that must be applied at an interface face.
+   * @param fi Face information
+   * @param time Evaluation time/state
+   * @param downwind_is_elem Returns true if the element (as opposed to the neighbor)
+   *        is downstream of the flow direction
+   * @return Pair of (is jump active, pressure jump value)
+   */
+  std::pair<bool, Real>
+  computePressureJump(const FaceInfo & fi, const Moose::StateArg & time, bool & downwind_is_elem)
+      const;
+
+  /// Helper to evaluate porosity on a face argument
+  Real evaluatePorosity(const Moose::FaceArg & face, const Moose::StateArg & time) const;
+
+  /// Helper to evaluate porosity on an element
+  Real evaluatePorosity(const Moose::ElemArg & elem_arg, const Moose::StateArg & time) const;
+
+  /// Interpolate porosity to a face according to the selected method
+  Real interpolateFacePorosity(const FaceInfo & fi,
+                               const Moose::StateArg & time,
+                               Real elem_eps,
+                               Real neighbor_eps) const;
+
+  /// Multiply the supplied value by porosity using the appropriate interpolation
+  void applyPorosityWeighting(const FaceInfo & fi,
+                              const Moose::StateArg & time,
+                              Real elem_eps,
+                              Real neighbor_eps,
+                              Real elem_value,
+                              Real neighbor_value,
+                              Real & face_value) const;
+
+  /// Functor describing the porosity
+  const Moose::Functor<Real> & _eps;
+
+  /// Interpolation scheme used to evaluate porosity-dependent quantities on faces
+  MooseEnum _eps_face_interp_method;
+
+  /// Whether automatic block-to-block porosity jumps should trigger interface logic
+  const bool _detect_block_interfaces;
+
+  /// Map from boundary id to user-supplied pressure jump options
+  std::unordered_map<BoundaryID, JumpOptions> _interface_jump_data;
+};

modules/navier_stokes/include/userobjects/RhieChowMassFlux.h

@@ -44,7 +44,7 @@ class RhieChowMassFlux : public RhieChowFaceFluxProvider, public NonADFunctorInt
   Real getMassFlux(const FaceInfo & fi) const;
 
   /// Get the volumetric face flux (used in advection terms)
-  Real getVolumetricFaceFlux(const FaceInfo & fi) const;
+  virtual Real getVolumetricFaceFlux(const FaceInfo & fi) const;
 
   virtual Real getVolumetricFaceFlux(const Moose::FV::InterpMethod m,
                                      const FaceInfo & fi,
@@ -53,11 +53,11 @@ class RhieChowMassFlux : public RhieChowFaceFluxProvider, public NonADFunctorInt
                                      bool subtract_mesh_velocity) const override;
 
   /// Initialize the container for face velocities
-  void initFaceMassFlux();
+  virtual void initFaceMassFlux();
   /// Initialize the coupling fields (HbyA and Ainv)
   void initCouplingField();
   /// Update the values of the face velocities in the containers
-  void computeFaceMassFlux();
+  virtual void computeFaceMassFlux();
   /// Update the cell values of the velocity variables
   void computeCellVelocity();
 
@@ -93,7 +93,7 @@ class RhieChowMassFlux : public RhieChowFaceFluxProvider, public NonADFunctorInt
   void setupMeshInformation();
 
   /// Populate the face values of the H/A and 1/A fields
-  void
+  virtual void
   populateCouplingFunctors(const std::vector<std::unique_ptr<NumericVector<Number>>> & raw_hbya,
                            const std::vector<std::unique_ptr<NumericVector<Number>>> & raw_Ainv);
 
@@ -105,6 +105,9 @@ class RhieChowMassFlux : public RhieChowFaceFluxProvider, public NonADFunctorInt
 
   virtual bool supportMeshVelocity() const override { return false; }
 
+  /// Accessor for cached face info (to be used by derived classes)
+  const std::vector<const FaceInfo *> & flowFaceInfo() const { return _flow_face_info; }
+
   /// The \p MooseMesh that this user object operates on
   const MooseMesh & _moose_mesh;
 

modules/navier_stokes/src/linearfvkernels/LinearPNSFVMomentumPressure.C

@@ -0,0 +1,38 @@
+//* This file is part of the MOOSE framework
+//* https://mooseframework.inl.gov
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#include "LinearPNSFVMomentumPressure.h"
+#include "NS.h"
+
+registerMooseObject("NavierStokesApp", LinearPNSFVMomentumPressure);
+
+InputParameters
+LinearPNSFVMomentumPressure::validParams()
+{
+  InputParameters params = LinearFVMomentumPressure::validParams();
+  params.addClassDescription("Porous momentum pressure gradient kernel that multiplies the "
+                             "contribution by the local porosity.");
+  params.addRequiredParam<MooseFunctorName>(NS::porosity, "Porosity functor.");
+  return params;
+}
+
+LinearPNSFVMomentumPressure::LinearPNSFVMomentumPressure(const InputParameters & params)
+  : LinearFVMomentumPressure(params), _eps(getFunctor<Real>(NS::porosity))
+{
+}
+
+Real
+LinearPNSFVMomentumPressure::computeRightHandSideContribution()
+{
+  const auto dof_value = _current_elem_info->dofIndices()[_pressure_sys_num][_pressure_var_num];
+  const auto elem_arg = makeElemArg(_current_elem_info->elem());
+  const auto state = determineState();
+  const Real porosity = _eps(elem_arg, state);
+  return -porosity * (*_pressure_gradient[_index])(dof_value) * _current_elem_volume;
+}