PorousRhieChowMassFlux

Rhie-Chow mass flux object specialized for porous flow/baffle cases.

Overview

PorousRhieChowMassFlux extends RhieChowMassFlux for porous-medium and porous-baffle calculations in the linear finite volume SIMPLE workflow.

Compared with the base object, it adds:

cell-porosity scaling of the pressure-coupling fields,
porous-baffle pressure jumps stored and relaxed on coupled interface faces,
compatibility with FVPressureJumpGreenGaussGradient, which reconstructs the pressure separately on each side of every baffle, and
a pressure-velocity coupling gradient that remains separate from the solution gradient used by the baffle diffusion correction.

Pressure jumps are supplied by the models listed in "pressure_jump_models", such as BernoulliFormLossPressureJump. The jump is updated from the current face mass flux and under-relaxed with "pressure_jump_relaxation". See Porous Media Weakly Compressible Linear Finite Volume for the discrete coupling.

Before each pressure corrector, the object updates the jump from the preceding corrected face flux. With PIMPLE, this happens before each of the $var element = document.getElementById("moose-equation-c9a4817e-ff98-4f2a-a0aa-a49f66cecdf5");katex.render("n+1", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ pressure correctors of an outer iteration, so the effective relaxation grows with "num_piso_iterations". Without a reconstructed pressure gradient, the face flux is only recomputed after the last corrector, so a relaxation factor $var element = document.getElementById("moose-equation-816110b0-beed-4083-8a4d-4148caebd57f");katex.render("\\omega", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ acts per outer iteration as $var element = document.getElementById("moose-equation-862072c5-6483-4c61-94ee-bbabf4c8bbbf");katex.render("1 - (1 - \\omega)^{n+1}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ . The pressure variable's FVPressureJumpGreenGaussGradient then removes that jump when interpolating pressure across each baffle. A zero jump recovers the standard Green-Gauss face value exactly.

The tangential superficial velocity also jumps across a baffle wherever porosity does. Density and viscosity must vary continuously. Every velocity variable must therefore use FVVelocityJumpGreenGaussGradient as its default gradient method. On baffle faces the $var element = document.getElementById("moose-equation-c1e352d4-60c8-47a5-aebb-71950c7ccc53");katex.render("\\rho\\mathbf{H}/\\mathbf{A}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ flux of each side is weighted with the same half-cell conductances as the pressure (see LinearFVPressureCorrectionDiffusionJump), so neither the velocity jump nor a body force balanced by the pressure on each side is read as a mass-flux error.

This is the user object expected by LinearPWCNSFVMomentumFlux, LinearFVPressureCorrectionDiffusionJump, and porous uses of LinearFVEnergyAdvection.

In this example, PorousRhieChowMassFlux uses a separately configured pressure-jump model on two porous interfaces.

[UserObjects<<<{"href": "../../syntax/UserObjects/index.html"}>>>]
  [pressure_jump]
    type = BernoulliFormLossPressureJump<<<{"description": "Provides Bernoulli and form-loss pressure jumps to PorousRhieChowMassFlux.", "href": "BernoulliFormLossPressureJump.html"}>>>
    boundary<<<{"description": "The boundaries on which this pressure jump model applies."}>>> = 'baffle baffle2'
    porosity<<<{"description": "The porosity functor. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number."}>>> = porosity
    density<<<{"description": "The density functor. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number."}>>> = ${rho}
  []
  [rc]
    type = PorousRhieChowMassFlux<<<{"description": "Rhie-Chow mass flux object specialized for porous flow/baffle cases.", "href": "PorousRhieChowMassFlux.html"}>>>
    u<<<{"description": "The x-component of velocity"}>>> = superficial_u
    pressure<<<{"description": "The pressure variable."}>>> = pressure
    rho<<<{"description": "Density functor. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number."}>>> = ${rho}
    porosity<<<{"description": "Porosity functor (defaults to 1 for non-porous calculations). A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number."}>>> = porosity
    p_diffusion_kernel<<<{"description": "The LinearFVPressureCorrectionDiffusion kernel acting on the pressure."}>>> = p_diffusion
    pressure_jump_models<<<{"description": "Pressure jump models applied by this object."}>>> = pressure_jump
    pressure_jump_relaxation<<<{"description": "Under-relaxation factor for pressure jump updates."}>>> = 0.1
  []
[]

Input Parameters

p_diffusion_kernelThe LinearFVPressureCorrectionDiffusion kernel acting on the pressure.
C++ Type:std::string
Controllable:No
Description:The LinearFVPressureCorrectionDiffusion kernel acting on the pressure.
pressureThe pressure variable.
C++ Type:VariableName
Unit:(no unit assumed)
Controllable:No
Description:The pressure variable.
rhoDensity functor. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number.
C++ Type:MooseFunctorName
Unit:(no unit assumed)
Controllable:No
Description:Density functor. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number.
uThe x-component of velocity
C++ Type:VariableName
Unit:(no unit assumed)
Controllable:No
Description:The x-component of velocity

Required Parameters

blockThe list of blocks (ids or names) that this object will be applied
C++ Type:std::vector<SubdomainName>
Controllable:No
Description:The list of blocks (ids or names) that this object will be applied
porosity1Porosity functor (defaults to 1 for non-porous calculations). A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number.
Default:1
C++ Type:MooseFunctorName
Unit:(no unit assumed)
Controllable:No
Description:Porosity functor (defaults to 1 for non-porous calculations). A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number.
pressure_diffusion_interpolationaverageThe face interpolation method for Ainv in the pressure correction diffusion term.
Default:average
C++ Type:MooseEnum
Options:average, harmonic
Controllable:No
Description:The face interpolation method for Ainv in the pressure correction diffusion term.
pressure_jump_modelsPressure jump models applied by this object.
C++ Type:std::vector<UserObjectName>
Controllable:No
Description:Pressure jump models applied by this object.
pressure_jump_relaxation1Under-relaxation factor for pressure jump updates.
Default:1
C++ Type:Real
Unit:(no unit assumed)
Range:0.0
Controllable:No
Description:Under-relaxation factor for pressure jump updates.
pressure_projection_methodstandardThe method to use in the pressure projection for Ainv - standard (SIMPLE) or consistent (SIMPLEC)
Default:standard
C++ Type:MooseEnum
Options:standard, consistent
Controllable:No
Description:The method to use in the pressure projection for Ainv - standard (SIMPLE) or consistent (SIMPLEC)
vThe y-component of velocity
C++ Type:VariableName
Unit:(no unit assumed)
Controllable:No
Description:The y-component of velocity
wThe z-component of velocity
C++ Type:VariableName
Unit:(no unit assumed)
Controllable:No
Description:The z-component of velocity

Optional Parameters

allow_duplicate_execution_on_initialFalseIn the case where this UserObject is depended upon by an initial condition, allow it to be executed twice during the initial setup (once before the IC and again after mesh adaptivity (if applicable).
Default:False
C++ Type:bool
Controllable:No
Description:In the case where this UserObject is depended upon by an initial condition, allow it to be executed twice during the initial setup (once before the IC and again after mesh adaptivity (if applicable).
execution_order_group0Execution order groups are executed in increasing order (e.g., the lowest number is executed first). Note that negative group numbers may be used to execute groups before the default (0) group. Please refer to the user object documentation for ordering of user object execution within a group.
Default:0
C++ Type:int
Controllable:No
Description:Execution order groups are executed in increasing order (e.g., the lowest number is executed first). Note that negative group numbers may be used to execute groups before the default (0) group. Please refer to the user object documentation for ordering of user object execution within a group.
force_postauxFalseForces the UserObject to be executed in POSTAUX
Default:False
C++ Type:bool
Controllable:No
Description:Forces the UserObject to be executed in POSTAUX
force_preauxFalseForces the UserObject to be executed in PREAUX
Default:False
C++ Type:bool
Controllable:No
Description:Forces the UserObject to be executed in PREAUX
force_preicFalseForces the UserObject to be executed in PREIC during initial setup
Default:False
C++ Type:bool
Controllable:No
Description:Forces the UserObject to be executed in PREIC during initial setup

Execution Scheduling Parameters

control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Controllable:No
Description:Adds user-defined labels for accessing object parameters via control logic.
enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Controllable:Yes
Description:Set the enabled status of the MooseObject.
use_displaced_meshFalseWhether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Default:False
C++ Type:bool
Controllable:No
Description:Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.

Advanced Parameters

prop_getter_suffixAn optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:An optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
use_interpolated_stateFalseFor the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.
Default:False
C++ Type:bool
Controllable:No
Description:For the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.

Material Property Retrieval Parameters

Input Files

(modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/linear-segregated/rz-pipe-channel/porous-baffle-rz.i)
(modules/navier_stokes/test/tests/finite_volume/pins/mms/linear-segregated/smooth-properties/smooth-properties.i)
(modules/navier_stokes/test/tests/finite_volume/ins/mms/linear-segregated/porosity-jump-traction/porosity-jump-traction.i)
(modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/linear-segregated/2d-simple-channel-enthalpy/porous-baffle-2d-enthalpy.i)
(modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/linear-segregated/errors/pressure-jump-errors.i)
(modules/navier_stokes/test/tests/finite_volume/pins/mms/linear-segregated/interface-reconstruction/interface-reconstruction.i)
(modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/linear-segregated/2d-diagonal-baffle/porous-baffle-straight-diagonal.i)
(modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/linear-segregated/2d-variable-density/porous-baffle-2d-variable-density-bernoulli.i)
(modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/linear-segregated/1d-simple-channel/porous-baffle-1d.i)
(modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/linear-segregated/radial-baffle-rz/radial-baffle-rz.i)
(modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/linear-segregated/1d-simple-channel/form-loss-orientation.i)
(modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/linear-segregated/2d-simple-channel-enthalpy/porous-baffle-2d-enthalpy-fluid-only.i)
(modules/navier_stokes/test/tests/finite_volume/pins/mms/linear-segregated/hydrostatic-jump/hydrostatic-jump.i)

(modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/linear-segregated/1d-simple-channel/porous-baffle-1d.i)

# One-dimensional verification of Bernoulli jumps, porous drag, and form loss.
# Mass conservation gives u_block_1 = u_block_2 = 0.1 and u_block_jump = 0.
# Without drag or form loss, the two pressure jumps (right minus left) are 0.03 and -0.03,
# so p_left = p_right = p_jump = 0. With mu = 1 and Forchheimer = 50, the porous block
# adds a pressure drop of 0.25. A form-loss coefficient of 10 adds 0.4 at each baffle,
# giving p_jump = p_left - p_right = 1.05 and baffle jumps of -0.37 and -0.43.
# Referencing the higher-porosity velocity 0.1 instead of the lower-porosity velocity 0.2 makes the
# form loss 0.1 at that baffle. With reference_velocity_side = 'lower_porosity higher_porosity',
# the baffle jumps are -0.37 and -0.13; with 'higher_porosity lower_porosity', they are -0.07 and
# -0.43. Both give p_jump = 0.75.

mu = 0.0 # 1e-2
rho = 2.0
advected_interp_method = 'average'
forchheimer_value = 50
inlet_u = 0.1

[FVGradientMethods]
  [jump_aware]
    type = FVPressureJumpGreenGaussGradient
  []
  [velocity_jump_aware]
    type = FVVelocityJumpGreenGaussGradient
  []
  [reconstructed]
    type = FVReconstructedPressureGradient
    base_gradient_method = jump_aware
    gradient_relaxation = 1.0
  []
[]

[Mesh]
  [mesh]
    type = CartesianMeshGenerator
    dim = 1
    dx = '0.5 0.5 0.5'
    ix = '200 200 200'
    subdomain_id = '1 2 3'
  []
  [baffle]
    type = SideSetsBetweenSubdomainsGenerator
    input = mesh
    primary_block = '1'
    paired_block = '2'
    new_boundary = 'baffle'
  []
  [baffle2]
    type = SideSetsBetweenSubdomainsGenerator
    input = baffle
    primary_block = '2'
    paired_block = '3'
    new_boundary = 'baffle2'
  []
[]

[Problem]
  linear_sys_names = 'u_system pressure_system'
  previous_nl_solution_required = true
[]

[UserObjects]
  [pressure_jump]
    type = BernoulliFormLossPressureJump
    boundary = 'baffle baffle2'
    porosity = porosity
    density = ${rho}
  []
  [rc]
    type = PorousRhieChowMassFlux
    u = superficial_u
    pressure = pressure
    rho = ${rho}
    porosity = porosity
    p_diffusion_kernel = p_diffusion
    pressure_jump_models = pressure_jump
    pressure_jump_relaxation = 0.1
  []
[]

[Variables]
  [superficial_u]
    type = MooseLinearVariableFVReal
    solver_sys = u_system
    initial_condition = ${inlet_u}
    gradient_method = velocity_jump_aware
  []
  [pressure]
    type = MooseLinearVariableFVReal
    solver_sys = pressure_system
    initial_condition = 0.0
    gradient_method = jump_aware
  []
[]

[FVInterpolationMethods]
  [average]
    type = FVGeometricAverage
  []
[]

[LinearFVKernels]
  [u_advection]
    type = LinearPWCNSFVMomentumFlux
    variable = superficial_u
    advected_interp_method_name = ${advected_interp_method}
    mu = ${mu}
    u = superficial_u
    momentum_component = 'x'
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = false
  []
  [u_pressure]
    type = LinearPWCNSFVMomentumPressure
    variable = superficial_u
    momentum_component = 'x'
    porosity = porosity
    pressure = pressure
    gradient_method = reconstructed
  []
  [u_friction]
    type = LinearFVMomentumFriction
    variable = superficial_u
    Forchheimer_name = forch
    porosity = porosity
    rho = ${rho}
    u = superficial_u
    momentum_component = 'x'
    block = 2
  []
  [p_diffusion]
    type = LinearFVPressureCorrectionDiffusionJump
    variable = pressure
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = false
  []
  [HbyA_divergence]
    type = LinearFVDivergence
    variable = pressure
    face_flux = HbyA
    force_boundary_execution = true
  []
[]

[LinearFVBCs]
  [inlet_u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = left
    variable = superficial_u
    functor = ${inlet_u}
  []
  [outlet_u]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = right
    variable = superficial_u
    use_two_term_expansion = false
  []
  [outlet_p]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = right
    variable = pressure
    functor = 0.0
  []
[]

[FunctorMaterials]
  [forch]
    type = GenericVectorFunctorMaterial
    prop_names = forch
    prop_values = '${forchheimer_value} ${forchheimer_value} ${forchheimer_value}'
  []
  [porosity]
    type = PiecewiseByBlockFunctorMaterial
    prop_name = porosity
    subdomain_to_prop_value = '1 0.5 2 1.0 3 0.5'
  []
[]

[Postprocessors]
  [p_left]
    type = SideAverageValue
    variable = pressure
    boundary = left
  []
  [p_right]
    type = SideAverageValue
    variable = pressure
    boundary = right
  []
  [p_jump]
    type = ParsedPostprocessor
    expression = 'p_left - p_right'
    pp_names = 'p_left p_right'
  []
  # These points are the centers of the cells adjacent to x = 0.5 and x = 1.0. The mesh spacing
  # is 0.5 / 200 = 0.0025, so the adjacent centers are offset by half that distance.
  [p_baffle_1_left]
    type = PointValue
    variable = pressure
    point = '0.49875 0 0'
  []
  [p_baffle_1_right]
    type = PointValue
    variable = pressure
    point = '0.50125 0 0'
  []
  [p_baffle_2_left]
    type = PointValue
    variable = pressure
    point = '0.99875 0 0'
  []
  [p_baffle_2_right]
    type = PointValue
    variable = pressure
    point = '1.00125 0 0'
  []
  [p_baffle_1_jump]
    type = ParsedPostprocessor
    expression = 'p_baffle_1_right - p_baffle_1_left'
    pp_names = 'p_baffle_1_right p_baffle_1_left'
  []
  [p_baffle_2_jump]
    type = ParsedPostprocessor
    expression = 'p_baffle_2_right - p_baffle_2_left'
    pp_names = 'p_baffle_2_right p_baffle_2_left'
  []
  [u_block_1]
    type = ElementAverageValue
    variable = superficial_u
    block = 1
  []
  [u_block_2]
    type = ElementAverageValue
    variable = superficial_u
    block = 2
  []
  [u_block_jump]
    type = ParsedPostprocessor
    expression = 'u_block_1 - u_block_2'
    pp_names = 'u_block_1 u_block_2'
  []
[]

[AuxVariables]
  [porosity_aux]
    type = MooseLinearVariableFVReal
  []
[]

[AuxKernels]
  [por]
    type = FunctorAux
    variable = porosity_aux
    functor = porosity
    execute_on = 'timestep_end'
  []
[]

[Executioner]
  type = SIMPLE
  momentum_l_abs_tol = 1e-12
  pressure_l_abs_tol = 1e-12
  momentum_l_tol = 0
  pressure_l_tol = 0
  rhie_chow_user_object = rc
  momentum_systems = 'u_system'
  pressure_system = pressure_system
  momentum_equation_relaxation = 0.3
  pressure_variable_relaxation = 0.1
  num_iterations = 1000
  pressure_absolute_tolerance = 1e-8
  # The exact superficial velocity is uniform, so the momentum normalization factor
  # sum(|Ax - Ax_avg| + |b - Ax_avg|) nearly vanishes and the normalized momentum residual measures
  # round-off. It levels off at a value that depends on the partitioning (near 1e-2 on one process,
  # near 1e-1 on three), so only the pressure residual is used as a convergence criterion.
  momentum_absolute_tolerance = 1
  momentum_petsc_options_iname = '-pc_type -pc_hypre_type'
  momentum_petsc_options_value = 'hypre boomeramg'
  pressure_petsc_options_iname = '-pc_type -pc_hypre_type'
  pressure_petsc_options_value = 'hypre boomeramg'
[]

[Outputs]
  csv = true
  execute_on = 'timestep_end'
[]

(modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/linear-segregated/rz-pipe-channel/porous-baffle-rz.i)

# Axisymmetric pipe version of ../1d-simple-channel/porous-baffle-1d.i: the same axial flow in an
# RZ pipe of radius 0.5 with symmetry conditions on the axis and on the wall (slip), so the
# solution is independent of r and must reproduce the one-dimensional results. The comments below
# are those of the one-dimensional input.
# One-dimensional verification of Bernoulli jumps, porous drag, and form loss.
# Mass conservation gives u_block_1 = u_block_2 = 0.1 and u_block_jump = 0.
# Without drag or form loss, the two pressure jumps (right minus left) are 0.03 and -0.03,
# so p_left = p_right = p_jump = 0. With mu = 1 and Forchheimer = 50, the porous block
# adds a pressure drop of 0.25. A form-loss coefficient of 10 adds 0.4 at each baffle,
# giving p_jump = p_left - p_right = 1.05 and baffle jumps of -0.37 and -0.43.

mu = 0.0 # 1e-2
rho = 2.0
advected_interp_method = 'average'
forchheimer_value = 50
inlet_u = 0.1

[FVGradientMethods]
  [jump_aware]
    type = FVPressureJumpGreenGaussGradient
  []
  [velocity_jump_aware]
    type = FVVelocityJumpGreenGaussGradient
  []
  [reconstructed]
    type = FVReconstructedPressureGradient
    base_gradient_method = jump_aware
    gradient_relaxation = 1.0
  []
[]

[Mesh]
  [mesh]
    type = CartesianMeshGenerator
    dim = 2
    dx = '0.5 0.5 0.5'
    ix = '200 200 200'
    dy = '0.5'
    iy = '2'
    subdomain_id = '1 2 3'
  []
  [baffle]
    type = SideSetsBetweenSubdomainsGenerator
    input = mesh
    primary_block = '1'
    paired_block = '2'
    new_boundary = 'baffle'
  []
  [baffle2]
    type = SideSetsBetweenSubdomainsGenerator
    input = baffle
    primary_block = '2'
    paired_block = '3'
    new_boundary = 'baffle2'
  []
  coord_type = RZ
  rz_coord_axis = X
[]

[Problem]
  linear_sys_names = 'u_system v_system pressure_system'
  previous_nl_solution_required = true
[]

[UserObjects]
  [pressure_jump]
    type = BernoulliFormLossPressureJump
    boundary = 'baffle baffle2'
    porosity = porosity
    density = ${rho}
  []
  [rc]
    type = PorousRhieChowMassFlux
    u = superficial_u
    v = superficial_v
    pressure = pressure
    rho = ${rho}
    porosity = porosity
    p_diffusion_kernel = p_diffusion
    pressure_jump_models = pressure_jump
    pressure_jump_relaxation = 0.1
  []
[]

[Variables]
  [superficial_u]
    type = MooseLinearVariableFVReal
    solver_sys = u_system
    initial_condition = ${inlet_u}
    gradient_method = velocity_jump_aware
  []
  [superficial_v]
    type = MooseLinearVariableFVReal
    solver_sys = v_system
    initial_condition = 0
    gradient_method = velocity_jump_aware
  []
  [pressure]
    type = MooseLinearVariableFVReal
    solver_sys = pressure_system
    initial_condition = 0.0
    gradient_method = jump_aware
  []
[]

[FVInterpolationMethods]
  [average]
    type = FVGeometricAverage
  []
[]

[LinearFVKernels]
  [u_advection]
    type = LinearPWCNSFVMomentumFlux
    variable = superficial_u
    advected_interp_method_name = ${advected_interp_method}
    mu = ${mu}
    u = superficial_u
    v = superficial_v
    momentum_component = 'x'
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = false
  []
  [v_advection]
    type = LinearPWCNSFVMomentumFlux
    variable = superficial_v
    advected_interp_method_name = ${advected_interp_method}
    mu = ${mu}
    u = superficial_u
    v = superficial_v
    momentum_component = 'y'
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = false
  []
  [u_pressure]
    type = LinearPWCNSFVMomentumPressure
    variable = superficial_u
    momentum_component = 'x'
    porosity = porosity
    pressure = pressure
    gradient_method = reconstructed
  []
  [v_pressure]
    type = LinearPWCNSFVMomentumPressure
    variable = superficial_v
    momentum_component = 'y'
    porosity = porosity
    pressure = pressure
    gradient_method = reconstructed
  []
  [u_friction]
    type = LinearFVMomentumFriction
    variable = superficial_u
    Forchheimer_name = forch
    porosity = porosity
    rho = ${rho}
    u = superficial_u
    v = superficial_v
    momentum_component = 'x'
    block = 2
  []
  [v_friction]
    type = LinearFVMomentumFriction
    variable = superficial_v
    Forchheimer_name = forch
    porosity = porosity
    rho = ${rho}
    u = superficial_u
    v = superficial_v
    momentum_component = 'y'
    block = 2
  []
  [p_diffusion]
    type = LinearFVPressureCorrectionDiffusionJump
    variable = pressure
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = false
  []
  [HbyA_divergence]
    type = LinearFVDivergence
    variable = pressure
    face_flux = HbyA
    force_boundary_execution = true
  []
[]

[LinearFVBCs]
  [inlet_u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = left
    variable = superficial_u
    functor = ${inlet_u}
  []
  [inlet_v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = left
    variable = superficial_v
    functor = 0
  []
  [outlet_u]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = right
    variable = superficial_u
    use_two_term_expansion = false
  []
  [outlet_v]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = right
    variable = superficial_v
    use_two_term_expansion = false
  []
  # The axis (bottom) and the slip wall (top) are both symmetry planes for this flow.
  [symmetry_u]
    type = LinearFVVelocitySymmetryBC
    boundary = 'bottom top'
    variable = superficial_u
    u = superficial_u
    v = superficial_v
    momentum_component = 'x'
  []
  [symmetry_v]
    type = LinearFVVelocitySymmetryBC
    boundary = 'bottom top'
    variable = superficial_v
    u = superficial_u
    v = superficial_v
    momentum_component = 'y'
  []
  [symmetry_p]
    type = LinearFVPressureSymmetryBC
    boundary = 'bottom top'
    variable = pressure
    HbyA_flux = HbyA
  []
  [outlet_p]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = right
    variable = pressure
    functor = 0.0
  []
[]

[FunctorMaterials]
  [forch]
    type = GenericVectorFunctorMaterial
    prop_names = forch
    prop_values = '${forchheimer_value} ${forchheimer_value} ${forchheimer_value}'
  []
  [porosity]
    type = PiecewiseByBlockFunctorMaterial
    prop_name = porosity
    subdomain_to_prop_value = '1 0.5 2 1.0 3 0.5'
  []
[]

[Postprocessors]
  [p_left]
    type = SideAverageValue
    variable = pressure
    boundary = left
  []
  [p_right]
    type = SideAverageValue
    variable = pressure
    boundary = right
  []
  [p_jump]
    type = ParsedPostprocessor
    expression = 'p_left - p_right'
    pp_names = 'p_left p_right'
  []
  # These points are the centers of the cells adjacent to x = 0.5 and x = 1.0. The mesh spacing
  # is 0.5 / 200 = 0.0025, so the adjacent centers are offset by half that distance. The radial
  # coordinate 0.375 is the center of the outer row of cells.
  [p_baffle_1_left]
    type = PointValue
    variable = pressure
    point = '0.49875 0.375 0'
  []
  [p_baffle_1_right]
    type = PointValue
    variable = pressure
    point = '0.50125 0.375 0'
  []
  [p_baffle_2_left]
    type = PointValue
    variable = pressure
    point = '0.99875 0.375 0'
  []
  [p_baffle_2_right]
    type = PointValue
    variable = pressure
    point = '1.00125 0.375 0'
  []
  [p_baffle_1_jump]
    type = ParsedPostprocessor
    expression = 'p_baffle_1_right - p_baffle_1_left'
    pp_names = 'p_baffle_1_right p_baffle_1_left'
  []
  [p_baffle_2_jump]
    type = ParsedPostprocessor
    expression = 'p_baffle_2_right - p_baffle_2_left'
    pp_names = 'p_baffle_2_right p_baffle_2_left'
  []
  [u_block_1]
    type = ElementAverageValue
    variable = superficial_u
    block = 1
  []
  [u_block_2]
    type = ElementAverageValue
    variable = superficial_u
    block = 2
  []
  [u_block_jump]
    type = ParsedPostprocessor
    expression = 'u_block_1 - u_block_2'
    pp_names = 'u_block_1 u_block_2'
  []
[]

[AuxVariables]
  [porosity_aux]
    type = MooseLinearVariableFVReal
  []
[]

[AuxKernels]
  [por]
    type = FunctorAux
    variable = porosity_aux
    functor = porosity
    execute_on = 'timestep_end'
  []
[]

[Executioner]
  type = SIMPLE
  momentum_l_abs_tol = 1e-12
  pressure_l_abs_tol = 1e-12
  momentum_l_tol = 0
  pressure_l_tol = 0
  rhie_chow_user_object = rc
  momentum_systems = 'u_system v_system'
  pressure_system = pressure_system
  momentum_equation_relaxation = 0.3
  pressure_variable_relaxation = 0.1
  num_iterations = 1000
  pressure_absolute_tolerance = 1e-8
  # The Forchheimer case plateaus near 1e-2 after pressure and velocity stabilize; 2e-2 stops
  # with the total pressure drop within 2e-6 relative to the value after 1000 iterations.
  momentum_absolute_tolerance = 2e-2
  momentum_petsc_options_iname = '-pc_type -pc_hypre_type'
  momentum_petsc_options_value = 'hypre boomeramg'
  pressure_petsc_options_iname = '-pc_type -pc_hypre_type'
  pressure_petsc_options_value = 'hypre boomeramg'
[]

[Outputs]
  csv = true
  execute_on = 'timestep_end'
[]

(modules/navier_stokes/test/tests/finite_volume/pins/mms/linear-segregated/smooth-properties/smooth-properties.i)

# Manufactured solution with continuously varying density and viscosity, optionally crossed by a
# pressure baffle with a porosity jump. In coordinates normal (xn) and tangent (xt) to the line
# x - skew * y = 0.5:
#  - rho = 1 + xn/2 + xt/4 + xn xt/5 and mu = mu0 (m0 + 4 xn^2)(1 + xt^2) are smooth. With m0 = 0.05
#    the viscosity varies by a factor of 40 across the domain and by a factor of several between
#    neighboring cells on the coarse meshes.
#  - rho U = curl(psi), psi = xt/2 + xt^2/10 + eps (xn (xt/10 - 3/10) + xn^2/5), so the mass flux is
#    divergence-free. U.n is continuous across xn = 0 and U_t jumps with eps, while the interstitial
#    tangential velocity U_t / eps is continuous.
#  - The pressure is smooth on each side. With a baffle, the right side adds the
#    BernoulliFormLossPressureJump evaluated with the interface density and normal velocity.
# The forcing contains the porous advection div(rho U U / eps), the full stress
# div(mu (grad U + grad U^T - 2/3 div U I)) of superficial velocity, eps grad p, and a reaction k U.
# Without a baffle (baffle = 0), eps_left = eps_right must hold.
# The expressions were generated with sympy.

eps_left = 1
eps_right = 1
form_loss = 0
mu0 = 0.01
m0 = 0.05
reaction = 1
skew = 0
nx = 8
ny = 8
# Set baffle = 1 together with UserObjects/rc/pressure_jump_models=pressure_jump to make the
# interface a pressure baffle; otherwise it is an ordinary internal face.
baffle = 0
deviatoric = true

# Porosity of the side whose interstitial velocity is the form-loss reference: the lower one, the
# BernoulliFormLossPressureJump default.
eps_ref = '${fparse min(eps_left, eps_right)}'

[Problem]
  linear_sys_names = 'u_system v_system pressure_system'
  previous_nl_solution_required = true
[]

[Mesh]
  [generated]
    type = GeneratedMeshGenerator
    dim = 2
    nx = ${nx}
    ny = ${ny}
  []
  [right]
    type = ParsedSubdomainMeshGenerator
    input = generated
    combinatorial_geometry = 'x > 0.5'
    block_id = 2
    block_name = right
  []
  [left]
    type = ParsedSubdomainMeshGenerator
    input = right
    combinatorial_geometry = 'x < 0.5'
    block_id = 1
    block_name = left
  []
  [skew]
    type = ParsedNodeTransformGenerator
    input = left
    x_function = 'x + ${skew} * y'
  []
  [baffle]
    type = SideSetsBetweenSubdomainsGenerator
    input = skew
    primary_block = left
    paired_block = right
    new_boundary = baffle
  []
[]

[FVGradientMethods]
  [jump_aware]
    type = FVPressureJumpGreenGaussGradient
  []
  [velocity_jump_aware]
    type = FVVelocityJumpGreenGaussGradient
  []
  # The reconstructed pressure gradient keeps the cell pressure gradient consistent with the face
  # fluxes and suppresses odd-even pressure modes.
  [reconstructed]
    type = FVReconstructedPressureGradient
    base_gradient_method = jump_aware
    gradient_relaxation = 1
  []
[]

[FVInterpolationMethods]
  [average]
    type = FVGeometricAverage
  []
[]

[UserObjects]
  [pressure_jump]
    type = BernoulliFormLossPressureJump
    boundary = baffle
    porosity = porosity
    density = density
    form_loss = ${form_loss}
  []
  [rc]
    type = PorousRhieChowMassFlux
    u = superficial_u
    v = superficial_v
    pressure = pressure
    rho = density
    porosity = porosity
    p_diffusion_kernel = p_diffusion
    pressure_projection_method = consistent
  []
[]

[Variables]
  [superficial_u]
    type = MooseLinearVariableFVReal
    solver_sys = u_system
    initial_condition = 0
    gradient_method = velocity_jump_aware
  []
  [superficial_v]
    type = MooseLinearVariableFVReal
    solver_sys = v_system
    initial_condition = 0
    gradient_method = velocity_jump_aware
  []
  [pressure]
    type = MooseLinearVariableFVReal
    solver_sys = pressure_system
    initial_condition = 0
    gradient_method = jump_aware
  []
[]

[LinearFVKernels]
  [u_advection_stress]
    type = LinearPWCNSFVMomentumFlux
    variable = superficial_u
    advected_interp_method_name = average
    mu = viscosity
    u = superficial_u
    v = superficial_v
    momentum_component = x
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = true
    use_deviatoric_terms = ${deviatoric}
  []
  [v_advection_stress]
    type = LinearPWCNSFVMomentumFlux
    variable = superficial_v
    advected_interp_method_name = average
    mu = viscosity
    u = superficial_u
    v = superficial_v
    momentum_component = y
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = true
    use_deviatoric_terms = ${deviatoric}
  []
  [u_pressure]
    type = LinearPWCNSFVMomentumPressure
    variable = superficial_u
    pressure = pressure
    porosity = porosity
    momentum_component = x
    gradient_method = reconstructed
  []
  [v_pressure]
    type = LinearPWCNSFVMomentumPressure
    variable = superficial_v
    pressure = pressure
    porosity = porosity
    momentum_component = y
    gradient_method = reconstructed
  []
  [u_source]
    type = LinearFVSource
    variable = superficial_u
    source_density = forcing_u
  []
  [v_source]
    type = LinearFVSource
    variable = superficial_v
    source_density = forcing_v
  []
  [u_reaction]
    type = LinearFVReaction
    variable = superficial_u
    coeff = ${reaction}
  []
  [v_reaction]
    type = LinearFVReaction
    variable = superficial_v
    coeff = ${reaction}
  []
  [p_diffusion]
    type = LinearFVPressureCorrectionDiffusionJump
    variable = pressure
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = true
  []
  [HbyA_divergence]
    type = LinearFVDivergence
    variable = pressure
    face_flux = HbyA
    force_boundary_execution = true
  []
[]

[LinearFVBCs]
  [velocity_x]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left right top bottom'
    variable = superficial_u
    functor = exact_u
  []
  [velocity_y]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left right top bottom'
    variable = superficial_v
    functor = exact_v
  []
  [pressure_flux]
    type = LinearFVPressureFluxBC
    boundary = 'left right top bottom'
    variable = pressure
    HbyA_flux = HbyA
    Ainv = Ainv
    rho = density
    u = superficial_u
    v = superficial_v
    # The two-term boundary pressure expansion is inaccurate on this nonorthogonal mesh.
    use_two_term_expansion = false
  []
[]

[FunctorMaterials]
  [density]
    type = ParsedFunctorMaterial
    property_name = density
    functor_names = 'xn xt'
    expression = '1+xn/2+xt/4+xn*xt/5'
  []
  [viscosity]
    type = ParsedFunctorMaterial
    property_name = viscosity
    functor_names = 'xn xt'
    expression = '${mu0}*(${m0}+4*xn^2)*(1+xt^2)'
  []
  [porosity]
    type = ParsedFunctorMaterial
    property_name = porosity
    functor_names = 'exact_porosity'
    expression = 'exact_porosity'
  []
  [pressure_error]
    type = ParsedFunctorMaterial
    property_name = pressure_error
    functor_names = 'pressure exact_pressure'
    expression = 'pressure-exact_pressure'
  []
  [u_error]
    type = ParsedFunctorMaterial
    property_name = u_error
    functor_names = 'superficial_u exact_u'
    expression = 'superficial_u-exact_u'
  []
  [v_error]
    type = ParsedFunctorMaterial
    property_name = v_error
    functor_names = 'superficial_v exact_v'
    expression = 'superficial_v-exact_v'
  []
[]

[Functions]
  [xn]
    type = ParsedFunction
    expression = '(x-${skew}*y-0.5)/sqrt(1+${skew}^2)'
  []
  [xt]
    type = ParsedFunction
    expression = '(${skew}*x+y)/sqrt(1+${skew}^2)'
  []
  # Side indicator: 1 on the right of the interface when it is a baffle, otherwise 0.
  [rs]
    type = ParsedFunction
    expression = 'if(x-${skew}*y>0.5,${baffle},0)'
  []
  [exact_porosity]
    type = ParsedFunction
    expression = 'if(x-${skew}*y<0.5,${eps_left},${eps_right})'
  []
  [exact_normal_velocity]
    type = ParsedFunction
    symbol_names = 'xn xt eps'
    symbol_values = 'xn xt exact_porosity'
    expression = '2*(eps*xn + 2*xt + 5)/(4*xn*xt + 10*xn + 5*xt + 20)'
  []
  [exact_tangential_velocity]
    type = ParsedFunction
    symbol_names = 'xn xt eps'
    symbol_values = 'xn xt exact_porosity'
    expression = '2*eps*(-4*xn - xt + 3)/(4*xn*xt + 10*xn + 5*xt + 20)'
  []
  [exact_pressure]
    type = ParsedFunction
    symbol_names = 'xn xt rs el er eref K'
    symbol_values = 'xn xt rs ${eps_left} ${eps_right} ${eps_ref} ${form_loss}'
    expression = '(el^2*er^2*eref^2*(xt + 4)*((4*xn + 1)^2 + 2*(4*xn + 1)*(2*xt - 1) + 4*(2*xt - 1)^2)/160 - '
                 'rs*(2*xt + 5)^2*(K*el^2*er^2 + eref^2*(el^2 - er^2))/50)/(el^2*er^2*eref^2*(xt + 4))'
  []
  [forcing_normal]
    type = ParsedFunction
    symbol_names = 'xn xt eps mu0 m0 k'
    symbol_values = 'xn xt exact_porosity ${mu0} ${m0} ${reaction}'
    expression = '(3*eps^2*(2*xn + xt)*(4*xn*xt + 10*xn + 5*xt + 20)^3 + 60*eps*k*(eps*xn + 2*xt + '
                 '5)*(4*xn*xt + 10*xn + 5*xt + 20)^2 + 320*eps*mu0*xn*(xt^2 + 1)*(eps*(4*xn + 5)*(4*xn + xt '
                 '- 3) - 2*(2*xt + 5)*(eps*xn + 2*xt + 5))*(4*xn*xt + 10*xn + 5*xt + 20) + 40*eps*mu0*(m0 + '
                 '4*xn^2)*(xt^2 + 1)*(2*eps*(16*xn - 1)*(4*xn*xt + 10*xn + 5*xt + 20) - 3*(4*xn + '
                 '5)*(-2*eps*(2*xt + 5)*(4*xn + xt - 3) + (4*xn + 5)*(eps*xn + 2*xt + 5)) + 12*(2*xt + '
                 '5)*(eps*(4*xn*xt + 10*xn + 5*xt + 20) - 2*(2*xt + 5)*(eps*xn + 2*xt + 5)) + 4*(2*xt + '
                 '5)*(-eps*(4*xn + 5)*(4*xn + xt - 3) + 2*(2*xt + 5)*(eps*xn + 2*xt + 5)) - 3*(4*xn*xt + '
                 '10*xn + 5*xt + 20)*(2*eps*(4*xn + 5) + eps*(2*xt + 5) + 2*eps*(4*xn + xt - 3) - 8*xn - '
                 '10)) - 120*eps*mu0*(8*xn*(xt^2 + 1)*(eps*(4*xn*xt + 10*xn + 5*xt + 20) - 2*(2*xt + '
                 '5)*(eps*xn + 2*xt + 5)) + xt*(m0 + 4*xn^2)*(2*eps*(2*xt + 5)*(4*xn + xt - 3) + 2*(1 - '
                 '2*eps)*(4*xn*xt + 10*xn + 5*xt + 20) - (4*xn + 5)*(eps*xn + 2*xt + 5)))*(4*xn*xt + 10*xn + '
                 '5*xt + 20) + 6*(eps*(eps*xn - 8*xn + 11)*(4*xn*xt + 10*xn + 5*xt + 20) + (eps*(4*xn + '
                 '5)*(4*xn + xt - 3) - 2*(2*xt + 5)*(eps*xn + 2*xt + 5))*(eps*xn + 2*xt + 5))*(4*xn*xt + '
                 '10*xn + 5*xt + 20))/(30*eps*(4*xn*xt + 10*xn + 5*xt + 20)^3)'
  []
  [forcing_tangential]
    type = ParsedFunction
    symbol_names = 'xn xt eps rs el er eref K mu0 m0 k'
    symbol_values = 'xn xt exact_porosity rs ${eps_left} ${eps_right} ${eps_ref} ${form_loss} ${mu0} ${m0} ${reaction}'
    expression = '(-1200*el^2*eps*er^2*eref^2*k*(xt + 4)^2*(4*xn + xt - 3)*(4*xn*xt + 10*xn + 5*xt + 20)^2 + '
                 '1600*el^2*er^2*eref^2*mu0*xt*(m0 + 4*xn^2)*(xt + 4)^2*(eps*(4*xn + 5)*(4*xn + xt - 3) - '
                 '2*(2*xt + 5)*(eps*xn + 2*xt + 5))*(4*xn*xt + 10*xn + 5*xt + 20) + '
                 '400*el^2*er^2*eref^2*mu0*(m0 + 4*xn^2)*(xt + 4)^2*(xt^2 + 1)*(12*eps*(4*xn + 5)*(-4*xn*xt '
                 '- 10*xn - 5*xt + (4*xn + 5)*(4*xn + xt - 3) - 20) - 4*(4*xn + 5)*(eps*(4*xn + 5)*(4*xn + '
                 'xt - 3) - 2*(2*xt + 5)*(eps*xn + 2*xt + 5)) + 12*(2*xt + 5)*(2*eps*(2*xt + 5)*(4*xn + xt - '
                 '3) - (4*xn + 5)*(eps*xn + 2*xt + 5)) - 2*(4*eps*xn - eps*(4*xn + 5) + 16*xt + 40)*(4*xn*xt '
                 '+ 10*xn + 5*xt + 20) + 3*(4*xn*xt + 10*xn + 5*xt + 20)*(4*eps*xn + eps*(4*xn + 5) - '
                 '16*eps*(2*xt + 5) + 16*xt + 40)) - 4800*el^2*er^2*eref^2*mu0*(xt + 4)^2*(eps*xt*(m0 + '
                 '4*xn^2)*(-4*xn*xt - 10*xn - 5*xt + (4*xn + 5)*(4*xn + xt - 3) - 20) + 2*xn*(xt^2 + '
                 '1)*(2*eps*(2*xt + 5)*(4*xn + xt - 3) + 2*(1 - 2*eps)*(4*xn*xt + 10*xn + 5*xt + 20) - (4*xn '
                 '+ 5)*(eps*xn + 2*xt + 5)))*(4*xn*xt + 10*xn + 5*xt + 20) + 120*el^2*er^2*eref^2*(xt + '
                 '4)^2*(-(eps*(4*xn + 5)*(4*xn + xt - 3) - 2*(2*xt + 5)*(eps*xn + 2*xt + 5))*(4*xn + xt - 3) '
                 '+ (-4*eps*xn + eps*(4*xn + xt - 3) - 8*xt - 20)*(4*xn*xt + 10*xn + 5*xt + 20))*(4*xn*xt + '
                 '10*xn + 5*xt + 20) - 3*eps*(5*el^2*er^2*eref^2*(xt + 4)^2*(-4*xn - 8*xt + 3) + 4*rs*(2*xt '
                 '+ 5)*(4*K*el^2*er^2*(xt + 4) - K*el^2*er^2*(2*xt + 5) + eref^2*(-el^2 + er^2)*(2*xt + 5) + '
                 '4*eref^2*(el^2 - er^2)*(xt + 4)))*(4*xn*xt + 10*xn + 5*xt + '
                 '20)^3)/(600*el^2*er^2*eref^2*(xt + 4)^2*(4*xn*xt + 10*xn + 5*xt + 20)^3)'
  []
  [exact_u]
    type = ParsedFunction
    expression = '(un+${skew}*ut)/sqrt(1+${skew}^2)'
    symbol_names = 'un ut'
    symbol_values = 'exact_normal_velocity exact_tangential_velocity'
  []
  [exact_v]
    type = ParsedFunction
    expression = '(-${skew}*un+ut)/sqrt(1+${skew}^2)'
    symbol_names = 'un ut'
    symbol_values = 'exact_normal_velocity exact_tangential_velocity'
  []
  [forcing_u]
    type = ParsedFunction
    expression = '(fn+${skew}*ft)/sqrt(1+${skew}^2)'
    symbol_names = 'fn ft'
    symbol_values = 'forcing_normal forcing_tangential'
  []
  [forcing_v]
    type = ParsedFunction
    expression = '(-${skew}*fn+ft)/sqrt(1+${skew}^2)'
    symbol_names = 'fn ft'
    symbol_values = 'forcing_normal forcing_tangential'
  []
[]

[Executioner]
  type = SIMPLE
  momentum_l_abs_tol = 1e-12
  pressure_l_abs_tol = 1e-12
  momentum_l_tol = 0
  pressure_l_tol = 0
  momentum_l_max_its = 100
  pressure_l_max_its = 100
  rhie_chow_user_object = rc
  momentum_systems = 'u_system v_system'
  pressure_system = pressure_system
  momentum_equation_relaxation = 0.5
  pressure_variable_relaxation = 0.2
  num_iterations = 3000
  pressure_absolute_tolerance = 1e-9
  momentum_absolute_tolerance = 1e-9
  momentum_petsc_options_iname = '-pc_type -pc_hypre_type'
  momentum_petsc_options_value = 'hypre boomeramg'
  pressure_petsc_options_iname = '-pc_type -pc_hypre_type'
  pressure_petsc_options_value = 'hypre boomeramg'
  print_fields = false

  # The exact pressure is zero at (xn, xt) = (-1/4, 1/2), on the left side.
  pin_pressure = true
  pressure_pin_value = 0
  pressure_pin_point = '${fparse 0.5 + (-0.25 + skew * (0.5 - 0.5 * skew / sqrt(1 + skew^2))) / sqrt(1 + skew^2)} '
                     '${fparse (0.25 * skew + 0.5 - 0.5 * skew / sqrt(1 + skew^2)) / sqrt(1 + skew^2)} 0'
[]

[Outputs]
  [csv]
    type = CSV
    execute_on = FINAL
  []
[]

[Postprocessors]
  [h]
    type = AverageElementSize
    execute_on = FINAL
  []
  [L2u]
    type = ElementL2FunctorError
    approximate = superficial_u
    exact = exact_u
    execute_on = FINAL
  []
  [L2v]
    type = ElementL2FunctorError
    approximate = superficial_v
    exact = exact_v
    execute_on = FINAL
  []
  [L2p]
    type = ElementL2FunctorError
    approximate = pressure
    exact = exact_pressure
    execute_on = FINAL
  []
  # The maximum pointwise errors detect cell-to-cell oscillations that the L2 norms average out.
  [max_u_error]
    type = ElementExtremeFunctorValue
    functor = u_error
    value_type = max_abs
    execute_on = FINAL
  []
  [max_v_error]
    type = ElementExtremeFunctorValue
    functor = v_error
    value_type = max_abs
    execute_on = FINAL
  []
  [max_p_error]
    type = ElementExtremeFunctorValue
    functor = pressure_error
    value_type = max_abs
    execute_on = FINAL
  []
[]

(modules/navier_stokes/test/tests/finite_volume/ins/mms/linear-segregated/porosity-jump-traction/porosity-jump-traction.i)

rho = 1
mu = 1
eps_left = 0.4
eps_right = 0.8
interface = 0.5
skew = 0.35
normal_velocity = 0
tangential_velocity = 1
tangential_gradient = 0
tangential_curvature = 0.5
shear_traction = 2

# In interface-normal coordinates (r,t), the manufactured velocity is
# v_n = normal_velocity + tangential_gradient*t and
# v_t = tangential_velocity + (shear_traction/mu - tangential_gradient)*r
#       + tangential_curvature*r^2. It is divergence-free and
# mu*(d(v_t)/dr + d(v_n)/dt) = shear_traction on both sides of the porosity interface.

[Problem]
  linear_sys_names = 'u_system v_system pressure_system'
  previous_nl_solution_required = true
[]

[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 4
    ny = 4
  []
  [right]
    type = ParsedSubdomainMeshGenerator
    input = gen
    combinatorial_geometry = 'x > ${interface}'
    block_id = 1
    block_name = right
  []
  [left]
    type = ParsedSubdomainMeshGenerator
    input = right
    combinatorial_geometry = 'x < ${interface}'
    block_id = 2
    block_name = left
  []
  [skew]
    type = ParsedNodeTransformGenerator
    input = left
    x_function = 'x + ${skew} * y'
  []
[]

[UserObjects]
  [rc]
    type = PorousRhieChowMassFlux
    u = vel_x
    v = vel_y
    pressure = pressure
    rho = ${rho}
    porosity = porosity
    p_diffusion_kernel = p_diffusion
    pressure_projection_method = consistent
  []
[]

[Variables]
  [vel_x]
    type = MooseLinearVariableFVReal
    solver_sys = u_system
    initial_condition = 0
  []
  [vel_y]
    type = MooseLinearVariableFVReal
    solver_sys = v_system
    initial_condition = 0
  []
  [pressure]
    type = MooseLinearVariableFVReal
    solver_sys = pressure_system
    initial_condition = 0
  []
[]

[FVInterpolationMethods]
  [average]
    type = FVGeometricAverage
  []
[]

[LinearFVKernels]
  [u_advection_stress]
    type = LinearPWCNSFVMomentumFlux
    variable = vel_x
    advected_interp_method_name = average
    mu = dynamic_viscosity
    u = vel_x
    v = vel_y
    momentum_component = x
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = true
    use_deviatoric_terms = true
  []
  [v_advection_stress]
    type = LinearPWCNSFVMomentumFlux
    variable = vel_y
    advected_interp_method_name = average
    mu = dynamic_viscosity
    u = vel_x
    v = vel_y
    momentum_component = y
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = true
    use_deviatoric_terms = true
  []
  [u_pressure]
    type = LinearPWCNSFVMomentumPressure
    variable = vel_x
    pressure = pressure
    porosity = porosity
    momentum_component = x
  []
  [v_pressure]
    type = LinearPWCNSFVMomentumPressure
    variable = vel_y
    pressure = pressure
    porosity = porosity
    momentum_component = y
  []
  [u_forcing]
    type = LinearFVSource
    variable = vel_x
    source_density = forcing_u
  []
  [v_forcing]
    type = LinearFVSource
    variable = vel_y
    source_density = forcing_v
  []
  [p_diffusion]
    type = LinearFVPressureCorrectionDiffusion
    variable = pressure
    diffusion_tensor = Ainv
    use_nonorthogonal_correction = true
    use_nonorthogonal_correction_on_boundary = true
  []
  [HbyA_divergence]
    type = LinearFVDivergence
    variable = pressure
    face_flux = HbyA
    force_boundary_execution = true
  []
[]

[LinearFVBCs]
  [velocity_x]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left right top bottom'
    variable = vel_x
    functor = exact_u
  []
  [velocity_y]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left right top bottom'
    variable = vel_y
    functor = exact_v
  []
  [pressure-extrapoaltion]
    type = LinearFVPressureFluxBC
    boundary = 'left right top bottom'
    variable = pressure
    HbyA_flux = HbyA
    Ainv = Ainv
    rho = ${rho}
    u = vel_x
    v = vel_y
    # The two-term boundary pressure expansion leaves an O(1) pressure error in the boundary
    # cells on this nonorthogonal mesh, which limits the L2 pressure slope to about 0.5.
    use_two_term_expansion = false
  []
[]

[FunctorMaterials]
  [viscosity]
    type = GenericFunctorMaterial
    prop_names = dynamic_viscosity
    prop_values = ${mu}
  []
  [porosity_left]
    type = GenericFunctorMaterial
    block = left
    prop_names = porosity
    prop_values = ${eps_left}
  []
  [porosity_right]
    type = GenericFunctorMaterial
    block = right
    prop_names = porosity
    prop_values = ${eps_right}
  []
[]

[Functions]
  [exact_u]
    type = ParsedFunction
    expression = '(vn+skew*vt)/sqrt(1+skew^2)'
    symbol_names = 'skew vn vt'
    symbol_values = '${skew} exact_normal_velocity exact_tangential_velocity'
  []
  [exact_v]
    type = ParsedFunction
    expression = '(-skew*vn+vt)/sqrt(1+skew^2)'
    symbol_names = 'skew vn vt'
    symbol_values = '${skew} exact_normal_velocity exact_tangential_velocity'
  []
  [exact_normal_velocity]
    type = ParsedFunction
    expression = 'U+b*(skew*x+y)/sqrt(1+skew^2)'
    symbol_names = 'skew U b'
    symbol_values = '${skew} ${normal_velocity} ${tangential_gradient}'
  []
  [exact_tangential_velocity]
    type = ParsedFunction
    expression = 'Vt+(q/mu-b)*(x-skew*y-interface)/sqrt(1+skew^2)'
                 '+c*((x-skew*y-interface)/sqrt(1+skew^2))^2'
    symbol_names = 'interface skew mu q b Vt c'
    symbol_values = '${interface} ${skew} ${mu} ${shear_traction} ${tangential_gradient} ${tangential_velocity} ${tangential_curvature}'
  []
  [tangential_velocity_derivative]
    type = ParsedFunction
    expression = 'q/mu-b+2*c*(x-skew*y-interface)/sqrt(1+skew^2)'
    symbol_names = 'interface skew mu q b c'
    symbol_values = '${interface} ${skew} ${mu} ${shear_traction} ${tangential_gradient} ${tangential_curvature}'
  []
  [forcing_u]
    type = ParsedFunction
    expression = '(b*vt+skew*(dvtdr*vn-2*mu*c))/sqrt(1+skew^2)'
    symbol_names = 'skew mu vn vt dvtdr b c'
    symbol_values = '${skew} ${mu} exact_normal_velocity exact_tangential_velocity tangential_velocity_derivative ${tangential_gradient} ${tangential_curvature}'
  []
  [forcing_v]
    type = ParsedFunction
    expression = '(-skew*b*vt+dvtdr*vn-2*mu*c)/sqrt(1+skew^2)'
    symbol_names = 'skew mu vn vt dvtdr b c'
    symbol_values = '${skew} ${mu} exact_normal_velocity exact_tangential_velocity tangential_velocity_derivative ${tangential_gradient} ${tangential_curvature}'
  []
[]

[Executioner]
  type = SIMPLE
  momentum_l_abs_tol = 1e-11
  pressure_l_abs_tol = 1e-11
  momentum_l_max_its = 50
  pressure_l_max_its = 50
  momentum_l_tol = 0
  pressure_l_tol = 0
  rhie_chow_user_object = rc
  momentum_systems = 'u_system v_system'
  pressure_system = pressure_system
  momentum_equation_relaxation = 0.7
  pressure_variable_relaxation = 0.5
  num_iterations = 5000
  pressure_absolute_tolerance = 1e-10
  momentum_absolute_tolerance = 1e-10
  momentum_petsc_options_iname = '-pc_type -pc_hypre_type'
  momentum_petsc_options_value = 'hypre boomeramg'
  pressure_petsc_options_iname = '-pc_type -pc_hypre_type'
  pressure_petsc_options_value = 'hypre boomeramg'
  print_fields = false

  pin_pressure = true
  pressure_pin_value = 0
  pressure_pin_point = '0.425 0.5 0'
[]

[Outputs]
  exodus = false
  [csv]
    type = CSV
    execute_on = FINAL
  []
[]

[Postprocessors]
  [h]
    type = AverageElementSize
    outputs = csv
    execute_on = FINAL
  []
  [L2u]
    type = ElementL2FunctorError
    approximate = vel_x
    exact = exact_u
    outputs = csv
    execute_on = FINAL
  []
  [L2v]
    type = ElementL2FunctorError
    approximate = vel_y
    exact = exact_v
    outputs = csv
    execute_on = FINAL
  []
  [L2p]
    type = ElementL2FunctorError
    approximate = pressure
    exact = 0
    outputs = csv
    execute_on = FINAL
  []
[]

(modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/linear-segregated/2d-simple-channel-enthalpy/porous-baffle-2d-enthalpy.i)

# Steady conjugate heat transfer from a heated solid to fluid in a 0.5-m porous block.
# The solid fraction scales the 20 MW/m^3 source to 10 MW/m^3. Conservation therefore gives
# u_block_2 = 2, T_fluid_out = 300.59523809524, and h_fluid_out = 1262500.
# For 10 upwind cells, T_fluid_block_2 = 300.32738095238. The local balance gives
# T_solid_block_2 - T_fluid_block_2 = 10e6 / 30000 = 333.33333333333,
# hence T_solid_block_2 = 633.66071428571.

mu = 2e-3
rho = 1000
advected_interp_method = 'upwind'
u_in = 2
forchheimer = 0
bf = '0 0'

cp_fluid = 4200
cp_solid = 500
k_fluid = 0.4
k_solid = 0.2
h_vol = 30000
q_solid = 20000000
T_inlet = 300
T_initial = 300

[FVGradientMethods]
  [jump_aware]
    type = FVPressureJumpGreenGaussGradient
  []
  [velocity_jump_aware]
    type = FVVelocityJumpGreenGaussGradient
  []
  [reconstructed]
    type = FVReconstructedPressureGradient
    base_gradient_method = jump_aware
    gradient_relaxation = 1.0
  []
[]

[Mesh]
  [mesh]
    type = CartesianMeshGenerator
    dim = 2
    dx = '0.5 0.5 0.5'
    dy = '0.5'
    ix = '10 10 10'
    iy = '3'
    subdomain_id = '1 2 3'
  []
  [baffle]
    type = SideSetsBetweenSubdomainsGenerator
    input = mesh
    primary_block = '1'
    paired_block = '2'
    new_boundary = 'baffle'
  []
  [baffle2]
    type = SideSetsBetweenSubdomainsGenerator
    input = baffle
    primary_block = '2'
    paired_block = '3'
    new_boundary = 'baffle2'
  []
  [split_top_bottom]
    type = BreakBoundaryOnSubdomainGenerator
    input = baffle2
    boundaries = 'top bottom'
  []
  [delete]
    type = BoundaryDeletionGenerator
    boundary_names = 'top bottom'
    input = split_top_bottom
  []
[]

[FVInterpolationMethods]
  [upwind]
    type = FVAdvectedUpwind
  []
[]

[Problem]
  linear_sys_names = 'u_system v_system pressure_system energy_system solid_energy_system'
  previous_nl_solution_required = true
[]

[UserObjects]
  [pressure_jump]
    type = BernoulliFormLossPressureJump
    boundary = 'baffle baffle2'
    porosity = porosity
    density = ${rho}
    form_loss = ${bf}
    reference_velocity_side = 'lower_porosity higher_porosity'
  []
  [rc]
    type = PorousRhieChowMassFlux
    u = superficial_u
    v = superficial_v
    pressure = pressure
    rho = ${rho}
    porosity = porosity
    p_diffusion_kernel = p_diffusion
    pressure_jump_models = pressure_jump
    pressure_jump_relaxation = 0.1
  []
[]

[Variables]
  [superficial_u]
    type = MooseLinearVariableFVReal
    solver_sys = u_system
    initial_condition = ${u_in}
    gradient_method = velocity_jump_aware
  []
  [superficial_v]
    type = MooseLinearVariableFVReal
    solver_sys = v_system
    initial_condition = 0.0
    gradient_method = velocity_jump_aware
  []
  [pressure]
    type = MooseLinearVariableFVReal
    solver_sys = pressure_system
    initial_condition = 0.0
    gradient_method = jump_aware
  []
  [h_fluid]
    type = MooseLinearVariableFVReal
    solver_sys = energy_system
    initial_condition = '${fparse cp_fluid * T_initial}'
  []
  [T_solid]
    type = MooseLinearVariableFVReal
    solver_sys = solid_energy_system
    initial_condition = ${T_initial}
    block = 2
  []
[]

[LinearFVKernels]
  [u_advection]
    type = LinearPWCNSFVMomentumFlux
    variable = superficial_u
    advected_interp_method_name = ${advected_interp_method}
    mu = ${mu}
    u = superficial_u
    v = superficial_v
    momentum_component = 'x'
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = false
  []
  [v_advection]
    type = LinearPWCNSFVMomentumFlux
    variable = superficial_v
    advected_interp_method_name = ${advected_interp_method}
    mu = ${mu}
    u = superficial_u
    v = superficial_v
    momentum_component = 'y'
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = false
  []
  [u_pressure]
    type = LinearPWCNSFVMomentumPressure
    variable = superficial_u
    momentum_component = 'x'
    porosity = porosity
    pressure = pressure
    gradient_method = reconstructed
  []
  [v_pressure]
    type = LinearPWCNSFVMomentumPressure
    variable = superficial_v
    momentum_component = 'y'
    porosity = porosity
    pressure = pressure
    gradient_method = reconstructed
  []
  [u_friction]
    type = LinearFVMomentumFriction
    variable = superficial_u
    Forchheimer_name = forch
    porosity = porosity
    rho = ${rho}
    u = superficial_u
    v = superficial_v
    momentum_component = 'x'
    block = 2
  []
  [v_friction]
    type = LinearFVMomentumFriction
    variable = superficial_v
    Forchheimer_name = forch
    porosity = porosity
    rho = ${rho}
    u = superficial_u
    v = superficial_v
    momentum_component = 'y'
    block = 2
  []
  [p_diffusion]
    type = LinearFVPressureCorrectionDiffusionJump
    variable = pressure
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = false
  []
  [HbyA_divergence]
    type = LinearFVDivergence
    variable = pressure
    face_flux = HbyA
    force_boundary_execution = true
  []

  [fluid_advection]
    type = LinearFVEnergyAdvection
    variable = h_fluid
    advected_quantity = enthalpy
    advected_interp_method = ${advected_interp_method}
    rhie_chow_user_object = rc
  []
  [fluid_diffusion]
    type = LinearFVDiffusion
    variable = h_fluid
    diffusion_coeff = kappa_fluid
    use_nonorthogonal_correction = false
  []
  [fluid_exchange]
    type = LinearFVEnthalpyVolumetricHeatTransfer
    variable = h_fluid
    h_solid_fluid = h_vol
    cp = cp_fluid
    T_fluid = T_fluid
    T_solid = T_solid
    is_solid = false
    block = 2
  []
  [fluid_source]
    type = LinearFVSource
    variable = h_fluid
    source_density = q_fluid_scaled
  []

  [solid_diffusion]
    type = LinearFVDiffusion
    variable = T_solid
    diffusion_coeff = ${k_solid}
    use_nonorthogonal_correction = false
    block = 2
  []
  [solid_exchange]
    type = LinearFVEnthalpyVolumetricHeatTransfer
    variable = T_solid
    h_solid_fluid = h_vol
    cp = 1.0
    T_fluid = T_fluid
    T_solid = T_solid
    is_solid = true
    block = 2
  []
  [solid_source]
    type = LinearFVSource
    variable = T_solid
    source_density = q_solid_scaled
    block = 2
  []
[]

[LinearFVBCs]
  [left_u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = left
    variable = superficial_u
    functor = ${u_in}
  []
  [outlet_u]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = right
    variable = superficial_u
    use_two_term_expansion = false
  []

  [left_v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = left
    variable = superficial_v
    functor = 0.0
  []
  [outlet_v]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = right
    variable = superficial_v
    use_two_term_expansion = false
  []

  [symmetry-u]
    type = LinearFVVelocitySymmetryBC
    boundary = 'top_to_1 top_to_3 bottom_to_1 bottom_to_3 top_to_2 bottom_to_2'
    variable = superficial_u
    u = superficial_u
    v = superficial_v
    momentum_component = x
  []
  [symmetry-v]
    type = LinearFVVelocitySymmetryBC
    boundary = 'top_to_1 top_to_3 bottom_to_1 bottom_to_3 top_to_2 bottom_to_2'
    variable = superficial_v
    u = superficial_u
    v = superficial_v
    momentum_component = y
  []

  [pressure-extrapolation]
    type = LinearFVPressureFluxBC
    boundary = 'bottom_to_1 bottom_to_2 bottom_to_3 top_to_1 top_to_2 top_to_3'
    variable = pressure
    HbyA_flux = HbyA
    Ainv = Ainv
    rho = ${rho}
    u = superficial_u
    v = superficial_v
  []

  [outlet_p]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = right
    variable = pressure
    functor = 0.0
  []

  [left_T_fluid]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = left
    variable = T_fluid
    functor = ${T_inlet}
  []
  [left_h_fluid]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = left
    variable = h_fluid
    functor = h_from_p_T
  []
  [walls_h_fluid]
    type = LinearFVAdvectionDiffusionFunctorNeumannBC
    boundary = 'top_to_1 top_to_2 top_to_3 bottom_to_1 bottom_to_2 bottom_to_3'
    variable = h_fluid
    functor = 0.0
    diffusion_coeff = kappa_fluid
  []
  [outlet_h_fluid]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = right
    variable = h_fluid
    use_two_term_expansion = false
  []

  [solid_zero_flux]
    type = LinearFVAdvectionDiffusionFunctorNeumannBC
    boundary = 'baffle baffle2 top_to_2 bottom_to_2'
    variable = T_solid
    functor = 0.0
    diffusion_coeff = kappa_solid
  []
[]

[FunctorMaterials]
  [forch]
    type = GenericVectorFunctorMaterial
    prop_names = forch
    prop_values = '${forchheimer} ${forchheimer} ${forchheimer}'
  []
  [porosity]
    type = PiecewiseByBlockFunctorMaterial
    prop_name = porosity
    subdomain_to_prop_value = '1 1.0 2 0.5 3 1.0'
  []
  [fluid_constants]
    type = GenericFunctorMaterial
    prop_names = 'cp_fluid kappa_fluid h_vol q_fluid_scaled'
    prop_values = '${cp_fluid} ${fparse k_fluid/cp_fluid} ${h_vol} 0'
  []
  [solid_constants]
    type = GenericFunctorMaterial
    prop_names = 'cp_solid kappa_solid'
    prop_values = '${cp_solid} ${fparse k_solid/cp_solid}'
  []
  [q_solid]
    type = GenericFunctorMaterial
    prop_names = 'q_solid'
    prop_values = ${q_solid}
    block = 2
  []
  [q_solid_scaled]
    type = ParsedFunctorMaterial
    property_name = q_solid_scaled
    functor_names = 'porosity q_solid'
    expression = '(1 - porosity) * q_solid'
    block = 2
  []
  [fluid_enthalpy_material]
    type = LinearFVEnthalpyFunctorMaterial
    pressure = pressure
    T_fluid = T_fluid
    h = h_fluid
    h_from_p_T_functor = h_from_p_T_functor
    T_from_p_h_functor = T_from_p_h_functor
  []
  [h_from_p_T_functor]
    type = ParsedFunctorMaterial
    property_name = h_from_p_T_functor
    functor_names = 'T_fluid'
    expression = '${cp_fluid} * T_fluid'
  []
  [T_from_p_h_functor]
    type = ParsedFunctorMaterial
    property_name = T_from_p_h_functor
    functor_names = 'h_fluid'
    expression = 'h_fluid / ${cp_fluid}'
  []
[]

[Postprocessors]
  [u_block_2]
    type = ElementAverageValue
    variable = superficial_u
    block = 2
  []
  [T_fluid_out]
    type = SideAverageValue
    variable = T_fluid
    boundary = right
  []
  [T_fluid_block_2]
    type = ElementAverageValue
    variable = T_fluid
    block = 2
  []
  [T_solid_block_2]
    type = ElementAverageValue
    variable = T_solid
    block = 2
  []
  [h_fluid_out]
    type = SideAverageValue
    variable = h_fluid
    boundary = right
  []
[]

[AuxVariables]
  [porosity_aux]
    type = MooseLinearVariableFVReal
  []
  [T_fluid]
    type = MooseLinearVariableFVReal
    initial_condition = ${T_initial}
  []
[]

[AuxKernels]
  [por]
    type = FunctorAux
    variable = porosity_aux
    functor = porosity
    execute_on = 'NONLINEAR'
  []
  [fluid_temperature]
    type = FunctorAux
    variable = T_fluid
    functor = T_from_p_h
    execute_on = 'NONLINEAR'
  []
[]

[VectorPostprocessors]
  # Sample each x-cell once at its center.
  [centerline]
    type = LineValueSampler
    variable = 'superficial_u T_fluid T_solid'
    start_point = '0.0119047619047619 0.25 0'
    end_point = '1.4880952380952381 0.25 0'
    num_points = 30
    sort_by = x
  []
[]

[Executioner]
  type = SIMPLE
  momentum_l_abs_tol = 1e-14
  pressure_l_abs_tol = 1e-14
  energy_l_abs_tol = 1e-14
  solid_energy_l_abs_tol = 1e-14
  momentum_l_tol = 0
  pressure_l_tol = 0
  energy_l_tol = 0
  solid_energy_l_tol = 0
  rhie_chow_user_object = rc
  momentum_systems = 'u_system v_system'
  pressure_system = pressure_system
  energy_system = energy_system
  solid_energy_system = solid_energy_system
  momentum_equation_relaxation = 0.2
  pressure_variable_relaxation = 0.1
  energy_equation_relaxation = 0.8
  num_iterations = 1000
  # These outer tolerances stop once the velocity and zero-pressure fields have stabilized to
  # substantially below the accuracy exercised by the temperature and velocity profile checks.
  pressure_absolute_tolerance = 1e-6
  momentum_absolute_tolerance = 1e-2
  energy_absolute_tolerance = 1e-8
  solid_energy_absolute_tolerance = 1e-8
  momentum_petsc_options_iname = '-pc_type -pc_hypre_type'
  momentum_petsc_options_value = 'hypre boomeramg'
  pressure_petsc_options_iname = '-pc_type -pc_hypre_type'
  pressure_petsc_options_value = 'hypre boomeramg'
  energy_petsc_options_iname = '-pc_type -pc_hypre_type'
  energy_petsc_options_value = 'hypre boomeramg'
  solid_energy_petsc_options_iname = '-pc_type -pc_hypre_type'
  solid_energy_petsc_options_value = 'hypre boomeramg'
[]

[Outputs]
  csv = true
  execute_on = FINAL
[]

(modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/linear-segregated/errors/pressure-jump-errors.i)

# Compact fixture for constructor and setup errors on pressure-jump objects.

mu = 0.0
rho = 2.0
advected_interp_method = 'average'
inlet_u = 0.1

[FVGradientMethods]
  [jump_aware]
    type = FVPressureJumpGreenGaussGradient
  []
  [velocity_jump_aware]
    type = FVVelocityJumpGreenGaussGradient
  []
  [reconstructed]
    type = FVReconstructedPressureGradient
    base_gradient_method = jump_aware
    gradient_relaxation = 1.0
  []
[]

[Mesh]
  [mesh]
    type = CartesianMeshGenerator
    dim = 1
    dx = '0.5 0.5 0.5'
    ix = '2 2 2'
    subdomain_id = '1 2 3'
  []
  [baffle]
    type = SideSetsBetweenSubdomainsGenerator
    input = mesh
    primary_block = '1'
    paired_block = '2'
    new_boundary = 'baffle'
  []
  [baffle2]
    type = SideSetsBetweenSubdomainsGenerator
    input = baffle
    primary_block = '2'
    paired_block = '3'
    new_boundary = 'baffle2'
  []
[]

[Problem]
  linear_sys_names = 'u_system pressure_system'
  previous_nl_solution_required = true
[]

[UserObjects]
  active = 'pressure_jump rc'
  [pressure_jump]
    type = BernoulliFormLossPressureJump
    boundary = 'baffle baffle2'
    porosity = porosity
    density = ${rho}
  []
  [duplicate_jump]
    type = BernoulliFormLossPressureJump
    boundary = 'baffle'
    porosity = porosity
    density = ${rho}
  []
  [rc]
    type = PorousRhieChowMassFlux
    u = superficial_u
    pressure = pressure
    rho = ${rho}
    porosity = porosity
    p_diffusion_kernel = p_diffusion
    pressure_jump_models = pressure_jump
    pressure_jump_relaxation = 0.1
  []
[]

[Variables]
  [superficial_u]
    type = MooseLinearVariableFVReal
    solver_sys = u_system
    initial_condition = ${inlet_u}
    gradient_method = velocity_jump_aware
  []
  [pressure]
    type = MooseLinearVariableFVReal
    solver_sys = pressure_system
    initial_condition = 0.0
    gradient_method = jump_aware
  []
[]

[FVInterpolationMethods]
  [average]
    type = FVGeometricAverage
  []
[]

[LinearFVKernels]
  active = 'u_advection u_pressure p_diffusion HbyA_divergence'
  [u_advection]
    type = LinearPWCNSFVMomentumFlux
    variable = superficial_u
    advected_interp_method_name = ${advected_interp_method}
    mu = ${mu}
    u = superficial_u
    momentum_component = 'x'
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = false
  []
  [u_pressure]
    type = LinearPWCNSFVMomentumPressure
    variable = superficial_u
    momentum_component = 'x'
    porosity = porosity
    pressure = pressure
    gradient_method = reconstructed
  []
  [p_diffusion]
    type = LinearFVPressureCorrectionDiffusionJump
    variable = pressure
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = false
  []
  [ordinary_p_diffusion]
    type = LinearFVPressureCorrectionDiffusion
    variable = pressure
    diffusion_tensor = Ainv
  []
  [HbyA_divergence]
    type = LinearFVDivergence
    variable = pressure
    face_flux = HbyA
    force_boundary_execution = true
  []
[]

[LinearFVBCs]
  [inlet_u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = left
    variable = superficial_u
    functor = ${inlet_u}
  []
  [outlet_u]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = right
    variable = superficial_u
    use_two_term_expansion = false
  []
  [outlet_p]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = right
    variable = pressure
    functor = 0.0
  []
[]

[FunctorMaterials]
  [porosity]
    type = PiecewiseByBlockFunctorMaterial
    prop_name = porosity
    subdomain_to_prop_value = '1 0.5 2 1.0 3 0.5'
  []
[]

[Executioner]
  type = SIMPLE
  rhie_chow_user_object = rc
  momentum_systems = 'u_system'
  pressure_system = pressure_system
  num_iterations = 1
  continue_on_max_its = true
[]

[Outputs]
  csv = false
  exodus = false
[]

(modules/navier_stokes/test/tests/finite_volume/pins/mms/linear-segregated/interface-reconstruction/interface-reconstruction.i)

# Manufactured pressure-baffle problem whose superficial velocity is linear on each side of the
# baffle. In coordinates normal (r) and tangent (t) to the interface:
#  - rho * U_r = mass_flux, so mass is conserved;
#  - U_t = eps_s * tangential_velocity + tangential_slope_s * r, so the interstitial tangential
#    velocity U_t / eps is continuous at r = 0 while the one-sided slopes may differ.
# The forcing balances advection and reaction exactly, so the pressure is constant on each side.
# Its jump is the BernoulliFormLossPressureJump with side velocities mass_flux / (rho eps_s):
#   p_right - p_left = mass_flux^2 / (2 rho) * (1 / eps_left^2 - 1 / eps_right^2).
# With viscosity, the baffle transmits only a shear traction, and only from a clear side (unit
# porosity). Between two clear sides it is continuous, so the slopes must be equal. Two porous
# sides transmit none, so viscous cases with porosity below one on both sides need zero slopes. A
# clear side shears against the superficial velocity of a porous side at the interface, so with
# the clear side on the left the traction balance on both rows requires equal slopes and a
# continuous superficial tangential velocity, eps_left * tangential_velocity =
# eps_right * tangential_velocity, which with a porosity jump means tangential_velocity = 0.
# The normal velocity is constant and carries no viscous stress. A consistent scheme reproduces
# this piecewise-linear solution to round-off on uniform meshes.

rho = 1
eps_left = 1
eps_right = 1
mass_flux = 0.02
mu = 0
reaction = 10
tangential_velocity = 0.1
tangential_slope_left = 0.2
tangential_slope_right = 0.2
interface = 0.5
skew = 0
nx = 8
ny = 8
pressure_gradient_method = reconstructed
nonorthogonal = false

[Problem]
  linear_sys_names = 'u_system v_system pressure_system'
  previous_nl_solution_required = true
[]

[Mesh]
  [generated]
    type = GeneratedMeshGenerator
    dim = 2
    nx = ${nx}
    ny = ${ny}
  []
  [right]
    type = ParsedSubdomainMeshGenerator
    input = generated
    combinatorial_geometry = 'x > ${interface}'
    block_id = 2
    block_name = right
  []
  [left]
    type = ParsedSubdomainMeshGenerator
    input = right
    combinatorial_geometry = 'x < ${interface}'
    block_id = 1
    block_name = left
  []
  [skew]
    type = ParsedNodeTransformGenerator
    input = left
    x_function = 'x + ${skew} * y'
  []
  [baffle]
    type = SideSetsBetweenSubdomainsGenerator
    input = skew
    primary_block = left
    paired_block = right
    new_boundary = baffle
  []
[]

[FVGradientMethods]
  [jump_aware]
    type = FVPressureJumpGreenGaussGradient
  []
  [velocity_jump_aware]
    type = FVVelocityJumpGreenGaussGradient
  []
  [reconstructed]
    type = FVReconstructedPressureGradient
    base_gradient_method = jump_aware
    gradient_relaxation = 1
  []
[]

[FVInterpolationMethods]
  [average]
    type = FVGeometricAverage
  []
[]

[UserObjects]
  [pressure_jump]
    type = BernoulliFormLossPressureJump
    boundary = baffle
    porosity = porosity
    density = density
  []
  [rc]
    type = PorousRhieChowMassFlux
    u = superficial_u
    v = superficial_v
    pressure = pressure
    rho = density
    porosity = porosity
    p_diffusion_kernel = p_diffusion
    pressure_jump_models = pressure_jump
    pressure_projection_method = consistent
  []
[]

[Variables]
  [superficial_u]
    type = MooseLinearVariableFVReal
    solver_sys = u_system
    initial_condition = ${mass_flux}
    gradient_method = velocity_jump_aware
  []
  [superficial_v]
    type = MooseLinearVariableFVReal
    solver_sys = v_system
    initial_condition = ${tangential_velocity}
    gradient_method = velocity_jump_aware
  []
  [pressure]
    type = MooseLinearVariableFVReal
    solver_sys = pressure_system
    initial_condition = 0
    gradient_method = jump_aware
  []
[]

[LinearFVKernels]
  [u_advection]
    type = LinearPWCNSFVMomentumFlux
    variable = superficial_u
    advected_interp_method_name = average
    mu = viscosity
    u = superficial_u
    v = superficial_v
    momentum_component = x
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = ${nonorthogonal}
  []
  [v_advection]
    type = LinearPWCNSFVMomentumFlux
    variable = superficial_v
    advected_interp_method_name = average
    mu = viscosity
    u = superficial_u
    v = superficial_v
    momentum_component = y
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = ${nonorthogonal}
  []
  [u_pressure]
    type = LinearPWCNSFVMomentumPressure
    variable = superficial_u
    pressure = pressure
    porosity = porosity
    momentum_component = x
    gradient_method = ${pressure_gradient_method}
  []
  [v_pressure]
    type = LinearPWCNSFVMomentumPressure
    variable = superficial_v
    pressure = pressure
    porosity = porosity
    momentum_component = y
    gradient_method = ${pressure_gradient_method}
  []
  [u_source]
    type = LinearFVSource
    variable = superficial_u
    source_density = forcing_u
  []
  [u_reaction]
    type = LinearFVReaction
    variable = superficial_u
    coeff = ${reaction}
  []
  [v_source]
    type = LinearFVSource
    variable = superficial_v
    source_density = forcing_v
  []
  [v_reaction]
    type = LinearFVReaction
    variable = superficial_v
    coeff = ${reaction}
  []
  [p_diffusion]
    type = LinearFVPressureCorrectionDiffusionJump
    variable = pressure
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = ${nonorthogonal}
  []
  [HbyA_divergence]
    type = LinearFVDivergence
    variable = pressure
    face_flux = HbyA
    force_boundary_execution = true
  []
[]

[LinearFVBCs]
  [velocity_x]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left right top bottom'
    variable = superficial_u
    functor = exact_u
  []
  [velocity_y]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left right top bottom'
    variable = superficial_v
    functor = exact_v
  []
  [pressure_flux]
    type = LinearFVPressureFluxBC
    boundary = 'left right top bottom'
    variable = pressure
    HbyA_flux = HbyA
    Ainv = Ainv
    rho = density
    u = superficial_u
    v = superficial_v
  []
[]

[FunctorMaterials]
  [density]
    type = ParsedFunctorMaterial
    property_name = density
    expression = '${rho}'
  []
  [porosity]
    type = ParsedFunctorMaterial
    property_name = porosity
    expression = 'if(x-${skew}*y<${interface},${eps_left},${eps_right})'
  []
  [viscosity]
    type = ParsedFunctorMaterial
    property_name = viscosity
    expression = '${mu}'
  []
  [pressure_error]
    type = ParsedFunctorMaterial
    property_name = pressure_error
    functor_names = 'pressure exact_pressure'
    expression = 'pressure-exact_pressure'
  []
[]

[AuxVariables]
  [reconstructed_grad_p]
    family = MONOMIAL_VEC
    order = CONSTANT
  []
[]

[AuxKernels]
  [reconstructed_grad_p]
    type = LinearFVGradientAux
    variable = reconstructed_grad_p
    linear_fv_variable = pressure
    gradient_method = reconstructed
    execute_on = FINAL
  []
[]

[Functions]
  [normal_coordinate]
    type = ParsedFunction
    expression = '(x-skew*y-interface)/sqrt(1+skew^2)'
    symbol_names = 'skew interface'
    symbol_values = '${skew} ${interface}'
  []
  [exact_normal_velocity]
    type = ParsedFunction
    expression = '${mass_flux}/${rho}'
  []
  [exact_tangential_velocity]
    type = ParsedFunction
    expression = 'velocity*porosity+if(x-skew*y<interface,slope_left,slope_right)*r'
    symbol_names = 'skew interface velocity porosity slope_left slope_right r'
    symbol_values = '${skew} ${interface} ${tangential_velocity} exact_porosity ${tangential_slope_left} ${tangential_slope_right} normal_coordinate'
  []
  [exact_pressure]
    type = ParsedFunction
    expression = 'if(x-skew*y<interface,0,mass_flux^2/(2*rho)*(1/eps_left^2-1/eps_right^2))'
    symbol_names = 'skew interface mass_flux rho eps_left eps_right'
    symbol_values = '${skew} ${interface} ${mass_flux} ${rho} ${eps_left} ${eps_right}'
  []
  [exact_u]
    type = ParsedFunction
    expression = '(normal+skew*tangential)/sqrt(1+skew^2)'
    symbol_names = 'skew normal tangential'
    symbol_values = '${skew} exact_normal_velocity exact_tangential_velocity'
  []
  [exact_v]
    type = ParsedFunction
    expression = '(-skew*normal+tangential)/sqrt(1+skew^2)'
    symbol_names = 'skew normal tangential'
    symbol_values = '${skew} exact_normal_velocity exact_tangential_velocity'
  []
  [forcing_u]
    type = ParsedFunction
    expression = 'skew*mass_flux*slope/(porosity*sqrt(1+skew^2))+reaction*u'
    symbol_names = 'skew interface mass_flux slope porosity reaction u'
    symbol_values = '${skew} ${interface} ${mass_flux} tangential_slope exact_porosity ${reaction} exact_u'
  []
  [forcing_v]
    type = ParsedFunction
    expression = 'mass_flux*slope/(porosity*sqrt(1+skew^2))+reaction*v'
    symbol_names = 'skew interface mass_flux slope porosity reaction v'
    symbol_values = '${skew} ${interface} ${mass_flux} tangential_slope exact_porosity ${reaction} exact_v'
  []
  [tangential_slope]
    type = ParsedFunction
    expression = 'if(x-skew*y<interface,slope_left,slope_right)'
    symbol_names = 'skew interface slope_left slope_right'
    symbol_values = '${skew} ${interface} ${tangential_slope_left} ${tangential_slope_right}'
  []
  [exact_porosity]
    type = ParsedFunction
    expression = 'if(x-skew*y<interface,eps_left,eps_right)'
    symbol_names = 'skew interface eps_left eps_right'
    symbol_values = '${skew} ${interface} ${eps_left} ${eps_right}'
  []
[]

[Executioner]
  type = SIMPLE
  momentum_l_abs_tol = 1e-12
  pressure_l_abs_tol = 1e-12
  momentum_l_tol = 0
  pressure_l_tol = 0
  momentum_l_max_its = 100
  pressure_l_max_its = 100
  rhie_chow_user_object = rc
  momentum_systems = 'u_system v_system'
  pressure_system = pressure_system
  momentum_equation_relaxation = 0.7
  pressure_variable_relaxation = 0.5
  # A solution that is exact to round-off stalls the normalized residuals near 1e-6 instead of
  # meeting the absolute tolerances, so run a fixed number of iterations instead.
  num_iterations = 200
  continue_on_max_its = true
  pressure_absolute_tolerance = 1e-11
  momentum_absolute_tolerance = 1e-11
  momentum_petsc_options_iname = '-pc_type -pc_hypre_type'
  momentum_petsc_options_value = 'hypre boomeramg'
  pressure_petsc_options_iname = '-pc_type -pc_hypre_type'
  pressure_petsc_options_value = 'hypre boomeramg'
  print_fields = false

  pin_pressure = true
  pressure_pin_value = 0
  pressure_pin_point = '0.25 0.5 0'
[]

[Postprocessors]
  [h]
    type = AverageElementSize
    execute_on = FINAL
  []
  [L2u]
    type = ElementL2FunctorError
    approximate = superficial_u
    exact = exact_u
    execute_on = FINAL
  []
  [L2v]
    type = ElementL2FunctorError
    approximate = superficial_v
    exact = exact_v
    execute_on = FINAL
  []
  [L2p]
    type = ElementL2FunctorError
    approximate = pressure
    exact = exact_pressure
    execute_on = FINAL
  []
  [L2gradp]
    type = ElementVectorL2Error
    variable = reconstructed_grad_p
    function_x = 0
    function_y = 0
    execute_on = FINAL
  []
  [boundary_mass_balance]
    type = VolumetricFlowRate
    boundary = 'left right top bottom'
    vel_x = superficial_u
    vel_y = superficial_v
    advected_quantity = density
    rhie_chow_user_object = rc
    execute_on = FINAL
  []
  [max_pressure_error]
    type = ElementExtremeFunctorValue
    functor = pressure_error
    value_type = max_abs
    execute_on = FINAL
  []
[]

[Outputs]
  exodus = false
  [csv]
    type = CSV
    execute_on = FINAL
  []
[]

(modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/linear-segregated/2d-diagonal-baffle/porous-baffle-straight-diagonal.i)

# Straight-channel hydraulic-only porous-baffle showcase.
# Flow enters from the left, crosses a diagonal porous baffle in a straight section,
# and exits to the right through a channel with slip walls.
# The solution should conserve the imposed unit superficial flow, produce positive pressure
# drops across the clean-to-porous, diagonal, and porous-to-clean interfaces, and keep the
# transverse velocity smooth on both sides of the diagonal. The skew geometry and porous drag
# preclude a useful closed-form pressure distribution, so the profiles provide the reference.

inlet_velocity = 1.0
rho = 998.2
mu = 2e-3
mu_zone_2 = ${mu}
eps_zone_1 = 0.5
eps_zone_2 = 0.7
forch_zone_1 = 10
forch_zone_2 = 10
entry_form_loss = 20
corner_form_loss = 20
exit_form_loss = 20
advected_interp_method = 'upwind'

[FVGradientMethods]
  [jump_aware]
    type = FVPressureJumpGreenGaussGradient
  []
  [velocity_jump_aware]
    type = FVVelocityJumpGreenGaussGradient
  []
  [reconstructed]
    type = FVReconstructedPressureGradient
    base_gradient_method = jump_aware
    gradient_relaxation = 1.0
  []
[]

[Mesh]
  [clean_inlet_gmg]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 7
    ny = 7
    xmin = 0
    xmax = 1
    ymin = 0
    ymax = 1
  []
  [clean_inlet_id]
    type = SubdomainIDGenerator
    input = clean_inlet_gmg
    subdomain_id = 1
  []
  [clean_inlet_block]
    type = RenameBlockGenerator
    input = clean_inlet_id
    old_block = '1'
    new_block = 'clean_inlet'
  []
  [clean_inlet_boundary]
    type = RenameBoundaryGenerator
    input = clean_inlet_block
    old_boundary = 'bottom right top left'
    new_boundary = 'bottom_clean clean_to_porous_tmp top_clean inlet'
  []

  [porous_1_patch]
    type = TransfiniteMeshGenerator
    corners = '1 0 0
               1.7 0 0
               2.3 1 0
               1 1 0'
    nx = 8
    ny = 8
    bottom_type = LINE
    top_type = LINE
    right_type = LINE
    left_type = LINE
  []
  [porous_1_id]
    type = SubdomainIDGenerator
    input = porous_1_patch
    subdomain_id = 2
  []
  [porous_1_block]
    type = RenameBlockGenerator
    input = porous_1_id
    old_block = '2'
    new_block = 'porous_zone_1'
  []
  [porous_1_boundary]
    type = RenameBoundaryGenerator
    input = porous_1_block
    old_boundary = 'bottom right top left'
    new_boundary = 'bottom_porous diagonal_p_tmp top_porous clean_to_porous_tmp'
  []

  [porous_2_patch]
    type = TransfiniteMeshGenerator
    corners = '2.3 1 0
               1.7 0 0
               3 0 0
               3 1 0'
    nx = 8
    ny = 8
    bottom_type = LINE
    top_type = LINE
    right_type = LINE
    left_type = LINE
  []
  [porous_2_id]
    type = SubdomainIDGenerator
    input = porous_2_patch
    subdomain_id = 3
  []
  [porous_2_block]
    type = RenameBlockGenerator
    input = porous_2_id
    old_block = '3'
    new_block = 'porous_zone_2'
  []
  [porous_2_boundary]
    type = RenameBoundaryGenerator
    input = porous_2_block
    old_boundary = 'bottom right top left'
    new_boundary = 'diagonal_p_tmp bottom_porous porous_to_clean_tmp top_porous'
  []

  [clean_outlet_gmg]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 14
    ny = 7
    xmin = 3
    xmax = 5
    ymin = 0
    ymax = 1
  []
  [clean_outlet_id]
    type = SubdomainIDGenerator
    input = clean_outlet_gmg
    subdomain_id = 4
  []
  [clean_outlet_block]
    type = RenameBlockGenerator
    input = clean_outlet_id
    old_block = '4'
    new_block = 'clean_outlet'
  []
  [clean_outlet_boundary]
    type = RenameBoundaryGenerator
    input = clean_outlet_block
    old_boundary = 'bottom right top left'
    new_boundary = 'bottom_clean outlet top_clean porous_to_clean_tmp'
  []

  [stitched_mesh]
    type = StitchMeshGenerator
    inputs = 'clean_inlet_boundary
              porous_1_boundary
              porous_2_boundary
              clean_outlet_boundary'
    stitch_boundaries_pairs = 'clean_to_porous_tmp clean_to_porous_tmp;
                               diagonal_p_tmp diagonal_p_tmp;
                               porous_to_clean_tmp porous_to_clean_tmp'
    clear_stitched_boundary_ids = true
    merge_boundaries_with_same_name = true
  []

  [clean_to_porous]
    type = SideSetsBetweenSubdomainsGenerator
    input = stitched_mesh
    primary_block = 'clean_inlet'
    paired_block = 'porous_zone_1'
    new_boundary = 'clean_to_porous'
  []
  [diagonal_baffle]
    type = SideSetsBetweenSubdomainsGenerator
    input = clean_to_porous
    primary_block = 'porous_zone_1'
    paired_block = 'porous_zone_2'
    new_boundary = 'diagonal_baffle'
  []
  [porous_to_clean]
    type = SideSetsBetweenSubdomainsGenerator
    input = diagonal_baffle
    primary_block = 'porous_zone_2'
    paired_block = 'clean_outlet'
    new_boundary = 'porous_to_clean'
  []
[]

[Problem]
  linear_sys_names = 'u_system v_system pressure_system'
  previous_nl_solution_required = true
[]

[UserObjects]
  [pressure_jump]
    type = BernoulliFormLossPressureJump
    boundary = 'clean_to_porous diagonal_baffle porous_to_clean'
    porosity = porosity
    density = ${rho}
    form_loss = '${entry_form_loss} ${corner_form_loss} ${exit_form_loss}'
    reference_velocity_side = 'lower_porosity lower_porosity lower_porosity'
  []
  [rc]
    type = PorousRhieChowMassFlux
    u = superficial_u
    v = superficial_v
    pressure = pressure
    rho = ${rho}
    porosity = porosity
    p_diffusion_kernel = p_diffusion
    pressure_projection_method = consistent
    pressure_jump_models = pressure_jump
    pressure_jump_relaxation = 0.1
  []
[]

[FVInterpolationMethods]
  [upwind]
    type = FVAdvectedUpwind
  []
[]

[Variables]
  [superficial_u]
    type = MooseLinearVariableFVReal
    solver_sys = u_system
    initial_condition = ${inlet_velocity}
    gradient_method = velocity_jump_aware
  []
  [superficial_v]
    type = MooseLinearVariableFVReal
    solver_sys = v_system
    initial_condition = 0.0
    gradient_method = velocity_jump_aware
  []
  [pressure]
    type = MooseLinearVariableFVReal
    solver_sys = pressure_system
    initial_condition = 0.0
    gradient_method = jump_aware
  []
[]

[LinearFVKernels]
  [u_advection]
    type = LinearPWCNSFVMomentumFlux
    variable = superficial_u
    advected_interp_method_name = ${advected_interp_method}
    mu = dynamic_viscosity
    u = superficial_u
    v = superficial_v
    momentum_component = 'x'
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = false
  []
  [v_advection]
    type = LinearPWCNSFVMomentumFlux
    variable = superficial_v
    advected_interp_method_name = ${advected_interp_method}
    mu = dynamic_viscosity
    u = superficial_u
    v = superficial_v
    momentum_component = 'y'
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = false
  []
  [u_pressure]
    type = LinearPWCNSFVMomentumPressure
    variable = superficial_u
    momentum_component = 'x'
    porosity = porosity
    pressure = pressure
    gradient_method = reconstructed
  []
  [v_pressure]
    type = LinearPWCNSFVMomentumPressure
    variable = superficial_v
    momentum_component = 'y'
    porosity = porosity
    pressure = pressure
    gradient_method = reconstructed
  []
  [u_friction]
    type = LinearFVMomentumFriction
    variable = superficial_u
    Forchheimer_name = forch
    porosity = porosity
    rho = ${rho}
    u = superficial_u
    v = superficial_v
    momentum_component = 'x'
    block = 'porous_zone_1 porous_zone_2'
  []
  [v_friction]
    type = LinearFVMomentumFriction
    variable = superficial_v
    Forchheimer_name = forch
    porosity = porosity
    rho = ${rho}
    u = superficial_u
    v = superficial_v
    momentum_component = 'y'
    block = 'porous_zone_1 porous_zone_2'
  []
  [p_diffusion]
    type = LinearFVPressureCorrectionDiffusionJump
    variable = pressure
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = true
  []
  [HbyA_divergence]
    type = LinearFVDivergence
    variable = pressure
    face_flux = HbyA
    force_boundary_execution = true
  []
[]

[LinearFVBCs]
  [inlet_u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = inlet
    variable = superficial_u
    functor = ${inlet_velocity}
  []
  [inlet_v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = inlet
    variable = superficial_v
    functor = 0.0
  []

  [pressure-extrapolation]
    type = LinearFVExtrapolatedPressureBC
    boundary = 'inlet'
    variable = pressure
    use_two_term_expansion = true
  []

  [outlet_u]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = outlet
    variable = superficial_u
    use_two_term_expansion = true
  []
  [outlet_v]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = outlet
    variable = superficial_v
    use_two_term_expansion = true
  []
  [outlet_p]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = outlet
    variable = pressure
    functor = 0.0
  []

  [slip_u]
    type = LinearFVVelocitySymmetryBC
    boundary = 'bottom_clean top_clean bottom_porous top_porous'
    variable = superficial_u
    u = superficial_u
    v = superficial_v
    momentum_component = x
  []
  [slip_v]
    type = LinearFVVelocitySymmetryBC
    boundary = 'bottom_clean top_clean bottom_porous top_porous'
    variable = superficial_v
    u = superficial_u
    v = superficial_v
    momentum_component = y
  []
  [pressure-split]
    type = LinearFVPressureFluxBC
    boundary = 'bottom_clean top_clean bottom_porous top_porous'
    variable = pressure
    HbyA_flux = HbyA
    Ainv = Ainv
    rho = ${rho}
    u = superficial_u
    v = superficial_v
  []
[]

[FunctorMaterials]
  [viscosity_clean_and_zone_1]
    type = GenericFunctorMaterial
    prop_names = 'dynamic_viscosity'
    prop_values = '${mu}'
    block = 'clean_inlet porous_zone_1 clean_outlet'
  []
  [viscosity_zone_2]
    type = GenericFunctorMaterial
    prop_names = 'dynamic_viscosity'
    prop_values = '${mu_zone_2}'
    block = 'porous_zone_2'
  []
  [porosity_clean]
    type = GenericFunctorMaterial
    prop_names = 'porosity'
    prop_values = '1.0'
    block = 'clean_inlet clean_outlet'
  []
  [porosity_zone_1]
    type = GenericFunctorMaterial
    prop_names = 'porosity'
    prop_values = '${eps_zone_1}'
    block = 'porous_zone_1'
  []
  [porosity_zone_2]
    type = GenericFunctorMaterial
    prop_names = 'porosity'
    prop_values = '${eps_zone_2}'
    block = 'porous_zone_2'
  []
  [interstitial_tangential_velocity]
    # The diagonal baffle runs from (1.7, 0) to (2.3, 1), so its unit tangent is (0.6, 1) / |(0.6, 1)|.
    type = ParsedFunctorMaterial
    property_name = interstitial_tangential_velocity_functor
    functor_names = 'superficial_u superficial_v porosity'
    expression = '(0.6 * superficial_u + superficial_v) / (sqrt(1.36) * porosity)'
  []
  [forch_zone_1]
    type = GenericVectorFunctorMaterial
    prop_names = 'forch'
    prop_values = '${forch_zone_1} ${fparse forch_zone_1} ${forch_zone_1}'
    block = 'porous_zone_1 clean_outlet'
  []
  [forch_zone_2]
    type = GenericVectorFunctorMaterial
    prop_names = 'forch'
    prop_values = '${fparse forch_zone_2} ${forch_zone_2} ${forch_zone_2}'
    block = 'porous_zone_2 clean_inlet'
  []
[]

[Postprocessors]
  [p_inlet]
    type = SideAverageValue
    variable = pressure
    boundary = inlet
  []
  [p_outlet]
    type = SideAverageValue
    variable = pressure
    boundary = outlet
  []
  [p_total_drop]
    type = ParsedPostprocessor
    expression = 'p_inlet - p_outlet'
    pp_names = 'p_inlet p_outlet'
  []
  [p_clean_in]
    type = ElementAverageValue
    variable = pressure
    block = 'clean_inlet'
  []
  [p_porous_1]
    type = ElementAverageValue
    variable = pressure
    block = 'porous_zone_1'
  []
  [p_porous_2]
    type = ElementAverageValue
    variable = pressure
    block = 'porous_zone_2'
  []
  [p_clean_out]
    type = ElementAverageValue
    variable = pressure
    block = 'clean_outlet'
  []
  [p_entry_drop]
    type = ParsedPostprocessor
    expression = 'p_clean_in - p_porous_1'
    pp_names = 'p_clean_in p_porous_1'
  []
  [p_baffle_drop]
    type = ParsedPostprocessor
    expression = 'p_porous_1 - p_porous_2'
    pp_names = 'p_porous_1 p_porous_2'
  []
  [p_exit_drop]
    type = ParsedPostprocessor
    expression = 'p_porous_2 - p_clean_out'
    pp_names = 'p_porous_2 p_clean_out'
  []
[]

[VectorPostprocessors]
  # The cells intersected by this horizontal line have uniform width; sample each once.
  [horizontal_profile]
    type = LineValueSampler
    variable = 'pressure superficial_u'
    warn_discontinuous_face_values = false
    start_point = '0.0714285714285714 0.5 0'
    end_point = '3.9285714285714286 0.5 0'
    num_points = 28
    sort_by = x
  []

  # Sample the two cell-center rows adjacent and parallel to the diagonal baffle.
  [v_y_left_of_interface]
    type = LineValueSampler
    variable = superficial_v
    warn_discontinuous_face_values = false
    start_point = '1.689795918367347 0.0714285714285714 0'
    end_point = '2.167346938775510 0.9285714285714286 0'
    num_points = 7
    sort_by = y
  []
  [v_y_right_of_interface]
    type = LineValueSampler
    variable = superficial_v
    warn_discontinuous_face_values = false
    start_point = '1.832653061224490 0.0714285714285714 0'
    end_point = '2.310204081632653 0.9285714285714286 0'
    num_points = 7
    sort_by = y
  []
  [ut_left_of_interface]
    type = LineValueSampler
    variable = interstitial_tangential_velocity
    warn_discontinuous_face_values = false
    start_point = '1.689795918367347 0.0714285714285714 0'
    end_point = '2.167346938775510 0.9285714285714286 0'
    num_points = 7
    sort_by = y
  []
  [ut_right_of_interface]
    type = LineValueSampler
    variable = interstitial_tangential_velocity
    warn_discontinuous_face_values = false
    start_point = '1.832653061224490 0.0714285714285714 0'
    end_point = '2.310204081632653 0.9285714285714286 0'
    num_points = 7
    sort_by = y
  []

  [skew_pressure_profile]
    # This skew line intersects cells over unequal distances, so list one interior point per cell.
    type = PointValueSampler
    variable = pressure
    warn_discontinuous_face_values = false
    points = '2.0282 0.9437 0
              2.0639 0.8722 0
              2.1039 0.7922 0
              2.1396 0.7208 0
              2.1700 0.6600 0
              2.2057 0.5885 0
              2.2418 0.5165 0
              2.2775 0.4451 0
              2.3208 0.3584 0
              2.3565 0.2869 0
              2.3929 0.2143 0
              2.4454 0.1092 0
              2.4811 0.0377 0'
    sort_by = x
  []
[]

[AuxVariables]
  [porosity_aux]
    type = MooseLinearVariableFVReal
  []
  [interstitial_tangential_velocity]
    type = MooseLinearVariableFVReal
  []
[]

[AuxKernels]
  [porosity_display]
    type = FunctorAux
    variable = porosity_aux
    functor = porosity
    execute_on = timestep_end
  []
  [interstitial_tangential_velocity]
    type = FunctorAux
    variable = interstitial_tangential_velocity
    functor = interstitial_tangential_velocity_functor
    execute_on = timestep_end
  []
[]

[Executioner]
  type = SIMPLE
  momentum_l_abs_tol = 1e-14
  pressure_l_abs_tol = 1e-14
  momentum_l_tol = 0
  pressure_l_tol = 0
  rhie_chow_user_object = rc
  momentum_systems = 'u_system v_system'
  pressure_system = pressure_system
  momentum_equation_relaxation = 0.5
  pressure_variable_relaxation = 0.2
  num_iterations = 1000
  pressure_absolute_tolerance = 1e-8
  momentum_absolute_tolerance = 1e-8
  momentum_petsc_options_iname = '-pc_type -pc_hypre_type'
  momentum_petsc_options_value = 'hypre boomeramg'
  pressure_petsc_options_iname = '-pc_type -pc_hypre_type'
  pressure_petsc_options_value = 'hypre boomeramg'
[]

[Outputs]
  file_base = porous-baffle-straight-diagonal_out
  csv = true
  execute_on = timestep_end
[]

(modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/linear-segregated/2d-variable-density/porous-baffle-2d-variable-density-bernoulli.i)

# Manufactured variable-density porous-baffle solution. The exact fields are
# rho = 1 - 0.05 x, rho * superficial_u = 1, porosity = 1 for x < 5 and 0.5
# otherwise, and p_right = 0. The exact pressure solution gives p_left = delta_p =
# 4.99973888889. Its block averages are p_block_1 = 4.84902726575 and
# p_block_2 = 1.51242493694, so p_block_jump = 3.33660232881. The exact velocity
# averages are u_block_1 = 1.15072828981 and u_block_2 = 1.62186043243.
# L2_rho and L2_porosity should be zero; the velocity and pressure errors measure
# discretization of their nonconstant exact profiles.

length_left = 5.0
length_right = 5.0
height = 1.0

nx_left = 30
nx_right = 30
ny = 5

epsilon_left = 1.0
epsilon_right = 0.5
u_in = 1.0
p_out = 0.0
mu = 1e-3

density_factor = 1e-3
density_slope_exact = '${fparse -0.05e3 * density_factor}'
rho_in_exact = '${fparse 1e3 * density_factor}'
mass_flux_exact = '${fparse rho_in_exact * u_in}'
rho_baffle_exact = '${fparse (1e3 - 0.05e3 * length_left) * density_factor}'
rho_out_exact = '${fparse (1e3 - 0.05e3 * (length_left + length_right)) * density_factor}'
u_prime_baffle_exact = '${fparse -mass_flux_exact * density_slope_exact / (rho_baffle_exact * rho_baffle_exact)}'
u_prime_out_exact = '${fparse -mass_flux_exact * density_slope_exact / (rho_out_exact * rho_out_exact)}'
transport_baffle_left_exact = '${fparse mass_flux_exact * mass_flux_exact / (rho_baffle_exact * epsilon_left * epsilon_left) - mu * u_prime_baffle_exact / epsilon_left}'
transport_baffle_right_exact = '${fparse mass_flux_exact * mass_flux_exact / (rho_baffle_exact * epsilon_right * epsilon_right) - mu * u_prime_baffle_exact / epsilon_right}'
transport_out_right_exact = '${fparse mass_flux_exact * mass_flux_exact / (rho_out_exact * epsilon_right * epsilon_right) - mu * u_prime_out_exact / epsilon_right}'
bernoulli_jump_exact = '${fparse 0.5 * mass_flux_exact * mass_flux_exact * (1 / (rho_baffle_exact * epsilon_right * epsilon_right) - 1 / (rho_baffle_exact * epsilon_left * epsilon_left))}'
pressure_right_constant_exact = '${fparse p_out + transport_out_right_exact}'
pressure_left_constant_exact = '${fparse pressure_right_constant_exact - transport_baffle_right_exact + bernoulli_jump_exact + transport_baffle_left_exact}'

advected_interp_method = 'average'

[FVInterpolationMethods]
  [average]
    type = FVGeometricAverage
  []
[]

[FVGradientMethods]
  [jump_aware]
    type = FVPressureJumpGreenGaussGradient
  []
  [velocity_jump_aware]
    type = FVVelocityJumpGreenGaussGradient
  []
  [reconstructed]
    type = FVReconstructedPressureGradient
    base_gradient_method = jump_aware
    gradient_relaxation = 1.0
  []
[]

[Mesh]
  [mesh]
    type = CartesianMeshGenerator
    dim = 2

    # Two blocks are enough for this verification problem: one block on each side of the
    # internal porous-baffle interface.
    dx = '${length_left} ${length_right}'
    dy = '${height}'
    ix = '${nx_left} ${nx_right}'
    iy = '${ny}'
    subdomain_id = '1 2'
  []
  [baffle]
    type = SideSetsBetweenSubdomainsGenerator
    input = mesh
    primary_block = '1'
    paired_block = '2'
    new_boundary = 'baffle'
  []
  [split_top_bottom]
    type = BreakBoundaryOnSubdomainGenerator
    input = baffle
    boundaries = 'top bottom'
  []
  [delete]
    type = BoundaryDeletionGenerator
    boundary_names = 'top bottom'
    input = split_top_bottom
  []
[]

[Problem]
  linear_sys_names = 'u_system v_system pressure_system'
  previous_nl_solution_required = true
[]

[UserObjects]
  [pressure_jump]
    type = BernoulliFormLossPressureJump
    boundary = baffle
    porosity = porosity
    density = rho_aux
  []
  [rc]
    type = PorousRhieChowMassFlux
    u = superficial_u
    v = superficial_v
    pressure = pressure
    rho = 'rho_aux'
    porosity = 'porosity'
    p_diffusion_kernel = p_diffusion

    pressure_jump_models = pressure_jump
    pressure_jump_relaxation = 0.1
  []
[]

[Variables]
  [superficial_u]
    type = MooseLinearVariableFVReal
    solver_sys = u_system
    initial_condition = ${u_in}
    gradient_method = velocity_jump_aware
  []
  [superficial_v]
    type = MooseLinearVariableFVReal
    solver_sys = v_system
    initial_condition = 0.0
    gradient_method = velocity_jump_aware
  []
  [pressure]
    type = MooseLinearVariableFVReal
    solver_sys = pressure_system
    initial_condition = ${p_out}
    gradient_method = jump_aware
  []
[]

[LinearFVKernels]
  [u_advection]
    type = LinearPWCNSFVMomentumFlux
    variable = superficial_u
    advected_interp_method_name = ${advected_interp_method}
    mu = 'mu'
    u = superficial_u
    v = superficial_v
    momentum_component = 'x'
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = false
  []
  [v_advection]
    type = LinearPWCNSFVMomentumFlux
    variable = superficial_v
    advected_interp_method_name = ${advected_interp_method}
    mu = 'mu'
    u = superficial_u
    v = superficial_v
    momentum_component = 'y'
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = false
  []
  [u_pressure]
    type = LinearPWCNSFVMomentumPressure
    variable = superficial_u
    momentum_component = 'x'
    porosity = 'porosity'
    pressure = pressure
    gradient_method = reconstructed
  []
  [v_pressure]
    type = LinearPWCNSFVMomentumPressure
    variable = superficial_v
    momentum_component = 'y'
    porosity = 'porosity'
    pressure = pressure
    gradient_method = reconstructed
  []
  [p_diffusion]
    type = LinearFVPressureCorrectionDiffusionJump
    variable = pressure
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = false
  []
  [HbyA_divergence]
    type = LinearFVDivergence
    variable = pressure
    face_flux = HbyA
    force_boundary_execution = true
  []
[]

[LinearFVBCs]
  [left_u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = left
    variable = superficial_u
    functor = ${u_in}
  []
  [left_v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = left
    variable = superficial_v
    functor = 0.0
  []

  [pressure-extrapolation]
    type = LinearFVExtrapolatedPressureBC
    boundary = 'left'
    variable = pressure
    use_two_term_expansion = true
  []

  # Fix the outlet pressure and leave the outlet velocity free.
  [outlet_u]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = right
    variable = superficial_u
    use_two_term_expansion = true
  []
  [outlet_v]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = right
    variable = superficial_v
    use_two_term_expansion = true
  []
  [outlet_p]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = right
    variable = pressure
    functor = ${p_out}
  []

  # Symmetry removes any wall losses and keeps the exact solution one-dimensional.
  [symmetry_u]
    type = LinearFVVelocitySymmetryBC
    boundary = 'top_to_1 top_to_2 bottom_to_1 bottom_to_2'
    variable = superficial_u
    u = superficial_u
    v = superficial_v
    momentum_component = x
  []
  [symmetry_v]
    type = LinearFVVelocitySymmetryBC
    boundary = 'top_to_1 top_to_2 bottom_to_1 bottom_to_2'
    variable = superficial_v
    u = superficial_u
    v = superficial_v
    momentum_component = y
  []
  [pressure_symmetric]
    type = LinearFVPressureSymmetryBC
    boundary = 'top_to_1 top_to_2 bottom_to_1 bottom_to_2'
    variable = pressure
    HbyA_flux = 'HbyA'
  []
[]

[Functions]
  [rho]
    type = ParsedFunction
    expression = '(1e3 - 0.05e3 * x)*${density_factor}'
  []
[]

[FunctorMaterials]
  [rho_exact]
    type = ParsedFunctorMaterial
    property_name = rho_exact
    expression = '(1e3 - 0.05e3 * x)*${density_factor}'
  []
  [porosity_exact]
    type = PiecewiseByBlockFunctorMaterial
    prop_name = porosity_exact
    subdomain_to_prop_value = '1 ${epsilon_left} 2 ${epsilon_right}'
  []
  [porosity]
    type = PiecewiseByBlockFunctorMaterial
    prop_name = porosity
    subdomain_to_prop_value = '1 ${epsilon_left} 2 ${epsilon_right}'
  []
  [mu]
    type = ParsedFunctorMaterial
    property_name = mu
    expression = '${mu}'
  []
  [u_superficial_exact]
    type = ParsedFunctorMaterial
    property_name = u_superficial_exact
    functor_names = 'rho_exact'
    expression = '${mass_flux_exact} / rho_exact'
  []
  [u_superficial_prime_exact]
    type = ParsedFunctorMaterial
    property_name = u_superficial_prime_exact
    functor_names = 'rho_exact'
    expression = '-${mass_flux_exact} * ${density_slope_exact} / (rho_exact * rho_exact)'
  []
  [p_exact]
    type = ParsedFunctorMaterial
    property_name = p_exact
    functor_names = 'rho_exact u_superficial_prime_exact'
    expression = 'if(x < ${length_left}, ${pressure_left_constant_exact} - ${mass_flux_exact} * ${mass_flux_exact} / (rho_exact * ${epsilon_left} * ${epsilon_left}) + ${mu} * u_superficial_prime_exact / ${epsilon_left}, ${pressure_right_constant_exact} - ${mass_flux_exact} * ${mass_flux_exact} / (rho_exact * ${epsilon_right} * ${epsilon_right}) + ${mu} * u_superficial_prime_exact / ${epsilon_right})'
  []
[]

[AuxVariables]
  [rho_aux]
    type = MooseLinearVariableFVReal
  []
  [porosity_aux]
    type = MooseLinearVariableFVReal
  []
[]

[AuxKernels]
  [assign_rho_aux]
    type = FunctorAux
    variable = rho_aux
    functor = 'rho'
    execute_on = 'initial timestep_end'
  []
  [assign_porosity_aux]
    type = FunctorAux
    variable = porosity_aux
    functor = 'porosity'
    execute_on = 'initial timestep_end'
  []
[]

[Postprocessors]
  [L2_rho]
    type = ElementL2FunctorError
    approximate = rho_aux
    exact = rho_exact
  []
  [L2_porosity]
    type = ElementL2FunctorError
    approximate = porosity_aux
    exact = porosity_exact
  []
  [L2_superficial_u]
    type = ElementL2FunctorError
    approximate = superficial_u
    exact = u_superficial_exact
  []
  [L2_pressure]
    type = ElementL2FunctorError
    approximate = pressure
    exact = p_exact
  []

  # Overall pressure drop from the solved field.
  [p_left]
    type = SideAverageValue
    variable = pressure
    boundary = left
  []
  [p_right]
    type = SideAverageValue
    variable = pressure
    boundary = right
  []
  [delta_p]
    type = ParsedPostprocessor
    expression = 'p_left - p_right'
    pp_names = 'p_left p_right'
  []

  # The block averages give a more direct measure of the interface jump because the exact
  # solution is piecewise constant in each block.
  [p_block_1]
    type = ElementAverageValue
    variable = pressure
    block = 1
  []
  [p_block_2]
    type = ElementAverageValue
    variable = pressure
    block = 2
  []
  [p_block_jump]
    type = ParsedPostprocessor
    expression = 'p_block_1 - p_block_2'
    pp_names = 'p_block_1 p_block_2'
  []

  # Check the superficial velocities on both sides of the jump.
  [u_block_1]
    type = ElementAverageValue
    variable = superficial_u
    block = 1
  []
  [u_block_2]
    type = ElementAverageValue
    variable = superficial_u
    block = 2
  []
[]

[VectorPostprocessors]
  [centerline_solution]
    type = LineValueSampler
    variable = 'rho_aux porosity_aux superficial_u pressure'
    start_point = '0.0 0.5 0'
    end_point = '${fparse length_left + length_right} 0.5 0'
    num_points = 50
    sort_by = x
    execute_on = timestep_end
  []
[]

[Executioner]
  type = SIMPLE
  momentum_l_abs_tol = 1e-14
  pressure_l_abs_tol = 1e-16
  momentum_l_tol = 0
  pressure_l_tol = 0
  rhie_chow_user_object = rc
  momentum_systems = 'u_system v_system'
  pressure_system = pressure_system
  momentum_equation_relaxation = 0.4
  pressure_variable_relaxation = 0.2
  num_iterations = 500
  pressure_absolute_tolerance = 1e-10
  # The exact transverse velocity is zero, so the y-momentum system has a vanishing normalization
  # factor and its normalized residual converges slowly without changing the solution. Its
  # tolerance is therefore not used as a convergence criterion.
  momentum_absolute_tolerance = '1e-10 1'
  momentum_petsc_options_iname = '-pc_type -pc_hypre_type'
  momentum_petsc_options_value = 'hypre boomeramg'
  pressure_petsc_options_iname = '-pc_type -pc_hypre_type'
  pressure_petsc_options_value = 'hypre boomeramg'
[]

[Outputs]
  csv = true
  execute_on = timestep_end
[]

(modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/linear-segregated/1d-simple-channel/porous-baffle-1d.i)

# One-dimensional verification of Bernoulli jumps, porous drag, and form loss.
# Mass conservation gives u_block_1 = u_block_2 = 0.1 and u_block_jump = 0.
# Without drag or form loss, the two pressure jumps (right minus left) are 0.03 and -0.03,
# so p_left = p_right = p_jump = 0. With mu = 1 and Forchheimer = 50, the porous block
# adds a pressure drop of 0.25. A form-loss coefficient of 10 adds 0.4 at each baffle,
# giving p_jump = p_left - p_right = 1.05 and baffle jumps of -0.37 and -0.43.
# Referencing the higher-porosity velocity 0.1 instead of the lower-porosity velocity 0.2 makes the
# form loss 0.1 at that baffle. With reference_velocity_side = 'lower_porosity higher_porosity',
# the baffle jumps are -0.37 and -0.13; with 'higher_porosity lower_porosity', they are -0.07 and
# -0.43. Both give p_jump = 0.75.

mu = 0.0 # 1e-2
rho = 2.0
advected_interp_method = 'average'
forchheimer_value = 50
inlet_u = 0.1

[FVGradientMethods]
  [jump_aware]
    type = FVPressureJumpGreenGaussGradient
  []
  [velocity_jump_aware]
    type = FVVelocityJumpGreenGaussGradient
  []
  [reconstructed]
    type = FVReconstructedPressureGradient
    base_gradient_method = jump_aware
    gradient_relaxation = 1.0
  []
[]

[Mesh]
  [mesh]
    type = CartesianMeshGenerator
    dim = 1
    dx = '0.5 0.5 0.5'
    ix = '200 200 200'
    subdomain_id = '1 2 3'
  []
  [baffle]
    type = SideSetsBetweenSubdomainsGenerator
    input = mesh
    primary_block = '1'
    paired_block = '2'
    new_boundary = 'baffle'
  []
  [baffle2]
    type = SideSetsBetweenSubdomainsGenerator
    input = baffle
    primary_block = '2'
    paired_block = '3'
    new_boundary = 'baffle2'
  []
[]

[Problem]
  linear_sys_names = 'u_system pressure_system'
  previous_nl_solution_required = true
[]

[UserObjects]
  [pressure_jump]
    type = BernoulliFormLossPressureJump
    boundary = 'baffle baffle2'
    porosity = porosity
    density = ${rho}
  []
  [rc]
    type = PorousRhieChowMassFlux
    u = superficial_u
    pressure = pressure
    rho = ${rho}
    porosity = porosity
    p_diffusion_kernel = p_diffusion
    pressure_jump_models = pressure_jump
    pressure_jump_relaxation = 0.1
  []
[]

[Variables]
  [superficial_u]
    type = MooseLinearVariableFVReal
    solver_sys = u_system
    initial_condition = ${inlet_u}
    gradient_method = velocity_jump_aware
  []
  [pressure]
    type = MooseLinearVariableFVReal
    solver_sys = pressure_system
    initial_condition = 0.0
    gradient_method = jump_aware
  []
[]

[FVInterpolationMethods]
  [average]
    type = FVGeometricAverage
  []
[]

[LinearFVKernels]
  [u_advection]
    type = LinearPWCNSFVMomentumFlux
    variable = superficial_u
    advected_interp_method_name = ${advected_interp_method}
    mu = ${mu}
    u = superficial_u
    momentum_component = 'x'
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = false
  []
  [u_pressure]
    type = LinearPWCNSFVMomentumPressure
    variable = superficial_u
    momentum_component = 'x'
    porosity = porosity
    pressure = pressure
    gradient_method = reconstructed
  []
  [u_friction]
    type = LinearFVMomentumFriction
    variable = superficial_u
    Forchheimer_name = forch
    porosity = porosity
    rho = ${rho}
    u = superficial_u
    momentum_component = 'x'
    block = 2
  []
  [p_diffusion]
    type = LinearFVPressureCorrectionDiffusionJump
    variable = pressure
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = false
  []
  [HbyA_divergence]
    type = LinearFVDivergence
    variable = pressure
    face_flux = HbyA
    force_boundary_execution = true
  []
[]

[LinearFVBCs]
  [inlet_u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = left
    variable = superficial_u
    functor = ${inlet_u}
  []
  [outlet_u]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = right
    variable = superficial_u
    use_two_term_expansion = false
  []
  [outlet_p]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = right
    variable = pressure
    functor = 0.0
  []
[]

[FunctorMaterials]
  [forch]
    type = GenericVectorFunctorMaterial
    prop_names = forch
    prop_values = '${forchheimer_value} ${forchheimer_value} ${forchheimer_value}'
  []
  [porosity]
    type = PiecewiseByBlockFunctorMaterial
    prop_name = porosity
    subdomain_to_prop_value = '1 0.5 2 1.0 3 0.5'
  []
[]

[Postprocessors]
  [p_left]
    type = SideAverageValue
    variable = pressure
    boundary = left
  []
  [p_right]
    type = SideAverageValue
    variable = pressure
    boundary = right
  []
  [p_jump]
    type = ParsedPostprocessor
    expression = 'p_left - p_right'
    pp_names = 'p_left p_right'
  []
  # These points are the centers of the cells adjacent to x = 0.5 and x = 1.0. The mesh spacing
  # is 0.5 / 200 = 0.0025, so the adjacent centers are offset by half that distance.
  [p_baffle_1_left]
    type = PointValue
    variable = pressure
    point = '0.49875 0 0'
  []
  [p_baffle_1_right]
    type = PointValue
    variable = pressure
    point = '0.50125 0 0'
  []
  [p_baffle_2_left]
    type = PointValue
    variable = pressure
    point = '0.99875 0 0'
  []
  [p_baffle_2_right]
    type = PointValue
    variable = pressure
    point = '1.00125 0 0'
  []
  [p_baffle_1_jump]
    type = ParsedPostprocessor
    expression = 'p_baffle_1_right - p_baffle_1_left'
    pp_names = 'p_baffle_1_right p_baffle_1_left'
  []
  [p_baffle_2_jump]
    type = ParsedPostprocessor
    expression = 'p_baffle_2_right - p_baffle_2_left'
    pp_names = 'p_baffle_2_right p_baffle_2_left'
  []
  [u_block_1]
    type = ElementAverageValue
    variable = superficial_u
    block = 1
  []
  [u_block_2]
    type = ElementAverageValue
    variable = superficial_u
    block = 2
  []
  [u_block_jump]
    type = ParsedPostprocessor
    expression = 'u_block_1 - u_block_2'
    pp_names = 'u_block_1 u_block_2'
  []
[]

[AuxVariables]
  [porosity_aux]
    type = MooseLinearVariableFVReal
  []
[]

[AuxKernels]
  [por]
    type = FunctorAux
    variable = porosity_aux
    functor = porosity
    execute_on = 'timestep_end'
  []
[]

[Executioner]
  type = SIMPLE
  momentum_l_abs_tol = 1e-12
  pressure_l_abs_tol = 1e-12
  momentum_l_tol = 0
  pressure_l_tol = 0
  rhie_chow_user_object = rc
  momentum_systems = 'u_system'
  pressure_system = pressure_system
  momentum_equation_relaxation = 0.3
  pressure_variable_relaxation = 0.1
  num_iterations = 1000
  pressure_absolute_tolerance = 1e-8
  # The exact superficial velocity is uniform, so the momentum normalization factor
  # sum(|Ax - Ax_avg| + |b - Ax_avg|) nearly vanishes and the normalized momentum residual measures
  # round-off. It levels off at a value that depends on the partitioning (near 1e-2 on one process,
  # near 1e-1 on three), so only the pressure residual is used as a convergence criterion.
  momentum_absolute_tolerance = 1
  momentum_petsc_options_iname = '-pc_type -pc_hypre_type'
  momentum_petsc_options_value = 'hypre boomeramg'
  pressure_petsc_options_iname = '-pc_type -pc_hypre_type'
  pressure_petsc_options_value = 'hypre boomeramg'
[]

[Outputs]
  csv = true
  execute_on = 'timestep_end'
[]

(modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/linear-segregated/radial-baffle-rz/radial-baffle-rz.i)

# Steady inviscid radial outflow through a porous annulus (r < r0) into a region of unit porosity
# (r > r0), separated by a pressure baffle on the cylinder r = r0, in RZ coordinates with x as the
# radial direction. With mass flow per radian m and the friction
# f = mu D U + 0.5 rho F |U| / eps U of LinearFVMomentumFriction, the exact solution is
#   U_r = m / (rho r), U_z = 0,
#   p = -m^2 / (2 rho eps^2 r^2) - mu D m / (rho eps) ln(r) + F m^2 / (2 rho eps^2 r) + c
# in each region. At the baffle the pressure jumps by the Bernoulli and form-loss terms of
# BernoulliFormLossPressureJump, with the porous side as the lower-porosity reference. The flow is
# inviscid because baffle faces carry no normal viscous stress, which a viscous solution of this
# form would need.

rho = 1
m = 1
eps_bed = 0.4
# Darcy coefficients mu D of each region. Central advection of a divergence-free flow adds almost
# nothing to the momentum diagonal, so without viscosity the friction keeps the momentum systems
# nonsingular, and the plenum needs it too.
darcy_bed = 1
darcy_plenum = 0.5
forchheimer = 2
form_loss = 0.5
r_inner = 1
r_baffle = 1.5
r_outer = 2
height = 0.25
nx = 8
ny = 2

[Problem]
  linear_sys_names = 'u_system v_system pressure_system'
  previous_nl_solution_required = true
[]

[Mesh]
  [generated]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = ${r_inner}
    xmax = ${r_outer}
    ymax = ${height}
    nx = ${nx}
    ny = ${ny}
  []
  [bed]
    type = ParsedSubdomainMeshGenerator
    input = generated
    combinatorial_geometry = 'x < ${r_baffle}'
    block_id = 1
    block_name = bed
  []
  [plenum]
    type = ParsedSubdomainMeshGenerator
    input = bed
    combinatorial_geometry = 'x > ${r_baffle}'
    block_id = 2
    block_name = plenum
  []
  [baffle]
    type = SideSetsBetweenSubdomainsGenerator
    input = plenum
    primary_block = bed
    paired_block = plenum
    new_boundary = baffle
  []
  coord_type = RZ
  rz_coord_axis = Y
[]

[FVGradientMethods]
  [jump_aware]
    type = FVPressureJumpGreenGaussGradient
  []
  [velocity_jump_aware]
    type = FVVelocityJumpGreenGaussGradient
  []
  [reconstructed]
    type = FVReconstructedPressureGradient
    base_gradient_method = jump_aware
    gradient_relaxation = 1
  []
[]

[FVInterpolationMethods]
  [average]
    type = FVGeometricAverage
  []
[]

[UserObjects]
  [pressure_jump]
    type = BernoulliFormLossPressureJump
    boundary = baffle
    porosity = porosity
    density = ${rho}
    form_loss = ${form_loss}
  []
  [rc]
    type = PorousRhieChowMassFlux
    u = superficial_u
    v = superficial_v
    pressure = pressure
    rho = ${rho}
    porosity = porosity
    p_diffusion_kernel = p_diffusion
    pressure_projection_method = consistent
    pressure_jump_models = pressure_jump
  []
[]

[Variables]
  [superficial_u]
    type = MooseLinearVariableFVReal
    solver_sys = u_system
    initial_condition = 0
    gradient_method = velocity_jump_aware
  []
  [superficial_v]
    type = MooseLinearVariableFVReal
    solver_sys = v_system
    initial_condition = 0
    gradient_method = velocity_jump_aware
  []
  [pressure]
    type = MooseLinearVariableFVReal
    solver_sys = pressure_system
    initial_condition = 0
    gradient_method = jump_aware
  []
[]

[LinearFVKernels]
  [u_advection_stress]
    type = LinearPWCNSFVMomentumFlux
    variable = superficial_u
    advected_interp_method_name = average
    mu = 0
    u = superficial_u
    v = superficial_v
    momentum_component = x
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = false
  []
  [v_advection_stress]
    type = LinearPWCNSFVMomentumFlux
    variable = superficial_v
    advected_interp_method_name = average
    mu = 0
    u = superficial_u
    v = superficial_v
    momentum_component = y
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = false
  []
  [u_pressure]
    type = LinearPWCNSFVMomentumPressure
    variable = superficial_u
    pressure = pressure
    porosity = porosity
    momentum_component = x
    gradient_method = reconstructed
  []
  [v_pressure]
    type = LinearPWCNSFVMomentumPressure
    variable = superficial_v
    pressure = pressure
    porosity = porosity
    momentum_component = y
    gradient_method = reconstructed
  []
  [u_friction]
    type = LinearFVMomentumFriction
    variable = superficial_u
    Darcy_name = darcy
    Forchheimer_name = forchheimer
    mu = 1
    rho = ${rho}
    porosity = porosity
    u = superficial_u
    v = superficial_v
    momentum_component = x
  []
  [v_friction]
    type = LinearFVMomentumFriction
    variable = superficial_v
    Darcy_name = darcy
    Forchheimer_name = forchheimer
    mu = 1
    rho = ${rho}
    porosity = porosity
    u = superficial_u
    v = superficial_v
    momentum_component = y
  []
  [p_diffusion]
    type = LinearFVPressureCorrectionDiffusionJump
    variable = pressure
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = false
  []
  [HbyA_divergence]
    type = LinearFVDivergence
    variable = pressure
    face_flux = HbyA
    force_boundary_execution = true
  []
[]

[LinearFVBCs]
  [velocity_x]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left right top bottom'
    variable = superficial_u
    functor = exact_u
  []
  [velocity_y]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left right top bottom'
    variable = superficial_v
    functor = 0
  []
  [pressure_flux]
    type = LinearFVPressureFluxBC
    boundary = 'left right top bottom'
    variable = pressure
    HbyA_flux = HbyA
    Ainv = Ainv
    rho = ${rho}
    u = superficial_u
    v = superficial_v
  []
[]

[FunctorMaterials]
  [porosity]
    type = PiecewiseByBlockFunctorMaterial
    prop_name = porosity
    subdomain_to_prop_value = 'bed ${eps_bed} plenum 1'
  []
  [friction_bed]
    type = GenericVectorFunctorMaterial
    prop_names = 'darcy forchheimer'
    prop_values = '${darcy_bed} ${darcy_bed} ${darcy_bed} ${forchheimer} ${forchheimer} ${forchheimer}'
    block = bed
  []
  [friction_plenum]
    type = GenericVectorFunctorMaterial
    prop_names = 'darcy forchheimer'
    prop_values = '${darcy_plenum} ${darcy_plenum} ${darcy_plenum} 0 0 0'
    block = plenum
  []
  [pressure_difference]
    type = ParsedFunctorMaterial
    property_name = pressure_difference
    functor_names = 'pressure exact_pressure'
    expression = 'pressure-exact_pressure'
  []
  # The pin fixes the pressure of the cell containing the pin point, so the computed pressure is
  # compared after removing the mean difference to the exact pressure.
  [shifted_pressure]
    type = ParsedFunctorMaterial
    property_name = shifted_pressure
    functor_names = 'pressure mean_pressure_difference'
    expression = 'pressure-mean_pressure_difference'
  []
  [pressure_error]
    type = ParsedFunctorMaterial
    property_name = pressure_error
    functor_names = 'shifted_pressure exact_pressure'
    expression = 'shifted_pressure-exact_pressure'
  []
  [u_error]
    type = ParsedFunctorMaterial
    property_name = u_error
    functor_names = 'superficial_u exact_u'
    expression = 'superficial_u-exact_u'
  []
[]

[Functions]
  [exact_u]
    type = ParsedFunction
    expression = '${m}/(${rho}*x)'
  []
  [bed_pressure]
    type = ParsedFunction
    expression = '-${m}^2/(2*${rho}*${eps_bed}^2*x^2) - ${darcy_bed}*${m}/(${rho}*${eps_bed})*log(x) '
                 '+ ${forchheimer}*${m}^2/(2*${rho}*${eps_bed}^2*x)'
  []
  [plenum_pressure]
    type = ParsedFunction
    expression = '-${m}^2/(2*${rho}*x^2) - ${darcy_plenum}*${m}/${rho}*log(x)'
  []
  # p_plenum(r0) - p_bed(r0) = 0.5 rho u0^2 (1/eps_bed^2 - 1) - 0.5 K rho (u0/eps_bed)^2 with
  # u0 = m / (rho r0).
  [exact_pressure]
    type = ParsedFunction
    symbol_names = 'pb pp u0'
    symbol_values = 'bed_pressure plenum_pressure ${fparse m / (rho * r_baffle)}'
    expression = 'if(x < ${r_baffle}, pb, pp + '
                 '(-${m}^2/(2*${rho}*${eps_bed}^2*${r_baffle}^2) - ${darcy_bed}*${m}/(${rho}*${eps_bed})*log(${r_baffle}) '
                 '+ ${forchheimer}*${m}^2/(2*${rho}*${eps_bed}^2*${r_baffle}) + ${m}^2/(2*${rho}*${r_baffle}^2) '
                 '+ ${darcy_plenum}*${m}/${rho}*log(${r_baffle})) '
                 '+ 0.5*${rho}*u0^2*(1/${eps_bed}^2 - 1) - 0.5*${form_loss}*${rho}*(u0/${eps_bed})^2)'
  []
[]

[Executioner]
  type = SIMPLE
  momentum_l_abs_tol = 1e-12
  pressure_l_abs_tol = 1e-12
  momentum_l_tol = 0
  pressure_l_tol = 0
  momentum_l_max_its = 100
  pressure_l_max_its = 100
  rhie_chow_user_object = rc
  momentum_systems = 'u_system v_system'
  pressure_system = pressure_system
  # Without viscosity nothing damps round-off in the axial velocity, and at a momentum relaxation
  # of 0.5 the iteration amplifies it by about 1% per iteration on the twice-refined mesh.
  momentum_equation_relaxation = 0.3
  pressure_variable_relaxation = 0.2
  num_iterations = 3000
  pressure_absolute_tolerance = 1e-9
  # The exact axial velocity is zero, so the axial momentum system has a vanishing normalization
  # factor and its normalized residual stays O(1e-2) at round-off. Its tolerance is therefore not
  # used as a convergence criterion; L2v checks the axial velocity instead.
  momentum_absolute_tolerance = '1e-9 1'
  momentum_petsc_options_iname = '-pc_type -pc_hypre_type'
  momentum_petsc_options_value = 'hypre boomeramg'
  pressure_petsc_options_iname = '-pc_type -pc_hypre_type'
  pressure_petsc_options_value = 'hypre boomeramg'
  print_fields = false

  pin_pressure = true
  pressure_pin_value = 0
  pressure_pin_point = '${fparse 0.5 * (r_baffle + r_outer)} ${fparse 0.5 * height} 0'
[]

[Outputs]
  [csv]
    type = CSV
    execute_on = FINAL
  []
[]

[Postprocessors]
  [mean_pressure_difference]
    type = ElementAverageFunctorPostprocessor
    functor = pressure_difference
    execute_on = FINAL
    # Must run before the pressure errors that subtract it.
    execution_order_group = -1
    outputs = none
  []
  [h]
    type = AverageElementSize
    execute_on = FINAL
  []
  [L2u]
    type = ElementL2FunctorError
    approximate = superficial_u
    exact = exact_u
    execute_on = FINAL
  []
  [L2v]
    type = ElementL2FunctorError
    approximate = superficial_v
    exact = 0
    execute_on = FINAL
  []
  [L2p]
    type = ElementL2FunctorError
    approximate = shifted_pressure
    exact = exact_pressure
    execute_on = FINAL
  []
  # The maximum pointwise errors detect cell-to-cell oscillations that the L2 norms average out.
  [max_u_error]
    type = ElementExtremeFunctorValue
    functor = u_error
    value_type = max_abs
    execute_on = FINAL
  []
  [max_v_error]
    type = ElementExtremeFunctorValue
    functor = superficial_v
    value_type = max_abs
    execute_on = FINAL
  []
  [max_p_error]
    type = ElementExtremeFunctorValue
    functor = pressure_error
    value_type = max_abs
    execute_on = FINAL
  []
[]

(modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/linear-segregated/1d-simple-channel/form-loss-orientation.i)

# Isolated form-loss orientation check. Equal porosity removes the reversible Bernoulli jump,
# so the remaining pressure drop is the irreversible form loss of 0.1.

mu = 0.0
rho = 2.0
advected_interp_method = 'average'
inlet_u = 0.1
form_loss = 10
primary_block = 1
paired_block = 2
velocity_inlet_boundary = left
pressure_outlet_boundary = right
direction = 1

[FVGradientMethods]
  [jump_aware]
    type = FVPressureJumpGreenGaussGradient
  []
  [velocity_jump_aware]
    type = FVVelocityJumpGreenGaussGradient
  []
  [reconstructed]
    type = FVReconstructedPressureGradient
    base_gradient_method = jump_aware
    gradient_relaxation = 1.0
  []
[]

[Mesh]
  [mesh]
    type = CartesianMeshGenerator
    dim = 1
    dx = '0.5 0.5'
    ix = '8 8'
    subdomain_id = '1 2'
  []
  [baffle]
    type = SideSetsBetweenSubdomainsGenerator
    input = mesh
    primary_block = ${primary_block}
    paired_block = ${paired_block}
    new_boundary = 'baffle'
  []
[]

[Problem]
  linear_sys_names = 'u_system pressure_system'
  previous_nl_solution_required = true
[]

[UserObjects]
  [pressure_jump]
    type = BernoulliFormLossPressureJump
    boundary = 'baffle'
    porosity = porosity
    density = ${rho}
    form_loss = ${form_loss}
  []
  [rc]
    type = PorousRhieChowMassFlux
    u = superficial_u
    pressure = pressure
    rho = ${rho}
    porosity = porosity
    p_diffusion_kernel = p_diffusion
    pressure_jump_models = pressure_jump
    # Dampen the quadratic form-loss feedback during the SIMPLE outer iteration.
    pressure_jump_relaxation = 0.1
  []
[]

[Variables]
  [superficial_u]
    type = MooseLinearVariableFVReal
    solver_sys = u_system
    initial_condition = ${inlet_u}
    gradient_method = velocity_jump_aware
  []
  [pressure]
    type = MooseLinearVariableFVReal
    solver_sys = pressure_system
    initial_condition = 0.0
    gradient_method = jump_aware
  []
[]

[FVInterpolationMethods]
  [average]
    type = FVGeometricAverage
  []
[]

[LinearFVKernels]
  [u_advection]
    type = LinearPWCNSFVMomentumFlux
    variable = superficial_u
    advected_interp_method_name = ${advected_interp_method}
    mu = ${mu}
    u = superficial_u
    momentum_component = 'x'
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = false
  []
  [u_pressure]
    type = LinearPWCNSFVMomentumPressure
    variable = superficial_u
    momentum_component = 'x'
    porosity = porosity
    pressure = pressure
    gradient_method = reconstructed
  []
  [p_diffusion]
    type = LinearFVPressureCorrectionDiffusionJump
    variable = pressure
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = false
  []
  [HbyA_divergence]
    type = LinearFVDivergence
    variable = pressure
    face_flux = HbyA
    force_boundary_execution = true
  []
[]

[LinearFVBCs]
  [inlet_u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = ${velocity_inlet_boundary}
    variable = superficial_u
    functor = ${inlet_u}
  []
  [outlet_u]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = ${pressure_outlet_boundary}
    variable = superficial_u
    use_two_term_expansion = false
  []
  [outlet_p]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = ${pressure_outlet_boundary}
    variable = pressure
    functor = 0.0
  []
[]

[FunctorMaterials]
  [porosity]
    type = GenericFunctorMaterial
    prop_names = porosity
    prop_values = 1
  []
[]

[Postprocessors]
  [p_left]
    type = SideAverageValue
    variable = pressure
    boundary = left
  []
  [p_right]
    type = SideAverageValue
    variable = pressure
    boundary = right
  []
  [directed_pressure_loss]
    type = ParsedPostprocessor
    expression = '${direction} * (p_left - p_right)'
    pp_names = 'p_left p_right'
  []
  [u_avg]
    type = ElementAverageValue
    variable = superficial_u
  []
[]

[Executioner]
  type = SIMPLE
  momentum_l_abs_tol = 1e-12
  pressure_l_abs_tol = 1e-12
  momentum_l_tol = 0
  pressure_l_tol = 0
  rhie_chow_user_object = rc
  momentum_systems = 'u_system'
  pressure_system = pressure_system
  momentum_equation_relaxation = 0.3
  pressure_variable_relaxation = 0.1
  # All four orientations recover the analytic solution within 500 iterations.
  num_iterations = 500
  pressure_absolute_tolerance = 1e-10
  # The exact superficial velocity is uniform, so the momentum normalization factor nearly vanishes
  # and the normalized momentum residual measures round-off, which cycles at a level that depends
  # on the partitioning. Only the pressure residual is used as a convergence criterion.
  momentum_absolute_tolerance = 1
  momentum_petsc_options_iname = '-pc_type -pc_hypre_type'
  momentum_petsc_options_value = 'hypre boomeramg'
  pressure_petsc_options_iname = '-pc_type -pc_hypre_type'
  pressure_petsc_options_value = 'hypre boomeramg'
[]

[Outputs]
  csv = true
  execute_on = 'timestep_end'
[]

(modules/navier_stokes/test/tests/finite_volume/pins/channel-flow/linear-segregated/2d-simple-channel-enthalpy/porous-baffle-2d-enthalpy-fluid-only.i)

# Steady fluid-only enthalpy transport through a 0.5-m heated porous block.
# With rho = 1000, u = 2, cp = 4200, and q = 20 MW/m^3, conservation gives
# u_block_2 = 2, T_fluid_out = 301.19047619048, and h_fluid_out = 1265000.
# For 10 upwind cells in the heated block, T_fluid_block_2 = 300.65476190476.

mu = 2e-3
rho = 1000
advected_interp_method = 'upwind'
u_in = 2
forchheimer = 0
bf = '0 0'

cp_fluid = 4200
k_fluid = 0.4
h_vol = 300
T_inlet = 300
T_initial = 300
q = 20000000

[FVGradientMethods]
  [jump_aware]
    type = FVPressureJumpGreenGaussGradient
  []
  [velocity_jump_aware]
    type = FVVelocityJumpGreenGaussGradient
  []
  [reconstructed]
    type = FVReconstructedPressureGradient
    base_gradient_method = jump_aware
    gradient_relaxation = 1.0
  []
[]

[Mesh]
  [mesh]
    type = CartesianMeshGenerator
    dim = 2
    dx = '0.5 0.5 0.5'
    dy = '0.5'
    ix = '10 10 10'
    iy = '3'
    subdomain_id = '1 2 3'
  []
  [baffle]
    type = SideSetsBetweenSubdomainsGenerator
    input = mesh
    primary_block = '1'
    paired_block = '2'
    new_boundary = 'baffle'
  []
  [baffle2]
    type = SideSetsBetweenSubdomainsGenerator
    input = baffle
    primary_block = '2'
    paired_block = '3'
    new_boundary = 'baffle2'
  []
  [split_top_bottom]
    type = BreakBoundaryOnSubdomainGenerator
    input = baffle2
    boundaries = 'top bottom'
  []
  [delete]
    type = BoundaryDeletionGenerator
    boundary_names = 'top bottom'
    input = split_top_bottom
  []
[]

[FVInterpolationMethods]
  [upwind]
    type = FVAdvectedUpwind
  []
[]

[Problem]
  linear_sys_names = 'u_system v_system pressure_system energy_system'
  previous_nl_solution_required = true
[]

[UserObjects]
  [pressure_jump]
    type = BernoulliFormLossPressureJump
    boundary = 'baffle baffle2'
    porosity = porosity
    density = ${rho}
    form_loss = ${bf}
    reference_velocity_side = 'lower_porosity higher_porosity'
  []
  [rc]
    type = PorousRhieChowMassFlux
    u = superficial_u
    v = superficial_v
    pressure = pressure
    rho = ${rho}
    porosity = porosity
    p_diffusion_kernel = p_diffusion
    pressure_jump_models = pressure_jump
    pressure_jump_relaxation = 0.1
  []
[]

[Variables]
  [superficial_u]
    type = MooseLinearVariableFVReal
    solver_sys = u_system
    initial_condition = ${u_in}
    gradient_method = velocity_jump_aware
  []
  [superficial_v]
    type = MooseLinearVariableFVReal
    solver_sys = v_system
    initial_condition = 0.0
    gradient_method = velocity_jump_aware
  []
  [pressure]
    type = MooseLinearVariableFVReal
    solver_sys = pressure_system
    initial_condition = 0.0
    gradient_method = jump_aware
  []
  [h_fluid]
    type = MooseLinearVariableFVReal
    solver_sys = energy_system
    initial_condition = '${fparse cp_fluid * T_initial}'
  []
[]

[LinearFVKernels]
  [u_advection]
    type = LinearPWCNSFVMomentumFlux
    variable = superficial_u
    advected_interp_method_name = ${advected_interp_method}
    mu = ${mu}
    u = superficial_u
    v = superficial_v
    momentum_component = 'x'
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = false
  []
  [v_advection]
    type = LinearPWCNSFVMomentumFlux
    variable = superficial_v
    advected_interp_method_name = ${advected_interp_method}
    mu = ${mu}
    u = superficial_u
    v = superficial_v
    momentum_component = 'y'
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = false
  []
  [u_pressure]
    type = LinearPWCNSFVMomentumPressure
    variable = superficial_u
    momentum_component = 'x'
    porosity = porosity
    pressure = pressure
    gradient_method = reconstructed
  []
  [v_pressure]
    type = LinearPWCNSFVMomentumPressure
    variable = superficial_v
    momentum_component = 'y'
    porosity = porosity
    pressure = pressure
    gradient_method = reconstructed
  []
  [u_friction]
    type = LinearFVMomentumFriction
    variable = superficial_u
    Forchheimer_name = forch
    porosity = porosity
    rho = ${rho}
    u = superficial_u
    v = superficial_v
    momentum_component = 'x'
    block = 2
  []
  [v_friction]
    type = LinearFVMomentumFriction
    variable = superficial_v
    Forchheimer_name = forch
    porosity = porosity
    rho = ${rho}
    u = superficial_u
    v = superficial_v
    momentum_component = 'y'
    block = 2
  []
  [p_diffusion]
    type = LinearFVPressureCorrectionDiffusionJump
    variable = pressure
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = false
  []
  [HbyA_divergence]
    type = LinearFVDivergence
    variable = pressure
    face_flux = HbyA
    force_boundary_execution = true
  []

  [fluid_advection]
    type = LinearFVEnergyAdvection
    variable = h_fluid
    advected_quantity = enthalpy
    advected_interp_method = ${advected_interp_method}
    rhie_chow_user_object = rc
  []
  [fluid_diffusion]
    type = LinearFVDiffusion
    variable = h_fluid
    diffusion_coeff = kappa_fluid
    use_nonorthogonal_correction = false
  []
  [fluid_source]
    type = LinearFVSource
    variable = h_fluid
    source_density = q_fluid_scaled
    block = 2
  []
[]

[LinearFVBCs]
  [left_u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = left
    variable = superficial_u
    functor = ${u_in}
  []
  [outlet_u]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = right
    variable = superficial_u
    use_two_term_expansion = false
  []

  [left_v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = left
    variable = superficial_v
    functor = 0.0
  []
  [outlet_v]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = right
    variable = superficial_v
    use_two_term_expansion = false
  []

  [symmetry-u]
    type = LinearFVVelocitySymmetryBC
    boundary = 'top_to_2 bottom_to_2 top_to_1 top_to_3 bottom_to_1 bottom_to_3'
    variable = superficial_u
    u = superficial_u
    v = superficial_v
    momentum_component = x
  []
  [symmetry-v]
    type = LinearFVVelocitySymmetryBC
    boundary = 'top_to_2 bottom_to_2 top_to_1 top_to_3 bottom_to_1 bottom_to_3'
    variable = superficial_v
    u = superficial_u
    v = superficial_v
    momentum_component = y
  []

  [outlet_p]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = right
    variable = pressure
    functor = 0.0
  []

  [pressure-extrapolation]
    type = LinearFVPressureFluxBC
    boundary = 'bottom_to_1 bottom_to_2 bottom_to_3 top_to_1 top_to_2 top_to_3'
    variable = pressure
    HbyA_flux = HbyA
    Ainv = Ainv
    rho = ${rho}
    u = superficial_u
    v = superficial_v
  []

  [left_T_fluid]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = left
    variable = T_fluid
    functor = ${T_inlet}
  []
  [left_h_fluid]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = left
    variable = h_fluid
    functor = h_from_p_T
  []
  [walls_h_fluid]
    type = LinearFVAdvectionDiffusionFunctorNeumannBC
    boundary = 'top_to_1 top_to_2 top_to_3 bottom_to_1 bottom_to_2 bottom_to_3'
    variable = h_fluid
    functor = 0.0
    diffusion_coeff = kappa_fluid
  []
  [outlet_h_fluid]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = right
    variable = h_fluid
    use_two_term_expansion = false
  []
[]

[FunctorMaterials]
  [forch]
    type = GenericVectorFunctorMaterial
    prop_names = forch
    prop_values = '${forchheimer} ${forchheimer} ${forchheimer}'
  []
  [porosity]
    type = PiecewiseByBlockFunctorMaterial
    prop_name = porosity
    subdomain_to_prop_value = '1 1.0 2 0.5 3 1.0'
  []
  [fluid_constants]
    type = GenericFunctorMaterial
    prop_names = 'cp_fluid kappa_fluid h_vol q_fluid_scaled'
    prop_values = '${cp_fluid} ${fparse k_fluid/cp_fluid} ${h_vol} ${q}'
  []
  [fluid_enthalpy_material]
    type = LinearFVEnthalpyFunctorMaterial
    pressure = pressure
    T_fluid = T_fluid
    h = h_fluid
    h_from_p_T_functor = h_from_p_T_functor
    T_from_p_h_functor = T_from_p_h_functor
  []
  [h_from_p_T_functor]
    type = ParsedFunctorMaterial
    property_name = h_from_p_T_functor
    functor_names = 'T_fluid'
    expression = '${cp_fluid} * T_fluid'
  []
  [T_from_p_h_functor]
    type = ParsedFunctorMaterial
    property_name = T_from_p_h_functor
    functor_names = 'h_fluid'
    expression = 'h_fluid / ${cp_fluid}'
  []
[]

[Postprocessors]
  [u_block_2]
    type = ElementAverageValue
    variable = superficial_u
    block = 2
  []
  [T_fluid_out]
    type = SideAverageValue
    variable = T_fluid
    boundary = right
  []
  [T_fluid_block_2]
    type = ElementAverageValue
    variable = T_fluid
    block = 2
  []
  [h_fluid_out]
    type = SideAverageValue
    variable = h_fluid
    boundary = right
  []
[]

[AuxVariables]
  [porosity_aux]
    type = MooseLinearVariableFVReal
  []
  [T_fluid]
    type = MooseLinearVariableFVReal
    initial_condition = ${T_initial}
  []
[]

[AuxKernels]
  [por]
    type = FunctorAux
    variable = porosity_aux
    functor = porosity
    execute_on = 'NONLINEAR'
  []
  [fluid_temperature]
    type = FunctorAux
    variable = T_fluid
    functor = T_from_p_h
    execute_on = 'NONLINEAR'
  []
[]

[VectorPostprocessors]
  # Sample each x-cell once at its center.
  [centerline]
    type = LineValueSampler
    variable = 'superficial_u T_fluid'
    start_point = '0.0119047619047619 0.25 0'
    end_point = '1.4880952380952381 0.25 0'
    num_points = 30
    sort_by = x
  []
[]

[Executioner]
  type = SIMPLE
  momentum_l_abs_tol = 1e-14
  pressure_l_abs_tol = 1e-14
  energy_l_abs_tol = 1e-14
  momentum_l_tol = 0
  pressure_l_tol = 0
  energy_l_tol = 0
  rhie_chow_user_object = rc
  momentum_systems = 'u_system v_system'
  pressure_system = pressure_system
  energy_system = energy_system
  momentum_equation_relaxation = 0.2
  pressure_variable_relaxation = 0.2
  energy_equation_relaxation = 0.7
  num_iterations = 1000
  # These outer tolerances stop once the velocity and zero-pressure fields have stabilized to
  # substantially below the accuracy exercised by the temperature and velocity profile checks.
  pressure_absolute_tolerance = 1e-6
  momentum_absolute_tolerance = 1e-2
  energy_absolute_tolerance = 1e-8
  momentum_petsc_options_iname = '-pc_type -pc_hypre_type'
  momentum_petsc_options_value = 'hypre boomeramg'
  pressure_petsc_options_iname = '-pc_type -pc_hypre_type'
  pressure_petsc_options_value = 'hypre boomeramg'
  energy_petsc_options_iname = '-pc_type -pc_hypre_type'
  energy_petsc_options_value = 'hypre boomeramg'
[]

[Outputs]
  csv = true
  execute_on = FINAL
[]

(modules/navier_stokes/test/tests/finite_volume/pins/mms/linear-segregated/hydrostatic-jump/hydrostatic-jump.i)

# Fluid at rest under gravity in a closed box split by a pressure baffle across which porosity
# jumps. The interface x - skew * y = interface is a mesh line.
# The superficial momentum balance at rest is 0 = -eps grad(p) + eps rho g, so
#   U = 0 and grad(p) = rho g on both sides,
# independent of porosity. With density rho = rho_0 + strat * z linear in the height
# z = -g . x / |g|, the exact pressure is
#   p = -|g| (rho_0 z + strat z^2 / 2) + const.
# The body force is the porosity-weighted buoyancy eps rho g, applied either by
# LinearFVMomentumBuoyancy or, with LinearFVKernels/inactive = 'u_buoyancy v_buoyancy', by
# LinearFVMomentumBoussinesq with a temperature chosen so that -alpha rho_ref (T - T_ref) = rho.

rho_0 = 1
strat = -0.1
eps_left = 1
eps_right = 1
mu = 0.01
reaction = 1
gx = 0
gy = -9.81
interface = 0.5
skew = 0
nx = 8
ny = 8
nonorthogonal = false
alpha = 1e-3
T_ref = 300
pressure_gradient_method = reconstructed

[Problem]
  linear_sys_names = 'u_system v_system pressure_system'
  previous_nl_solution_required = true
[]

[Mesh]
  [generated]
    type = GeneratedMeshGenerator
    dim = 2
    nx = ${nx}
    ny = ${ny}
  []
  [right]
    type = ParsedSubdomainMeshGenerator
    input = generated
    combinatorial_geometry = 'x > ${interface}'
    block_id = 2
    block_name = right
  []
  [left]
    type = ParsedSubdomainMeshGenerator
    input = right
    combinatorial_geometry = 'x < ${interface}'
    block_id = 1
    block_name = left
  []
  [skew]
    type = ParsedNodeTransformGenerator
    input = left
    x_function = 'x + ${skew} * y'
  []
  [baffle]
    type = SideSetsBetweenSubdomainsGenerator
    input = skew
    primary_block = left
    paired_block = right
    new_boundary = baffle
  []
[]

[FVGradientMethods]
  [jump_aware]
    type = FVPressureJumpGreenGaussGradient
  []
  [velocity_jump_aware]
    type = FVVelocityJumpGreenGaussGradient
  []
  [reconstructed]
    type = FVReconstructedPressureGradient
    base_gradient_method = jump_aware
    gradient_relaxation = 1
  []
[]

[FVInterpolationMethods]
  [average]
    type = FVGeometricAverage
  []
[]

[UserObjects]
  [pressure_jump]
    type = BernoulliFormLossPressureJump
    boundary = baffle
    porosity = porosity
    density = density
  []
  [rc]
    type = PorousRhieChowMassFlux
    u = superficial_u
    v = superficial_v
    pressure = pressure
    rho = density
    porosity = porosity
    p_diffusion_kernel = p_diffusion
    pressure_jump_models = pressure_jump
    pressure_projection_method = consistent
    # The body force enters H/A through face averages, so the 1/A diffusion must be averaged the
    # same way for the hydrostatic balance to hold on the face when the viscous A varies.
    pressure_diffusion_interpolation = average
  []
[]

[Variables]
  [superficial_u]
    type = MooseLinearVariableFVReal
    solver_sys = u_system
    initial_condition = 0
    gradient_method = velocity_jump_aware
  []
  [superficial_v]
    type = MooseLinearVariableFVReal
    solver_sys = v_system
    initial_condition = 0
    gradient_method = velocity_jump_aware
  []
  [pressure]
    type = MooseLinearVariableFVReal
    solver_sys = pressure_system
    initial_condition = 0
    gradient_method = jump_aware
  []
[]

[LinearFVKernels]
  inactive = 'u_boussinesq v_boussinesq'
  [u_advection]
    type = LinearPWCNSFVMomentumFlux
    variable = superficial_u
    advected_interp_method_name = average
    mu = ${mu}
    u = superficial_u
    v = superficial_v
    momentum_component = x
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = ${nonorthogonal}
  []
  [v_advection]
    type = LinearPWCNSFVMomentumFlux
    variable = superficial_v
    advected_interp_method_name = average
    mu = ${mu}
    u = superficial_u
    v = superficial_v
    momentum_component = y
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = ${nonorthogonal}
  []
  [u_pressure]
    type = LinearPWCNSFVMomentumPressure
    variable = superficial_u
    pressure = pressure
    porosity = porosity
    momentum_component = x
    gradient_method = ${pressure_gradient_method}
  []
  [v_pressure]
    type = LinearPWCNSFVMomentumPressure
    variable = superficial_v
    pressure = pressure
    porosity = porosity
    momentum_component = y
    gradient_method = ${pressure_gradient_method}
  []
  [u_buoyancy]
    type = LinearFVMomentumBuoyancy
    variable = superficial_u
    rho = density
    porosity = porosity
    reference_rho = 0
    gravity = '${gx} ${gy} 0'
    momentum_component = x
  []
  [v_buoyancy]
    type = LinearFVMomentumBuoyancy
    variable = superficial_v
    rho = density
    porosity = porosity
    reference_rho = 0
    gravity = '${gx} ${gy} 0'
    momentum_component = y
  []
  [u_boussinesq]
    type = LinearFVMomentumBoussinesq
    variable = superficial_u
    T_fluid = T_fluid
    rho = 1
    alpha_name = ${alpha}
    ref_temperature = ${T_ref}
    porosity = porosity
    gravity = '${gx} ${gy} 0'
    momentum_component = x
  []
  [v_boussinesq]
    type = LinearFVMomentumBoussinesq
    variable = superficial_v
    T_fluid = T_fluid
    rho = 1
    alpha_name = ${alpha}
    ref_temperature = ${T_ref}
    porosity = porosity
    gravity = '${gx} ${gy} 0'
    momentum_component = y
  []
  [u_reaction]
    type = LinearFVReaction
    variable = superficial_u
    coeff = ${reaction}
  []
  [v_reaction]
    type = LinearFVReaction
    variable = superficial_v
    coeff = ${reaction}
  []
  [p_diffusion]
    type = LinearFVPressureCorrectionDiffusionJump
    variable = pressure
    rhie_chow_user_object = rc
    use_nonorthogonal_correction = ${nonorthogonal}
  []
  [HbyA_divergence]
    type = LinearFVDivergence
    variable = pressure
    face_flux = HbyA
    force_boundary_execution = true
  []
[]

[LinearFVBCs]
  [velocity_x]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left right top bottom'
    variable = superficial_u
    functor = 0
  []
  [velocity_y]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left right top bottom'
    variable = superficial_v
    functor = 0
  []
  [pressure_flux]
    type = LinearFVPressureFluxBC
    boundary = 'left right top bottom'
    variable = pressure
    HbyA_flux = HbyA
    Ainv = Ainv
    rho = density
    u = superficial_u
    v = superficial_v
  []
[]

[FunctorMaterials]
  [density]
    type = ParsedFunctorMaterial
    property_name = density
    functor_names = 'exact_density'
    expression = 'exact_density'
  []
  [porosity]
    type = ParsedFunctorMaterial
    property_name = porosity
    expression = 'if(x-${skew}*y<${interface},${eps_left},${eps_right})'
  []
  [pressure_error]
    type = ParsedFunctorMaterial
    property_name = pressure_error
    functor_names = 'pressure exact_pressure'
    expression = 'pressure-exact_pressure'
  []
  [speed]
    type = ParsedFunctorMaterial
    property_name = speed
    functor_names = 'superficial_u superficial_v'
    expression = 'sqrt(superficial_u^2+superficial_v^2)'
  []
[]

[AuxVariables]
  [T_fluid]
    type = MooseLinearVariableFVReal
  []
[]

[AuxKernels]
  [T_fluid]
    type = FunctorAux
    variable = T_fluid
    functor = boussinesq_temperature
    execute_on = INITIAL
  []
[]

[Functions]
  [boussinesq_temperature]
    type = ParsedFunction
    expression = 'T_ref-rho/alpha'
    symbol_names = 'T_ref rho alpha'
    symbol_values = '${T_ref} exact_density ${alpha}'
  []
  [height]
    type = ParsedFunction
    expression = '-(gx*x+gy*y)/sqrt(gx^2+gy^2)'
    symbol_names = 'gx gy'
    symbol_values = '${gx} ${gy}'
  []
  [exact_density]
    type = ParsedFunction
    expression = 'rho_0+strat*hz'
    symbol_names = 'rho_0 strat hz'
    symbol_values = '${rho_0} ${strat} height'
  []
  # Hydrostatic pressure with the pin at (0.3, 0.3), inside a cell on 2^n meshes, set to zero.
  [exact_pressure]
    type = ParsedFunction
    expression = '-sqrt(gx^2+gy^2)*(rho_0*(hz-hp)+strat*(hz^2-hp^2)/2)'
    symbol_names = 'gx gy rho_0 strat hz hp'
    symbol_values = '${gx} ${gy} ${rho_0} ${strat} height '
                    '${fparse -(gx*(0.3+skew*0.3)+gy*0.3)/sqrt(gx^2+gy^2)}'
  []
[]

[Executioner]
  type = SIMPLE
  momentum_l_abs_tol = 1e-13
  pressure_l_abs_tol = 1e-13
  momentum_l_tol = 0
  pressure_l_tol = 0
  momentum_l_max_its = 100
  pressure_l_max_its = 100
  rhie_chow_user_object = rc
  momentum_systems = 'u_system v_system'
  pressure_system = pressure_system
  momentum_equation_relaxation = 0.7
  pressure_variable_relaxation = 0.5
  # At rest the normalized momentum residual has no flow to scale against and stalls near 1e-4,
  # so the absolute tolerances below are not reached. The solution is fixed well before 50
  # iterations.
  num_iterations = 50
  continue_on_max_its = true
  pressure_absolute_tolerance = 1e-12
  momentum_absolute_tolerance = 1e-12
  momentum_petsc_options_iname = '-pc_type -pc_hypre_type'
  momentum_petsc_options_value = 'hypre boomeramg'
  pressure_petsc_options_iname = '-pc_type -pc_hypre_type'
  pressure_petsc_options_value = 'hypre boomeramg'
  print_fields = false

  pin_pressure = true
  pressure_pin_value = 0
  pressure_pin_point = '${fparse 0.3 + skew * 0.3} 0.3 0'
[]

[Postprocessors]
  [max_speed]
    type = ElementExtremeFunctorValue
    functor = speed
    value_type = max
    execute_on = FINAL
  []
  # The pin sits at a point, not a cell centroid, so the error is compared up to a constant.
  [max_pressure_error]
    type = ElementExtremeFunctorValue
    functor = pressure_error
    value_type = max
    execute_on = FINAL
    outputs = none
  []
  [min_pressure_error]
    type = ElementExtremeFunctorValue
    functor = pressure_error
    value_type = min
    execute_on = FINAL
    outputs = none
  []
  [pressure_error_range]
    type = ParsedPostprocessor
    pp_names = 'max_pressure_error min_pressure_error'
    expression = 'max_pressure_error-min_pressure_error'
    execute_on = FINAL
  []
[]

[Outputs]
  [csv]
    type = CSV
    execute_on = FINAL
  []
[]

Overview
Input Parameters
Input Files