LinearFVDiffusion

Description

This kernel contributes to the system matrix and the right hand side of a system which is solved for a linear finite volume variable MooseVariableLinearFVReal. The contributions can be derived using the integral of the diffusion term in the following form:

var element = document.getElementById("moose-equation-5f287c1c-1f6a-434b-8107-cd8a8e6b8fd2");katex.render("-\\int\\limits_{V_C} \\nabla \\cdot D \\nabla u dV = \\sum\\limits_f \\int\\limits_{S_f} -D \\nabla u \\cdot \\vec{n} dS,", element, {displayMode:true,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});

where we used the divergence theorem to transform a volumetric integral over cell $var element = document.getElementById("moose-equation-907c2094-8d62-458d-b688-c0673a886ab2");katex.render("V_C", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ of a vector field to a sum of surface integrals over the faces of the cell. Furthermore, $var element = document.getElementById("moose-equation-bbe1e67a-6a0b-45ae-babc-f491b1c62fc7");katex.render("D", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ denotes a space dependent diffusion coefficient. With this, we can use the finite volume approximation for the face integrals in the following way:

var element = document.getElementById("moose-equation-c5f3307a-6ec7-44af-a4e1-071a446c0865");katex.render("-\\int\\limits_{S_f} D\\nabla u \\cdot \\vec{n} dS \\approx -D_f\\left(|\\vec{\\Delta}|[u_N - u_C] + \\overline{\\nabla u} \\cdot \\vec{k}\\right) |S_f|,", element, {displayMode:true,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});

where $var element = document.getElementById("moose-equation-12885388-3605-48d7-82c0-ff497465298d");katex.render("|S_f|", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ denotes the surface area. Vectors $var element = document.getElementById("moose-equation-fcdbd386-eecd-478d-82b5-acfbb391ff5f");katex.render("\\vec{\\Delta}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ and $var element = document.getElementById("moose-equation-beffccca-674a-41cb-8042-b088806825cf");katex.render("\\vec{k}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ are determined to respect $var element = document.getElementById("moose-equation-a10535ed-4bf2-4d21-96e1-8cdc19c5fbaa");katex.render("\\vec{n} = \\vec{\\Delta} + \\vec{k}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ , where $var element = document.getElementById("moose-equation-bbf8b35e-d440-4421-bb16-3d901ea34423");katex.render("\\vec{\\Delta}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is always parallel to the line connecting the current and neighbor cell centroids. We use the over-relaxed approach for the split of the normal vector, described in Moukalled et al. (2016) and Jasak (1996) in detail. As shown above, using these two vectors, the approximate form of the normal-gradient is typically split into two terms:

$var element = document.getElementById("moose-equation-9e4e4c7b-9143-4cd0-beb3-6e40cc4d57ea");katex.render("|\\vec{\\Delta}|[u_N - u_C]", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ which describes a contribution that comes from a finite difference approximation of the gradient on orthogonal grids. Hence, it is referred to as an orthogonal contribution. For orthogonal meshes, $var element = document.getElementById("moose-equation-4028744b-9835-42de-a583-682513d5ab29");katex.render("|\\Delta|", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is just $var element = document.getElementById("moose-equation-d83aae05-db40-47bc-a7b9-2584eae6702f");katex.render("\\frac{1}{|\\vec{r}_{CN}|}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ where $var element = document.getElementById("moose-equation-42eb19c0-08cd-4e56-a53b-86d82b303975");katex.render("|\\vec{r}_{CN}|", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the distance between the current and neighbor cell centroids. This term contributes a ( $var element = document.getElementById("moose-equation-d5cb3068-7dbf-44ee-abe3-77c585bc8471");katex.render("|\\vec{\\Delta}||S_f|", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ ) to the diagonal and off-diagonal entries of the system matrix with different signs.
On non-orthogonal meshes, besides $var element = document.getElementById("moose-equation-bf824163-6a4c-45df-a61f-7b3dc9df666f");katex.render("|\\vec{\\Delta}| = \\frac{1}{|\\vec{r}_{CN}*\\vec{n}|}", 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 following correction term is needed: $var element = document.getElementById("moose-equation-a8422253-902a-431b-8eb0-46f346cf12f3");katex.render("\\overline{\\nabla u} \\cdot \\vec{k}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ , where $var element = document.getElementById("moose-equation-06341e65-b7c0-4a0b-bd09-f77b07944f87");katex.render("\\overline{\\nabla u}", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ denotes the interpolated gradient at the face center computed using the cell gradients on the current and neighbor cells. This term is treated in an explicit manner meaning that it is added to the right hand side vector of the system.

For more information on the numerical representation of the diffusion term and the different techniques used for applying boundary conditions through this kernel, see Moukalled et al. (2016) and Jasak (1996).

The diffusion coefficient parameter ("diffusion_coeff") accepts anything that supports functor-based evaluations. For more information on functors in MOOSE, see Functor system.

Example input syntax

The input file below shows a pure diffusion problem on a two-dimensional domain.

[LinearFVKernels<<<{"href": "../../syntax/LinearFVKernels/index.html"}>>>]
  [diffusion]
    type = LinearFVDiffusion<<<{"description": "Represents the matrix and right hand side contributions of a diffusion term in a partial differential equation.", "href": "LinearFVDiffusion.html"}>>>
    variable<<<{"description": "The name of the variable whose linear system this object contributes to"}>>> = u
    diffusion_coeff<<<{"description": "The diffusion coefficient. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number."}>>> = coeff_func
    use_nonorthogonal_correction<<<{"description": "If the nonorthogonal correction should be used when computing the normal gradient."}>>> = false
  []
  [source]
    type = LinearFVSource<<<{"description": "Represents the matrix and right hand side contributions of a solution-independent source term in a partial differential equation.", "href": "LinearFVSource.html"}>>>
    variable<<<{"description": "The name of the variable whose linear system this object contributes to"}>>> = u
    source_density<<<{"description": "The source density. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number."}>>> = source_func
  []
[]

Input Parameters

variableThe name of the variable whose linear system this object contributes to
C++ Type:LinearVariableName
Unit:(no unit assumed)
Controllable:No
Description:The name of the variable whose linear system this object contributes to

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
diffusion_coeff1The diffusion coefficient. 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:The diffusion coefficient. A functor is any of the following: a variable, a functor material property, a function, a postprocessor or a number.
force_boundary_executionFalseWhether to force execution of this object on all external boundaries.
Default:False
C++ Type:bool
Controllable:No
Description:Whether to force execution of this object on all external boundaries.
matrix_onlyFalseWhether this object is only doing assembly to matrices (no vectors)
Default:False
C++ Type:bool
Controllable:No
Description:Whether this object is only doing assembly to matrices (no vectors)
use_nonorthogonal_correctionTrueIf the nonorthogonal correction should be used when computing the normal gradient.
Default:True
C++ Type:bool
Controllable:No
Description:If the nonorthogonal correction should be used when computing the normal gradient.

Optional Parameters

absolute_value_vector_tagsThe tags for the vectors this residual object should fill with the absolute value of the residual contribution
C++ Type:std::vector<TagName>
Controllable:No
Description:The tags for the vectors this residual object should fill with the absolute value of the residual contribution
extra_matrix_tagsThe extra tags for the matrices this Kernel should fill
C++ Type:std::vector<TagName>
Controllable:No
Description:The extra tags for the matrices this Kernel should fill
extra_vector_tagsThe extra tags for the vectors this Kernel should fill
C++ Type:std::vector<TagName>
Controllable:No
Description:The extra tags for the vectors this Kernel should fill
matrix_tagssystemThe tag for the matrices this Kernel should fill
Default:system
C++ Type:MultiMooseEnum
Options:nontime, system
Controllable:No
Description:The tag for the matrices this Kernel should fill
vector_tagsrhsThe tag for the vectors this Kernel should fill
Default:rhs
C++ Type:MultiMooseEnum
Options:rhs, time
Controllable:No
Description:The tag for the vectors this Kernel should fill

Contribution To Tagged Field Data 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.
implicitTrueDetermines whether this object is calculated using an implicit or explicit form
Default:True
C++ Type:bool
Controllable:No
Description:Determines whether this object is calculated using an implicit or explicit form
seed0The seed for the master random number generator
Default:0
C++ Type:unsigned int
Controllable:No
Description:The seed for the master random number generator
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

ghost_layers1The number of layers of elements to ghost.
Default:1
C++ Type:unsigned short
Controllable:No
Description:The number of layers of elements to ghost.
use_point_neighborsFalseWhether to use point neighbors, which introduces additional ghosting to that used for simple face neighbors.
Default:False
C++ Type:bool
Controllable:No
Description:Whether to use point neighbors, which introduces additional ghosting to that used for simple face neighbors.

Parallel Ghosting Parameters

Input Files

(test/tests/postprocessors/side_diffusive_flux_integral/side_diffusive_flux_integral_linear_fv.i)
(modules/navier_stokes/test/tests/finite_volume/ins/cht/conjugate_heat_transfer/cht_rob-rob.i)
(test/tests/linearfvkernels/diffusion/diffusion-2d.i)
(modules/navier_stokes/test/tests/finite_volume/ins/cht/bulk_heat_transfer/flow-around-square-linear.i)
(modules/navier_stokes/test/tests/finite_volume/ins/turbulence/channel/linear-segregated-transient/channel_ERCOFTAC.i)
(test/tests/variables/linearfv/diffusion-1d-pp.i)
(test/tests/linearfvbcs/robin/diffusion-1d-robin.i)
(modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/steady-transient-compare/common-blocks.i)
(test/tests/linearfvkernels/diffusion/diffusion-2d-rz.i)
(test/tests/multiapps/linearfv_nonlinearfv/linearfv.i)
(modules/navier_stokes/test/tests/finite_volume/wcns/enthalpy_equation/enthalpy_equation.i)
(test/tests/variables/linearfv/diffusion-1d-aux.i)
(test/tests/linearfvkernels/diffusion/diffusion-2d_neumann.i)
(modules/navier_stokes/test/tests/finite_volume/ins/turbulence/bfs/linear-segregated/BFS_ERCOFTAC.i)
(modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d-scalar/channel.i)
(modules/navier_stokes/test/tests/finite_volume/ins/turbulence/channel/linear-segregated/channel_ERCOFTAC.i)
(test/tests/linearfvkernels/block-restriction/block-restricted-adr.i)
(modules/navier_stokes/test/tests/finite_volume/ins/natural_convection/linear_segregated/2d/diff_heated_cavity_linear_segregated.i)
(modules/navier_stokes/test/tests/finite_volume/ins/cht/bulk_heat_transfer/flow-around-square-linear-fluidonly.i)
(test/tests/multisystem/restore_multiapp/sub.i)
(modules/navier_stokes/test/tests/finite_volume/ins/turbulence/block-restricted/block-ke.i)
(test/tests/time_steppers/iteration_adaptive/adapt_linear_systems.i)
(modules/navier_stokes/test/tests/finite_volume/ins/natural_convection/linear_segregated/2d/diff_heated_cavity_linear_buoyancy.i)
(modules/navier_stokes/test/tests/finite_volume/ins/turbulence/lid-driven/linear-segregated/lid-driven-turb-std-wall.i)
(test/tests/linearfvkernels/block-restriction/block-restricted-diffusion.i)
(modules/navier_stokes/test/tests/finite_volume/ins/turbulence/lid-driven/linear-segregated/lid-driven-turb-non-eq-wall.i)
(test/tests/tag/tag-linearfv.i)
(test/tests/time_integrators/implicit-euler/ie-linearfv.i)
(modules/navier_stokes/test/tests/finite_volume/ins/turbulence/bfs/linear-segregated-transient/BFS_ERCOFTAC.i)
(modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d-heated/fluid.i)
(modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d/2d-boussinesq-transient.i)
(test/tests/linearfvkernels/diffusion-reaction-advection/advection-diffusion-reaction-1d.i)
(test/tests/linearfvkernels/diffusion/diffusion-1d_neumann.i)
(modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d/2d-boussinesq.i)
(test/tests/linearfvkernels/diffusion-reaction-advection/advection-diffusion-reaction-2d.i)
(test/tests/linearfvkernels/diffusion/diffusion-1d.i)
(test/tests/multisystem/picard/linearfv_nonlinearfv/same_input.i)
(test/tests/linearfvbcs/robin/diffusion-2d-robin.i)
(test/tests/multisystem/restore_multiapp/parent.i)
(test/tests/outputs/debug/show_execution_linear_fv_flux.i)
(test/tests/vectorpostprocessors/side_value_sampler/side_value_sampler_fv.i)
(modules/navier_stokes/test/tests/finite_volume/ins/cht/conjugate_heat_transfer/cht_neu-dir.i)
(test/tests/linearfvkernels/block-restriction/block-restricted-diffusion-react.i)
(test/tests/transfers/multiapp_copy_transfer/linear_sys_to_aux/linear_sub.i)

Child Objects

(modules/navier_stokes/include/linearfvkernels/LinearFVTurbulentDiffusion.h)

References

Hrvoje Jasak. Error analysis and estimation for the finite volume method with applications to fluid flows. PhD thesis, Imperial College London (University of London), 1996.

@phdthesis{jasak1996error,
    author = "Jasak, Hrvoje",
    title = "Error analysis and estimation for the finite volume method with applications to fluid flows.",
    year = "1996",
    school = "Imperial College London (University of London)"
}

Fadl Moukalled, L Mangani, Marwan Darwish, and others. The finite volume method in computational fluid dynamics. Volume 6. Springer, 2016.

@book{moukalled2016finite,
    author = "Moukalled, Fadl and Mangani, L and Darwish, Marwan and others",
    title = "The finite volume method in computational fluid dynamics",
    volume = "6",
    year = "2016",
    publisher = "Springer"
}

(test/tests/linearfvkernels/diffusion/diffusion-2d.i)

[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 2
    ny = 1
    ymax = 0.5
  []
[]

[Problem]
  linear_sys_names = 'u_sys'
[]

[Variables]
  [u]
    type = MooseLinearVariableFVReal
    solver_sys = 'u_sys'
    initial_condition = 1.0
  []
[]

[LinearFVKernels]
  [diffusion]
    type = LinearFVDiffusion
    variable = u
    diffusion_coeff = coeff_func
    use_nonorthogonal_correction = false
  []
  [source]
    type = LinearFVSource
    variable = u
    source_density = source_func
  []
[]

[LinearFVBCs]
  [dir]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = u
    boundary = "left right top bottom"
    functor = analytic_solution
  []
[]

[Functions]
  [coeff_func]
    type = ParsedFunction
    expression = '1+0.5*x*y'
  []
  [source_func]
    type = ParsedFunction
    expression = '2*(1.5-y*y)+2*x*y*(1.5-y*y)+2*(1.5-x*x)+2*x*y*(1.5-x*x)'
  []
  [analytic_solution]
    type = ParsedFunction
    expression = '(1.5-x*x)*(1.5-y*y)'
  []
[]

[Postprocessors]
  [h]
    type = AverageElementSize
    execute_on = FINAL
  []
  [error]
    type = ElementL2FunctorError
    approximate = u
    exact = analytic_solution
    execute_on = FINAL
  []
[]

[Convergence]
  [linear]
    type = IterationCountConvergence
    max_iterations = 1
    converge_at_max_iterations = true
  []
[]

[Executioner]
  type = Steady
  system_names = u_sys
  l_tol = 1e-10
  multi_system_fixed_point=true
  multi_system_fixed_point_convergence=linear
  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre boomeramg'
[]

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

(test/tests/postprocessors/side_diffusive_flux_integral/side_diffusive_flux_integral_linear_fv.i)

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 10
  ny = 10
[]

[Problem]
  linear_sys_names = 'u_sys'
[]

[Variables]
  [u]
    solver_sys = 'u_sys'
    type = MooseLinearVariableFVReal
  []
[]

[LinearFVKernels]
  [diff]
    type = LinearFVDiffusion
    variable = u
    diffusion_coeff = diffusivity
  []
[]

[LinearFVBCs]
  [left]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = u
    boundary = left
    functor = 0
  []
  [right]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = u
    boundary = right
    functor = 1
  []
[]

[FunctorMaterials]
  [mat_props]
    type = GenericFunctorMaterial
    prop_names = diffusivity
    prop_values = 1
  []
[]

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

[AuxKernels]
  [grad]
    variable = grad
    functor = u
    type = FunctorElementalGradientAux
  []
[]

[Postprocessors]
  [avg_left]
    type = SideAverageValue
    variable = u
    boundary = left
  []
  [avg_right]
    type = SideAverageValue
    variable = u
    boundary = right
  []
  [avg_flux_left]
    type = SideDiffusiveFluxAverage
    variable = u
    boundary = left
    functor_diffusivity = diffusivity
  []
  [avg_flux_right]
    type = SideDiffusiveFluxAverage
    variable = u
    boundary = right
    functor_diffusivity = diffusivity
  []
[]

[Executioner]
  type = Steady
  system_names = u_sys
  l_tol = 1e-10
  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre boomeramg'
[]

[Outputs]
  exodus = true
[]

(modules/navier_stokes/test/tests/finite_volume/ins/cht/conjugate_heat_transfer/cht_rob-rob.i)

### benchmark sources:
### https://doi.org/10.1016/j.compfluid.2018.06.016
### https://doi.org/10.1016/0017-9310(74)90087-8

b = 0.01 # plate thickness
l = 0.2 # plate length

nxi = 24 # nx in the inlet/entrance region
nyf = 18 # ny in fluid
nxf = 24 # nx in the main fluid region
nys = 8 # ny in the solid domain

fx1_bias = 1.00 # bdry layer bias - fluid
fx2_bias = '${fparse 1.0/1.00}' # bdry layer bias - solid
fy_bias = 1.20 # bdry layer bias - fluid
sy_bias = '${fparse 1.0/1.05}' # bdry layer bias - solid

k_s = 0.2876
rho = 0.3525
mu = 3.95e-5
k = 0.06808
cp = 1142.6

vin = 12.0
Tin = 1000.0
T_s_bottom = 600.0
P_out = 1.03e5

h_s = 1.0
h_f = 1.0

advected_interp_method = 'upwind'

[Mesh]
  [fluid_channel]
    type = GeneratedMeshGenerator
    dim = 2
    nx = ${nxf}
    ny = ${nyf}
    xmin = 0
    xmax = ${l}
    ymin = 0
    ymax = '${fparse 10.0*b}'
    subdomain_ids = '1'
    subdomain_name = 'fluid'
    bias_x = '${fx1_bias}'
    bias_y = '${fparse fy_bias}'
    boundary_name_prefix = 'fluid'
  []
  [solid_base]
    type = GeneratedMeshGenerator
    dim = 2
    nx = ${nxf}
    ny = ${nys}
    xmin = 0
    xmax = ${l}
    ymin = '${fparse -b}'
    ymax = 0
    subdomain_ids = '2'
    subdomain_name = 'solid'
    bias_x = ${fx1_bias}
    bias_y = '${fparse sy_bias}'
    boundary_id_offset = 10
    boundary_name_prefix = 'solid'
  []
  [entrance]
    type = GeneratedMeshGenerator
    dim = 2
    nx = '${fparse 2.0*nxi}'
    ny = ${nyf}
    xmin = '${fparse -2.0*l}'
    xmax = 0
    ymin = 0
    ymax = '${fparse 10.0*b}'
    subdomain_ids = '0'
    subdomain_name = 'entrance'
    bias_x = ${fx2_bias}
    bias_y = '${fparse fy_bias}'
    boundary_id_offset = 20
    boundary_name_prefix = 'ent'
  []
  [smg]
    type = StitchedMeshGenerator
    inputs = 'entrance fluid_channel solid_base'
    stitch_boundaries_pairs = 'ent_right fluid_left;
                              fluid_bottom solid_top'
    prevent_boundary_ids_overlap = false
  []
  [interface]
    type = SideSetsBetweenSubdomainsGenerator
    input = 'smg'
    primary_block = 'fluid'
    paired_block = 'solid'
    new_boundary = interface
  []
  [symmetry_transform]
    type = SymmetryTransformGenerator
    input = interface
    mirror_point = '0 0 0'
    mirror_normal_vector = '0 1 0'
  []
  inactive = 'symmetry_transform'
[]

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

[UserObjects]
  [rc]
    type = RhieChowMassFlux
    u = vel_x
    v = vel_y
    pressure = pressure
    rho = ${rho}
    p_diffusion_kernel = p_diffusion
    block = '0 1'
  []
[]

[Variables]
  [vel_x]
    type = MooseLinearVariableFVReal
    initial_condition = ${vin}
    solver_sys = u_system
    block = '0 1'
  []
  [vel_y]
    type = MooseLinearVariableFVReal
    solver_sys = v_system
    initial_condition = 0.0
    block = '0 1'
  []
  [pressure]
    type = MooseLinearVariableFVReal
    solver_sys = pressure_system
    initial_condition = ${P_out}
    block = '0 1'
  []
  [T_fluid]
    type = MooseLinearVariableFVReal
    solver_sys = energy_system
    initial_condition = ${Tin}
    block = '0 1'
  []
  [T_solid]
    type = MooseLinearVariableFVReal
    solver_sys = solid_energy_system
    initial_condition = ${T_s_bottom}
    block = 2
  []
[]

[LinearFVKernels]
  [u_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_x
    advected_interp_method = ${advected_interp_method}
    mu = ${mu}
    u = vel_x
    v = vel_y
    momentum_component = 'x'
    rhie_chow_user_object = 'rc'
    use_nonorthogonal_correction = false
  []
  [v_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_y
    advected_interp_method = ${advected_interp_method}
    mu = ${mu}
    u = vel_x
    v = vel_y
    momentum_component = 'y'
    rhie_chow_user_object = 'rc'
    use_nonorthogonal_correction = false
  []
  [u_pressure]
    type = LinearFVMomentumPressure
    variable = vel_x
    pressure = pressure
    momentum_component = 'x'
  []
  [v_pressure]
    type = LinearFVMomentumPressure
    variable = vel_y
    pressure = pressure
    momentum_component = 'y'
  []

  [p_diffusion]
    type = LinearFVAnisotropicDiffusion
    variable = pressure
    diffusion_tensor = Ainv
    use_nonorthogonal_correction = false
  []
  [HbyA_divergence]
    type = LinearFVDivergence
    variable = pressure
    face_flux = HbyA
    force_boundary_execution = true
  []

  [h_advection]
    type = LinearFVEnergyAdvection
    variable = T_fluid
    advected_quantity = temperature
    cp = ${cp}
    advected_interp_method = ${advected_interp_method}
    rhie_chow_user_object = 'rc'
  []
  [conduction]
    type = LinearFVDiffusion
    variable = T_fluid
    diffusion_coeff = ${k}
    use_nonorthogonal_correction = false
  []

  [solid-conduction]
    type = LinearFVDiffusion
    variable = T_solid
    diffusion_coeff = ${k_s}
    use_nonorthogonal_correction = false
  []
[]

[LinearFVBCs]
  # velocity BCs
  [inlet-u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'ent_left'
    variable = vel_x
    functor = ${vin}
  []
  [inlet-v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'ent_left'
    variable = vel_y
    functor = '0.000'
  []
  [walls-u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'ent_bottom interface'
    variable = vel_x
    functor = 0.0
  []
  [walls-v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'ent_bottom interface'
    variable = vel_y
    functor = 0.0
  []
  [outlet_p]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'fluid_right'
    variable = pressure
    functor = ${P_out}
  []
  [outlet_u]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = 'fluid_right'
    variable = vel_x
    use_two_term_expansion = false
  []
  [outlet_v]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = 'fluid_right'
    variable = vel_y
    use_two_term_expansion = false
  []

  # freestream BCs for top of fluid domain
  [freestream_u]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = 'fluid_top ent_top'
    variable = vel_x
    use_two_term_expansion = false
  []
  [freestream_v]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = 'fluid_top ent_top'
    variable = vel_y
    use_two_term_expansion = false
  []
  [freestream_p]
    type = LinearFVAdvectionDiffusionFunctorNeumannBC
    boundary = 'fluid_top ent_top'
    variable = pressure
    functor = 0
  []

  # temperature BCs
  [inlet_T]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = T_fluid
    functor = ${Tin}
    boundary = 'ent_left'
  []
  [heated_wall_solid]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = T_solid
    functor = ${T_s_bottom}
    boundary = 'solid_bottom'
  []
  [insulated_fluid]
    type = LinearFVAdvectionDiffusionFunctorNeumannBC
    variable = T_fluid
    functor = 0
    boundary = 'ent_top ent_bottom fluid_top'
  []
  [insulated_solid]
    type = LinearFVAdvectionDiffusionFunctorNeumannBC
    variable = T_solid
    functor = 0
    boundary = 'solid_left solid_right'
  []
  [outlet_T]
    type = LinearFVAdvectionDiffusionOutflowBC
    variable = T_fluid
    use_two_term_expansion = false
    boundary = 'fluid_right'
  []

  [fluid_solid]
    type = LinearFVRobinCHTBC
    variable = T_fluid
    boundary = interface
    h = ${h_f}
    incoming_flux = heat_flux_to_fluid_interface
    surface_temperature = interface_temperature_solid_interface
    thermal_conductivity = ${k}
  []
  [solid_fluid]
    type = LinearFVRobinCHTBC
    variable = T_solid
    boundary = interface
    h = ${h_s}
    incoming_flux = heat_flux_to_solid_interface
    surface_temperature = interface_temperature_fluid_interface
    thermal_conductivity = ${k_s}
  []
[]

[FunctorMaterials]
  [rhocpT]
    property_name = 'rhocpT'
    type = ParsedFunctorMaterial
    functor_names = 'T_fluid'
    expression = '${rho}*${cp}*T_fluid'
  []
[]

[Postprocessors]
  [h_in]
    type = VolumetricFlowRate
    boundary = 'ent_left'
    vel_x = vel_x
    vel_y = vel_y
    rhie_chow_user_object = rc
    advected_quantity = 'rhocpT'
    subtract_mesh_velocity = false
  []
  [h_out]
    type = VolumetricFlowRate
    boundary = 'fluid_right fluid_top ent_top interface'
    vel_x = vel_x
    vel_y = vel_y
    rhie_chow_user_object = rc
    advected_quantity = 'rhocpT'
    advected_interp_method = upwind
    subtract_mesh_velocity = false
  []
[]

[VectorPostprocessors]
  [y_vs_ts]
    type = LineValueSampler
    variable = 'T_solid'
    start_point = '0.05 -1e-9 0' # making sure we are always in the domain
    end_point = '0.05 ${fparse -b+1e-9} 0'
    num_points = 8
    sort_by = id
    warn_discontinuous_face_values = false
  []
  [y_vs_tf]
    type = LineValueSampler
    variable = 'T_fluid'
    start_point = '0.05 1e-9 0' # making sure we are always in the domain
    end_point = '0.05 ${fparse b-1e-9} 0'
    num_points = 12
    sort_by = id
    warn_discontinuous_face_values = false
  []
[]

[Executioner]
  type = SIMPLE

  num_iterations = 1000

  momentum_systems = 'u_system v_system'
  pressure_system = 'pressure_system'
  rhie_chow_user_object = 'rc'
  momentum_l_abs_tol = 1e-10
  pressure_l_abs_tol = 1e-10
  momentum_l_tol = 0
  pressure_l_tol = 0
  momentum_equation_relaxation = 0.9
  pressure_variable_relaxation = 0.3
  momentum_absolute_tolerance = 1e-7
  pressure_absolute_tolerance = 1e-7
  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_system = 'energy_system'
  solid_energy_system = 'solid_energy_system'

  energy_l_abs_tol = 1e-10
  solid_energy_l_abs_tol = 1e-10
  energy_l_tol = 0
  solid_energy_l_tol = 0
  energy_equation_relaxation = 1.0
  energy_absolute_tolerance = 1e-7
  solid_energy_absolute_tolerance = 1e-7
  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'

  cht_interfaces = 'interface'
  cht_solid_flux_relaxation = 1.0
  cht_fluid_flux_relaxation = 1.0
  cht_solid_temperature_relaxation = 1.0
  cht_fluid_temperature_relaxation = 1.0
  max_cht_fpi = 2

  print_fields = false
[]

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

(test/tests/linearfvkernels/diffusion/diffusion-2d.i)

[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 2
    ny = 1
    ymax = 0.5
  []
[]

[Problem]
  linear_sys_names = 'u_sys'
[]

[Variables]
  [u]
    type = MooseLinearVariableFVReal
    solver_sys = 'u_sys'
    initial_condition = 1.0
  []
[]

[LinearFVKernels]
  [diffusion]
    type = LinearFVDiffusion
    variable = u
    diffusion_coeff = coeff_func
    use_nonorthogonal_correction = false
  []
  [source]
    type = LinearFVSource
    variable = u
    source_density = source_func
  []
[]

[LinearFVBCs]
  [dir]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = u
    boundary = "left right top bottom"
    functor = analytic_solution
  []
[]

[Functions]
  [coeff_func]
    type = ParsedFunction
    expression = '1+0.5*x*y'
  []
  [source_func]
    type = ParsedFunction
    expression = '2*(1.5-y*y)+2*x*y*(1.5-y*y)+2*(1.5-x*x)+2*x*y*(1.5-x*x)'
  []
  [analytic_solution]
    type = ParsedFunction
    expression = '(1.5-x*x)*(1.5-y*y)'
  []
[]

[Postprocessors]
  [h]
    type = AverageElementSize
    execute_on = FINAL
  []
  [error]
    type = ElementL2FunctorError
    approximate = u
    exact = analytic_solution
    execute_on = FINAL
  []
[]

[Convergence]
  [linear]
    type = IterationCountConvergence
    max_iterations = 1
    converge_at_max_iterations = true
  []
[]

[Executioner]
  type = Steady
  system_names = u_sys
  l_tol = 1e-10
  multi_system_fixed_point=true
  multi_system_fixed_point_convergence=linear
  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre boomeramg'
[]

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

(modules/navier_stokes/test/tests/finite_volume/ins/cht/bulk_heat_transfer/flow-around-square-linear.i)

mu = 0.01
rho = 1.1
k = 0.0005
cp = 10
k_s = 3.0
h_conv = 5

power_density = 10000

advected_interp_method = 'upwind'

[Mesh]
  [generated_mesh]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 10
    ny = 10
    xmin = 0
    ymin = 0
    ymax = 0.1
    xmax = 0.1
  []
  [subdomain1]
    type = SubdomainBoundingBoxGenerator
    input = generated_mesh
    block_name = subdomain1
    bottom_left = '0.04 0.04 0'
    block_id = 1
    top_right = '0.06 0.06 0'
  []
  [interface]
    type = SideSetsBetweenSubdomainsGenerator
    input = subdomain1
    primary_block = 0
    paired_block = 1
    new_boundary = interface
  []
[]

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

[UserObjects]
  [rc]
    type = RhieChowMassFlux
    u = vel_x
    v = vel_y
    pressure = pressure
    rho = ${rho}
    p_diffusion_kernel = p_diffusion
    block = 0
  []
[]

[Variables]
  [vel_x]
    type = MooseLinearVariableFVReal
    initial_condition = 0.1
    solver_sys = u_system
    block = 0
  []
  [vel_y]
    type = MooseLinearVariableFVReal
    solver_sys = v_system
    initial_condition = 0.01
    block = 0
  []
  [pressure]
    type = MooseLinearVariableFVReal
    solver_sys = pressure_system
    initial_condition = 0.2
    block = 0
  []
  [T_fluid]
    type = MooseLinearVariableFVReal
    solver_sys = energy_system
    initial_condition = 300
    block = 0
  []
  [T_solid]
    type = MooseLinearVariableFVReal
    solver_sys = solid_energy_system
    initial_condition = 500
    block = 1
  []
[]

[LinearFVKernels]
  [u_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_x
    advected_interp_method = ${advected_interp_method}
    mu = ${mu}
    u = vel_x
    v = vel_y
    momentum_component = 'x'
    rhie_chow_user_object = 'rc'
    use_nonorthogonal_correction = true
  []
  [v_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_y
    advected_interp_method = ${advected_interp_method}
    mu = ${mu}
    u = vel_x
    v = vel_y
    momentum_component = 'y'
    rhie_chow_user_object = 'rc'
    use_nonorthogonal_correction = true
  []
  [u_pressure]
    type = LinearFVMomentumPressure
    variable = vel_x
    pressure = pressure
    momentum_component = 'x'
  []
  [v_pressure]
    type = LinearFVMomentumPressure
    variable = vel_y
    pressure = pressure
    momentum_component = 'y'
  []

  [p_diffusion]
    type = LinearFVAnisotropicDiffusion
    variable = pressure
    diffusion_tensor = Ainv
    use_nonorthogonal_correction = true
  []
  [HbyA_divergence]
    type = LinearFVDivergence
    variable = pressure
    face_flux = HbyA
    force_boundary_execution = true
  []

  [h_advection]
    type = LinearFVEnergyAdvection
    variable = T_fluid
    advected_quantity = temperature
    cp = ${cp}
    advected_interp_method = ${advected_interp_method}
    rhie_chow_user_object = 'rc'
  []
  [conduction]
    type = LinearFVDiffusion
    variable = T_fluid
    diffusion_coeff = ${k}
    use_nonorthogonal_correction = true
  []

  [solid-conduction]
    type = LinearFVDiffusion
    variable = T_solid
    diffusion_coeff = ${k_s}
    use_nonorthogonal_correction = true
  []
  [solid-source]
    type = LinearFVSource
    variable = T_solid
    source_density = ${power_density}
  []
[]

[LinearFVBCs]
  [inlet-u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left'
    variable = vel_x
    functor = '0.1'
  []
  [inlet-v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left'
    variable = vel_y
    functor = '0.0'
  []
  [walls-u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'top bottom interface'
    variable = vel_x
    functor = 0.0
  []
  [walls-v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'top bottom interface'
    variable = vel_y
    functor = 0.0
  []
  [outlet_p]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'right'
    variable = pressure
    functor = 1.4
  []
  [outlet_u]
    type = LinearFVAdvectionDiffusionOutflowBC
    variable = vel_x
    use_two_term_expansion = false
    boundary = right
  []
  [outlet_v]
    type = LinearFVAdvectionDiffusionOutflowBC
    variable = vel_y
    use_two_term_expansion = false
    boundary = right
  []

  [inlet_T]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = T_fluid
    functor = 300
    boundary = 'left'
  []
  [walls_T]
    type = LinearFVAdvectionDiffusionFunctorNeumannBC
    variable = T_fluid
    functor = 0.0
    boundary = 'top bottom'
  []
  [outlet_T]
    type = LinearFVAdvectionDiffusionOutflowBC
    variable = T_fluid
    use_two_term_expansion = false
    boundary = right
  []

  [fluid_solid]
    type = LinearFVConvectiveHeatTransferBC
    variable = T_fluid
    T_solid = T_solid
    T_fluid = T_fluid
    boundary = interface
    h = ${h_conv}
  []
  [solid_fluid]
    type = LinearFVConvectiveHeatTransferBC
    variable = T_solid
    T_solid = T_solid
    T_fluid = T_fluid
    boundary = interface
    h = ${h_conv}
  []
[]

[FunctorMaterials]
  [rhocpT]
    property_name = 'rhocpT'
    type = ParsedFunctorMaterial
    functor_names = 'T_fluid'
    expression = '${rho}*${cp}*T_fluid'
  []
[]

[Postprocessors]
  [h_in]
    type = VolumetricFlowRate
    boundary = left
    vel_x = vel_x
    vel_y = vel_y
    rhie_chow_user_object = rc
    advected_quantity = 'rhocpT'
    subtract_mesh_velocity = false
  []
  [h_out]
    type = VolumetricFlowRate
    boundary = right
    vel_x = vel_x
    vel_y = vel_y
    rhie_chow_user_object = rc
    advected_quantity = 'rhocpT'
    advected_interp_method = upwind
    subtract_mesh_velocity = false
  []
  [power]
    type = ElementIntegralFunctorPostprocessor
    functor = ${power_density}
    block = 1
  []
[]

[Executioner]
  type = SIMPLE
  momentum_l_abs_tol = 1e-13
  pressure_l_abs_tol = 1e-13
  energy_l_abs_tol = 1e-13
  solid_energy_l_abs_tol = 1e-13
  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.8
  energy_equation_relaxation = 1.0
  pressure_variable_relaxation = 0.3
  num_iterations = 1000
  pressure_absolute_tolerance = 1e-10
  momentum_absolute_tolerance = 1e-10
  energy_absolute_tolerance = 1e-10
  solid_energy_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'
  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'
  print_fields = false
  continue_on_max_its = true
[]

[Outputs]
  exodus = true
  execute_on = timestep_end
[]

(modules/navier_stokes/test/tests/finite_volume/ins/turbulence/channel/linear-segregated-transient/channel_ERCOFTAC.i)

H = 1 #halfwidth of the channel
L = 100

Re = 13700

rho = 1
bulk_u = 1
mu = '${fparse rho * bulk_u * 2 * H / Re}'

advected_interp_method = 'upwind'

### k-epsilon Closure Parameters ###
sigma_k =1.0
sigma_eps = 1.3
C1_eps = 1.44
C2_eps = 1.92
C_mu = 0.09

### Initial and Boundary Conditions ###
intensity = 0.01
k_init = '${fparse 1.5*(intensity * bulk_u)^2}'
eps_init = '${fparse C_mu^0.75 * k_init^1.5 / (2*H)}'

### Modeling parameters ###
bulk_wall_treatment = false
walls = 'top bottom'
wall_treatment = 'eq_newton'  # Options: eq_newton, eq_incremental, eq_linearized, neq

[Mesh]
  [block_1]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0
    xmax = ${L}
    ymin = 0
    ymax = ${H}
    nx = 4
    ny = 3
    bias_y = 0.7
  []
  [block_2]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0
    xmax = ${L}
    ymin = ${fparse -H}
    ymax = 0
    nx = 4
    ny = 3
    bias_y = ${fparse 1/0.7}
  []
  [smg]
    type = StitchedMeshGenerator
    inputs = 'block_1 block_2'
    clear_stitched_boundary_ids = true
    stitch_boundaries_pairs = 'bottom top'
    merge_boundaries_with_same_name = true
  []
  # Prevent test diffing on distributed parallel element numbering
  allow_renumbering = false
[]

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

[GlobalParams]
  rhie_chow_user_object = 'rc'
  advected_interp_method = ${advected_interp_method}
[]

[UserObjects]
  [rc]
    type = RhieChowMassFlux
    u = vel_x
    v = vel_y
    pressure = pressure
    rho = ${rho}
    p_diffusion_kernel = p_diffusion
  []
[]

[Variables]
  [vel_x]
    type = MooseLinearVariableFVReal
    initial_condition = ${bulk_u}
    solver_sys = u_system
  []
  [vel_y]
    type = MooseLinearVariableFVReal
    initial_condition = 0
    solver_sys = v_system
  []
  [pressure]
    type = MooseLinearVariableFVReal
    initial_condition = 1e-8
    solver_sys = pressure_system
  []
  [TKE]
    type = MooseLinearVariableFVReal
    solver_sys = TKE_system
    initial_condition = ${k_init}
  []
  [TKED]
    type = MooseLinearVariableFVReal
    solver_sys = TKED_system
    initial_condition = ${eps_init}
  []
[]

[LinearFVKernels]
  [u_time]
    type = LinearFVTimeDerivative
    variable = vel_x
    factor = ${rho}
  []
  [u_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_x
    advected_interp_method = ${advected_interp_method}
    mu = 'mu_t'
    u = vel_x
    v = vel_y
    momentum_component = 'x'
    rhie_chow_user_object = 'rc'
    use_nonorthogonal_correction = false
    use_deviatoric_terms = yes
  []
  [u_diffusion]
    type = LinearFVDiffusion
    variable = vel_x
    diffusion_coeff = ${mu}
    use_nonorthogonal_correction = false
  []
  [u_pressure]
    type = LinearFVMomentumPressure
    variable = vel_x
    pressure = pressure
    momentum_component = 'x'
  []

  [v_time]
    type = LinearFVTimeDerivative
    variable = vel_y
    factor = ${rho}
  []
  [v_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_y
    advected_interp_method = ${advected_interp_method}
    mu = 'mu_t'
    u = vel_x
    v = vel_y
    momentum_component = 'y'
    rhie_chow_user_object = 'rc'
    use_nonorthogonal_correction = false
    use_deviatoric_terms = yes
  []
  [v_diffusion]
    type = LinearFVDiffusion
    variable = vel_y
    diffusion_coeff = ${mu}
    use_nonorthogonal_correction = false
  []
  [v_pressure]
    type = LinearFVMomentumPressure
    variable = vel_y
    pressure = pressure
    momentum_component = 'y'
  []

  [p_diffusion]
    type = LinearFVAnisotropicDiffusion
    variable = pressure
    diffusion_tensor = Ainv
    use_nonorthogonal_correction = false
  []
  [HbyA_divergence]
    type = LinearFVDivergence
    variable = pressure
    face_flux = HbyA
    force_boundary_execution = true
  []

  [TKE_time]
    type = LinearFVTimeDerivative
    variable = TKE
    factor = ${rho}
  []
  [TKE_advection]
    type = LinearFVTurbulentAdvection
    variable = TKE
  []
  [TKE_diffusion]
    type = LinearFVTurbulentDiffusion
    variable = TKE
    diffusion_coeff = ${mu}
    use_nonorthogonal_correction = false
  []
  [TKE_turb_diffusion]
    type = LinearFVTurbulentDiffusion
    variable = TKE
    diffusion_coeff = 'mu_t'
    scaling_coeff = ${sigma_k}
    use_nonorthogonal_correction = false
  []
  [TKE_source_sink]
    type = LinearFVTKESourceSink
    variable = TKE
    u = vel_x
    v = vel_y
    epsilon = TKED
    rho = ${rho}
    mu = ${mu}
    mu_t = 'mu_t'
    walls = ${walls}
    wall_treatment = ${wall_treatment}
    C_pl = 1e10
  []

  [TKED_time]
    type = LinearFVTimeDerivative
    variable = TKED
    factor = ${rho}
  []
  [TKED_advection]
    type = LinearFVTurbulentAdvection
    variable = TKED
    walls = ${walls}
  []
  [TKED_diffusion]
    type = LinearFVTurbulentDiffusion
    variable = TKED
    diffusion_coeff = ${mu}
    use_nonorthogonal_correction = false
    walls = ${walls}
  []
  [TKED_turb_diffusion]
    type = LinearFVTurbulentDiffusion
    variable = TKED
    diffusion_coeff = 'mu_t'
    scaling_coeff = ${sigma_eps}
    use_nonorthogonal_correction = false
    walls = ${walls}
  []
  [TKED_source_sink]
    type = LinearFVTKEDSourceSink
    variable = TKED
    u = vel_x
    v = vel_y
    tke = TKE
    rho = ${rho}
    mu = ${mu}
    mu_t = 'mu_t'
    C1_eps = ${C1_eps}
    C2_eps = ${C2_eps}
    walls = ${walls}
    wall_treatment = ${wall_treatment}
    C_pl = 1e10
  []
[]

[LinearFVBCs]
  [inlet-u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left'
    variable = vel_x
    functor = '${bulk_u}'
  []
  [inlet-v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left'
    variable = vel_y
    functor = '0.0'
  []
  [walls-u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'top bottom'
    variable = vel_x
    functor = 0.0
  []
  [walls-v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'top bottom'
    variable = vel_y
    functor = 0.0
  []
  [outlet_u]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = 'right'
    variable = vel_x
    use_two_term_expansion = false
  []
  [outlet_v]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = 'right'
    variable = vel_y
    use_two_term_expansion = false
  []
  [outlet_p]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'right'
    variable = pressure
    functor = 0.0
  []

  [inlet_TKE]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left'
    variable = TKE
    functor = '${k_init}'
  []
  [outlet_TKE]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = 'right'
    variable = TKE
    use_two_term_expansion = false
  []
  [inlet_TKED]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left'
    variable = TKED
    functor = '${eps_init}'
  []
  [outlet_TKED]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = 'right'
    variable = TKED
    use_two_term_expansion = false
  []
  [walls_mu_t]
    type = LinearFVTurbulentViscosityWallFunctionBC
    boundary = 'bottom top'
    variable = 'mu_t'
    u = vel_x
    v = vel_y
    rho = ${rho}
    mu = ${mu}
    tke = TKE
    wall_treatment = ${wall_treatment}
  []
[]

[AuxVariables]
  [mu_t]
    type = MooseLinearVariableFVReal
    initial_condition = '${fparse rho * C_mu * ${k_init}^2 / eps_init}'
  []
  [yplus]
    type = MooseLinearVariableFVReal
  []
  [mu_eff]
    type = MooseLinearVariableFVReal
    initial_condition = '${fparse rho * C_mu * ${k_init}^2 / eps_init}'
  []
[]

[AuxKernels]
  [compute_mu_t]
    type = kEpsilonViscosityAux
    variable = mu_t
    C_mu = ${C_mu}
    tke = TKE
    epsilon = TKED
    mu = ${mu}
    rho = ${rho}
    u = vel_x
    v = vel_y
    bulk_wall_treatment = ${bulk_wall_treatment}
    walls = ${walls}
    wall_treatment = ${wall_treatment}
    execute_on = 'NONLINEAR'
    mu_t_ratio_max = 1e20
  []
  [compute_y_plus]
    type = RANSYPlusAux
    variable = yplus
    tke = TKE
    mu = ${mu}
    rho = ${rho}
    u = vel_x
    v = vel_y
    walls = ${walls}
    wall_treatment = ${wall_treatment}
    execute_on = 'NONLINEAR'
  []
  [compute_mu_eff]
    type = ParsedAux
    variable = 'mu_eff'
    coupled_variables = 'mu_t'
    expression = 'mu_t + ${mu}'
    execute_on = 'NONLINEAR'
  []
[]

[Executioner]
  type = PIMPLE

  rhie_chow_user_object = 'rc'
  momentum_systems = 'u_system v_system'
  pressure_system = 'pressure_system'
  turbulence_systems = 'TKE_system TKED_system'

  momentum_l_abs_tol = 1e-14
  pressure_l_abs_tol = 1e-14
  turbulence_l_abs_tol = 1e-14
  momentum_l_tol = 1e-14
  pressure_l_tol = 1e-14
  turbulence_l_tol = 1e-14

  momentum_equation_relaxation = 0.7
  pressure_variable_relaxation = 0.3
  turbulence_equation_relaxation = '0.25 0.25'
  num_iterations = 1000
  pressure_absolute_tolerance = 1e-12
  momentum_absolute_tolerance = 1e-12
  turbulence_absolute_tolerance = '1e-12 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'
  turbulence_petsc_options_iname = '-pc_type -pc_hypre_type'
  turbulence_petsc_options_value = 'hypre boomeramg'

  print_fields = false
  continue_on_max_its = true
  dt = 1.0
  num_steps = 2
[]

[Outputs]
  csv = true
[]

[VectorPostprocessors]
  [side_bottom]
    type = SideValueSampler
    boundary = 'bottom'
    variable = 'vel_x vel_y pressure TKE TKED'
    sort_by = 'x'
    execute_on = 'timestep_end'
  []
  [side_top]
    type = SideValueSampler
    boundary = 'top'
    variable = 'vel_x vel_y pressure TKE TKED'
    sort_by = 'x'
    execute_on = 'timestep_end'
  []
  [line_center_channel]
    type = LineValueSampler
    start_point = '${fparse 0.125 * L} ${fparse 0.0001} 0'
    end_point = '${fparse 0.875 * L} ${fparse 0.0001} 0'
    num_points = ${Mesh/block_1/nx}
    variable = 'vel_x vel_y pressure TKE TKED'
    sort_by = 'x'
    execute_on = 'timestep_end'
  []
  [line_quarter_radius_channel]
    type = LineValueSampler
    start_point = '${fparse 0.125 * L} ${fparse 0.5 * H} 0'
    end_point = '${fparse 0.875 * L} ${fparse 0.5 * H} 0'
    num_points =  ${Mesh/block_1/nx}
    variable = 'vel_x vel_y pressure TKE TKED'
    sort_by = 'x'
    execute_on = 'timestep_end'
  []
[]

(test/tests/variables/linearfv/diffusion-1d-pp.i)

[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim = 1
    nx = 50
  []
[]

[Problem]
  linear_sys_names = 'u_sys'
[]

[Variables]
  [u]
    type = MooseLinearVariableFVReal
    solver_sys = 'u_sys'
    initial_condition = 1.0
  []
[]

[LinearFVKernels]
  [diffusion]
    type = LinearFVDiffusion
    variable = u
    diffusion_coeff = coeff_func
  []
  [source]
    type = LinearFVSource
    variable = u
    source_density = source_func
  []
[]

[LinearFVBCs]
  [dir]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = u
    boundary = "left right"
    functor = analytic_solution
  []
[]

[Functions]
  [coeff_func]
    type = ParsedFunction
    expression = '0.5*x'
  []
  [source_func]
    type = ParsedFunction
    expression = '2*x'
  []
  [analytic_solution]
    type = ParsedFunction
    expression = '1-x*x'
  []
[]

[Postprocessors]
  [average]
    type = ElementAverageValue
    variable = u
    execute_on = TIMESTEP_END
    outputs = csv
  []
  [min]
    type = ElementExtremeValue
    variable = u
    value_type = min
    execute_on = TIMESTEP_END
    outputs = csv
  []
  [max]
    type = ElementExtremeValue
    variable = u
    value_type = max
    execute_on = TIMESTEP_END
    outputs = csv
  []
  [num_dofs]
    type = NumDOFs
    execute_on = TIMESTEP_END
    outputs = csv
  []
  [elem_value]
    type = ElementalVariableValue
    variable = u
    elementid = 10
    execute_on = TIMESTEP_END
    outputs = csv
  []
  [point_value]
    type = PointValue
    variable = u
    point = '0.33333 0 0'
    execute_on = TIMESTEP_END
    outputs = csv
  []
[]

[VectorPostprocessors]
  [line-sample]
    type = LineValueSampler
    variable = u
    start_point = '0.13333 0 0'
    end_point = '0.766666 0 0'
    num_points = 9
    sort_by = x
    execute_on = TIMESTEP_END
    outputs = vpp_csv
  []
[]

[Executioner]
  type = Steady
  system_names = u_sys
  l_tol = 1e-10
  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre boomeramg'
[]

[Outputs]
  [csv]
    type = CSV
    execute_on = TIMESTEP_END
  []
  [vpp_csv]
    type = CSV
    execute_on = TIMESTEP_END
  []
[]

(test/tests/linearfvbcs/robin/diffusion-1d-robin.i)

##################################################################
k     = 7.0 # diffusion coeff.
amp   = 3.6 # sinusoid amplitude, for u_exact

x_l   = ${fparse 0.0*pi} # domain bound (left)
x_r   = ${fparse 0.9*pi}     # domain bound (right)

alpha = 5.000 # robin BC coeff for gradient term
beta  = 2.000 # robin BC coeff for variable term
gamma = ${fparse (alpha*amp*cos(x_r) ) + (beta*amp*sin(x_r))} # RHS of Robin BC, applied at right boundary

npts = 2
##################################################################

[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim  = 1
    nx   = ${npts}
    xmin = ${x_l}
    xmax = ${x_r}
  []
[]

[Problem]
  linear_sys_names = 'u_sys'
[]

[Variables]
  [u]
    type = MooseLinearVariableFVReal
    solver_sys = 'u_sys'
    initial_condition = 0.0
  []
[]

[Functions]
  [u_exact]
    type = ParsedFunction
    expression = '${amp}*sin(x)'
  []
  [source_fn]
    type = ParsedFunction
    expression = '${fparse k*amp}*sin(x)'
  []
[]

[LinearFVKernels]
  [diffusion]
    type = LinearFVDiffusion
    variable = u
    diffusion_coeff = ${k}
    use_nonorthogonal_correction = False
  []
  [source]
    type = LinearFVSource
    variable = u
    source_density = source_fn
  []
[]

[LinearFVBCs]
  [dir_r]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = u
    boundary = "left"
    functor = 0
  []
  [rob_l]
    type = LinearFVAdvectionDiffusionFunctorRobinBC
    variable = u
    boundary = "right"
    alpha = ${alpha}
    beta  = ${beta}
    gamma = ${gamma}
  []
[]

[Postprocessors]
  [h]
    type = AverageElementSize
    execute_on = FINAL
  []
  [error]
    type = ElementL2FunctorError
    approximate = u
    exact = u_exact
    execute_on = FINAL
  []
[]

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

[Convergence]
  [linear]
    type = IterationCountConvergence
    max_iterations = 4
    converge_at_max_iterations = true
  []
[]

[Executioner]
  type = Steady
  system_names = u_sys
  l_tol = 1e-7
  linear_convergence = linear
  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre boomeramg'
[]

(modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/steady-transient-compare/common-blocks.i)

mu = 2.6
rho = 1.0
advected_interp_method = 'average'
cp = 300
k = 20
alpha_b = 1e-4

[Mesh]
  [mesh]
    type = CartesianMeshGenerator
    dim = 2
    dx = '1.'
    dy = '0.2'
    ix = '10'
    iy = '5'
  []
[]

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

[UserObjects]
  [rc]
    type = RhieChowMassFlux
    u = vel_x
    v = vel_y
    pressure = pressure
    rho = ${rho}
    p_diffusion_kernel = p_diffusion
  []
[]

[Variables]
  [vel_x]
    type = MooseLinearVariableFVReal
    initial_condition = 0.5
    solver_sys = u_system
  []
  [vel_y]
    type = MooseLinearVariableFVReal
    solver_sys = v_system
    initial_condition = 0.0
  []
  [pressure]
    type = MooseLinearVariableFVReal
    solver_sys = pressure_system
    initial_condition = 0.2
  []
  [T]
    type = MooseLinearVariableFVReal
    solver_sys = energy_system
    initial_condition = 300
  []
[]

[LinearFVKernels]
  [u_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_x
    advected_interp_method = ${advected_interp_method}
    mu = ${mu}
    u = vel_x
    v = vel_y
    momentum_component = 'x'
    rhie_chow_user_object = 'rc'
    use_nonorthogonal_correction = false
  []
  [v_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_y
    advected_interp_method = ${advected_interp_method}
    mu = ${mu}
    u = vel_x
    v = vel_y
    momentum_component = 'y'
    rhie_chow_user_object = 'rc'
    use_nonorthogonal_correction = false
  []
  [u_pressure]
    type = LinearFVMomentumPressure
    variable = vel_x
    pressure = pressure
    momentum_component = 'x'
  []
  [v_pressure]
    type = LinearFVMomentumPressure
    variable = vel_y
    pressure = pressure
    momentum_component = 'y'
  []
  [u_boussinesq]
    type = LinearFVMomentumBoussinesq
    variable = vel_x
    rho = ${rho}
    gravity = '0 -9.81 0'
    alpha_name = ${alpha_b}
    ref_temperature = 300.0
    T_fluid = T
    momentum_component = 'x'
  []
  [v_boussinesq]
    type = LinearFVMomentumBoussinesq
    variable = vel_y
    rho = ${rho}
    gravity = '0 -9.81 0'
    alpha_name = ${alpha_b}
    ref_temperature = 300.0
    T_fluid = T
    momentum_component = 'y'
  []

  [p_diffusion]
    type = LinearFVAnisotropicDiffusion
    variable = pressure
    diffusion_tensor = Ainv
    use_nonorthogonal_correction = false
  []
  [HbyA_divergence]
    type = LinearFVDivergence
    variable = pressure
    face_flux = HbyA
    force_boundary_execution = true
  []

  [h_advection]
    type = LinearFVEnergyAdvection
    variable = T
    advected_quantity = temperature
    cp = ${cp}
    advected_interp_method = ${advected_interp_method}
    rhie_chow_user_object = 'rc'
  []
  [conduction]
    type = LinearFVDiffusion
    variable = T
    diffusion_coeff = ${k}
    use_nonorthogonal_correction = false
  []
[]

[FunctorMaterials]
  [constant_functors]
    type = GenericFunctorMaterial
    prop_names = 'cp alpha_b'
    prop_values = '${cp} ${alpha_b}'
  []
[]

[LinearFVBCs]
  [inlet-u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left'
    variable = vel_x
    functor = '1.1'
  []
  [inlet-v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left'
    variable = vel_y
    functor = '0.0'
  []
  [walls-u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'top bottom'
    variable = vel_x
    functor = 0.0
  []
  [walls-v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'top bottom'
    variable = vel_y
    functor = 0.0
  []
  [outlet_p]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'right'
    variable = pressure
    functor = 1.4
  []
  [outlet_u]
    type = LinearFVAdvectionDiffusionOutflowBC
    variable = vel_x
    use_two_term_expansion = false
    boundary = right
  []
  [outlet_v]
    type = LinearFVAdvectionDiffusionOutflowBC
    variable = vel_y
    use_two_term_expansion = false
    boundary = right
  []
  [inlet_top_T]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = T
    functor = ${fparse 300.0}
    boundary = 'left top'
  []
  [bottom_T]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = T
    functor = ${fparse 350.0}
    boundary = bottom
  []
  [outlet_T]
    type = LinearFVAdvectionDiffusionOutflowBC
    variable = T
    use_two_term_expansion = false
    boundary = right
  []
[]

[Outputs]
  exodus = true
  execute_on = FINAL
[]

(test/tests/linearfvkernels/diffusion/diffusion-2d-rz.i)

[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 2
    ny = 1
    ymax = 0.5
  []
  coord_type = RZ
  rz_coord_axis = Y
[]

[Problem]
  linear_sys_names = 'u_sys'
[]

[Variables]
  [u]
    type = MooseLinearVariableFVReal
    solver_sys = 'u_sys'
    initial_condition = 1.0
  []
[]

[LinearFVKernels]
  [diffusion]
    type = LinearFVDiffusion
    variable = u
    diffusion_coeff = coeff_func
    use_nonorthogonal_correction = true
  []
  [source]
    type = LinearFVSource
    variable = u
    source_density = source_func
  []
[]

[LinearFVBCs]
  [dir]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = u
    boundary = "left right top bottom"
    functor = analytic_solution
  []
[]

[Functions]
  [coeff_func]
    type = ParsedFunction
    expression = '1+0.5*x*y'
  []
  [source_func]
    type = ParsedFunction
    expression = '-(-1.0*x^2*y*(1.5 - x^2) + x*(1.5 - x^2)*(-1.0*x*y - 2))/x - (-1.0*x^2*y*(1.5 - y^2) - 4*x*(1.5 - y^2)*(0.5*x*y + 1))/x'
  []
  [analytic_solution]
    type = ParsedFunction
    expression = '(1.5-x*x)*(1.5-y*y)'
  []
[]

[Postprocessors]
  [h]
    type = AverageElementSize
    execute_on = FINAL
  []
  [error]
    type = ElementL2FunctorError
    approximate = u
    exact = analytic_solution
    execute_on = FINAL
  []
[]

[Convergence]
  [linear]
    type = IterationCountConvergence
    max_iterations = 1
    converge_at_max_iterations = true
  []
[]

[Executioner]
  type = Steady
  system_names = u_sys
  l_tol = 1e-10
  multi_system_fixed_point=true
  multi_system_fixed_point_convergence=linear
  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre boomeramg'
[]

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

(test/tests/multiapps/linearfv_nonlinearfv/linearfv.i)

[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim = 1
    nx = 6
  []
[]

[Problem]
  linear_sys_names = 'u_sys'
[]

[Variables]
  [u]
    type = MooseLinearVariableFVReal
    solver_sys = 'u_sys'
    initial_condition = 1.0
  []
[]

[AuxVariables]
  [diff_var]
    type = MooseVariableFVReal
    initial_condition = 2.0
  []
[]

[LinearFVKernels]
  [diffusion]
    type = LinearFVDiffusion
    variable = u
    diffusion_coeff = diff_var
  []
  [source]
    type = LinearFVSource
    variable = u
    source_density = 1
  []
[]

[LinearFVBCs]
  [dir]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = u
    boundary = "left right"
    functor = 1
  []
[]

[MultiApps]
  inactive = 'nonlinear'
  [nonlinear]
    type = FullSolveMultiApp
    input_files = nonlinearfv.i
    execute_on = timestep_begin
    no_restore = true
  []
[]

[Transfers]
  inactive = 'from_nonlinear to_nonlinear'
  [from_nonlinear]
    type = MultiAppCopyTransfer
    from_multi_app = nonlinear
    source_variable = 'v'
    variable = 'diff_var'
    execute_on = timestep_begin
  []
  [to_nonlinear]
    type = MultiAppCopyTransfer
    to_multi_app = nonlinear
    source_variable = 'u'
    variable = 'diff_var'
    execute_on = timestep_begin
  []
[]

[Executioner]
  type = Steady
  system_names = u_sys
  l_tol = 1e-10
  petsc_options_iname = '-pc_type -pc_factor_mat_solver_type'
  petsc_options_value = 'lu mumps'
  fixed_point_rel_tol = 1e-10
[]

[Outputs]
  exodus = true
  execute_on = timestep_end
[]

(modules/navier_stokes/test/tests/finite_volume/wcns/enthalpy_equation/enthalpy_equation.i)

H = 0.015 #halfwidth of the channel, 10 cm of channel height
L = 1
bulk_u = 0.01
p_ref = 101325.0

advected_interp_method = 'upwind'

[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0
    xmax = ${L}
    ymin = -${H}
    ymax = ${H}
    nx = 30
    ny = 15
  []
[]

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

[UserObjects]
  [rc]
    type = RhieChowMassFlux
    u = vel_x
    v = vel_y
    pressure = pressure
    rho = 'rho'
    p_diffusion_kernel = p_diffusion
  []
[]

[Variables]
  [vel_x]
    type = MooseLinearVariableFVReal
    solver_sys = u_system
    initial_condition = ${bulk_u}
  []
  [vel_y]
    type = MooseLinearVariableFVReal
    solver_sys = v_system
    initial_condition = 0
  []
  [pressure]
    type = MooseLinearVariableFVReal
    solver_sys = pressure_system
    initial_condition = ${p_ref}
  []
  [h]
    type = MooseLinearVariableFVReal
    solver_sys = energy_system
    initial_condition = 44000 # 1900 is an approx of cp(T)
  []
[]

[AuxVariables]
  [rho_var]
    type = MooseLinearVariableFVReal
  []
  [cp_var]
    type = MooseLinearVariableFVReal
  []
  [mu_var]
    type = MooseLinearVariableFVReal
  []
  [k_var]
    type = MooseLinearVariableFVReal
  []
  [T]
    type = MooseLinearVariableFVReal
    initial_condition = 777.
  []
[]

[LinearFVKernels]

  [u_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_x
    mu = 'mu'
    momentum_component = 'x'
    use_nonorthogonal_correction = false
    advected_interp_method = ${advected_interp_method}
    rhie_chow_user_object = 'rc'
    u = vel_x
    v = vel_y
  []
  [u_pressure]
    type = LinearFVMomentumPressure
    variable = vel_x
    pressure = pressure
    momentum_component = 'x'
  []

  [v_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_y
    mu = 'mu'
    momentum_component = 'y'
    use_nonorthogonal_correction = false
    advected_interp_method = ${advected_interp_method}
    rhie_chow_user_object = 'rc'
    u = vel_x
    v = vel_y
  []
  [v_pressure]
    type = LinearFVMomentumPressure
    variable = vel_y
    pressure = pressure
    momentum_component = 'y'
  []

  [p_diffusion]
    type = LinearFVAnisotropicDiffusion
    variable = pressure
    diffusion_tensor = Ainv
    use_nonorthogonal_correction = false
  []
  [HbyA_divergence]
    type = LinearFVDivergence
    variable = pressure
    face_flux = HbyA
    force_boundary_execution = true
  []

  [temp_conduction]
    type = LinearFVDiffusion
    diffusion_coeff = 'alpha'
    variable = h
  []
  [temp_advection]
    type = LinearFVEnergyAdvection
    variable = h
    advected_interp_method = ${advected_interp_method}
    rhie_chow_user_object = 'rc'
  []
[]

[LinearFVBCs]
  [inlet_u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left'
    variable = vel_x
    functor = ${bulk_u} #${bulk_u} #'fully_developed_velocity'
  []
  [inlet-v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left'
    variable = vel_y
    functor = 0
  []
  [inlet_h]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = h
    boundary = 'left'
    functor = h_from_p_T # ${fparse 1900.*860.}
  []
  [inlet_T]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = T
    boundary = 'left'
    functor = 860.
  []

  [walls-u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = vel_x
    boundary = 'top bottom'
    functor = 0.
  []
  [walls-v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = vel_y
    boundary = 'top bottom'
    functor = 0.
  []
  [walls_h]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = h
    boundary = 'top bottom'
    functor = h_from_p_T # ${fparse 1900. * 950}
  []
  [walls_T]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = T
    boundary = 'top bottom'
    functor = 950.
  []
  [walls_p]
    type = LinearFVExtrapolatedPressureBC
    boundary = 'top bottom'
    variable = pressure
    use_two_term_expansion = false
  []

  [outlet_p]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'right'
    variable = pressure
    functor = ${p_ref}
  []
  [outlet_h]
    type = LinearFVAdvectionDiffusionOutflowBC
    variable = h
    use_two_term_expansion = false
    boundary = 'right'
  []
  [outlet_u]
    type = LinearFVAdvectionDiffusionOutflowBC
    variable = vel_x
    use_two_term_expansion = false
    boundary = right
  []
  [outlet_v]
    type = LinearFVAdvectionDiffusionOutflowBC
    variable = vel_y
    use_two_term_expansion = false
    boundary = right
  []
[]

[FluidProperties]
  [lead]
    type = LeadFluidProperties
  []
[]

[FunctorMaterials]
  [fluid_props_to_mat_props]
    type = GeneralFunctorFluidProps
    fp = lead
    pressure = ${p_ref}
    T_fluid = 'T'
    speed = 1
    porosity = 1
    characteristic_length = 1
  []
  [alpha]
    type = ADParsedFunctorMaterial
    property_name = 'alpha'
    functor_names = 'k cp'
    expression = 'k/cp'
  []
  [enthalpy_material]
    type = LinearFVEnthalpyFunctorMaterial
    pressure = ${p_ref}
    T_fluid = T
    h = h
    fp = lead
  []
[]

[AuxKernels]
  [rho_out]
    type = FunctorAux
    functor = 'rho'
    variable = 'rho_var'
    execute_on = 'NONLINEAR'
  []
  [cp_out]
    type = FunctorAux
    functor = 'cp'
    variable = 'cp_var'
    execute_on = 'NONLINEAR'
  []
  [mu_out]
    type = FunctorAux
    functor = 'mu'
    variable = 'mu_var'
    execute_on = 'NONLINEAR'
  []
  [k_out]
    type = FunctorAux
    functor = 'k'
    variable = 'k_var'
    execute_on = 'NONLINEAR'
  []
  [T_from_h_functor]
    type = FunctorAux
    functor = 'T_from_p_h'
    variable = 'T'
    execute_on = 'NONLINEAR'
  []
[]

[Executioner]
  type = SIMPLE
  momentum_l_abs_tol = 1e-6
  pressure_l_abs_tol = 1e-6
  energy_l_abs_tol = 1e-8
  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.7
  pressure_variable_relaxation = 0.3
  energy_equation_relaxation = 0.9
  num_iterations = 200
  pressure_absolute_tolerance = 1e-6
  momentum_absolute_tolerance = 1e-6
  energy_absolute_tolerance = 1e-6
  print_fields = false
  momentum_l_max_its = 1000

  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'
  continue_on_max_its = true
[]

[Outputs]
  exodus = true
  execute_on = 'TIMESTEP_BEGIN FINAL'
[]

(test/tests/variables/linearfv/diffusion-1d-aux.i)

[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim = 1
    nx = 10
  []
[]

[Problem]
  linear_sys_names = 'u_sys'
[]

[Variables]
  [u]
    type = MooseLinearVariableFVReal
    solver_sys = 'u_sys'
    initial_condition = 1.0
  []
[]

[AuxVariables]
  [v_volume]
    type = MooseLinearVariableFVReal
    initial_condition = 50
  []
  [v_functor]
    type = MooseLinearVariableFVReal
    initial_condition = 25
  []
  [v_parsed]
    type = MooseLinearVariableFVReal
    initial_condition = 12.5
  []
[]

[AuxKernels]
  [volume]
    type = VolumeAux
    variable = v_volume
  []
  [functor]
    type = FunctorAux
    variable = v_functor
    functor = u
  []
  [parsed]
    type = ParsedAux
    variable = v_parsed
    coupled_variables = 'v_volume v_functor'
    expression = '0.5*v_volume+0.5*v_functor'
  []
[]

[LinearFVKernels]
  [diffusion]
    type = LinearFVDiffusion
    variable = u
    diffusion_coeff = coeff_func
  []
  [source]
    type = LinearFVSource
    variable = u
    source_density = source_func
  []
[]

[LinearFVBCs]
  [dir]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = u
    boundary = "left right"
    functor = analytic_solution
  []
[]

[Functions]
  [coeff_func]
    type = ParsedFunction
    expression = '0.5*x'
  []
  [source_func]
    type = ParsedFunction
    expression = '2*x'
  []
  [analytic_solution]
    type = ParsedFunction
    expression = '1-x*x'
  []
[]

[Executioner]
  type = Steady
  system_names = u_sys
  l_tol = 1e-10
  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre boomeramg'
[]

[Outputs]
  exodus = true
  execute_on = TIMESTEP_END
[]

(test/tests/linearfvkernels/diffusion/diffusion-2d_neumann.i)

[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 2
    ny = 1
    ymax = 0.5
  []
[]

[Problem]
  linear_sys_names = 'u_sys'
[]

[Variables]
  [u]
    type = MooseLinearVariableFVReal
    solver_sys = 'u_sys'
    initial_condition = 1.0
  []
[]

[LinearFVKernels]
  [diffusion]
    type = LinearFVDiffusion
    variable = u
    diffusion_coeff = coeff_func
    use_nonorthogonal_correction = false
  []
  [source]
    type = LinearFVSource
    variable = u
    source_density = source_func
  []
[]

[LinearFVBCs]
  [dir]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = u
    boundary = "left right"
    functor = analytic_solution
  []
  [neu_bottom]
    type = LinearFVAdvectionDiffusionFunctorNeumannBC
    variable = u
    boundary = "bottom"
    functor = analytic_solution_neumann_bottom
    diffusion_coeff = coeff_func
  []
  [neu_top]
    type = LinearFVAdvectionDiffusionFunctorNeumannBC
    variable = u
    boundary = "top"
    functor = analytic_solution_neumann_top
    diffusion_coeff = coeff_func
  []
[]

[Functions]
  [coeff_func]
    type = ParsedFunction
    expression = '1+0.5*x*y'
  []
  [source_func]
    type = ParsedFunction
    expression = '2*(1.5-y*y)+2*x*y*(1.5-y*y)+2*(1.5-x*x)+2*x*y*(1.5-x*x)'
  []
  [analytic_solution]
    type = ParsedFunction
    expression = '(1.5-x*x)*(1.5-y*y)'
  []
  [analytic_solution_neumann_bottom]
    type = ParsedFunction
    expression = '-(1+0.5*x*y)*(1.5-x*x)*(-2*y)'
  []
  [analytic_solution_neumann_top]
    type = ParsedFunction
    expression = '(1+0.5*x*y)*(1.5-x*x)*(-2*y)'
  []
[]

[Postprocessors]
  [h]
    type = AverageElementSize
    execute_on = FINAL
  []
  [error]
    type = ElementL2FunctorError
    approximate = u
    exact = analytic_solution
    execute_on = FINAL
  []
[]

[Convergence]
  [linear]
    type = IterationCountConvergence
    max_iterations = 1
    converge_at_max_iterations = true
  []
[]

[Executioner]
  type = Steady
  system_names = u_sys
  l_tol = 1e-10
  multi_system_fixed_point=true
  multi_system_fixed_point_convergence=linear
  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre boomeramg'
[]

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

(modules/navier_stokes/test/tests/finite_volume/ins/turbulence/bfs/linear-segregated/BFS_ERCOFTAC.i)

Re = 5100

rho = 1.0
bulk_u = 1.0
H = 1.0
mu = '${fparse rho * bulk_u * H / Re}'

advected_interp_method = 'upwind'

### k-epsilon Closure Parameters ###
sigma_k = 1.0
sigma_eps = 1.3
C1_eps = 1.44
C2_eps = 1.92
C_mu = 0.09

### Initial and Boundary Conditions ###
intensity = 0.01
k_init = '${fparse 1.5*(intensity * bulk_u)^2}'
eps_init = '${fparse C_mu^0.75 * k_init^1.5 / H}'

### Modeling parameters ###
bulk_wall_treatment = false
walls = 'bottom wall-side top'
wall_treatment = 'eq_incremental' # Options: eq_newton, eq_incremental, eq_linearized, neq

[Mesh]
  [gen]
    type = CartesianMeshGenerator
    dim = 2
    dx = '${fparse 10.0*H} ${fparse 20.0*H}'
    dy = '${H} ${fparse 5*H}'
    ix = '8 16'
    iy = '2 8'
    subdomain_id = '
                    2 1
                    1 1
                  '
  []
  [corner_walls]
    type = SideSetsBetweenSubdomainsGenerator
    input = gen
    primary_block = '1'
    paired_block = '2'
    new_boundary = 'wall-side'
  []
  [delete_bottom]
    type = BlockDeletionGenerator
    input = corner_walls
    block = '2'
  []
  # Prevent test diffing on distributed parallel element numbering
  allow_renumbering = false
[]

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

[GlobalParams]
  rhie_chow_user_object = 'rc'
  advected_interp_method = ${advected_interp_method}
[]

[UserObjects]
  [rc]
    type = RhieChowMassFlux
    u = vel_x
    v = vel_y
    pressure = pressure
    rho = ${rho}
    p_diffusion_kernel = p_diffusion
  []
[]

[Variables]
  [vel_x]
    type = MooseLinearVariableFVReal
    initial_condition = ${bulk_u}
    solver_sys = u_system
  []
  [vel_y]
    type = MooseLinearVariableFVReal
    initial_condition = 0
    solver_sys = v_system
  []
  [pressure]
    type = MooseLinearVariableFVReal
    initial_condition = 1e-8
    solver_sys = pressure_system
  []
  [TKE]
    type = MooseLinearVariableFVReal
    solver_sys = TKE_system
    initial_condition = ${k_init}
  []
  [TKED]
    type = MooseLinearVariableFVReal
    solver_sys = TKED_system
    initial_condition = ${eps_init}
  []
[]

[LinearFVKernels]
  [u_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_x
    advected_interp_method = ${advected_interp_method}
    mu = 'mu_t'
    u = vel_x
    v = vel_y
    momentum_component = 'x'
    rhie_chow_user_object = 'rc'
    use_nonorthogonal_correction = false
    use_deviatoric_terms = yes
  []
  [u_diffusion]
    type = LinearFVDiffusion
    variable = vel_x
    diffusion_coeff = ${mu}
    use_nonorthogonal_correction = false
  []
  [u_pressure]
    type = LinearFVMomentumPressure
    variable = vel_x
    pressure = pressure
    momentum_component = 'x'
  []

  [v_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_y
    advected_interp_method = ${advected_interp_method}
    mu = 'mu_t'
    u = vel_x
    v = vel_y
    momentum_component = 'y'
    rhie_chow_user_object = 'rc'
    use_nonorthogonal_correction = false
    use_deviatoric_terms = yes
  []
  [v_diffusion]
    type = LinearFVDiffusion
    variable = vel_y
    diffusion_coeff = ${mu}
    use_nonorthogonal_correction = false
  []
  [v_pressure]
    type = LinearFVMomentumPressure
    variable = vel_y
    pressure = pressure
    momentum_component = 'y'
  []

  [p_diffusion]
    type = LinearFVAnisotropicDiffusion
    variable = pressure
    diffusion_tensor = Ainv
    use_nonorthogonal_correction = false
  []
  [HbyA_divergence]
    type = LinearFVDivergence
    variable = pressure
    face_flux = HbyA
    force_boundary_execution = true
  []

  [TKE_advection]
    type = LinearFVTurbulentAdvection
    variable = TKE
  []
  [TKE_diffusion]
    type = LinearFVTurbulentDiffusion
    variable = TKE
    diffusion_coeff = ${mu}
    use_nonorthogonal_correction = false
  []
  [TKE_turb_diffusion]
    type = LinearFVTurbulentDiffusion
    variable = TKE
    diffusion_coeff = 'mu_t'
    scaling_coeff = ${sigma_k}
    use_nonorthogonal_correction = false
  []
  [TKE_source_sink]
    type = LinearFVTKESourceSink
    variable = TKE
    u = vel_x
    v = vel_y
    epsilon = TKED
    rho = ${rho}
    mu = ${mu}
    mu_t = 'mu_t'
    walls = ${walls}
    wall_treatment = ${wall_treatment}
    C_pl = 1e10
  []

  [TKED_advection]
    type = LinearFVTurbulentAdvection
    variable = TKED
    walls = ${walls}
  []
  [TKED_diffusion]
    type = LinearFVTurbulentDiffusion
    variable = TKED
    diffusion_coeff = ${mu}
    use_nonorthogonal_correction = false
    walls = ${walls}
  []
  [TKED_turb_diffusion]
    type = LinearFVTurbulentDiffusion
    variable = TKED
    diffusion_coeff = 'mu_t'
    scaling_coeff = ${sigma_eps}
    use_nonorthogonal_correction = false
    walls = ${walls}
  []
  [TKED_source_sink]
    type = LinearFVTKEDSourceSink
    variable = TKED
    u = vel_x
    v = vel_y
    tke = TKE
    rho = ${rho}
    mu = ${mu}
    mu_t = 'mu_t'
    C1_eps = ${C1_eps}
    C2_eps = ${C2_eps}
    walls = ${walls}
    wall_treatment = ${wall_treatment}
    C_pl = 1e10
  []
[]

[LinearFVBCs]
  [inlet-u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left'
    variable = vel_x
    functor = '${bulk_u}'
  []
  [inlet-v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left'
    variable = vel_y
    functor = '0.0'
  []
  [inlet_TKE]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left'
    variable = TKE
    functor = '${k_init}'
  []
  [inlet_TKED]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left'
    variable = TKED
    functor = '${eps_init}'
  []

  [outlet_p]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'right'
    variable = pressure
    functor = 0.0
  []
  [outlet_u]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = 'right'
    variable = vel_x
    use_two_term_expansion = false
  []
  [outlet_v]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = 'right'
    variable = vel_y
    use_two_term_expansion = false
  []
  [outlet_TKE]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = 'right'
    variable = TKE
    use_two_term_expansion = false
  []
  [outlet_TKED]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = 'right'
    variable = TKED
    use_two_term_expansion = false
  []

  [walls-u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = ${walls}
    variable = vel_x
    functor = 0.0
  []
  [walls-v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = ${walls}
    variable = vel_y
    functor = 0.0
  []
  [walls_mu_t]
    type = LinearFVTurbulentViscosityWallFunctionBC
    boundary = ${walls}
    variable = 'mu_t'
    u = vel_x
    v = vel_y
    rho = ${rho}
    mu = ${mu}
    tke = TKE
    wall_treatment = ${wall_treatment}
  []
[]

[AuxVariables]
  [mu_t]
    type = MooseLinearVariableFVReal
    initial_condition = '${fparse rho * C_mu * ${k_init}^2 / eps_init}'
  []
  [yplus]
    type = MooseLinearVariableFVReal
  []
  [mu_eff]
    type = MooseLinearVariableFVReal
    initial_condition = '${fparse rho * C_mu * ${k_init}^2 / eps_init + mu}'
  []
[]

[AuxKernels]
  [compute_mu_t]
    type = kEpsilonViscosityAux
    variable = mu_t
    C_mu = ${C_mu}
    tke = TKE
    epsilon = TKED
    mu = ${mu}
    rho = ${rho}
    u = vel_x
    v = vel_y
    bulk_wall_treatment = ${bulk_wall_treatment}
    walls = ${walls}
    wall_treatment = ${wall_treatment}
    execute_on = 'NONLINEAR'
    mu_t_ratio_max = 1e20
  []
  [compute_y_plus]
    type = RANSYPlusAux
    variable = yplus
    tke = TKE
    mu = ${mu}
    rho = ${rho}
    u = vel_x
    v = vel_y
    walls = ${walls}
    wall_treatment = ${wall_treatment}
    execute_on = 'NONLINEAR'
  []
  [compute_mu_eff]
    type = ParsedAux
    variable = 'mu_eff'
    coupled_variables = 'mu_t'
    expression = 'mu_t + ${mu}'
    execute_on = 'NONLINEAR'
  []
[]

[Executioner]
  type = SIMPLE

  rhie_chow_user_object = 'rc'
  momentum_systems = 'u_system v_system'
  pressure_system = 'pressure_system'
  turbulence_systems = 'TKE_system TKED_system'

  momentum_l_abs_tol = 1e-14
  pressure_l_abs_tol = 1e-14
  turbulence_l_abs_tol = 1e-14
  momentum_l_tol = 1e-14
  pressure_l_tol = 1e-14
  turbulence_l_tol = 1e-14

  momentum_equation_relaxation = 0.7
  pressure_variable_relaxation = 0.3
  turbulence_equation_relaxation = '0.5 0.5'
  num_iterations = 1000
  pressure_absolute_tolerance = 1e-12
  momentum_absolute_tolerance = 1e-12
  turbulence_absolute_tolerance = '1e-12 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'
  turbulence_petsc_options_iname = '-pc_type -pc_hypre_type'
  turbulence_petsc_options_value = 'hypre boomeramg'

  print_fields = false
  continue_on_max_its = true
[]

[Outputs]
  csv = true
  [console]
    type = Console
    outlier_variable_norms = false
  []
[]


[VectorPostprocessors]
  [side_bottom]
    type = SideValueSampler
    boundary = 'bottom'
    variable = 'vel_x vel_y pressure TKE TKED'
    sort_by = 'x'
    execute_on = 'timestep_end'
  []
  [side_top]
    type = SideValueSampler
    boundary = 'top'
    variable = 'vel_x vel_y pressure TKE TKED'
    sort_by = 'x'
    execute_on = 'timestep_end'
  []
  [line_entry_channel_wall]
    type = LineValueSampler
    start_point = '${fparse 0.5 * H} ${fparse 1.00001 * H} 0'
    end_point = '${fparse 29.5 * H} ${fparse 1.00001 * H} 0'
    num_points = 24
    variable = 'vel_x vel_y pressure TKE TKED'
    sort_by = 'x'
    execute_on = 'timestep_end'
  []
  [line_quarter_entry_channel]
    type = LineValueSampler
    start_point = '${fparse 0.5 * H} ${fparse 2.25001 * H} 0'
    end_point = '${fparse 29.5 * H} ${fparse 2.25001 * H} 0'
    num_points = 24
    variable = 'vel_x vel_y pressure TKE TKED'
    sort_by = 'x'
    execute_on = 'timestep_end'
  []
[]

(modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d-scalar/channel.i)

mu = 2.6
rho = 1.0
advected_interp_method = 'upwind'
k1 = 0.1
k2 = 0.2

[Mesh]
  [mesh]
    type = CartesianMeshGenerator
    dim = 2
    dx = '0.25 0.25'
    dy = '0.2'
    ix = '5 5'
    iy = '5'
    subdomain_id = '0 1'
  []
[]

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

[UserObjects]
  [rc]
    type = RhieChowMassFlux
    u = vel_x
    v = vel_y
    pressure = pressure
    rho = ${rho}
    p_diffusion_kernel = p_diffusion
  []
[]

[Variables]
  [vel_x]
    type = MooseLinearVariableFVReal
    initial_condition = 0.5
    solver_sys = u_system
  []
  [vel_y]
    type = MooseLinearVariableFVReal
    solver_sys = v_system
    initial_condition = 0.0
  []
  [pressure]
    type = MooseLinearVariableFVReal
    solver_sys = pressure_system
    initial_condition = 0.2
  []
  [scalar1]
    type = MooseLinearVariableFVReal
    solver_sys = s1_system
    initial_condition = 1.1
  []
  [scalar2]
    type = MooseLinearVariableFVReal
    solver_sys = s2_system
    initial_condition = 3
  []
[]

[LinearFVKernels]
  [u_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_x
    advected_interp_method = ${advected_interp_method}
    mu = ${mu}
    u = vel_x
    v = vel_y
    momentum_component = 'x'
    rhie_chow_user_object = 'rc'
    use_nonorthogonal_correction = false
  []
  [v_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_y
    advected_interp_method = ${advected_interp_method}
    mu = ${mu}
    u = vel_x
    v = vel_y
    momentum_component = 'y'
    rhie_chow_user_object = 'rc'
    use_nonorthogonal_correction = false
  []
  [u_pressure]
    type = LinearFVMomentumPressure
    variable = vel_x
    pressure = pressure
    momentum_component = 'x'
  []
  [v_pressure]
    type = LinearFVMomentumPressure
    variable = vel_y
    pressure = pressure
    momentum_component = 'y'
  []

  [p_diffusion]
    type = LinearFVAnisotropicDiffusion
    variable = pressure
    diffusion_tensor = Ainv
    use_nonorthogonal_correction = false
  []
  [HbyA_divergence]
    type = LinearFVDivergence
    variable = pressure
    face_flux = HbyA
    force_boundary_execution = true
  []

  [s1_advection]
    type = LinearFVScalarAdvection
    variable = scalar1
    advected_interp_method = ${advected_interp_method}
    rhie_chow_user_object = 'rc'
  []
  [s1_diffusion]
    type = LinearFVDiffusion
    variable = scalar1
    diffusion_coeff = ${k1}
    use_nonorthogonal_correction = false
  []
  [s2_diffusion]
    type = LinearFVScalarAdvection
    variable = scalar2
    advected_interp_method = ${advected_interp_method}
    rhie_chow_user_object = 'rc'
  []
  [s2_conduction]
    type = LinearFVDiffusion
    variable = scalar2
    diffusion_coeff = ${k2}
    use_nonorthogonal_correction = false
  []
[]

[LinearFVBCs]
  [inlet-u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left'
    variable = vel_x
    functor = '1.1'
  []
  [inlet-v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left'
    variable = vel_y
    functor = '0.0'
  []
  [walls-u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'top bottom'
    variable = vel_x
    functor = 0.0
  []
  [walls-v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'top bottom'
    variable = vel_y
    functor = 0.0
  []
  [outlet_p]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'right'
    variable = pressure
    functor = 1.4
  []
  [outlet_u]
    type = LinearFVAdvectionDiffusionOutflowBC
    variable = vel_x
    use_two_term_expansion = false
    boundary = right
  []
  [outlet_v]
    type = LinearFVAdvectionDiffusionOutflowBC
    variable = vel_y
    use_two_term_expansion = false
    boundary = right
  []
  [inlet_wall_scalar1]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = scalar1
    functor = 1
    boundary = 'left'
  []
  [outlet_scalar1]
    type = LinearFVAdvectionDiffusionOutflowBC
    variable = scalar1
    use_two_term_expansion = false
    boundary = right
  []
  [inlet_wall_scalar2]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = scalar2
    functor = 2
    boundary = 'left'
  []
  [outlet_scalar2]
    type = LinearFVAdvectionDiffusionOutflowBC
    variable = scalar2
    use_two_term_expansion = false
    boundary = right
  []
[]

[Executioner]
  type = SIMPLE
  momentum_l_abs_tol = 1e-13
  pressure_l_abs_tol = 1e-13
  passive_scalar_l_abs_tol = 1e-13
  momentum_l_tol = 0
  pressure_l_tol = 0
  passive_scalar_l_tol = 0
  rhie_chow_user_object = 'rc'
  momentum_systems = 'u_system v_system'
  pressure_system = 'pressure_system'
  passive_scalar_systems = 's1_system s2_system'
  momentum_equation_relaxation = 0.8
  passive_scalar_equation_relaxation = '0.9 0.9'
  pressure_variable_relaxation = 0.3
  # We need to converge the problem to show conservation
  num_iterations = 200
  pressure_absolute_tolerance = 1e-10
  momentum_absolute_tolerance = 1e-10
  # The solution being flat, the normalization factor based on fluxes goes to
  # 0. The convergence criteria being multiplied by said factor, we won't do any
  # better than this. For a non-flat solution, use tighter tolerances
  passive_scalar_absolute_tolerance = '1e-2 1e-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'
  passive_scalar_petsc_options_iname = '-pc_type -pc_hypre_type'
  passive_scalar_petsc_options_value = 'hypre boomeramg'
  print_fields = false
  continue_on_max_its = true
[]

[Outputs]
  exodus = true
  execute_on = timestep_end
  hide = 'pressure vel_x vel_y'
[]

(modules/navier_stokes/test/tests/finite_volume/ins/turbulence/channel/linear-segregated/channel_ERCOFTAC.i)

H = 1 #halfwidth of the channel
L = 100

Re = 13700

rho = 1
bulk_u = 1
mu = '${fparse rho * bulk_u * 2 * H / Re}'

advected_interp_method = 'upwind'

### k-epsilon Closure Parameters ###
sigma_k = 1.0
sigma_eps = 1.3
C1_eps = 1.44
C2_eps = 1.92
C_mu = 0.09

### Initial and Boundary Conditions ###
intensity = ${fparse 0.16*Re^(-1./8.)}
k_init = '${fparse 1.5*(intensity * bulk_u)^2}'
eps_init = '${fparse C_mu^0.75 * k_init^1.5 / (2*H)}'

### Modeling parameters ###
bulk_wall_treatment = false
walls = 'top bottom'
wall_treatment = 'eq_newton'  # Options: eq_newton, eq_incremental, eq_linearized, neq

[Mesh]
  [block_1]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0
    xmax = ${L}
    ymin = 0
    ymax = ${H}
    nx = 4
    ny = 4
    bias_y = 0.7
  []
  [block_2]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0
    xmax = ${L}
    ymin = ${fparse -H}
    ymax = 0
    nx = 4
    ny = 4
    bias_y = ${fparse 1/0.7}
  []
  [smg]
    type = StitchedMeshGenerator
    inputs = 'block_1 block_2'
    clear_stitched_boundary_ids = true
    stitch_boundaries_pairs = 'bottom top'
    merge_boundaries_with_same_name = true
  []
  # Prevent test diffing on distributed parallel element numbering
  allow_renumbering = false
[]

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

[GlobalParams]
  rhie_chow_user_object = 'rc'
  advected_interp_method = ${advected_interp_method}
[]

[UserObjects]
  [rc]
    type = RhieChowMassFlux
    u = vel_x
    v = vel_y
    pressure = pressure
    rho = ${rho}
    p_diffusion_kernel = p_diffusion
  []
[]

[Variables]
  [vel_x]
    type = MooseLinearVariableFVReal
    initial_condition = ${bulk_u}
    solver_sys = u_system
  []
  [vel_y]
    type = MooseLinearVariableFVReal
    initial_condition = 0
    solver_sys = v_system
  []
  [pressure]
    type = MooseLinearVariableFVReal
    initial_condition = 1e-8
    solver_sys = pressure_system
  []
  [TKE]
    type = MooseLinearVariableFVReal
    solver_sys = TKE_system
    initial_condition = ${k_init}
  []
  [TKED]
    type = MooseLinearVariableFVReal
    solver_sys = TKED_system
    initial_condition = ${eps_init}
  []
[]

[LinearFVKernels]
  [u_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_x
    advected_interp_method = ${advected_interp_method}
    mu = 'mu_t'
    u = vel_x
    v = vel_y
    momentum_component = 'x'
    rhie_chow_user_object = 'rc'
    use_nonorthogonal_correction = false
    use_deviatoric_terms = yes
  []
  [u_diffusion]
    type = LinearFVDiffusion
    variable = vel_x
    diffusion_coeff = ${mu}
    use_nonorthogonal_correction = false
  []
  [u_pressure]
    type = LinearFVMomentumPressure
    variable = vel_x
    pressure = pressure
    momentum_component = 'x'
  []
  [v_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_y
    advected_interp_method = ${advected_interp_method}
    mu = 'mu_t'
    u = vel_x
    v = vel_y
    momentum_component = 'y'
    rhie_chow_user_object = 'rc'
    use_nonorthogonal_correction = false
    use_deviatoric_terms = yes
  []
  [v_diffusion]
    type = LinearFVDiffusion
    variable = vel_y
    diffusion_coeff = ${mu}
    use_nonorthogonal_correction = false
  []
  [v_pressure]
    type = LinearFVMomentumPressure
    variable = vel_y
    pressure = pressure
    momentum_component = 'y'
  []

  [p_diffusion]
    type = LinearFVAnisotropicDiffusion
    variable = pressure
    diffusion_tensor = Ainv
    use_nonorthogonal_correction = false
  []
  [HbyA_divergence]
    type = LinearFVDivergence
    variable = pressure
    face_flux = HbyA
    force_boundary_execution = true
  []

  [TKE_advection]
    type = LinearFVTurbulentAdvection
    variable = TKE
  []
  [TKE_diffusion]
    type = LinearFVTurbulentDiffusion
    variable = TKE
    diffusion_coeff = ${mu}
    use_nonorthogonal_correction = false
  []
  [TKE_turb_diffusion]
    type = LinearFVTurbulentDiffusion
    variable = TKE
    diffusion_coeff = 'mu_t'
    scaling_coeff = ${sigma_k}
    use_nonorthogonal_correction = false
  []
  [TKE_source_sink]
    type = LinearFVTKESourceSink
    variable = TKE
    u = vel_x
    v = vel_y
    epsilon = TKED
    rho = ${rho}
    mu = ${mu}
    mu_t = 'mu_t'
    walls = ${walls}
    wall_treatment = ${wall_treatment}
    C_pl = 1e10
  []

  [TKED_advection]
    type = LinearFVTurbulentAdvection
    variable = TKED
    walls = ${walls}
  []
  [TKED_diffusion]
    type = LinearFVTurbulentDiffusion
    variable = TKED
    diffusion_coeff = ${mu}
    use_nonorthogonal_correction = false
    walls = ${walls}
  []
  [TKED_turb_diffusion]
    type = LinearFVTurbulentDiffusion
    variable = TKED
    diffusion_coeff = 'mu_t'
    scaling_coeff = ${sigma_eps}
    use_nonorthogonal_correction = false
    walls = ${walls}
  []
  [TKED_source_sink]
    type = LinearFVTKEDSourceSink
    variable = TKED
    u = vel_x
    v = vel_y
    tke = TKE
    rho = ${rho}
    mu = ${mu}
    mu_t = 'mu_t'
    C1_eps = ${C1_eps}
    C2_eps = ${C2_eps}
    walls = ${walls}
    wall_treatment = ${wall_treatment}
    C_pl = 1e10
  []
[]

[LinearFVBCs]
  [inlet-u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left'
    variable = vel_x
    functor = '${bulk_u}'
  []
  [inlet-v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left'
    variable = vel_y
    functor = '0.0'
  []
  [walls-u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'top bottom'
    variable = vel_x
    functor = 0.0
  []
  [walls-v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'top bottom'
    variable = vel_y
    functor = 0.0
  []
  [outlet_u]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = 'right'
    variable = vel_x
    use_two_term_expansion = false
  []
  [outlet_v]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = 'right'
    variable = vel_y
    use_two_term_expansion = false
  []
  [outlet_p]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'right'
    variable = pressure
    functor = 0.0
  []

  [inlet_TKE]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left'
    variable = TKE
    functor = '${k_init}'
  []
  [outlet_TKE]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = 'right'
    variable = TKE
    use_two_term_expansion = false
  []
  [inlet_TKED]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left'
    variable = TKED
    functor = '${eps_init}'
  []
  [outlet_TKED]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = 'right'
    variable = TKED
    use_two_term_expansion = false
  []
  [walls_mu_t]
    type = LinearFVTurbulentViscosityWallFunctionBC
    boundary = 'bottom top'
    variable = 'mu_t'
    u = vel_x
    v = vel_y
    rho = ${rho}
    mu = ${mu}
    tke = TKE
    wall_treatment = ${wall_treatment}
  []
[]

[AuxVariables]
  [mu_t]
    type = MooseLinearVariableFVReal
    initial_condition = '${fparse rho * C_mu * ${k_init}^2 / eps_init}'
  []
  [yplus]
    type = MooseLinearVariableFVReal
  []
[]

[AuxKernels]
  [compute_mu_t]
    type = kEpsilonViscosityAux
    variable = mu_t
    C_mu = ${C_mu}
    tke = TKE
    epsilon = TKED
    mu = ${mu}
    rho = ${rho}
    u = vel_x
    v = vel_y
    bulk_wall_treatment = ${bulk_wall_treatment}
    walls = ${walls}
    wall_treatment = ${wall_treatment}
    execute_on = 'NONLINEAR'
    mu_t_ratio_max = 1e20
  []
  [compute_y_plus]
    type = RANSYPlusAux
    variable = yplus
    tke = TKE
    mu = ${mu}
    rho = ${rho}
    u = vel_x
    v = vel_y
    walls = ${walls}
    wall_treatment = ${wall_treatment}
    execute_on = 'NONLINEAR'
  []
[]

[Executioner]
  type = SIMPLE

  rhie_chow_user_object = 'rc'
  momentum_systems = 'u_system v_system'
  pressure_system = 'pressure_system'
  turbulence_systems = 'TKE_system TKED_system'

  momentum_l_abs_tol = 1e-14
  pressure_l_abs_tol = 1e-14
  turbulence_l_abs_tol = 1e-14
  momentum_l_tol = 1e-14
  pressure_l_tol = 1e-14
  turbulence_l_tol = 1e-14

  momentum_equation_relaxation = 0.7
  pressure_variable_relaxation = 0.3
  turbulence_equation_relaxation = '0.2 0.2'
  turbulence_field_relaxation = '0.2 0.2'
  num_iterations = 1000
  pressure_absolute_tolerance = 1e-12
  momentum_absolute_tolerance = 1e-12
  turbulence_absolute_tolerance = '1e-12 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'
  turbulence_petsc_options_iname = '-pc_type -pc_hypre_type'
  turbulence_petsc_options_value = 'hypre boomeramg'

  print_fields = false
  continue_on_max_its = true
[]

[Outputs]
  csv = true
[]

[AuxVariables]
  [pressure_over_density]
    type = MooseLinearVariableFVReal
    solver_sys = TKE_system
    initial_condition = ${k_init}
  []
[]

[AuxKernels]
  [compute_pressure_over_density]
    type = ParsedAux
    variable = pressure_over_density
    coupled_variables = 'pressure'
    expression = 'pressure/${rho}'
  []
[]

[VectorPostprocessors]
  [side_bottom]
    type = SideValueSampler
    boundary = 'bottom'
    variable = 'vel_x vel_y pressure_over_density TKE TKED'
    sort_by = 'x'
    execute_on = 'timestep_end'
  []
  [side_top]
    type = SideValueSampler
    boundary = 'top'
    variable = 'vel_x vel_y pressure_over_density TKE TKED'
    sort_by = 'x'
    execute_on = 'timestep_end'
  []
  [line_center_channel]
    type = LineValueSampler
    start_point = '${fparse 0.125 * L} ${fparse 0.0001} 0'
    end_point = '${fparse 0.875 * L} ${fparse 0.0001} 0'
    num_points = ${Mesh/block_1/nx}
    variable = 'vel_x vel_y pressure_over_density TKE TKED'
    sort_by = 'x'
    execute_on = 'timestep_end'
  []
  [line_quarter_radius_channel]
    type = LineValueSampler
    start_point = '${fparse 0.125 * L} ${fparse 0.5 * H} 0'
    end_point = '${fparse 0.875 * L} ${fparse 0.5 * H} 0'
    num_points =  ${Mesh/block_1/nx}
    variable = 'vel_x vel_y pressure_over_density TKE TKED'
    sort_by = 'x'
    execute_on = 'timestep_end'
  []
[]

(test/tests/linearfvkernels/block-restriction/block-restricted-adr.i)

[Mesh]
  [cmg]
    type = CartesianMeshGenerator
    dim = 2
    dx = '0.1 1 0.1'
    dy = '0.1 0.5 0.1'
    ix = '1 2 1'
    iy = '1 1 1'
    subdomain_id = '1 1 1 1 2 3 1 1 1'
  []
  [transform]
    type = TransformGenerator
    input = cmg
    transform = TRANSLATE
    vector_value = '-0.1 -0.1 0.0'
  []
  [create_sides]
    type = SideSetsBetweenSubdomainsGenerator
    input = transform
    new_boundary = sides
    primary_block = 2
    paired_block = 1
  []
  [create_outlet]
    type = SideSetsBetweenSubdomainsGenerator
    input = create_sides
    new_boundary = outlet
    primary_block = 2
    paired_block = 3
  []
[]

[Problem]
  linear_sys_names = 'u_sys'
[]

[Variables]
  [u]
    type = MooseLinearVariableFVReal
    solver_sys = 'u_sys'
    initial_condition = 1.0
    block = 2
  []
[]

[LinearFVKernels]
  [diffusion]
    type = LinearFVDiffusion
    variable = u
    diffusion_coeff = diff_coeff_func
    use_nonorthogonal_correction = false
  []
  [advection]
    type = LinearFVAdvection
    variable = u
    velocity = "0.5 0 0"
    advected_interp_method = average
  []
  [reaction]
    type = LinearFVReaction
    variable = u
    coeff = coeff_func
  []
  [source]
    type = LinearFVSource
    variable = u
    source_density = source_func
  []
[]

[LinearFVBCs]
  inactive = "outflow"
  [dir]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = u
    boundary = "sides outlet"
    functor = analytic_solution
  []
  [outflow]
    type = LinearFVAdvectionDiffusionOutflowBC
    variable = u
    boundary = "right"
    use_two_term_expansion = true
  []
[]

[Functions]
  [diff_coeff_func]
    type = ParsedFunction
    expression = '1.0+0.5*x*y'
  []
  [coeff_func]
    type = ParsedFunction
    expression = '1.0+1.0/(1+x*y)'
  []
  [source_func]
    type = ParsedFunction
    expression = '-1.0*x*pi*sin((1/2)*x*pi)*cos(2*y*pi) - 0.25*y*pi*sin(2*y*pi)*cos((1/2)*x*pi) + (1.0 + 1.0/(x*y + 1))*(sin((1/2)*x*pi)*sin(2*y*pi) + 1.5) + (17/4)*pi^2*(0.5*x*y + 1.0)*sin((1/2)*x*pi)*sin(2*y*pi) + 0.25*pi*sin(2*y*pi)*cos((1/2)*x*pi)'
  []
  [analytic_solution]
    type = ParsedFunction
    expression = 'sin((1/2)*x*pi)*sin(2*y*pi) + 1.5'
  []
[]

[Postprocessors]
  [h]
    type = AverageElementSize
    execute_on = FINAL
    block = 2
  []
  [error]
    type = ElementL2FunctorError
    approximate = u
    exact = analytic_solution
    execute_on = FINAL
    block = 2
  []
[]

[Convergence]
  [linear]
    type = IterationCountConvergence
    max_iterations = 1
    converge_at_max_iterations = true
  []
[]

[Executioner]
  type = Steady
  system_names = u_sys
  l_tol = 1e-10
  multi_system_fixed_point=true
  multi_system_fixed_point_convergence=linear
  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre boomeramg'
[]

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

(modules/navier_stokes/test/tests/finite_volume/ins/natural_convection/linear_segregated/2d/diff_heated_cavity_linear_segregated.i)

################################################################################
# MATERIAL PROPERTIES
################################################################################
rho = 3279.
T_0 = 875.0
mu = 1.
k_cond = 38.0
cp = 640.
alpha = 3.26e-5

walls = 'right left top bottom'

[GlobalParams]
  rhie_chow_user_object = 'ins_rhie_chow_interpolator'
  advected_interp_method = 'upwind'
  u = vel_x
  v = vel_y
[]

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

################################################################################
# GEOMETRY
################################################################################

[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0
    xmax = 1
    ymin = 0
    ymax = 1
    nx = 15
    ny = 15
  []
[]

################################################################################
# EQUATIONS: VARIABLES, KERNELS & BCS
################################################################################

[UserObjects]
  [ins_rhie_chow_interpolator]
    type = RhieChowMassFlux
    u = vel_x
    v = vel_y
    pressure = pressure
    rho = ${rho}
    p_diffusion_kernel = p_diffusion
    body_force_kernel_names = "u_buoyancy; v_buoyancy"
  []
[]

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

[LinearFVKernels]
  [u_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_x
    mu = ${mu}
    momentum_component = 'x'
    use_nonorthogonal_correction = false
  []
  [u_pressure]
    type = LinearFVMomentumPressure
    variable = vel_x
    pressure = pressure
    momentum_component = 'x'
  []
  [u_buoyancy]
    type = LinearFVMomentumBoussinesq
    variable = vel_x
    T_fluid = T_fluid
    gravity = '0 -9.8 0'
    rho = ${rho}
    ref_temperature = ${T_0}
    alpha_name = ${alpha}
    momentum_component = 'x'
  []

  [v_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_y
    mu = ${mu}
    momentum_component = 'y'
    use_nonorthogonal_correction = false
  []
  [v_pressure]
    type = LinearFVMomentumPressure
    variable = vel_y
    pressure = pressure
    momentum_component = 'y'
  []
  [v_buoyancy]
    type = LinearFVMomentumBoussinesq
    variable = vel_y
    T_fluid = T_fluid
    gravity = '0 -9.81 0'
    rho = ${rho}
    ref_temperature = ${T_0}
    alpha_name = ${alpha}
    momentum_component = 'y'
  []

  [p_diffusion]
    type = LinearFVAnisotropicDiffusion
    variable = pressure
    diffusion_tensor = Ainv
    use_nonorthogonal_correction = false
  []
  [HbyA_divergence]
    type = LinearFVDivergence
    variable = pressure
    face_flux = HbyA
    force_boundary_execution = false
  []

   ####### FUEL ENERGY EQUATION #######

  [heat_advection]
    type = LinearFVEnergyAdvection
    variable = T_fluid
    advected_quantity = temperature
    cp = ${cp}
  []
  [conduction]
    type = LinearFVDiffusion
    variable = T_fluid
    diffusion_coeff = ${fparse k_cond}
  []
[]


[LinearFVBCs]
  [no-slip-u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = vel_x
    boundary = ${walls}
    functor = 0
  []
  [no-slip-v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = vel_y
    boundary = ${walls}
    functor = 0
  []
  [T_cold]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = T_fluid
    boundary = 'right'
    functor = 870.0
  []
  [T_hot]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = T_fluid
    boundary = 'left'
    functor = 880.0
  []
  [T_all]
    type = LinearFVAdvectionDiffusionExtrapolatedBC
    variable = T_fluid
    boundary = 'top bottom'
    use_two_term_expansion = false
  []
  [pressure]
    type = LinearFVPressureFluxBC
    boundary = 'top bottom left right'
    variable = pressure
    HbyA_flux = HbyA
    Ainv = Ainv
  []
[]

[FunctorMaterials]
  [constant_functors]
    type = GenericFunctorMaterial
    prop_names = 'cp alpha_b'
    prop_values = '${cp} ${alpha}'
  []
[]

################################################################################
# EXECUTION / SOLVE
################################################################################

[Executioner]
  type = SIMPLE
  momentum_l_abs_tol = 1e-11
  pressure_l_abs_tol = 1e-11
  energy_l_abs_tol = 1e-11
  momentum_l_tol = 0
  pressure_l_tol = 0
  energy_l_tol = 0
  rhie_chow_user_object = 'ins_rhie_chow_interpolator'
  momentum_systems = 'u_system v_system'
  pressure_system = 'pressure_system'
  energy_system = 'energy_system'
  momentum_equation_relaxation = 0.7
  pressure_variable_relaxation = 0.3
  energy_equation_relaxation = 0.9
  num_iterations = 1500
  pressure_absolute_tolerance = 1e-8
  momentum_absolute_tolerance = 1e-8
  energy_absolute_tolerance = 1e-8
  print_fields = false
  momentum_l_max_its = 300

  pin_pressure = true
  pressure_pin_value = 0.0
  pressure_pin_point = '0.5 0.0 0.0'

  # momentum_petsc_options = '-ksp_monitor'
  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'

  continue_on_max_its = true
[]

################################################################################
# SIMULATION OUTPUTS
################################################################################

[Outputs]
  #exodus = true
  [out]
    type = Exodus
    file_base = 'diff_heated_cavity_linear_segregated_out'
  []
[]

(modules/navier_stokes/test/tests/finite_volume/ins/cht/bulk_heat_transfer/flow-around-square-linear-fluidonly.i)

mu = 0.01
rho = 1.1
k = 0.0005
cp = 10
h_conv = 5

advected_interp_method = 'upwind'

[Mesh]
  [generated_mesh]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 10
    ny = 10
    xmin = 0
    ymin = 0
    ymax = 0.1
    xmax = 0.1
  []
  [subdomain1]
    type = SubdomainBoundingBoxGenerator
    input = generated_mesh
    block_name = subdomain1
    bottom_left = '0.04 0.04 0'
    block_id = 1
    top_right = '0.06 0.06 0'
  []
  [interface]
    type = SideSetsBetweenSubdomainsGenerator
    input = subdomain1
    primary_block = 0
    paired_block = 1
    new_boundary = interface
  []
  [delete]
    type = BlockDeletionGenerator
    input = interface
    block = 1
  []
[]

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

[UserObjects]
  [rc]
    type = RhieChowMassFlux
    u = vel_x
    v = vel_y
    pressure = pressure
    rho = ${rho}
    p_diffusion_kernel = p_diffusion
    block = 0
  []
[]

[Variables]
  [vel_x]
    type = MooseLinearVariableFVReal
    initial_condition = 0.1
    solver_sys = u_system
    block = 0
  []
  [vel_y]
    type = MooseLinearVariableFVReal
    solver_sys = v_system
    initial_condition = 0.01
    block = 0
  []
  [pressure]
    type = MooseLinearVariableFVReal
    solver_sys = pressure_system
    initial_condition = 0.2
    block = 0
  []
  [T_fluid]
    type = MooseLinearVariableFVReal
    solver_sys = energy_system
    initial_condition = 300
    block = 0
  []
[]

[LinearFVKernels]
  [u_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_x
    advected_interp_method = ${advected_interp_method}
    mu = ${mu}
    u = vel_x
    v = vel_y
    momentum_component = 'x'
    rhie_chow_user_object = 'rc'
    use_nonorthogonal_correction = true
  []
  [v_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_y
    advected_interp_method = ${advected_interp_method}
    mu = ${mu}
    u = vel_x
    v = vel_y
    momentum_component = 'y'
    rhie_chow_user_object = 'rc'
    use_nonorthogonal_correction = true
  []
  [u_pressure]
    type = LinearFVMomentumPressure
    variable = vel_x
    pressure = pressure
    momentum_component = 'x'
  []
  [v_pressure]
    type = LinearFVMomentumPressure
    variable = vel_y
    pressure = pressure
    momentum_component = 'y'
  []

  [p_diffusion]
    type = LinearFVAnisotropicDiffusion
    variable = pressure
    diffusion_tensor = Ainv
    use_nonorthogonal_correction = true
  []
  [HbyA_divergence]
    type = LinearFVDivergence
    variable = pressure
    face_flux = HbyA
    force_boundary_execution = true
  []

  [h_advection]
    type = LinearFVEnergyAdvection
    variable = T_fluid
    advected_quantity = temperature
    cp = ${cp}
    advected_interp_method = ${advected_interp_method}
    rhie_chow_user_object = 'rc'
  []
  [conduction]
    type = LinearFVDiffusion
    variable = T_fluid
    diffusion_coeff = ${k}
    use_nonorthogonal_correction = true
  []
[]

[LinearFVBCs]
  [inlet-u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left'
    variable = vel_x
    functor = '0.1'
  []
  [inlet-v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left'
    variable = vel_y
    functor = '0.0'
  []
  [walls-u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'top bottom interface'
    variable = vel_x
    functor = 0.0
  []
  [walls-v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'top bottom interface'
    variable = vel_y
    functor = 0.0
  []
  [outlet_p]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'right'
    variable = pressure
    functor = 1.4
  []
  [outlet_u]
    type = LinearFVAdvectionDiffusionOutflowBC
    variable = vel_x
    use_two_term_expansion = false
    boundary = right
  []
  [outlet_v]
    type = LinearFVAdvectionDiffusionOutflowBC
    variable = vel_y
    use_two_term_expansion = false
    boundary = right
  []

  [inlet_T]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = T_fluid
    functor = 300
    boundary = 'left'
  []
  [walls_T]
    type = LinearFVAdvectionDiffusionFunctorNeumannBC
    variable = T_fluid
    functor = 0.0
    boundary = 'top bottom'
  []
  [outlet_T]
    type = LinearFVAdvectionDiffusionOutflowBC
    variable = T_fluid
    use_two_term_expansion = false
    boundary = right
  []

  [fluid_solid]
    type = LinearFVConvectiveHeatTransferBC
    variable = T_fluid
    T_solid = boundary_value
    T_fluid = T_fluid
    boundary = interface
    h = ${h_conv}
  []
[]

[FunctorMaterials]
  [rhocpT]
    property_name = 'rhocpT'
    type = ParsedFunctorMaterial
    functor_names = 'T_fluid'
    expression = '${rho}*${cp}*T_fluid'
  []
[]

[Functions]
  [boundary_value]
    type = ConstantFunction
    value = 350
  []
[]

[Executioner]
  type = SIMPLE
  momentum_l_abs_tol = 1e-13
  pressure_l_abs_tol = 1e-13
  energy_l_abs_tol = 1e-13
  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.8
  energy_equation_relaxation = 1.0
  pressure_variable_relaxation = 0.3
  num_iterations = 1000
  pressure_absolute_tolerance = 1e-10
  momentum_absolute_tolerance = 1e-10
  energy_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'
  energy_petsc_options_iname = '-pc_type -pc_hypre_type'
  energy_petsc_options_value = 'hypre boomeramg'
  print_fields = false
  continue_on_max_its = true
[]

[Outputs]
  exodus = true
  execute_on = timestep_end
[]

(test/tests/multisystem/restore_multiapp/sub.i)

[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim = 1
    nx = 6
  []
[]

[Problem]
  nl_sys_names = 'v_sys'
  linear_sys_names = 'u_sys'
[]

[Variables]
  [v]
    type = MooseVariableFVReal
    initial_condition = 2.0
    solver_sys = v_sys
  []
  [u]
    type = MooseLinearVariableFVReal
    solver_sys = 'u_sys'
    initial_condition = 1.0
  []
[]

[FVKernels]
  [diffusion]
    type = FVDiffusion
    variable = v
    coeff = u
  []
  [source]
    type = FVBodyForce
    variable = v
    function = 3
  []
[]

[LinearFVKernels]
  [diffusion]
    type = LinearFVDiffusion
    variable = u
    diffusion_coeff = v
  []
  [source]
    type = LinearFVSource
    variable = u
    source_density = 1
  []
[]

[FVBCs]
  [dir]
    type = FVFunctorDirichletBC
    variable = v
    boundary = "left right"
    functor = 2
  []
[]

[LinearFVBCs]
  [dir]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = u
    boundary = "left right"
    functor = 1
  []
[]

[Convergence]
  [linear]
    type = IterationCountConvergence
    max_iterations = 6
    converge_at_max_iterations = true
  []
[]

[Executioner]
  type = Steady
  #type = Transient
  #steady_state_detection = true
  system_names = 'v_sys u_sys'
  l_abs_tol = 1e-12
  l_tol = 1e-10
  nl_abs_tol = 1e-10
  multi_system_fixed_point=true
  multi_system_fixed_point_convergence=linear
  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre boomeramg'
[]

[Outputs]
  exodus = true
  #execute_on = timestep_end
[]

(modules/navier_stokes/test/tests/finite_volume/ins/turbulence/block-restricted/block-ke.i)

H = 1 #halfwidth of the channel
L = 100

Re = 13700

rho = 1
bulk_u = 1
mu = '${fparse rho * bulk_u * 2 * H / Re}'

advected_interp_method = 'upwind'

### k-epsilon Closure Parameters ###
sigma_k = 1.0
sigma_eps = 1.3
C1_eps = 1.44
C2_eps = 1.92
C_mu = 0.09

### Initial and Boundary Conditions ###
intensity = '${fparse 0.16*Re^(-1./8.)}'
k_init = '${fparse 1.5*(intensity * bulk_u)^2}'
eps_init = '${fparse C_mu^0.75 * k_init^1.5 / (2*H)}'

### Modeling parameters ###
bulk_wall_treatment = false
walls = 'walls'
wall_treatment = 'eq_newton' # Options: eq_newton, eq_incremental, eq_linearized, neq

[Mesh]
  [block_left]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = '${fparse -2*H}'
    xmax = '${fparse -H}'
    ymin = 0
    ymax = ${L}
    nx = 3
    ny = 4
  []
  [block_1]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = '${fparse -H}'
    xmax = 0
    ymin = 0
    ymax = ${L}
    nx = 4
    ny = 4
    bias_x = '${fparse 1/0.7}'
  []
  [block_2]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0
    xmax = ${H}
    ymin = 0
    ymax = ${L}
    nx = 4
    ny = 4
    bias_x = 0.7
  []
  [block_right]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = ${H}
    xmax = '${fparse 2*H}'
    ymin = 0
    ymax = ${L}
    nx = 3
    ny = 4
  []
  [smg]
    type = StitchedMeshGenerator
    inputs = 'block_left block_1 block_2 block_right'
    stitch_boundaries_pairs = 'right left; right left; right left'
  []
  [middle]
    input = smg
    type = SubdomainBoundingBoxGenerator
    block_id = 1
    bottom_left = '${fparse -H} 0 0'
    top_right = '${H} ${L} 50'
  []
  [walls]
    type = SideSetsBetweenSubdomainsGenerator
    primary_block = 1
    paired_block = 0
    new_boundary = walls
    input = middle
  []
  # Prevent test diffing on distributed parallel element numbering
  allow_renumbering = false
[]

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

[GlobalParams]
  rhie_chow_user_object = 'rc'
  advected_interp_method = ${advected_interp_method}
[]

[UserObjects]
  [rc]
    type = RhieChowMassFlux
    u = vel_x
    v = vel_y
    pressure = pressure
    rho = ${rho}
    p_diffusion_kernel = p_diffusion
    block = 1
  []
[]

[Variables]
  [vel_x]
    type = MooseLinearVariableFVReal
    initial_condition = 0.0
    solver_sys = u_system
    block = 1
  []
  [vel_y]
    type = MooseLinearVariableFVReal
    initial_condition = ${bulk_u}
    solver_sys = v_system
    block = 1
  []
  [pressure]
    type = MooseLinearVariableFVReal
    initial_condition = 1e-8
    solver_sys = pressure_system
    block = 1
  []
  [TKE]
    type = MooseLinearVariableFVReal
    solver_sys = TKE_system
    initial_condition = ${k_init}
    block = 1
  []
  [TKED]
    type = MooseLinearVariableFVReal
    solver_sys = TKED_system
    initial_condition = ${eps_init}
    block = 1
  []
[]

[LinearFVKernels]
  [u_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_x
    advected_interp_method = ${advected_interp_method}
    mu = 'mu_t'
    u = vel_x
    v = vel_y
    momentum_component = 'x'
    rhie_chow_user_object = 'rc'
    use_nonorthogonal_correction = false
    use_deviatoric_terms = yes
  []
  [u_diffusion]
    type = LinearFVDiffusion
    variable = vel_x
    diffusion_coeff = ${mu}
    use_nonorthogonal_correction = false
  []
  [u_pressure]
    type = LinearFVMomentumPressure
    variable = vel_x
    pressure = pressure
    momentum_component = 'x'
  []
  [v_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_y
    advected_interp_method = ${advected_interp_method}
    mu = 'mu_t'
    u = vel_x
    v = vel_y
    momentum_component = 'y'
    rhie_chow_user_object = 'rc'
    use_nonorthogonal_correction = false
    use_deviatoric_terms = yes
  []
  [v_diffusion]
    type = LinearFVDiffusion
    variable = vel_y
    diffusion_coeff = ${mu}
    use_nonorthogonal_correction = false
  []
  [v_pressure]
    type = LinearFVMomentumPressure
    variable = vel_y
    pressure = pressure
    momentum_component = 'y'
  []

  [p_diffusion]
    type = LinearFVAnisotropicDiffusion
    variable = pressure
    diffusion_tensor = Ainv
    use_nonorthogonal_correction = false
  []
  [HbyA_divergence]
    type = LinearFVDivergence
    variable = pressure
    face_flux = HbyA
    force_boundary_execution = true
  []

  [TKE_advection]
    type = LinearFVTurbulentAdvection
    variable = TKE
  []
  [TKE_diffusion]
    type = LinearFVTurbulentDiffusion
    variable = TKE
    diffusion_coeff = ${mu}
    use_nonorthogonal_correction = false
  []
  [TKE_turb_diffusion]
    type = LinearFVTurbulentDiffusion
    variable = TKE
    diffusion_coeff = 'mu_t'
    scaling_coeff = ${sigma_k}
    use_nonorthogonal_correction = false
  []
  [TKE_source_sink]
    type = LinearFVTKESourceSink
    variable = TKE
    u = vel_x
    v = vel_y
    epsilon = TKED
    rho = ${rho}
    mu = ${mu}
    mu_t = 'mu_t'
    walls = ${walls}
    wall_treatment = ${wall_treatment}
    C_pl = 1e10
  []

  [TKED_advection]
    type = LinearFVTurbulentAdvection
    variable = TKED
    walls = ${walls}
  []
  [TKED_diffusion]
    type = LinearFVTurbulentDiffusion
    variable = TKED
    diffusion_coeff = ${mu}
    use_nonorthogonal_correction = false
    walls = ${walls}
  []
  [TKED_turb_diffusion]
    type = LinearFVTurbulentDiffusion
    variable = TKED
    diffusion_coeff = 'mu_t'
    scaling_coeff = ${sigma_eps}
    use_nonorthogonal_correction = false
    walls = ${walls}
  []
  [TKED_source_sink]
    type = LinearFVTKEDSourceSink
    variable = TKED
    u = vel_x
    v = vel_y
    tke = TKE
    rho = ${rho}
    mu = ${mu}
    mu_t = 'mu_t'
    C1_eps = ${C1_eps}
    C2_eps = ${C2_eps}
    walls = ${walls}
    wall_treatment = ${wall_treatment}
    C_pl = 1e10
  []
[]

[LinearFVBCs]
  [inlet-u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'bottom'
    variable = vel_x
    functor = 0
  []
  [inlet-v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'bottom'
    variable = vel_y
    functor = '${bulk_u}'
  []
  [walls-u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'walls'
    variable = vel_x
    functor = 0.0
  []
  [walls-v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'walls'
    variable = vel_y
    functor = 0.0
  []
  [outlet_u]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = 'top'
    variable = vel_x
    use_two_term_expansion = false
  []
  [outlet_v]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = 'top'
    variable = vel_y
    use_two_term_expansion = false
  []
  [outlet_p]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'top'
    variable = pressure
    functor = 0.0
  []

  [inlet_TKE]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'bottom'
    variable = TKE
    functor = '${k_init}'
  []
  [outlet_TKE]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = 'top'
    variable = TKE
    use_two_term_expansion = false
  []
  [inlet_TKED]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'bottom'
    variable = TKED
    functor = '${eps_init}'
  []
  [outlet_TKED]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = 'top'
    variable = TKED
    use_two_term_expansion = false
  []
  [walls_mu_t]
    type = LinearFVTurbulentViscosityWallFunctionBC
    boundary = 'walls'
    variable = 'mu_t'
    u = vel_x
    v = vel_y
    rho = ${rho}
    mu = ${mu}
    tke = TKE
    wall_treatment = ${wall_treatment}
  []
[]

[AuxVariables]
  [mu_t]
    type = MooseLinearVariableFVReal
    initial_condition = '${fparse rho * C_mu * ${k_init}^2 / eps_init}'
    block = 1
  []
  [yplus]
    type = MooseLinearVariableFVReal
    block = 1
  []
  [mu_eff]
    type = MooseLinearVariableFVReal
    initial_condition = '${fparse rho * C_mu * ${k_init}^2 / eps_init}'
    block = 1
  []
[]

[AuxKernels]
  [compute_mu_t]
    type = kEpsilonViscosityAux
    variable = mu_t
    C_mu = ${C_mu}
    tke = TKE
    epsilon = TKED
    mu = ${mu}
    rho = ${rho}
    u = vel_x
    v = vel_y
    bulk_wall_treatment = ${bulk_wall_treatment}
    walls = ${walls}
    wall_treatment = ${wall_treatment}
    execute_on = 'NONLINEAR'
    mu_t_ratio_max = 1e20
  []
  [compute_y_plus]
    type = RANSYPlusAux
    variable = yplus
    tke = TKE
    mu = ${mu}
    rho = ${rho}
    u = vel_x
    v = vel_y
    walls = ${walls}
    wall_treatment = ${wall_treatment}
    execute_on = 'NONLINEAR'
  []
  [compute_mu_eff]
    type = ParsedAux
    variable = 'mu_eff'
    coupled_variables = 'mu_t'
    expression = 'mu_t + ${mu}'
    execute_on = 'NONLINEAR'
  []
[]

[Executioner]
  type = SIMPLE

  rhie_chow_user_object = 'rc'
  momentum_systems = 'u_system v_system'
  pressure_system = 'pressure_system'
  turbulence_systems = 'TKE_system TKED_system'

  momentum_l_abs_tol = 1e-14
  pressure_l_abs_tol = 1e-14
  turbulence_l_abs_tol = 1e-14
  momentum_l_tol = 1e-14
  pressure_l_tol = 1e-14
  turbulence_l_tol = 1e-14

  momentum_equation_relaxation = 0.7
  pressure_variable_relaxation = 0.3
  turbulence_equation_relaxation = '0.2 0.2'
  turbulence_field_relaxation = '0.2 0.2'
  num_iterations = 1000
  pressure_absolute_tolerance = 1e-12
  momentum_absolute_tolerance = 1e-12
  turbulence_absolute_tolerance = '1e-12 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'
  turbulence_petsc_options_iname = '-pc_type -pc_hypre_type'
  turbulence_petsc_options_value = 'hypre boomeramg'

  print_fields = false
  continue_on_max_its = true
[]

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

[VectorPostprocessors]
  [line_wall]
    type = LineValueSampler
    start_point = '${fparse 0.99 * H} ${fparse 0.125 * L} 0'
    end_point = '${fparse 0.99 * H} ${fparse 0.875 * L} 0'
    num_points = ${Mesh/block_1/nx}
    variable = 'vel_x vel_y pressure TKE TKED'
    sort_by = 'x'
    execute_on = 'timestep_end'
  []
  [line_center_channel]
    type = LineValueSampler
    start_point = ' ${fparse 0.0001 * H} ${fparse 0.125 * L} 0'
    end_point = '${fparse 0.0001 * H} ${fparse 0.875 * L} 0'
    num_points = ${Mesh/block_1/nx}
    variable = 'vel_x vel_y pressure TKE TKED'
    sort_by = 'x'
    execute_on = 'timestep_end'
  []
  [line_quarter_radius_channel]
    type = LineValueSampler
    start_point = '${fparse 0.51 * H} ${fparse 0.125 * L} 0'
    end_point = '${fparse 0.51 * H} ${fparse 0.875 * L} 0'
    num_points =  ${Mesh/block_1/nx}
    variable = 'vel_x vel_y pressure TKE TKED'
    sort_by = 'x'
    execute_on = 'timestep_end'
  []
[]

(test/tests/time_steppers/iteration_adaptive/adapt_linear_systems.i)

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 50
  ny = 2
  xmax = 5
[]

[Problem]
  linear_sys_names = 'u_sys'
[]

[Variables]
  [u]
    type = MooseLinearVariableFVReal
    solver_sys = 'u_sys'
    initial_condition = 1.0
  []
[]

[LinearFVKernels]
  [time]
    type = LinearFVTimeDerivative
    variable = 'u'
  []
  [diffusion]
    type = LinearFVDiffusion
    variable = u
    diffusion_coeff = 5
  []
  [source]
    type = LinearFVSource
    variable = u
    source_density = 2
  []
[]

[LinearFVBCs]
  [dir]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = u
    boundary = "left right"
    functor = 12
  []
[]

[Executioner]
  type = Transient
  system_names = u_sys
  start_time = 0.0
  end_time = 19
  n_startup_steps = 2
  dtmax = 6.0
  [TimeStepper]
    type = IterationAdaptiveDT
    optimal_iterations = 8
    dt = 1.0
  []
  verbose = true
  linear_convergence = much_logic
[]

[Convergence]
  [force_grow]
    type = IterationCountConvergence
    min_iterations = 0
    max_iterations = 4
    converge_at_max_iterations = true
  []
  [force_shrink]
    type = IterationCountConvergence
    min_iterations = 12
    max_iterations = 13
    converge_at_max_iterations = true
  []
  [much_logic]
    type = ParsedConvergence
    convergence_expression = 'if(time < 5, force_grow, force_shrink)'
    symbol_names = 'time force_grow force_shrink'
    symbol_values = 'time force_grow force_shrink'
  []
[]

[Postprocessors]
  [_dt]
    type = TimestepSize
  []
  [time]
    type = TimePostprocessor
  []
[]

[Outputs]
  csv = true
[]

(modules/navier_stokes/test/tests/finite_volume/ins/natural_convection/linear_segregated/2d/diff_heated_cavity_linear_buoyancy.i)

################################################################################
# MATERIAL PROPERTIES
################################################################################
rho_0 = 3279.
mu = 1.0
k_cond = 38.0
cp = ${fparse 640}
alpha_b = 3.26e-5
T_0 = 875.0

walls = 'right left top bottom'

[GlobalParams]
  rhie_chow_user_object = 'ins_rhie_chow_interpolator'
  advected_interp_method = 'upwind'
  u = vel_x
  v = vel_y
[]

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

################################################################################
# GEOMETRY
################################################################################

[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 1
    xmax = 2.0
    ymin = 1.0
    ymax = 2.0
    nx = 15
    ny = 15
  []
[]

################################################################################
# EQUATIONS: VARIABLES, KERNELS & BCS
################################################################################

[UserObjects]
  [ins_rhie_chow_interpolator]
    type = RhieChowMassFlux
    u = vel_x
    v = vel_y
    pressure = pressure
    rho = 'rho'
    p_diffusion_kernel = p_diffusion
    body_force_kernel_names = 'u_buoyancy; v_buoyancy'
  []
[]

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

[LinearFVKernels]
  [u_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_x
    mu = ${mu}
    momentum_component = 'x'
    use_nonorthogonal_correction = false
  []
  [u_pressure]
    type = LinearFVMomentumPressure
    variable = vel_x
    pressure = pressure
    momentum_component = 'x'
  []
  [u_buoyancy]
    type = LinearFVMomentumBuoyancy
    variable = vel_x
    rho = 'rho'
    reference_rho = ${rho_0}
    gravity = '0 -9.81 0'
    momentum_component = 'x'
  []

  [v_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_y
    mu = ${mu}
    momentum_component = 'y'
    use_nonorthogonal_correction = false
  []
  [v_pressure]
    type = LinearFVMomentumPressure
    variable = vel_y
    pressure = pressure
    momentum_component = 'y'
  []
  [v_buoyancy]
    type = LinearFVMomentumBuoyancy
    variable = vel_y
    rho = 'rho'
    reference_rho = ${rho_0}
    gravity = '0 -9.81 0'
    momentum_component = 'y'
  []

  [p_diffusion]
    type = LinearFVAnisotropicDiffusion
    variable = pressure
    diffusion_tensor = Ainv
    use_nonorthogonal_correction = false
  []
  [HbyA_divergence]
    type = LinearFVDivergence
    variable = pressure
    face_flux = HbyA
    force_boundary_execution = false
  []

   ####### FUEL ENERGY EQUATION #######

  [heat_advection]
    type = LinearFVEnergyAdvection
    variable = T_fluid
    advected_quantity = temperature
    cp = ${cp}
  []
  [conduction]
    type = LinearFVDiffusion
    variable = T_fluid
    diffusion_coeff = ${k_cond}
  []
[]

[LinearFVBCs]
  [no-slip-u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = vel_x
    boundary = ${walls}
    functor = 0
  []
  [no-slip-v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = vel_y
    boundary = ${walls}
    functor = 0
  []
  [T_cold]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = T_fluid
    boundary = 'right'
    functor = 870.0
  []
  [T_hot]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = T_fluid
    boundary = 'left'
    functor = 880.0
  []
  [T_all]
    type = LinearFVAdvectionDiffusionFunctorNeumannBC
    variable = T_fluid
    boundary = 'top bottom'
    functor = 0.0
  []
  [pressure]
    type = LinearFVPressureFluxBC
    boundary = 'top bottom left right'
    variable = pressure
    HbyA_flux = HbyA
    Ainv = Ainv
  []
[]

[FunctorMaterials]
  [rho_function]
    type = ParsedFunctorMaterial
    property_name = 'rho'
    functor_names = 'T_fluid'
    expression = '${rho_0}*(1-${alpha_b}*(T_fluid-${T_0})) '
  []
[]

################################################################################
# EXECUTION / SOLVE
################################################################################

[Executioner]
  type = SIMPLE
  momentum_l_abs_tol = 1e-11
  pressure_l_abs_tol = 1e-11
  energy_l_abs_tol = 1e-11
  momentum_l_tol = 0
  pressure_l_tol = 0
  energy_l_tol = 0
  rhie_chow_user_object = 'ins_rhie_chow_interpolator'
  momentum_systems = 'u_system v_system'
  pressure_system = 'pressure_system'
  energy_system = 'energy_system'
  momentum_equation_relaxation = 0.7
  pressure_variable_relaxation = 0.3
  energy_equation_relaxation = 0.9
  num_iterations = 1500
  pressure_absolute_tolerance = 1e-8
  momentum_absolute_tolerance = 1e-8
  energy_absolute_tolerance = 1e-8
  print_fields = false
  momentum_l_max_its = 300

  pin_pressure = true
  pressure_pin_value = 0.0
  pressure_pin_point = '1.5 1.5 0.0'

  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'

  continue_on_max_its = true
[]

################################################################################
# SIMULATION OUTPUTS
################################################################################

[Outputs]
  exodus = true
[]

(modules/navier_stokes/test/tests/finite_volume/ins/turbulence/lid-driven/linear-segregated/lid-driven-turb-std-wall.i)

### Thermophysical Properties ###
mu = 1e-3
rho = 1.0

### Operation Conditions ###
lid_velocity = 1.0
side_length = 0.1

### Initial Conditions ###
intensity = 0.01
k_init = '${fparse 1.5*(intensity * lid_velocity)^2}'
eps_init = '${fparse C_mu^0.75 * k_init^1.5 / side_length}'

### k-epsilon Closure Parameters ###
sigma_k = 1.0
sigma_eps = 1.3
C1_eps = 1.44
C2_eps = 1.92
C_mu = 0.09

### Modeling parameters ###
bulk_wall_treatment = false
walls = 'left top right bottom'
wall_treatment = 'eq_newton' # Options: eq_newton, eq_incremental, eq_linearized, neq

[GlobalParams]
  rhie_chow_user_object = 'rc'
  advected_interp_method = 'upwind'
[]

[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0
    xmax = ${side_length}
    ymin = 0
    ymax = ${side_length}
    nx = 12
    ny = 12
  []
  # Prevent test diffing on distributed parallel element numbering
  allow_renumbering = false
[]

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

[UserObjects]
  [rc]
    type = RhieChowMassFlux
    u = vel_x
    v = vel_y
    pressure = pressure
    rho = ${rho}
    p_diffusion_kernel = p_diffusion
  []
[]

[Variables]
  [vel_x]
    type = MooseLinearVariableFVReal
    initial_condition = ${lid_velocity}
    solver_sys = u_system
  []
  [vel_y]
    type = MooseLinearVariableFVReal
    initial_condition = 0
    solver_sys = v_system
  []
  [pressure]
    type = MooseLinearVariableFVReal
    initial_condition = 1e-8
    solver_sys = pressure_system
  []
  [TKE]
    type = MooseLinearVariableFVReal
    solver_sys = TKE_system
    initial_condition = ${k_init}
  []
  [TKED]
    type = MooseLinearVariableFVReal
    solver_sys = TKED_system
    initial_condition = ${eps_init}
  []
[]

[LinearFVKernels]
  [u_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_x
    mu = 'mu_t'
    u = vel_x
    v = vel_y
    momentum_component = 'x'
    rhie_chow_user_object = 'rc'
    use_nonorthogonal_correction = true
    use_deviatoric_terms = yes
  []
  [u_diffusion]
    type = LinearFVDiffusion
    variable = vel_x
    diffusion_coeff = ${mu}
    use_nonorthogonal_correction = true
  []
  [u_pressure]
    type = LinearFVMomentumPressure
    variable = vel_x
    pressure = pressure
    momentum_component = 'x'
  []

  [v_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_y
    mu = 'mu_t'
    u = vel_x
    v = vel_y
    momentum_component = 'y'
    rhie_chow_user_object = 'rc'
    use_nonorthogonal_correction = true
    use_deviatoric_terms = yes
  []
  [v_diffusion]
    type = LinearFVDiffusion
    variable = vel_y
    diffusion_coeff = ${mu}
    use_nonorthogonal_correction = true
  []
  [v_pressure]
    type = LinearFVMomentumPressure
    variable = vel_y
    pressure = pressure
    momentum_component = 'y'
  []

  [p_diffusion]
    type = LinearFVAnisotropicDiffusion
    variable = pressure
    diffusion_tensor = Ainv
    use_nonorthogonal_correction = true
  []
  [HbyA_divergence]
    type = LinearFVDivergence
    variable = pressure
    face_flux = HbyA
    force_boundary_execution = true
  []

  [TKE_advection]
    type = LinearFVTurbulentAdvection
    variable = TKE
  []
  [TKE_diffusion]
    type = LinearFVTurbulentDiffusion
    variable = TKE
    diffusion_coeff = ${mu}
    use_nonorthogonal_correction = true
  []
  [TKE_turb_diffusion]
    type = LinearFVTurbulentDiffusion
    variable = TKE
    diffusion_coeff = 'mu_t'
    scaling_coeff = ${sigma_k}
    use_nonorthogonal_correction = true
  []
  [TKE_source_sink]
    type = LinearFVTKESourceSink
    variable = TKE
    u = vel_x
    v = vel_y
    epsilon = TKED
    rho = ${rho}
    mu = ${mu}
    mu_t = 'mu_t'
    walls = ${walls}
    wall_treatment = ${wall_treatment}
    C_pl = 1e10
  []

  [TKED_advection]
    type = LinearFVTurbulentAdvection
    variable = TKED
    walls = ${walls}
  []
  [TKED_diffusion]
    type = LinearFVTurbulentDiffusion
    variable = TKED
    diffusion_coeff = ${mu}
    use_nonorthogonal_correction = true
    walls = ${walls}
  []
  [TKED_turb_diffusion]
    type = LinearFVTurbulentDiffusion
    variable = TKED
    diffusion_coeff = 'mu_t'
    scaling_coeff = ${sigma_eps}
    use_nonorthogonal_correction = true
    walls = ${walls}
  []
  [TKED_source_sink]
    type = LinearFVTKEDSourceSink
    variable = TKED
    u = vel_x
    v = vel_y
    tke = TKE
    rho = ${rho}
    mu = ${mu}
    mu_t = 'mu_t'
    C1_eps = ${C1_eps}
    C2_eps = ${C2_eps}
    walls = ${walls}
    wall_treatment = ${wall_treatment}
    C_pl = 1e10
  []
[]

[LinearFVBCs]
  [top_x]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = vel_x
    boundary = 'top'
    functor = 1
  []
  [no_slip_x]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = vel_x
    boundary = 'left right bottom'
    functor = 0
  []
  [no_slip_y]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = vel_y
    boundary = 'left right top bottom'
    functor = 0
  []
  [pressure-extrapolation]
    type = LinearFVExtrapolatedPressureBC
    boundary = 'left right top bottom'
    variable = pressure
    use_two_term_expansion = true
  []
  [walls_mu_t]
    type = LinearFVTurbulentViscosityWallFunctionBC
    boundary = ${walls}
    variable = 'mu_t'
    u = vel_x
    v = vel_y
    rho = ${rho}
    mu = ${mu}
    tke = TKE
    wall_treatment = ${wall_treatment}
  []
[]

[AuxVariables]
  [mu_t]
    type = MooseLinearVariableFVReal
    initial_condition = '${fparse rho * C_mu * ${k_init}^2 / eps_init}'
  []
  [yplus]
    type = MooseLinearVariableFVReal
  []
  [mu_eff]
    type = MooseLinearVariableFVReal
    initial_condition = '${fparse rho * C_mu * ${k_init}^2 / eps_init}'
  []
[]

[AuxKernels]
  [compute_mu_t]
    type = kEpsilonViscosityAux
    variable = mu_t
    C_mu = ${C_mu}
    tke = TKE
    epsilon = TKED
    mu = ${mu}
    rho = ${rho}
    u = vel_x
    v = vel_y
    bulk_wall_treatment = ${bulk_wall_treatment}
    walls = ${walls}
    wall_treatment = ${wall_treatment}
    execute_on = 'NONLINEAR'
    mu_t_ratio_max = 1e20
  []
  [compute_y_plus]
    type = RANSYPlusAux
    variable = yplus
    tke = TKE
    mu = ${mu}
    rho = ${rho}
    u = vel_x
    v = vel_y
    walls = ${walls}
    wall_treatment = ${wall_treatment}
    execute_on = 'NONLINEAR'
  []
  [compute_mu_eff]
    type = ParsedAux
    variable = 'mu_eff'
    coupled_variables = 'mu_t'
    expression = 'mu_t + ${mu}'
    execute_on = 'NONLINEAR'
  []
[]

[Executioner]
  type = SIMPLE

  rhie_chow_user_object = 'rc'
  momentum_systems = 'u_system v_system'
  pressure_system = 'pressure_system'
  turbulence_systems = 'TKE_system TKED_system'

  momentum_l_abs_tol = 1e-14
  pressure_l_abs_tol = 1e-14
  turbulence_l_abs_tol = 1e-14
  momentum_l_tol = 1e-14
  pressure_l_tol = 1e-14
  turbulence_l_tol = 1e-14

  momentum_equation_relaxation = 0.7
  pressure_variable_relaxation = 0.3
  turbulence_equation_relaxation = '0.5 0.5'
  num_iterations = 1000
  pressure_absolute_tolerance = 1e-12
  momentum_absolute_tolerance = 1e-12
  turbulence_absolute_tolerance = '1e-12 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'
  turbulence_petsc_options_iname = '-pc_type -pc_hypre_type'
  turbulence_petsc_options_value = 'hypre boomeramg'

  print_fields = false
  continue_on_max_its = true

  pin_pressure = true
  pressure_pin_value = 0.0
  pressure_pin_point = '0.01 0.099 0.0'
[]

[Outputs]
  csv = true
  perf_graph = false
  print_nonlinear_residuals = false
  print_linear_residuals = true
[]

[VectorPostprocessors]
  [side_bottom]
    type = SideValueSampler
    boundary = 'bottom'
    variable = 'vel_x vel_y pressure TKE TKED'
    sort_by = 'x'
    execute_on = 'timestep_end'
  []
  [side_top]
    type = SideValueSampler
    boundary = 'top'
    variable = 'vel_x vel_y pressure TKE TKED'
    sort_by = 'x'
    execute_on = 'timestep_end'
  []
  [side_left]
    type = SideValueSampler
    boundary = 'left'
    variable = 'vel_x vel_y pressure TKE TKED'
    sort_by = 'y'
    execute_on = 'timestep_end'
  []
  [side_right]
    type = SideValueSampler
    boundary = 'right'
    variable = 'vel_x vel_y pressure TKE TKED'
    sort_by = 'y'
    execute_on = 'timestep_end'
  []
  [horizontal_center]
    type = LineValueSampler
    start_point = '${fparse 0.01 * side_length} ${fparse 0.499 * side_length} 0'
    end_point = '${fparse 0.99 * side_length} ${fparse 0.499 * side_length} 0'
    num_points = ${Mesh/gen/nx}
    variable = 'vel_x vel_y pressure TKE TKED'
    sort_by = 'x'
    execute_on = 'timestep_end'
  []
  [vertical_center]
    type = LineValueSampler
    start_point = '${fparse 0.499 * side_length} ${fparse 0.01 * side_length} 0'
    end_point = '${fparse 0.499 * side_length} ${fparse 0.99 * side_length} 0'
    num_points =  ${Mesh/gen/ny}
    variable = 'vel_x vel_y pressure TKE TKED'
    sort_by = 'y'
    execute_on = 'timestep_end'
  []
[]

(test/tests/linearfvkernels/block-restriction/block-restricted-diffusion.i)

[Mesh]
  [cmg]
    type = CartesianMeshGenerator
    dim = 2
    dx = '0.1 1 0.1'
    dy = '0.1 0.5 0.1'
    ix = '1 2 1'
    iy = '1 1 1'
    subdomain_id = '1 1 1 1 2 3 1 1 1'
  []
  [transform]
    type = TransformGenerator
    input = cmg
    transform = TRANSLATE
    vector_value = '-0.1 -0.1 0.0'
  []
  [create_sides]
    type = SideSetsBetweenSubdomainsGenerator
    input = transform
    new_boundary = sides
    primary_block = 2
    paired_block = 1
  []
  [create_outlet]
    type = SideSetsBetweenSubdomainsGenerator
    input = create_sides
    new_boundary = outlet
    primary_block = 2
    paired_block = 3
  []
[]

[Problem]
  linear_sys_names = 'u_sys'
[]

[Variables]
  [u]
    type = MooseLinearVariableFVReal
    solver_sys = 'u_sys'
    initial_condition = 1.0
    block = 2
  []
[]

[LinearFVKernels]
  [diffusion]
    type = LinearFVDiffusion
    variable = u
    diffusion_coeff = diff_coeff_func
    use_nonorthogonal_correction = false
  []
  [source]
    type = LinearFVSource
    variable = u
    source_density = source_func
  []
[]

[LinearFVBCs]
  [dir]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = u
    boundary = "sides outlet"
    functor = analytic_solution
  []
[]

[Functions]
  [diff_coeff_func]
    type = ParsedFunction
    expression = '1.0+0.5*x*y'
  []
  [source_func]
    type = ParsedFunction
    expression = '-1.0*x*pi*sin(x*pi)*cos(2*y*pi) - 0.5*y*pi*sin(2*y*pi)*cos(x*pi) + 5*pi^2*(0.5*x*y + 1.0)*sin(x*pi)*sin(2*y*pi)'
  []
  [analytic_solution]
    type = ParsedFunction
    expression = 'sin(x*pi)*sin(2*y*pi) + 1.5'
  []
[]

[Postprocessors]
  [h]
    type = AverageElementSize
    execute_on = FINAL
    block = 2
  []
  [error]
    type = ElementL2FunctorError
    approximate = u
    exact = analytic_solution
    execute_on = FINAL
    block = 2
  []
[]

[Convergence]
  [linear]
    type = IterationCountConvergence
    max_iterations = 1
    converge_at_max_iterations = true
  []
[]

[Executioner]
  type = Steady
  system_names = u_sys
  l_tol = 1e-10
  multi_system_fixed_point=true
  multi_system_fixed_point_convergence=linear
  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre boomeramg'
[]

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

(modules/navier_stokes/test/tests/finite_volume/ins/turbulence/lid-driven/linear-segregated/lid-driven-turb-non-eq-wall.i)

### Thermophysical Properties ###
mu = 1e-3
rho = 1.0

### Operation Conditions ###
lid_velocity = 1.0
side_length = 0.1

### Initial Conditions ###
intensity = 0.01
k_init = '${fparse 1.5*(intensity * lid_velocity)^2}'
eps_init = '${fparse C_mu^0.75 * k_init^1.5 / side_length}'

### k-epsilon Closure Parameters ###
sigma_k = 1.0
sigma_eps = 1.3
C1_eps = 1.44
C2_eps = 1.92
C_mu = 0.09

### Modeling parameters ###
bulk_wall_treatment = false
walls = 'left top right bottom'
wall_treatment = 'neq' # Options: eq_newton, eq_incremental, eq_linearized, neq

[GlobalParams]
  rhie_chow_user_object = 'rc'
  advected_interp_method = 'upwind'
[]

[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0
    xmax = ${side_length}
    ymin = 0
    ymax = ${side_length}
    nx = 12
    ny = 12
  []
  # Prevent test diffing on distributed parallel element numbering
  allow_renumbering = false
[]

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

[UserObjects]
  [rc]
    type = RhieChowMassFlux
    u = vel_x
    v = vel_y
    pressure = pressure
    rho = ${rho}
    p_diffusion_kernel = p_diffusion
  []
[]

[Variables]
  [vel_x]
    type = MooseLinearVariableFVReal
    initial_condition = ${lid_velocity}
    solver_sys = u_system
  []
  [vel_y]
    type = MooseLinearVariableFVReal
    initial_condition = 0
    solver_sys = v_system
  []
  [pressure]
    type = MooseLinearVariableFVReal
    initial_condition = 1e-8
    solver_sys = pressure_system
  []
  [TKE]
    type = MooseLinearVariableFVReal
    solver_sys = TKE_system
    initial_condition = ${k_init}
  []
  [TKED]
    type = MooseLinearVariableFVReal
    solver_sys = TKED_system
    initial_condition = ${eps_init}
  []
[]

[LinearFVKernels]
  [u_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_x
    mu = 'mu_t'
    u = vel_x
    v = vel_y
    momentum_component = 'x'
    rhie_chow_user_object = 'rc'
    use_nonorthogonal_correction = false
    use_deviatoric_terms = yes
  []
  [u_diffusion]
    type = LinearFVDiffusion
    variable = vel_x
    diffusion_coeff = ${mu}
    use_nonorthogonal_correction = false
  []
  [u_pressure]
    type = LinearFVMomentumPressure
    variable = vel_x
    pressure = pressure
    momentum_component = 'x'
  []

  [v_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_y
    mu = 'mu_t'
    u = vel_x
    v = vel_y
    momentum_component = 'y'
    rhie_chow_user_object = 'rc'
    use_nonorthogonal_correction = false
    use_deviatoric_terms = yes
  []
  [v_diffusion]
    type = LinearFVDiffusion
    variable = vel_y
    diffusion_coeff = ${mu}
    use_nonorthogonal_correction = false
  []
  [v_pressure]
    type = LinearFVMomentumPressure
    variable = vel_y
    pressure = pressure
    momentum_component = 'y'
  []

  [p_diffusion]
    type = LinearFVAnisotropicDiffusion
    variable = pressure
    diffusion_tensor = Ainv
    use_nonorthogonal_correction = false
  []
  [HbyA_divergence]
    type = LinearFVDivergence
    variable = pressure
    face_flux = HbyA
    force_boundary_execution = true
  []

  [TKE_advection]
    type = LinearFVTurbulentAdvection
    variable = TKE
  []
  [TKE_diffusion]
    type = LinearFVTurbulentDiffusion
    variable = TKE
    diffusion_coeff = ${mu}
    use_nonorthogonal_correction = false
  []
  [TKE_turb_diffusion]
    type = LinearFVTurbulentDiffusion
    variable = TKE
    diffusion_coeff = 'mu_t'
    scaling_coeff = ${sigma_k}
    use_nonorthogonal_correction = false
  []
  [TKE_source_sink]
    type = LinearFVTKESourceSink
    variable = TKE
    u = vel_x
    v = vel_y
    epsilon = TKED
    rho = ${rho}
    mu = ${mu}
    mu_t = 'mu_t'
    walls = ${walls}
    wall_treatment = ${wall_treatment}
    C_pl = 1e10
  []

  [TKED_advection]
    type = LinearFVTurbulentAdvection
    variable = TKED
    walls = ${walls}
  []
  [TKED_diffusion]
    type = LinearFVTurbulentDiffusion
    variable = TKED
    diffusion_coeff = ${mu}
    use_nonorthogonal_correction = false
    walls = ${walls}
  []
  [TKED_turb_diffusion]
    type = LinearFVTurbulentDiffusion
    variable = TKED
    diffusion_coeff = 'mu_t'
    scaling_coeff = ${sigma_eps}
    use_nonorthogonal_correction = false
    walls = ${walls}
  []
  [TKED_source_sink]
    type = LinearFVTKEDSourceSink
    variable = TKED
    u = vel_x
    v = vel_y
    tke = TKE
    rho = ${rho}
    mu = ${mu}
    mu_t = 'mu_t'
    C1_eps = ${C1_eps}
    C2_eps = ${C2_eps}
    walls = ${walls}
    wall_treatment = ${wall_treatment}
    C_pl = 1e10
  []
[]

[LinearFVBCs]
  [top_x]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = vel_x
    boundary = 'top'
    functor = 1
  []
  [no_slip_x]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = vel_x
    boundary = 'left right bottom'
    functor = 0
  []
  [no_slip_y]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = vel_y
    boundary = 'left right top bottom'
    functor = 0
  []
  [pressure-extrapolation]
    type = LinearFVExtrapolatedPressureBC
    boundary = 'left right top bottom'
    variable = pressure
    use_two_term_expansion = true
  []
  [walls_mu_t]
    type = LinearFVTurbulentViscosityWallFunctionBC
    boundary = 'bottom top'
    variable = 'mu_t'
    u = vel_x
    v = vel_y
    rho = ${rho}
    mu = ${mu}
    tke = TKE
    wall_treatment = ${wall_treatment}
  []
[]

[AuxVariables]
  [mu_t]
    type = MooseLinearVariableFVReal
    initial_condition = '${fparse rho * C_mu * ${k_init}^2 / eps_init}'
  []
  [yplus]
    type = MooseLinearVariableFVReal
  []
  [mu_eff]
    type = MooseLinearVariableFVReal
    initial_condition = '${fparse rho * C_mu * ${k_init}^2 / eps_init}'
  []
[]

[AuxKernels]
  [compute_mu_t]
    type = kEpsilonViscosityAux
    variable = mu_t
    C_mu = ${C_mu}
    tke = TKE
    epsilon = TKED
    mu = ${mu}
    rho = ${rho}
    u = vel_x
    v = vel_y
    bulk_wall_treatment = ${bulk_wall_treatment}
    walls = ${walls}
    wall_treatment = ${wall_treatment}
    execute_on = 'NONLINEAR'
    mu_t_ratio_max = 1e20
  []
  [compute_y_plus]
    type = RANSYPlusAux
    variable = yplus
    tke = TKE
    mu = ${mu}
    rho = ${rho}
    u = vel_x
    v = vel_y
    walls = ${walls}
    wall_treatment = ${wall_treatment}
    execute_on = 'NONLINEAR'
  []
  [compute_mu_eff]
    type = ParsedAux
    variable = 'mu_eff'
    coupled_variables = 'mu_t'
    expression = 'mu_t + ${mu}'
    execute_on = 'NONLINEAR'
  []
[]

[Executioner]
  type = SIMPLE

  rhie_chow_user_object = 'rc'
  momentum_systems = 'u_system v_system'
  pressure_system = 'pressure_system'
  turbulence_systems = 'TKE_system TKED_system'

  momentum_l_abs_tol = 1e-14
  pressure_l_abs_tol = 1e-14
  turbulence_l_abs_tol = 1e-14
  momentum_l_tol = 1e-14
  pressure_l_tol = 1e-14
  turbulence_l_tol = 1e-14

  momentum_equation_relaxation = 0.7
  pressure_variable_relaxation = 0.3
  turbulence_equation_relaxation = '0.5 0.5'
  num_iterations = 1000
  pressure_absolute_tolerance = 1e-12
  momentum_absolute_tolerance = 1e-12
  turbulence_absolute_tolerance = '1e-12 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'
  turbulence_petsc_options_iname = '-pc_type -pc_hypre_type'
  turbulence_petsc_options_value = 'hypre boomeramg'

  print_fields = false
  continue_on_max_its = true

  pin_pressure = true
  pressure_pin_value = 0.0
  pressure_pin_point = '0.01 0.099 0.0'
[]

[Outputs]
  csv = true
  perf_graph = false
  print_nonlinear_residuals = false
  print_linear_residuals = true
[]

[VectorPostprocessors]
  [side_bottom]
    type = SideValueSampler
    boundary = 'bottom'
    variable = 'vel_x vel_y pressure TKE TKED'
    sort_by = 'x'
    execute_on = 'timestep_end'
  []
  [side_top]
    type = SideValueSampler
    boundary = 'top'
    variable = 'vel_x vel_y pressure TKE TKED'
    sort_by = 'x'
    execute_on = 'timestep_end'
  []
  [side_left]
    type = SideValueSampler
    boundary = 'left'
    variable = 'vel_x vel_y pressure TKE TKED'
    sort_by = 'y'
    execute_on = 'timestep_end'
  []
  [side_right]
    type = SideValueSampler
    boundary = 'right'
    variable = 'vel_x vel_y pressure TKE TKED'
    sort_by = 'y'
    execute_on = 'timestep_end'
  []
  [horizontal_center]
    type = LineValueSampler
    start_point = '${fparse 0.01 * side_length} ${fparse 0.499 * side_length} 0'
    end_point = '${fparse 0.99 * side_length} ${fparse 0.499 * side_length} 0'
    num_points = ${Mesh/gen/nx}
    variable = 'vel_x vel_y pressure TKE TKED'
    sort_by = 'x'
    execute_on = 'timestep_end'
  []
  [vertical_center]
    type = LineValueSampler
    start_point = '${fparse 0.499 * side_length} ${fparse 0.01 * side_length} 0'
    end_point = '${fparse 0.499 * side_length} ${fparse 0.99 * side_length} 0'
    num_points =  ${Mesh/gen/ny}
    variable = 'vel_x vel_y pressure TKE TKED'
    sort_by = 'y'
    execute_on = 'timestep_end'
  []
[]

(test/tests/tag/tag-linearfv.i)

[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim = 1
    nx = 3
  []
[]

[Problem]
  linear_sys_names = 'u_sys v_sys'
  extra_tag_matrices = 'mat_tag_u; mat_tag_v'
  extra_tag_vectors = 'vec_tag_u; vec_tag_v'
[]

[Variables]
  [u]
    type = MooseLinearVariableFVReal
    initial_condition = 1.0
    solver_sys = u_sys
  []
  [v]
    type = MooseLinearVariableFVReal
    initial_condition = 0.5
    solver_sys = v_sys
  []
[]

[LinearFVKernels]
  [diffusion]
    type = LinearFVDiffusion
    variable = u
    diffusion_coeff = 2.0
  []
  [reaction]
    type = LinearFVReaction
    variable = u
    coeff = 3.0
    matrix_tags = 'system mat_tag_u'
  []
  [source]
    type = LinearFVSource
    variable = u
    source_density = 60.0
    vector_tags = 'rhs vec_tag_u'
  []
  [diffusion_v]
    type = LinearFVDiffusion
    variable = v
    diffusion_coeff = 1.0
  []
  [reaction_v]
    type = LinearFVReaction
    variable = v
    coeff = 1.5
    matrix_tags = 'system mat_tag_v'
  []
  [source_v]
    type = LinearFVSource
    variable = v
    source_density = 20.0
    vector_tags = 'rhs vec_tag_v'
  []
[]

[LinearFVBCs]
  [left_u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = u
    boundary = left
    functor = 1.0
  []
  [right_u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = u
    boundary = right
    functor = 3.0
  []
  [left_v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = v
    boundary = left
    functor = 1.0
  []
  [right_v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = v
    boundary = right
    functor = 3.0
  []
[]

[AuxVariables]
  [soln_u_dof]
    type = MooseLinearVariableFVReal
  []
  [soln_u]
    type = MooseLinearVariableFVReal
  []
  [rhs_u_dof]
    type = MooseLinearVariableFVReal
  []
  [rhs_u]
    type = MooseLinearVariableFVReal
  []
  [vector_tag_u]
    type = MooseLinearVariableFVReal
  []
  [matrix_u_diag]
    type = MooseLinearVariableFVReal
  []
  [soln_v_dof]
    type = MooseLinearVariableFVReal
  []
  [soln_v]
    type = MooseLinearVariableFVReal
  []
  [rhs_v_dof]
    type = MooseLinearVariableFVReal
  []
  [rhs_v]
    type = MooseLinearVariableFVReal
  []
  [vector_tag_v]
    type = MooseLinearVariableFVReal
  []
  [matrix_v_diag]
    type = MooseLinearVariableFVReal
  []
[]

[AuxKernels]
  [soln_u_dof]
    type = TagVectorDofValueAux
    variable = soln_u_dof
    v = u
    vector_tag = 'solution'
  []
  [soln_u]
    type = TagVectorAux
    variable = soln_u
    v = u
    vector_tag = 'solution'
  []
  [rhs_u_dof]
    type = TagVectorDofValueAux
    variable = rhs_u_dof
    v = u
    vector_tag = 'rhs'
  []
  [rhs_u]
    type = TagVectorAux
    variable = rhs_u
    v = u
    vector_tag = 'rhs'
  []
  [extra_vector_u_dof]
    type = TagVectorDofValueAux
    variable = vector_tag_u
    v = u
    vector_tag = 'vec_tag_u'
  []
  [extra_vector_u]
    type = TagVectorAux
    variable = vector_tag_u
    v = u
    vector_tag = 'vec_tag_u'
  []
  [extra_matrix_u]
    type = TagMatrixAux
    variable = matrix_u_diag
    v = u
    matrix_tag = 'mat_tag_u'
  []
  [soln_v_dof]
    type = TagVectorDofValueAux
    variable = soln_v_dof
    v = v
    vector_tag = 'solution'
  []
  [soln_v]
    type = TagVectorAux
    variable = soln_v
    v = v
    vector_tag = 'solution'
  []
  [rhs_v_dof]
    type = TagVectorDofValueAux
    variable = rhs_v_dof
    v = v
    vector_tag = 'rhs'
  []
  [rhs_v]
    type = TagVectorAux
    variable = rhs_v
    v = v
    vector_tag = 'rhs'
  []
  [extra_vector_v_dof]
    type = TagVectorDofValueAux
    variable = vector_tag_v
    v = v
    vector_tag = 'vec_tag_v'
  []
  [extra_vector_v]
    type = TagVectorAux
    variable = vector_tag_v
    v = v
    vector_tag = 'vec_tag_v'
  []
  [extra_matrix_v]
    type = TagMatrixAux
    variable = matrix_v_diag
    v = v
    matrix_tag = 'mat_tag_v'
  []
[]

[Executioner]
  type = Steady
  solve_type = 'LINEAR'
  system_names = 'u_sys v_sys'
  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre boomeramg'
[]

[Outputs]
  exodus = true
  execute_on = timestep_end
[]

(test/tests/time_integrators/implicit-euler/ie-linearfv.i)

###########################################################
# This is a simple test with a time-dependent problem
# demonstrating the use of the TimeIntegrator system.
#
# Testing a solution that is second order in space
# and first order in time
#
# @Requirement F1.30
###########################################################

[Mesh]
  type = GeneratedMesh
  dim = 2
  xmin = -1
  xmax = 1
  ymin = -1
  ymax = 1
  nx = 10
  ny = 10
[]

[Problem]
  linear_sys_names = 'u_sys'
[]

[Variables]
  [u]
    type = MooseLinearVariableFVReal
    solver_sys = 'u_sys'
    initial_condition = 0.0
  []
[]

[Functions]
  [forcing_fn]
    type = ParsedFunction
    expression = ((x*x)+(y*y))-(4*t)
  []
  [exact_fn]
    type = ParsedFunction
    expression = t*((x*x)+(y*y))
  []
[]

[LinearFVKernels]
  [ie]
    type = LinearFVTimeDerivative
    variable = u
  []
  [diff]
    type = LinearFVDiffusion
    variable = u
  []
  [source]
    type = LinearFVSource
    variable = u
    source_density = forcing_fn
  []
[]

[LinearFVBCs]
  [all]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = u
    boundary = '0 1 2 3'
    functor = exact_fn
  []
[]

[Postprocessors]
  [l2_err]
    type = ElementL2Error
    variable = u
    function = exact_fn
  []
[]

[Executioner]
  type = Transient

  system_names = u_sys
  l_tol = 1e-10

  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre boomeramg'

  # Test of the TimeIntegrator System
  scheme = 'implicit-euler'

  start_time = 0.0
  num_steps = 5
  dt = 0.25
[]

[Outputs]
  exodus = true
[]

(modules/navier_stokes/test/tests/finite_volume/ins/turbulence/bfs/linear-segregated-transient/BFS_ERCOFTAC.i)

Re = 5100

rho = 1.0
bulk_u = 1.0
H = 1.0
mu = '${fparse rho * bulk_u * H / Re}'

advected_interp_method = 'upwind'

### k-epsilon Closure Parameters ###
sigma_k = 1.0
sigma_eps = 1.3
C1_eps = 1.44
C2_eps = 1.92
C_mu = 0.09

### Initial and Boundary Conditions ###
intensity = 0.01
k_init = '${fparse 1.5*(intensity * bulk_u)^2}'
eps_init = '${fparse C_mu^0.75 * k_init^1.5 / H}'

### Modeling parameters ###
bulk_wall_treatment = false
walls = 'bottom wall-side top'
wall_treatment = 'neq' # Options: eq_newton, eq_incremental, eq_linearized, neq

[Mesh]
  [gen]
    type = CartesianMeshGenerator
    dim = 2
    dx = '${fparse 10.0*H} ${fparse 20.0*H}'
    dy = '${H} ${fparse 5*H}'
    ix = '8 16'
    iy = '2 8'
    subdomain_id = '
                    2 1
                    1 1
                  '
  []
  [corner_walls]
    type = SideSetsBetweenSubdomainsGenerator
    input = gen
    primary_block ='1'
    paired_block ='2'
    new_boundary = 'wall-side'
  []
  [delete_bottom]
    type = BlockDeletionGenerator
    input = corner_walls
    block ='2'
  []
  # Prevent test diffing on distributed parallel element numbering
  allow_renumbering = false
[]

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

[GlobalParams]
  rhie_chow_user_object = 'rc'
  advected_interp_method = ${advected_interp_method}
[]

[UserObjects]
  [rc]
    type = RhieChowMassFlux
    u = vel_x
    v = vel_y
    pressure = pressure
    rho = ${rho}
    p_diffusion_kernel = p_diffusion
  []
[]

[Variables]
  [vel_x]
    type = MooseLinearVariableFVReal
    initial_condition = ${bulk_u}
    solver_sys = u_system
  []
  [vel_y]
    type = MooseLinearVariableFVReal
    initial_condition = 0
    solver_sys = v_system
  []
  [pressure]
    type = MooseLinearVariableFVReal
    initial_condition = 1e-8
    solver_sys = pressure_system
  []
  [TKE]
    type = MooseLinearVariableFVReal
    solver_sys = TKE_system
    initial_condition = ${k_init}
  []
  [TKED]
    type = MooseLinearVariableFVReal
    solver_sys = TKED_system
    initial_condition = ${eps_init}
  []
[]

[LinearFVKernels]
  [u_time]
    type = LinearFVTimeDerivative
    variable = vel_x
    factor = ${rho}
  []
  [u_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_x
    advected_interp_method = ${advected_interp_method}
    mu = 'mu_t'
    u = vel_x
    v = vel_y
    momentum_component = 'x'
    rhie_chow_user_object = 'rc'
    use_nonorthogonal_correction = false
    use_deviatoric_terms = yes
  []
  [u_diffusion]
    type = LinearFVDiffusion
    variable = vel_x
    diffusion_coeff = ${mu}
    use_nonorthogonal_correction = false
  []
  [u_pressure]
    type = LinearFVMomentumPressure
    variable = vel_x
    pressure = pressure
    momentum_component = 'x'
  []

  [v_time]
    type = LinearFVTimeDerivative
    variable = vel_y
    factor = ${rho}
  []
  [v_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_y
    advected_interp_method = ${advected_interp_method}
    mu = 'mu_t'
    u = vel_x
    v = vel_y
    momentum_component = 'y'
    rhie_chow_user_object = 'rc'
    use_nonorthogonal_correction = false
    use_deviatoric_terms = yes
  []
  [v_diffusion]
    type = LinearFVDiffusion
    variable = vel_y
    diffusion_coeff = ${mu}
    use_nonorthogonal_correction = false
  []
  [v_pressure]
    type = LinearFVMomentumPressure
    variable = vel_y
    pressure = pressure
    momentum_component = 'y'
  []

  [p_diffusion]
    type = LinearFVAnisotropicDiffusion
    variable = pressure
    diffusion_tensor = Ainv
    use_nonorthogonal_correction = false
  []
  [HbyA_divergence]
    type = LinearFVDivergence
    variable = pressure
    face_flux = HbyA
    force_boundary_execution = true
  []

  [TKE_time]
    type = LinearFVTimeDerivative
    variable = TKE
    factor = ${rho}
  []
  [TKE_advection]
    type = LinearFVTurbulentAdvection
    variable = TKE
  []
  [TKE_diffusion]
    type = LinearFVTurbulentDiffusion
    variable = TKE
    diffusion_coeff = ${mu}
    use_nonorthogonal_correction = false
  []
  [TKE_turb_diffusion]
    type = LinearFVTurbulentDiffusion
    variable = TKE
    diffusion_coeff = 'mu_t'
    scaling_coeff = ${sigma_k}
    use_nonorthogonal_correction = false
  []
  [TKE_source_sink]
    type = LinearFVTKESourceSink
    variable = TKE
    u = vel_x
    v = vel_y
    epsilon = TKED
    rho = ${rho}
    mu = ${mu}
    mu_t = 'mu_t'
    walls = ${walls}
    wall_treatment = ${wall_treatment}
    C_pl = 1e10
  []

  [TKED_time]
    type = LinearFVTimeDerivative
    variable = TKED
    factor = ${rho}
  []
  [TKED_advection]
    type = LinearFVTurbulentAdvection
    variable = TKED
    walls = ${walls}
  []
  [TKED_diffusion]
    type = LinearFVTurbulentDiffusion
    variable = TKED
    diffusion_coeff = ${mu}
    use_nonorthogonal_correction = false
    walls = ${walls}
  []
  [TKED_turb_diffusion]
    type = LinearFVTurbulentDiffusion
    variable = TKED
    diffusion_coeff = 'mu_t'
    scaling_coeff = ${sigma_eps}
    use_nonorthogonal_correction = false
    walls = ${walls}
  []
  [TKED_source_sink]
    type = LinearFVTKEDSourceSink
    variable = TKED
    u = vel_x
    v = vel_y
    tke = TKE
    rho = ${rho}
    mu = ${mu}
    mu_t = 'mu_t'
    C1_eps = ${C1_eps}
    C2_eps = ${C2_eps}
    walls = ${walls}
    wall_treatment = ${wall_treatment}
    C_pl = 1e10
  []
[]

[LinearFVBCs]
  [inlet-u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left'
    variable = vel_x
    functor = '${bulk_u}'
  []
  [inlet-v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left'
    variable = vel_y
    functor = '0.0'
  []
  [inlet_TKE]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left'
    variable = TKE
    functor = '${k_init}'
  []
  [inlet_TKED]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left'
    variable = TKED
    functor = '${eps_init}'
  []

  [outlet_p]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'right'
    variable = pressure
    functor = 0.0
  []
  [outlet_u]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = 'right'
    variable = vel_x
    use_two_term_expansion = false
  []
  [outlet_v]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = 'right'
    variable = vel_y
    use_two_term_expansion = false
  []
  [outlet_TKE]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = 'right'
    variable = TKE
    use_two_term_expansion = false
  []
  [outlet_TKED]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = 'right'
    variable = TKED
    use_two_term_expansion = false
  []

  [walls-u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = ${walls}
    variable = vel_x
    functor = 0.0
  []
  [walls-v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = ${walls}
    variable = vel_y
    functor = 0.0
  []
  [walls_mu_t]
    type = LinearFVTurbulentViscosityWallFunctionBC
    boundary = ${walls}
    variable = 'mu_t'
    u = vel_x
    v = vel_y
    rho = ${rho}
    mu = ${mu}
    tke = TKE
    wall_treatment = ${wall_treatment}
  []
[]

[AuxVariables]
  [mu_t]
    type = MooseLinearVariableFVReal
    initial_condition = '${fparse rho * C_mu * ${k_init}^2 / eps_init}'
  []
  [yplus]
    type = MooseLinearVariableFVReal
  []
  [mu_eff]
    type = MooseLinearVariableFVReal
    initial_condition = '${fparse rho * C_mu * ${k_init}^2 / eps_init + mu}'
  []
[]

[AuxKernels]
  [compute_mu_t]
    type = kEpsilonViscosityAux
    variable = mu_t
    C_mu = ${C_mu}
    tke = TKE
    epsilon = TKED
    mu = ${mu}
    rho = ${rho}
    u = vel_x
    v = vel_y
    bulk_wall_treatment = ${bulk_wall_treatment}
    walls = ${walls}
    wall_treatment = ${wall_treatment}
    execute_on = 'NONLINEAR'
    mu_t_ratio_max = 1e20
  []
  [compute_y_plus]
    type = RANSYPlusAux
    variable = yplus
    tke = TKE
    mu = ${mu}
    rho = ${rho}
    u = vel_x
    v = vel_y
    walls = ${walls}
    wall_treatment = ${wall_treatment}
    execute_on = 'NONLINEAR'
  []
  [compute_mu_eff]
    type = ParsedAux
    variable = 'mu_eff'
    coupled_variables = 'mu_t'
    expression = 'mu_t + ${mu}'
    execute_on = 'NONLINEAR'
  []
[]

[Executioner]
  type = PIMPLE

  rhie_chow_user_object = 'rc'
  momentum_systems = 'u_system v_system'
  pressure_system = 'pressure_system'
  turbulence_systems = 'TKE_system TKED_system'

  momentum_l_abs_tol = 1e-14
  pressure_l_abs_tol = 1e-14
  turbulence_l_abs_tol = 1e-14
  momentum_l_tol = 1e-14
  pressure_l_tol = 1e-14
  turbulence_l_tol = 1e-14

  momentum_equation_relaxation = 0.7
  pressure_variable_relaxation = 0.3
  turbulence_equation_relaxation = '0.5 0.5'
  num_iterations = 20
  pressure_absolute_tolerance = 1e-6
  momentum_absolute_tolerance = 1e-6
  turbulence_absolute_tolerance = '1e-6 1e-6'
  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'
  turbulence_petsc_options_iname = '-pc_type -pc_hypre_type'
  turbulence_petsc_options_value = 'hypre boomeramg'

  print_fields = false
  continue_on_max_its = true
  dt = 10.0
  num_steps = 10
  num_piso_iterations = 2
[]

[Outputs]
  csv = true
[]

[VectorPostprocessors]
  [side_bottom]
    type = SideValueSampler
    boundary = 'bottom'
    variable = 'vel_x vel_y pressure TKE TKED'
    sort_by = 'x'
    execute_on = 'timestep_end'
  []
  [side_top]
    type = SideValueSampler
    boundary = 'top'
    variable = 'vel_x vel_y pressure TKE TKED'
    sort_by = 'x'
    execute_on = 'timestep_end'
  []
  [line_entry_channel_wall]
    type = LineValueSampler
    start_point = '${fparse 0.5 * H} ${fparse 1.00001 * H} 0'
    end_point = '${fparse 29.5 * H} ${fparse 1.00001 * H} 0'
    num_points = 24
    variable = 'vel_x vel_y pressure TKE TKED'
    sort_by = 'x'
    execute_on = 'timestep_end'
  []
  [line_quarter_entry_channel]
    type = LineValueSampler
    start_point = '${fparse 0.5 * H} ${fparse 2.25001 * H} 0'
    end_point = '${fparse 29.5 * H} ${fparse 2.25001 * H} 0'
    num_points = 24
    variable = 'vel_x vel_y pressure TKE TKED'
    sort_by = 'x'
    execute_on = 'timestep_end'
  []
[]

(modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d-heated/fluid.i)

mu = 2.6
rho = 1.0
advected_interp_method = 'upwind'
cp = 1000
k = 1

[Mesh]
  [mesh]
    type = CartesianMeshGenerator
    dim = 2
    dx = '0.25 0.25'
    dy = '0.2'
    ix = '5 5'
    iy = '5'
    subdomain_id = '0 1'
  []
[]

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

[UserObjects]
  [rc]
    type = RhieChowMassFlux
    u = vel_x
    v = vel_y
    pressure = pressure
    rho = ${rho}
    p_diffusion_kernel = p_diffusion
  []
[]

[Variables]
  [vel_x]
    type = MooseLinearVariableFVReal
    initial_condition = 0.5
    solver_sys = u_system
  []
  [vel_y]
    type = MooseLinearVariableFVReal
    solver_sys = v_system
    initial_condition = 0.0
  []
  [pressure]
    type = MooseLinearVariableFVReal
    solver_sys = pressure_system
    initial_condition = 0.2
  []
  [T_fluid]
    type = MooseLinearVariableFVReal
    solver_sys = energy_system
    initial_condition = 300
  []
[]

[LinearFVKernels]
  [u_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_x
    advected_interp_method = ${advected_interp_method}
    mu = ${mu}
    u = vel_x
    v = vel_y
    momentum_component = 'x'
    rhie_chow_user_object = 'rc'
    use_nonorthogonal_correction = false
  []
  [v_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_y
    advected_interp_method = ${advected_interp_method}
    mu = ${mu}
    u = vel_x
    v = vel_y
    momentum_component = 'y'
    rhie_chow_user_object = 'rc'
    use_nonorthogonal_correction = false
  []
  [u_pressure]
    type = LinearFVMomentumPressure
    variable = vel_x
    pressure = pressure
    momentum_component = 'x'
  []
  [v_pressure]
    type = LinearFVMomentumPressure
    variable = vel_y
    pressure = pressure
    momentum_component = 'y'
  []

  [p_diffusion]
    type = LinearFVAnisotropicDiffusion
    variable = pressure
    diffusion_tensor = Ainv
    use_nonorthogonal_correction = false
  []
  [HbyA_divergence]
    type = LinearFVDivergence
    variable = pressure
    face_flux = HbyA
    force_boundary_execution = true
  []

  [h_advection]
    type = LinearFVEnergyAdvection
    variable = T_fluid
    advected_quantity = temperature
    cp = ${cp}
    advected_interp_method = ${advected_interp_method}
    rhie_chow_user_object = 'rc'
  []
  [conduction]
    type = LinearFVDiffusion
    variable = T_fluid
    diffusion_coeff = ${k}
    use_nonorthogonal_correction = false
  []
  [heat_exchange]
    type = LinearFVVolumetricHeatTransfer
    variable = T_fluid
    h_solid_fluid = 100
    T_fluid = T_fluid
    T_solid = T_solid
    is_solid = false
    block = 1
  []
[]

[FunctorMaterials]
  [constant_functors]
    type = GenericFunctorMaterial
    prop_names = 'cp'
    prop_values = '${cp}'
  []
[]

[LinearFVBCs]
  [inlet-u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left'
    variable = vel_x
    functor = '1.1'
  []
  [inlet-v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left'
    variable = vel_y
    functor = '0.0'
  []
  [walls-u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'top bottom'
    variable = vel_x
    functor = 0.0
  []
  [walls-v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'top bottom'
    variable = vel_y
    functor = 0.0
  []
  [outlet_p]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'right'
    variable = pressure
    functor = 1.4
  []
  [outlet_u]
    type = LinearFVAdvectionDiffusionOutflowBC
    variable = vel_x
    use_two_term_expansion = false
    boundary = right
  []
  [outlet_v]
    type = LinearFVAdvectionDiffusionOutflowBC
    variable = vel_y
    use_two_term_expansion = false
    boundary = right
  []
  [inlet_wall_T]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = T_fluid
    functor = 300
    boundary = 'left top bottom'
  []
  [outlet_T]
    type = LinearFVAdvectionDiffusionOutflowBC
    variable = T_fluid
    use_two_term_expansion = false
    boundary = right
  []
[]

[AuxVariables]
  [T_solid]
    type = MooseLinearVariableFVReal
    initial_condition = 300
    block = 1
  []
[]

[MultiApps]
  inactive = 'solid'
  [solid]
    type = FullSolveMultiApp
    input_files = solid.i
    execute_on = timestep_begin
    no_restore = true
  []
[]

[Transfers]
  inactive = 'from_solid to_solid'
  [from_solid]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    from_multi_app = solid
    source_variable = 'T_solid'
    variable = 'T_solid'
    execute_on = timestep_begin
    from_blocks = 1
  []
  [to_solid]
    type = MultiAppGeneralFieldShapeEvaluationTransfer
    to_multi_app = solid
    source_variable = 'T_fluid'
    variable = 'T_fluid'
    execute_on = timestep_begin
    to_blocks = 1
  []
[]

[Executioner]
  type = SIMPLE
  momentum_l_abs_tol = 1e-13
  pressure_l_abs_tol = 1e-13
  energy_l_abs_tol = 1e-13
  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.8
  energy_equation_relaxation = 0.9
  pressure_variable_relaxation = 0.3
  num_iterations = 20
  pressure_absolute_tolerance = 1e-10
  momentum_absolute_tolerance = 1e-10
  energy_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'
  energy_petsc_options_iname = '-pc_type -pc_hypre_type'
  energy_petsc_options_value = 'hypre boomeramg'
  print_fields = false
  continue_on_max_its = true
[]

[Outputs]
  exodus = true
  execute_on = timestep_end
[]

(modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d/2d-boussinesq-transient.i)

mu = 2.6
rho = 1.0
advected_interp_method = 'average'
cp = 300
k = 10
alpha_b = 1e-4

[Mesh]
  [mesh]
    type = CartesianMeshGenerator
    dim = 2
    dx = '1.'
    dy = '0.2'
    ix = '10'
    iy = '5'
  []
[]

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

[UserObjects]
  [rc]
    type = RhieChowMassFlux
    u = vel_x
    v = vel_y
    pressure = pressure
    rho = ${rho}
    p_diffusion_kernel = p_diffusion
  []
[]

[Variables]
  [vel_x]
    type = MooseLinearVariableFVReal
    initial_condition = 0.5
    solver_sys = u_system
  []
  [vel_y]
    type = MooseLinearVariableFVReal
    solver_sys = v_system
    initial_condition = 0.0
  []
  [pressure]
    type = MooseLinearVariableFVReal
    solver_sys = pressure_system
    initial_condition = 0.2
  []
  [T]
    type = MooseLinearVariableFVReal
    solver_sys = energy_system
    initial_condition = 300
  []
[]

[LinearFVKernels]
  [u_time]
    type = LinearFVTimeDerivative
    variable = vel_x
    factor = ${rho}
  []
  [v_time]
    type = LinearFVTimeDerivative
    variable = vel_y
    factor = ${rho}
  []
  [u_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_x
    advected_interp_method = ${advected_interp_method}
    mu = ${mu}
    u = vel_x
    v = vel_y
    momentum_component = 'x'
    rhie_chow_user_object = 'rc'
    use_nonorthogonal_correction = false
  []
  [v_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_y
    advected_interp_method = ${advected_interp_method}
    mu = ${mu}
    u = vel_x
    v = vel_y
    momentum_component = 'y'
    rhie_chow_user_object = 'rc'
    use_nonorthogonal_correction = false
  []
  [u_pressure]
    type = LinearFVMomentumPressure
    variable = vel_x
    pressure = pressure
    momentum_component = 'x'
  []
  [v_pressure]
    type = LinearFVMomentumPressure
    variable = vel_y
    pressure = pressure
    momentum_component = 'y'
  []
  [u_boussinesq]
    type = LinearFVMomentumBoussinesq
    variable = vel_x
    rho = ${rho}
    gravity = '0 -9.81 0'
    alpha_name = ${alpha_b}
    ref_temperature = 300.0
    T_fluid = T
    momentum_component = 'x'
  []
  [v_boussinesq]
    type = LinearFVMomentumBoussinesq
    variable = vel_y
    rho = ${rho}
    gravity = '0 -9.81 0'
    alpha_name = ${alpha_b}
    ref_temperature = 300.0
    T_fluid = T
    momentum_component = 'y'
  []

  [p_diffusion]
    type = LinearFVAnisotropicDiffusion
    variable = pressure
    diffusion_tensor = Ainv
    use_nonorthogonal_correction = false
  []
  [HbyA_divergence]
    type = LinearFVDivergence
    variable = pressure
    face_flux = HbyA
    force_boundary_execution = true
  []

  [h_time]
    type = LinearFVTimeDerivative
    variable = T
    factor = ${fparse rho*cp}
  []
  [h_advection]
    type = LinearFVEnergyAdvection
    variable = T
    advected_quantity = temperature
    cp = ${cp}
    advected_interp_method = ${advected_interp_method}
    rhie_chow_user_object = 'rc'
  []
  [conduction]
    type = LinearFVDiffusion
    variable = T
    diffusion_coeff = ${k}
    use_nonorthogonal_correction = false
  []
[]

[FunctorMaterials]
  [constant_functors]
    type = GenericFunctorMaterial
    prop_names = 'cp alpha_b'
    prop_values = '${cp} ${alpha_b}'
  []
[]

[LinearFVBCs]
  [inlet-u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left'
    variable = vel_x
    functor = '1.1'
  []
  [inlet-v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left'
    variable = vel_y
    functor = '0.0'
  []
  [walls-u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'top bottom'
    variable = vel_x
    functor = 0.0
  []
  [walls-v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'top bottom'
    variable = vel_y
    functor = 0.0
  []
  [outlet_p]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'right'
    variable = pressure
    functor = 1.4
  []
  [outlet_u]
    type = LinearFVAdvectionDiffusionOutflowBC
    variable = vel_x
    use_two_term_expansion = false
    boundary = right
  []
  [outlet_v]
    type = LinearFVAdvectionDiffusionOutflowBC
    variable = vel_y
    use_two_term_expansion = false
    boundary = right
  []
  [inlet_top_T]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = T
    functor = 300.0
    boundary = 'left top'
  []
  [bottom_T]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = T
    functor = wall-temperature
    boundary = bottom
  []
  [outlet_T]
    type = LinearFVAdvectionDiffusionOutflowBC
    variable = T
    use_two_term_expansion = false
    boundary = right
  []
[]

[Functions]
  [wall-temperature]
    type = ParsedFunction
    expression = '350 + 50 * sin(6.28*t)'
  []
[]

[Executioner]
  type = PIMPLE
  momentum_l_abs_tol = 1e-12
  pressure_l_abs_tol = 1e-12
  energy_l_abs_tol = 1e-12
  momentum_l_tol = 1e-12
  pressure_l_tol = 1e-12
  energy_l_tol = 1e-12
  rhie_chow_user_object = 'rc'
  momentum_systems = 'u_system v_system'
  pressure_system = 'pressure_system'
  energy_system = 'energy_system'
  momentum_equation_relaxation = 0.8
  pressure_variable_relaxation = 0.3
  energy_equation_relaxation = 0.9
  num_iterations = 100
  pressure_absolute_tolerance = 1e-11
  momentum_absolute_tolerance = 1e-11
  energy_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'
  energy_petsc_options_iname = '-pc_type -pc_hypre_type'
  energy_petsc_options_value = 'hypre boomeramg'
  print_fields = false
  continue_on_max_its = true
  dt = 0.01
  num_steps = 6
  num_piso_iterations = 0
[]


[Outputs]
  exodus = true
[]

(test/tests/linearfvkernels/diffusion-reaction-advection/advection-diffusion-reaction-1d.i)

[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim = 1
    nx = 2
  []
[]

[Problem]
  linear_sys_names = 'u_sys'
[]

[Variables]
  [u]
    type = MooseLinearVariableFVReal
    solver_sys = 'u_sys'
    initial_condition = 1.0
  []
[]

[LinearFVKernels]
  [diffusion]
    type = LinearFVDiffusion
    variable = u
    diffusion_coeff = diff_coeff_func
    use_nonorthogonal_correction = false
  []
  [advection]
    type = LinearFVAdvection
    variable = u
    velocity = "0.5 0 0"
    advected_interp_method = average
  []
  [reaction]
    type = LinearFVReaction
    variable = u
    coeff = coeff_func
  []
  [source]
    type = LinearFVSource
    variable = u
    source_density = source_func
  []
[]

[LinearFVBCs]
  inactive = "outflow neumann"
  [dir]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = u
    boundary = "left right"
    functor = analytic_solution
  []
  [outflow]
    type = LinearFVAdvectionDiffusionOutflowBC
    variable = u
    boundary = "right"
    use_two_term_expansion = true
  []
  [neumann]
    type = LinearFVAdvectionDiffusionFunctorNeumannBC
    variable = u
    boundary = "left"
    functor = analytic_solution_neumann_left
    diffusion_coeff = diff_coeff_func
  []
[]

[Functions]
  [diff_coeff_func]
    type = ParsedFunction
    expression = '1+0.5*x'
  []
  [coeff_func]
    type = ParsedFunction
    expression = '1+1/(1+x)'
  []
  [source_func]
    type = ParsedFunction
    expression = '(1+1/(x+1))*(sin(pi/2*x)+1.5)+0.25*pi*pi*(0.5*x+1)*sin(pi/2*x)'
  []
  [analytic_solution]
    type = ParsedFunction
    expression = 'sin(pi/2*x)+1.5'
  []
  [analytic_solution_neumann_left]
    type = ParsedFunction
    expression = '-(1+0.5*x)*cos(pi/2*x)*pi/2'
  []
[]

[Postprocessors]
  [h]
    type = AverageElementSize
    execute_on = FINAL
  []
  [error]
    type = ElementL2FunctorError
    approximate = u
    exact = analytic_solution
    execute_on = FINAL
  []
[]

[Convergence]
  [linear]
    type = IterationCountConvergence
    max_iterations = 1
    converge_at_max_iterations = true
  []
[]

[Executioner]
  type = Steady
  system_names = u_sys
  multi_system_fixed_point=true
  multi_system_fixed_point_convergence=linear
  petsc_options_iname = '-pc_type -pc_hypre_type -ksp_rtol'
  petsc_options_value = 'hypre boomeramg 1e-10'
[]

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

(test/tests/linearfvkernels/diffusion/diffusion-1d_neumann.i)

[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim = 1
    nx = 2
  []
[]

[Problem]
  linear_sys_names = 'u_sys'
[]

[Variables]
  [u]
    type = MooseLinearVariableFVReal
    solver_sys = 'u_sys'
    initial_condition = 1.0
  []
[]

[LinearFVKernels]
  [diffusion]
    type = LinearFVDiffusion
    variable = u
    diffusion_coeff = coeff_func
  []
  [source]
    type = LinearFVSource
    variable = u
    source_density = source_func
  []
[]

[LinearFVBCs]
  [dir]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = u
    boundary = "right"
    functor = analytic_solution
  []
  [neu]
    type = LinearFVAdvectionDiffusionFunctorNeumannBC
    variable = u
    boundary = "left"
    functor = analytic_solution_neumann
  []
[]

[Functions]
  [coeff_func]
    type = ParsedFunction
    expression = '0.5*x'
  []
  [source_func]
    type = ParsedFunction
    expression = '2*x'
  []
  [analytic_solution]
    type = ParsedFunction
    expression = '1-x*x'
  []
  [analytic_solution_neumann]
    type = ParsedFunction
    expression = '-(0.5*x)*(-2*x)'
  []
[]

[Postprocessors]
  [h]
    type = AverageElementSize
    execute_on = FINAL
  []
  [error]
    type = ElementL2FunctorError
    approximate = u
    exact = analytic_solution
    execute_on = FINAL
  []
[]

[Executioner]
  type = Steady
  system_names = u_sys
  l_tol = 1e-10
  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre boomeramg'
[]

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

(modules/navier_stokes/test/tests/finite_volume/ins/channel-flow/linear-segregated/2d/2d-boussinesq.i)

mu = 2.6
rho = 1.0
advected_interp_method = 'average'
cp = 300
k = 20
alpha_b = 1e-4

[Mesh]
  [mesh]
    type = CartesianMeshGenerator
    dim = 2
    dx = '1.'
    dy = '0.2'
    ix = '10'
    iy = '5'
  []
[]

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

[UserObjects]
  [rc]
    type = RhieChowMassFlux
    u = vel_x
    v = vel_y
    pressure = pressure
    rho = ${rho}
    p_diffusion_kernel = p_diffusion
  []
[]

[Variables]
  [vel_x]
    type = MooseLinearVariableFVReal
    initial_condition = 0.5
    solver_sys = u_system
  []
  [vel_y]
    type = MooseLinearVariableFVReal
    solver_sys = v_system
    initial_condition = 0.0
  []
  [pressure]
    type = MooseLinearVariableFVReal
    solver_sys = pressure_system
    initial_condition = 0.2
  []
  [T]
    type = MooseLinearVariableFVReal
    solver_sys = energy_system
    initial_condition = 300
  []
[]

[LinearFVKernels]
  [u_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_x
    advected_interp_method = ${advected_interp_method}
    mu = ${mu}
    u = vel_x
    v = vel_y
    momentum_component = 'x'
    rhie_chow_user_object = 'rc'
    use_nonorthogonal_correction = false
  []
  [v_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_y
    advected_interp_method = ${advected_interp_method}
    mu = ${mu}
    u = vel_x
    v = vel_y
    momentum_component = 'y'
    rhie_chow_user_object = 'rc'
    use_nonorthogonal_correction = false
  []
  [u_pressure]
    type = LinearFVMomentumPressure
    variable = vel_x
    pressure = pressure
    momentum_component = 'x'
  []
  [v_pressure]
    type = LinearFVMomentumPressure
    variable = vel_y
    pressure = pressure
    momentum_component = 'y'
  []
  [u_boussinesq]
    type = LinearFVMomentumBoussinesq
    variable = vel_x
    rho = ${rho}
    gravity = '0 -9.81 0'
    alpha_name = ${alpha_b}
    ref_temperature = 300.0
    T_fluid = T
    momentum_component = 'x'
  []
  [v_boussinesq]
    type = LinearFVMomentumBoussinesq
    variable = vel_y
    rho = ${rho}
    gravity = '0 -9.81 0'
    alpha_name = ${alpha_b}
    ref_temperature = 300.0
    T_fluid = T
    momentum_component = 'y'
  []

  [p_diffusion]
    type = LinearFVAnisotropicDiffusion
    variable = pressure
    diffusion_tensor = Ainv
    use_nonorthogonal_correction = false
  []
  [HbyA_divergence]
    type = LinearFVDivergence
    variable = pressure
    face_flux = HbyA
    force_boundary_execution = true
  []

  [h_advection]
    type = LinearFVEnergyAdvection
    variable = T
    advected_quantity = temperature
    cp = ${cp}
    advected_interp_method = ${advected_interp_method}
    rhie_chow_user_object = 'rc'
  []
  [conduction]
    type = LinearFVDiffusion
    variable = T
    diffusion_coeff = ${k}
    use_nonorthogonal_correction = false
  []
[]

[FunctorMaterials]
  [constant_functors]
    type = GenericFunctorMaterial
    prop_names = 'cp alpha_b'
    prop_values = '${cp} ${alpha_b}'
  []
[]

[LinearFVBCs]
  [inlet-u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left'
    variable = vel_x
    functor = '1.1'
  []
  [inlet-v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'left'
    variable = vel_y
    functor = '0.0'
  []
  [walls-u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'top bottom'
    variable = vel_x
    functor = 0.0
  []
  [walls-v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'top bottom'
    variable = vel_y
    functor = 0.0
  []
  [outlet_p]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'right'
    variable = pressure
    functor = 1.4
  []
  [outlet_u]
    type = LinearFVAdvectionDiffusionOutflowBC
    variable = vel_x
    use_two_term_expansion = false
    boundary = right
  []
  [outlet_v]
    type = LinearFVAdvectionDiffusionOutflowBC
    variable = vel_y
    use_two_term_expansion = false
    boundary = right
  []
  [inlet_top_T]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = T
    functor = ${fparse 300.0}
    boundary = 'left top'
  []
  [bottom_T]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = T
    functor = ${fparse 350.0}
    boundary = bottom
  []
  [outlet_T]
    type = LinearFVAdvectionDiffusionOutflowBC
    variable = T
    use_two_term_expansion = false
    boundary = right
  []
[]

[Executioner]
  type = SIMPLE
  momentum_l_abs_tol = 1e-10
  pressure_l_abs_tol = 1e-10
  energy_l_abs_tol = 1e-10
  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.8
  energy_equation_relaxation = 0.9
  pressure_variable_relaxation = 0.3
  num_iterations = 200
  pressure_absolute_tolerance = 1e-10
  momentum_absolute_tolerance = 1e-10
  energy_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'
  energy_petsc_options_iname = '-pc_type -pc_hypre_type'
  energy_petsc_options_value = 'hypre boomeramg'
  print_fields = false
[]

[Outputs]
  exodus = true
[]

(test/tests/linearfvkernels/diffusion-reaction-advection/advection-diffusion-reaction-2d.i)

[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 2
    ny = 1
    ymax = 0.5
  []
[]

[Problem]
  linear_sys_names = 'u_sys'
[]

[Variables]
  [u]
    type = MooseLinearVariableFVReal
    solver_sys = 'u_sys'
    initial_condition = 1.0
  []
[]

[LinearFVKernels]
  [diffusion]
    type = LinearFVDiffusion
    variable = u
    diffusion_coeff = diff_coeff_func
    use_nonorthogonal_correction = false
  []
  [advection]
    type = LinearFVAdvection
    variable = u
    velocity = "0.5 0 0"
    advected_interp_method = average
  []
  [reaction]
    type = LinearFVReaction
    variable = u
    coeff = coeff_func
  []
  [source]
    type = LinearFVSource
    variable = u
    source_density = source_func
  []
[]

[LinearFVBCs]
  inactive = "outflow neumann"
  [dir]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = u
    boundary = "left right top bottom"
    functor = analytic_solution
  []
  [outflow]
    type = LinearFVAdvectionDiffusionOutflowBC
    variable = u
    boundary = "right"
    use_two_term_expansion = true
  []
  [neumann]
    type = LinearFVAdvectionDiffusionFunctorNeumannBC
    variable = u
    boundary = "top"
    functor = analytic_solution_neumann_top
    diffusion_coeff = diff_coeff_func
  []
[]

[Functions]
  [diff_coeff_func]
    type = ParsedFunction
    expression = '1.0+0.5*x*y'
  []
  [coeff_func]
    type = ParsedFunction
    expression = '1.0+1.0/(1+x*y)'
  []
  [source_func]
    type = ParsedFunction
    expression = '-1.0*x*pi*sin((1/2)*x*pi)*cos(2*y*pi) - 0.25*y*pi*sin(2*y*pi)*cos((1/2)*x*pi) + (1.0 + 1.0/(x*y + 1))*(sin((1/2)*x*pi)*sin(2*y*pi) + 1.5) + (17/4)*pi^2*(0.5*x*y + 1.0)*sin((1/2)*x*pi)*sin(2*y*pi) + 0.25*pi*sin(2*y*pi)*cos((1/2)*x*pi)'
  []
  [analytic_solution]
    type = ParsedFunction
    expression = 'sin((1/2)*x*pi)*sin(2*y*pi) + 1.5'
  []
  [analytic_solution_neumann_top]
    type = ParsedFunction
    expression = '(1.0+0.5*x*y)*sin((1/2)*x*pi)*cos(2*y*pi)*2*pi'
  []
[]

[Postprocessors]
  [h]
    type = AverageElementSize
    execute_on = FINAL
  []
  [error]
    type = ElementL2FunctorError
    approximate = u
    exact = analytic_solution
    execute_on = FINAL
  []
[]

[Convergence]
  [linear]
    type = IterationCountConvergence
    max_iterations = 1
    converge_at_max_iterations = true
  []
[]

[Executioner]
  type = Steady
  system_names = u_sys
  l_tol = 1e-10
  multi_system_fixed_point=true
  multi_system_fixed_point_convergence=linear
  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre boomeramg'
[]

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

(test/tests/linearfvkernels/diffusion/diffusion-1d.i)

[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim = 1
    nx = 2
  []
[]

[Problem]
  linear_sys_names = 'u_sys'
[]

[Variables]
  [u]
    type = MooseLinearVariableFVReal
    solver_sys = 'u_sys'
    initial_condition = 1.0
  []
[]

[LinearFVKernels]
  [diffusion]
    type = LinearFVDiffusion
    variable = u
    diffusion_coeff = coeff_func
  []
  [source]
    type = LinearFVSource
    variable = u
    source_density = source_func
  []
[]

[LinearFVBCs]
  [dir]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = u
    boundary = "left right"
    functor = analytic_solution
  []
[]

[Functions]
  [coeff_func]
    type = ParsedFunction
    expression = '0.5*x'
  []
  [source_func]
    type = ParsedFunction
    expression = '2*x'
  []
  [analytic_solution]
    type = ParsedFunction
    expression = '1-x*x'
  []
[]

[Postprocessors]
  [h]
    type = AverageElementSize
    execute_on = FINAL
  []
  [error]
    type = ElementL2FunctorError
    approximate = u
    exact = analytic_solution
    execute_on = FINAL
  []
[]

[Executioner]
  type = Steady
  system_names = u_sys
  l_tol = 1e-10
  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre boomeramg'
[]

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

(test/tests/multisystem/picard/linearfv_nonlinearfv/same_input.i)

[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim = 1
    nx = 6
  []
[]

[Problem]
  nl_sys_names = 'v_sys'
  linear_sys_names = 'u_sys'
[]

[Variables]
  [v]
    type = MooseVariableFVReal
    initial_condition = 2.0
    solver_sys = v_sys
  []
  [u]
    type = MooseLinearVariableFVReal
    solver_sys = 'u_sys'
    initial_condition = 1.0
  []
[]

[FVKernels]
  [diffusion]
    type = FVDiffusion
    variable = v
    coeff = u
  []
  [source]
    type = FVBodyForce
    variable = v
    function = 3
  []
[]

[LinearFVKernels]
  [diffusion]
    type = LinearFVDiffusion
    variable = u
    diffusion_coeff = v
  []
  [source]
    type = LinearFVSource
    variable = u
    source_density = 1
  []
[]

[FVBCs]
  [dir]
    type = FVFunctorDirichletBC
    variable = v
    boundary = "left right"
    functor = 2
  []
[]

[LinearFVBCs]
  [dir]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = u
    boundary = "left right"
    functor = 1
  []
[]

[Convergence]
  [linear]
    type = IterationCountConvergence
    max_iterations = 6
    converge_at_max_iterations = true
  []
[]

[Executioner]
  type = Steady
  system_names = 'v_sys u_sys'
  l_abs_tol = 1e-12
  l_tol = 1e-10
  nl_abs_tol = 1e-10
  multi_system_fixed_point=true
  multi_system_fixed_point_convergence=linear
  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre boomeramg'
[]

[Outputs]
  exodus = true
  execute_on = timestep_end
[]

(test/tests/linearfvbcs/robin/diffusion-2d-robin.i)

#################################################################
k     = 7.0 # diffusion coeff.
amp   = 3.6 # sinusoid amplitude, for u_exact

x1   = ${fparse 0.1*pi}
x2   = ${fparse 1.0*pi}
y1   = ${fparse 0.0*pi}
y2   = ${fparse 1.0*pi}

alpha = 2.000 # robin BC coeff for gradient term
beta  = 5.000 # robin BC coeff for variable term

nx = 2
ny = 2
##################################################################

[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim  = 2
    nx   = ${nx}
    ny   = ${ny}
    xmin = ${x1}
    xmax = ${x2}
    ymin = ${y1}
    ymax = ${y2}
  []
[]

[Problem]
  linear_sys_names = 'u_sys'
[]

[Variables]
  [u]
    type = MooseLinearVariableFVReal
    solver_sys = 'u_sys'
    initial_condition = 0.0
  []
[]

[Functions]
  [u_exact]
    type = ParsedFunction
    expression = '(${amp}*sin(x)*sin(y))'
  []
  [source_fn]
    type = ParsedFunction
    expression = '${fparse k*amp}*2.0*sin(x)*sin(y)'
  []
  [gamma_fn]
    type = ParsedFunction
    expression = '${fparse -amp*alpha}*cos(x)*sin(y) + ${beta} * u_e'
    symbol_names = 'u_e'
    symbol_values = 'u_exact'
  []
[]

[LinearFVKernels]
  [diffusion]
    type = LinearFVDiffusion
    variable = u
    diffusion_coeff = ${k}
    use_nonorthogonal_correction = true
  []
  [source]
    type = LinearFVSource
    variable = u
    source_density = source_fn
  []
[]

[LinearFVBCs]
  [right]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = u
    boundary = "right"
    functor = u_exact
  []
  [top]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = u
    boundary = "top"
    functor = u_exact
  []
  [bottom]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = u
    boundary = "bottom"
    functor = u_exact
  []
  [robin_left]
    type = LinearFVAdvectionDiffusionFunctorRobinBC
    variable = u
    boundary = "left"
    alpha = ${alpha}
    beta  = ${beta}
    gamma = gamma_fn
  []
[]

[Postprocessors]
  [h]
    type = AverageElementSize
    execute_on = FINAL
  []
  [error]
    type = ElementL2FunctorError
    approximate = u
    exact = u_exact
    execute_on = FINAL
  []
[]

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

[Convergence]
  [linear]
    type = IterationCountConvergence
    max_iterations = 4
    converge_at_max_iterations = true
  []
[]

[Executioner]
  type = Steady
  system_names = u_sys
  l_tol = 1e-7
  multi_system_fixed_point = true
  multi_system_fixed_point_convergence = linear
  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre boomeramg'
  linear_convergence = linear
[]

(test/tests/multisystem/restore_multiapp/parent.i)

[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim = 1
    nx = 6
  []
[]

[Problem]
  nl_sys_names = 'v_sys'
  linear_sys_names = 'u_sys'
[]

[Variables]
  [v]
    type = MooseVariableFVReal
    initial_condition = 2.0
    solver_sys = v_sys
  []
  [u]
    type = MooseLinearVariableFVReal
    solver_sys = 'u_sys'
    initial_condition = 1.0
  []
[]

[FVKernels]
  [diffusion]
    type = FVDiffusion
    variable = v
    coeff = u
  []
  [source]
    type = FVBodyForce
    variable = v
    function = 3
  []
[]

[LinearFVKernels]
  [diffusion]
    type = LinearFVDiffusion
    variable = u
    diffusion_coeff = v
  []
  [source]
    type = LinearFVSource
    variable = u
    source_density = 1
  []
[]

[FVBCs]
  [dir]
    type = FVFunctorDirichletBC
    variable = v
    boundary = "left right"
    functor = 2
  []
[]

[LinearFVBCs]
  [dir]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = u
    boundary = "left right"
    functor = 1
  []
[]

[Convergence]
  [linear]
    type = IterationCountConvergence
    max_iterations = 6
    converge_at_max_iterations = true
  []
[]

[MultiApps]
  [sub]
    type = FullSolveMultiApp
    input_files = sub.i
    execute_on = TIMESTEP_END
    keep_solution_during_restore = true
  []
[]

[Executioner]
  type = Steady
  system_names = 'v_sys u_sys'
  l_abs_tol = 1e-12
  l_tol = 1e-10
  nl_abs_tol = 1e-10
  multi_system_fixed_point=true
  multi_system_fixed_point_convergence=linear
  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre boomeramg'
  fixed_point_max_its = 2
  fixed_point_min_its = 2
  #disable_fixed_point_residual_norm_check = true
  accept_on_max_fixed_point_iteration = true
[]

(test/tests/outputs/debug/show_execution_linear_fv_flux.i)

[Mesh]
  [gen_mesh]
    type = GeneratedMeshGenerator
    dim = 1
    xmin = 0
    xmax = 10
    nx = 50
  []
  [left]
    type = ParsedSubdomainMeshGenerator
    input = 'gen_mesh'
    combinatorial_geometry = 'x < 0.5'
    block_id = '1'
  []
[]

[Problem]
  linear_sys_names = 'u_sys'
[]

[Variables]
  [u]
    type = MooseLinearVariableFVReal
    solver_sys = 'u_sys'
  []
[]

[LinearFVKernels]
  [diffusion_1]
    type = LinearFVDiffusion
    variable = u
    diffusion_coeff = 1.5
    block = 0
  []
  [diffusion_2]
    type = LinearFVDiffusion
    variable = u
    diffusion_coeff = 2.5
    block = 1
  []
  [source_1]
    type = LinearFVSource
    variable = u
    source_density = 3.5
    block = 0
  []
  [source_2]
    type = LinearFVSource
    variable = u
    source_density = 4.5
    block = 1
  []
[]

[LinearFVBCs]
  [left_bc]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = u
    boundary = 'left right'
    functor = 0
  []
[]

[Executioner]
  type = Steady
  system_names = u_sys
  l_tol = 1e-10
[]

[Debug]
  show_execution_order = 'NONLINEAR'
[]

(test/tests/vectorpostprocessors/side_value_sampler/side_value_sampler_fv.i)

[Mesh]
  [mesh]
    type = CartesianMeshGenerator
    dim = 2
    dx = '0.5 0.5'
    dy = '1'
    ix = '5 5'
    iy = '10'
    subdomain_id = '1 1'
  []
  # Limited to 1 side to avoid inconsistencies in parallel
  [internal_sideset]
    type = ParsedGenerateSideset
    combinatorial_geometry = 'y<0.51 & y>0.49 & x<0.11'
    new_sideset_name = 'center'
    input = 'mesh'
  []
  # this keeps numbering continuous so tests dont fail on different ids in CSV
  allow_renumbering = false
[]

[Problem]
  nl_sys_names = 'u_sys'
  linear_sys_names = 'v_sys'
[]

[Variables]
  [u]
    type = MooseVariableFVReal
    solver_sys = u_sys
  []
  [v]
    type = MooseLinearVariableFVReal
    solver_sys = v_sys
  []
[]

[FVKernels]
  [diff]
    type = FVDiffusion
    variable = u
    coeff = 1.0
  []
[]

[LinearFVKernels]
  [diff]
    type = LinearFVDiffusion
    variable = v
  []
[]

[FVBCs]
  [all]
    type = FVFunctorDirichletBC
    variable = u
    boundary = 'top bottom left right'
    functor = linear_x
  []
[]

[LinearFVBCs]
  [all]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = v
    boundary = 'top bottom left right'
    functor = linear_x
  []
[]

[Functions]
  [linear_x]
    type = ParsedFunction
    expression = '5*x'
  []
[]

[VectorPostprocessors]
  [side_sample]
    type = SideValueSampler
    variable = 'u v'
    boundary = top
    sort_by = x
  []
  [internal_sample]
    type = SideValueSampler
    variable = 'u v'
    boundary = center
    sort_by = 'id'
  []
[]

[Executioner]
  type = Steady
  solve_type = PJFNK
  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre boomeramg'
  system_names = 'u_sys v_sys'
  l_tol = 1e-8
[]

[Outputs]
  execute_on = 'timestep_end'
  [vpp_csv]
    type = CSV
  []
[]

(modules/navier_stokes/test/tests/finite_volume/ins/cht/conjugate_heat_transfer/cht_neu-dir.i)

### benchmark sources:
### https://doi.org/10.1016/j.compfluid.2018.06.016
### https://doi.org/10.1016/0017-9310(74)90087-8

b = 0.01 # plate thickness
l = 0.2 # plate length

nxi = 24 # nx in the inlet/entrance region
nyf = 18 # ny in fluid
nxf = 24 # nx in the main fluid region
nys = 8 # ny in the solid domain

fx1_bias = 1.00 # bdry layer bias - fluid
fx2_bias = '${fparse 1.0/1.00}' # bdry layer bias - solid
fy_bias = 1.20 # bdry layer bias - fluid
sy_bias = '${fparse 1.0/1.05}' # bdry layer bias - solid

k_s = 0.2876
rho = 0.3525
mu = 3.95e-5
k = 0.06808
cp = 1142.6

vin = 12.0
Tin = 1000.0
T_s_bottom = 600.0
P_out = 1.03e5

h_s = 0.0

advected_interp_method = 'upwind'

[Mesh]
  [fluid_channel]
    type = GeneratedMeshGenerator
    dim = 2
    nx = ${nxf}
    ny = ${nyf}
    xmin = 0
    xmax = ${l}
    ymin = 0
    ymax = '${fparse 10.0*b}'
    subdomain_ids = '1'
    subdomain_name = 'fluid'
    bias_x = '${fx1_bias}'
    bias_y = '${fparse fy_bias}'
    boundary_name_prefix = 'fluid'
  []
  [solid_base]
    type = GeneratedMeshGenerator
    dim = 2
    nx = ${nxf}
    ny = ${nys}
    xmin = 0
    xmax = ${l}
    ymin = '${fparse -b}'
    ymax = 0
    subdomain_ids = '2'
    subdomain_name = 'solid'
    bias_x = ${fx1_bias}
    bias_y = '${fparse sy_bias}'
    boundary_id_offset = 10
    boundary_name_prefix = 'solid'
  []
  [entrance]
    type = GeneratedMeshGenerator
    dim = 2
    nx = '${fparse 2.0*nxi}'
    ny = ${nyf}
    xmin = '${fparse -2.0*l}'
    xmax = 0
    ymin = 0
    ymax = '${fparse 10.0*b}'
    subdomain_ids = '0'
    subdomain_name = 'entrance'
    bias_x = ${fx2_bias}
    bias_y = '${fparse fy_bias}'
    boundary_id_offset = 20
    boundary_name_prefix = 'ent'
  []
  [smg]
    type = StitchedMeshGenerator
    inputs = 'entrance fluid_channel solid_base'
    stitch_boundaries_pairs = 'ent_right fluid_left;
                              fluid_bottom solid_top'
    prevent_boundary_ids_overlap = false
  []
  [interface]
    type = SideSetsBetweenSubdomainsGenerator
    input = 'smg'
    primary_block = 'fluid'
    paired_block = 'solid'
    new_boundary = interface
  []
  [symmetry_transform]
    type = SymmetryTransformGenerator
    input = interface
    mirror_point = '0 0 0'
    mirror_normal_vector = '0 1 0'
  []
  inactive = 'symmetry_transform'
[]

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

[UserObjects]
  [rc]
    type = RhieChowMassFlux
    u = vel_x
    v = vel_y
    pressure = pressure
    rho = ${rho}
    p_diffusion_kernel = p_diffusion
    block = '0 1'
  []
[]

[Variables]
  [vel_x]
    type = MooseLinearVariableFVReal
    initial_condition = ${vin}
    solver_sys = u_system
    block = '0 1'
  []
  [vel_y]
    type = MooseLinearVariableFVReal
    solver_sys = v_system
    initial_condition = 0.0
    block = '0 1'
  []
  [pressure]
    type = MooseLinearVariableFVReal
    solver_sys = pressure_system
    initial_condition = ${P_out}
    block = '0 1'
  []
  [T_fluid]
    type = MooseLinearVariableFVReal
    solver_sys = energy_system
    initial_condition = ${Tin}
    block = '0 1'
  []
  [T_solid]
    type = MooseLinearVariableFVReal
    solver_sys = solid_energy_system
    initial_condition = ${T_s_bottom}
    block = 2
  []
[]

[LinearFVKernels]
  [u_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_x
    advected_interp_method = ${advected_interp_method}
    mu = ${mu}
    u = vel_x
    v = vel_y
    momentum_component = 'x'
    rhie_chow_user_object = 'rc'
    use_nonorthogonal_correction = false
  []
  [v_advection_stress]
    type = LinearWCNSFVMomentumFlux
    variable = vel_y
    advected_interp_method = ${advected_interp_method}
    mu = ${mu}
    u = vel_x
    v = vel_y
    momentum_component = 'y'
    rhie_chow_user_object = 'rc'
    use_nonorthogonal_correction = false
  []
  [u_pressure]
    type = LinearFVMomentumPressure
    variable = vel_x
    pressure = pressure
    momentum_component = 'x'
  []
  [v_pressure]
    type = LinearFVMomentumPressure
    variable = vel_y
    pressure = pressure
    momentum_component = 'y'
  []

  [p_diffusion]
    type = LinearFVAnisotropicDiffusion
    variable = pressure
    diffusion_tensor = Ainv
    use_nonorthogonal_correction = false
  []
  [HbyA_divergence]
    type = LinearFVDivergence
    variable = pressure
    face_flux = HbyA
    force_boundary_execution = true
  []

  [h_advection]
    type = LinearFVEnergyAdvection
    variable = T_fluid
    advected_quantity = temperature
    cp = ${cp}
    advected_interp_method = ${advected_interp_method}
    rhie_chow_user_object = 'rc'
  []
  [conduction]
    type = LinearFVDiffusion
    variable = T_fluid
    diffusion_coeff = ${k}
    use_nonorthogonal_correction = false
  []
  [solid-conduction]
    type = LinearFVDiffusion
    variable = T_solid
    diffusion_coeff = ${k_s}
    use_nonorthogonal_correction = false
  []
[]

[LinearFVBCs]
  # velocity BCs
  [inlet-u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'ent_left'
    variable = vel_x
    functor = ${vin}
  []
  [inlet-v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'ent_left'
    variable = vel_y
    functor = '0.000'
  []
  [walls-u]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'ent_bottom interface'
    variable = vel_x
    functor = 0.0
  []
  [walls-v]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'ent_bottom interface'
    variable = vel_y
    functor = 0.0
  []
  [outlet_p]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    boundary = 'fluid_right'
    variable = pressure
    functor = ${P_out}
  []
  [outlet_u]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = 'fluid_right'
    variable = vel_x
    use_two_term_expansion = false
  []
  [outlet_v]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = 'fluid_right'
    variable = vel_y
    use_two_term_expansion = false
  []

  # freestream BCs for top of fluid domain
  [freestream_u]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = 'fluid_top ent_top'
    variable = vel_x
    use_two_term_expansion = false
  []
  [freestream_v]
    type = LinearFVAdvectionDiffusionOutflowBC
    boundary = 'fluid_top ent_top'
    variable = vel_y
    use_two_term_expansion = false
  []
  [freestream_p]
    type = LinearFVAdvectionDiffusionFunctorNeumannBC
    boundary = 'fluid_top ent_top'
    variable = pressure
    functor = 0
  []

  # temperature BCs
  [inlet_T]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = T_fluid
    functor = ${Tin}
    boundary = 'ent_left'
  []
  [heated_wall_solid]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = T_solid
    functor = ${T_s_bottom}
    boundary = 'solid_bottom'
  []
  [insulated_fluid]
    type = LinearFVAdvectionDiffusionFunctorNeumannBC
    variable = T_fluid
    functor = 0
    boundary = 'ent_bottom ent_top fluid_top'
  []
  [insulated_solid]
    type = LinearFVAdvectionDiffusionFunctorNeumannBC
    variable = T_solid
    functor = 0
    boundary = 'solid_left solid_right'
  []
  [outlet_T]
    type = LinearFVAdvectionDiffusionOutflowBC
    variable = T_fluid
    use_two_term_expansion = false
    boundary = 'fluid_right'
  []

  [fluid_solid]
    type = LinearFVDirichletCHTBC
    variable = T_fluid
    boundary = interface
    functor = interface_temperature_solid_interface
  []
  [solid_fluid]
    type = LinearFVRobinCHTBC
    variable = T_solid
    boundary = interface
    h = ${h_s}
    thermal_conductivity = ${k_s}
    incoming_flux = heat_flux_to_solid_interface
    surface_temperature = interface_temperature_fluid_interface
  []
[]

[FunctorMaterials]
  [rhocpT]
    property_name = 'rhocpT'
    type = ParsedFunctorMaterial
    functor_names = 'T_fluid'
    expression = '${rho}*${cp}*T_fluid'
  []
[]

[Postprocessors]
  [h_in]
    type = VolumetricFlowRate
    boundary = 'ent_left'
    vel_x = vel_x
    vel_y = vel_y
    rhie_chow_user_object = rc
    advected_quantity = 'rhocpT'
    subtract_mesh_velocity = false
  []
  [h_out]
    type = VolumetricFlowRate
    boundary = 'fluid_right fluid_top ent_top interface'
    vel_x = vel_x
    vel_y = vel_y
    rhie_chow_user_object = rc
    advected_quantity = 'rhocpT'
    advected_interp_method = upwind
    subtract_mesh_velocity = false
  []
[]

[VectorPostprocessors]
  [y_vs_ts]
    type = LineValueSampler
    variable = 'T_solid'
    start_point = '0.05 -1e-9 0' # making sure we are always in the domain
    end_point = '0.05 ${fparse -b+1e-9} 0'
    num_points = 8
    sort_by = id
    warn_discontinuous_face_values = false
  []
  [y_vs_tf]
    type = LineValueSampler
    variable = 'T_fluid'
    start_point = '0.05 1e-9 0' # making sure we are always in the domain
    end_point = '0.05 ${fparse b-1e-9} 0'
    num_points = 12
    sort_by = id
    warn_discontinuous_face_values = false
  []
[]

[Executioner]
  type = SIMPLE

  num_iterations = 1000

  momentum_systems = 'u_system v_system'
  pressure_system = 'pressure_system'
  rhie_chow_user_object = 'rc'
  momentum_l_abs_tol = 1e-10
  pressure_l_abs_tol = 1e-10
  momentum_l_tol = 0
  pressure_l_tol = 0
  momentum_equation_relaxation = 0.9
  pressure_variable_relaxation = 0.3
  momentum_absolute_tolerance = 1e-7
  pressure_absolute_tolerance = 1e-7
  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_system = 'energy_system'
  solid_energy_system = 'solid_energy_system'
  energy_l_abs_tol = 1e-10
  solid_energy_l_abs_tol = 1e-10
  energy_l_tol = 0
  solid_energy_l_tol = 0
  energy_equation_relaxation = 1.0
  energy_absolute_tolerance = 1e-7
  solid_energy_absolute_tolerance = 1e-7
  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'

  cht_interfaces = 'interface'
  cht_solid_flux_relaxation = 0.4
  cht_fluid_flux_relaxation = 0.4
  cht_solid_temperature_relaxation = 0.4
  cht_fluid_temperature_relaxation = 0.4
  max_cht_fpi = 3

  print_fields = false
[]

[Outputs]
  csv = true
  execute_on = timestep_end
[]

(test/tests/linearfvkernels/block-restriction/block-restricted-diffusion-react.i)

source=1
diff_coeff=2
reac_coeff=3

[Mesh]
  [cmg]
    type = CartesianMeshGenerator
    dim = 1
    dx = '0.5 0.5'
    ix = '20 20'
    subdomain_id = '1 2'
  []
[]

[Problem]
  linear_sys_names = 'u_sys'
[]

[Variables]
  [u]
    type = MooseLinearVariableFVReal
    solver_sys = 'u_sys'
    initial_condition = 1.0
  []
[]

[LinearFVKernels]
  [diffusion]
    type = LinearFVDiffusion
    variable = u
    diffusion_coeff = ${diff_coeff}
    use_nonorthogonal_correction = false
    block = 1
  []
  [reaction]
    type = LinearFVReaction
    variable = u
    coeff = ${reac_coeff}
    block = 2
  []
  [source]
    type = LinearFVSource
    variable = u
    source_density = ${source}
  []
[]

[LinearFVBCs]
  [dir]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = u
    boundary = "left"
    functor = 0
  []
[]

[Functions]
  [analytic_solution]
    type = ParsedFunction
    expression = 'if(x<0.5, -x*x*S/2/D+(S/C+0.5*0.5/2/D*S)/0.5*x, S/C)'
    symbol_names = 'S D C'
    symbol_values = '${source} ${diff_coeff} ${reac_coeff}'
  []
[]

[Postprocessors]
  [h]
    type = AverageElementSize
    execute_on = TIMESTEP_END
    block = 2
  []
  [error]
    type = ElementL2FunctorError
    approximate = u
    exact = analytic_solution
    execute_on = TIMESTEP_END
    block = 2
  []
[]

[Convergence]
  [linear]
    type = IterationCountConvergence
    max_iterations = 1
    converge_at_max_iterations = true
  []
[]

[Executioner]
  type = Steady
  system_names = u_sys
  multi_system_fixed_point=true
  multi_system_fixed_point_convergence=linear
  petsc_options_iname = '-pc_type -pc_hypre_type -ksp_rtol'
  petsc_options_value = 'hypre boomeramg 1e-10'
[]

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

(test/tests/transfers/multiapp_copy_transfer/linear_sys_to_aux/linear_sub.i)

[Mesh]
  [gmg]
    type = GeneratedMeshGenerator
    dim = 1
    nx = 10
  []
[]

[Problem]
  linear_sys_names = 'u_sys'
[]

[Variables]
  [u]
    type = MooseLinearVariableFVReal
    solver_sys = 'u_sys'
    initial_condition = 1.0
  []
[]

[LinearFVKernels]
  [diffusion]
    type = LinearFVDiffusion
    variable = u
    diffusion_coeff = coeff_func
  []
  [source]
    type = LinearFVSource
    variable = u
    source_density = source_func
  []
[]

[LinearFVBCs]
  [dir]
    type = LinearFVAdvectionDiffusionFunctorDirichletBC
    variable = u
    boundary = "left right"
    functor = analytic_solution
  []
[]

[Functions]
  [coeff_func]
    type = ParsedFunction
    expression = '0.5*x'
  []
  [source_func]
    type = ParsedFunction
    expression = '2*x'
  []
  [analytic_solution]
    type = ParsedFunction
    expression = '1-x*x'
  []
[]

[Executioner]
  type = Steady
  system_names = u_sys
  l_tol = 1e-10
  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre boomeramg'
[]

[Outputs]
  exodus = true
  execute_on = TIMESTEP_END
[]

(modules/navier_stokes/include/linearfvkernels/LinearFVTurbulentDiffusion.h)

// This file is part of the MOOSE framework
// https://www.mooseframework.org
//
// All rights reserved, see COPYRIGHT for full restrictions
// https://github.com/idaholab/moose/blob/master/COPYRIGHT
//
// Licensed under LGPL 2.1, please see LICENSE for details
// https://www.gnu.org/licenses/lgpl-2.1.html

#pragma once

#include "LinearFVDiffusion.h"

/**
 * Kernel that adds contributions from a diffusion term of the turbulent variables,
 * limited in the near-wall regions,
 * discretized using the finite volume method to a linear system.
 */
class LinearFVTurbulentDiffusion : public LinearFVDiffusion
{
public:
  static InputParameters validParams();

  /**
   * Class constructor.
   * @param params The InputParameters for the kernel.
   */
  LinearFVTurbulentDiffusion(const InputParameters & params);

  virtual void initialSetup() override;

  virtual Real computeElemMatrixContribution() override;

  virtual Real computeNeighborMatrixContribution() override;

  virtual Real computeElemRightHandSideContribution() override;

  virtual Real computeNeighborRightHandSideContribution() override;

  virtual Real computeBoundaryMatrixContribution(const LinearFVBoundaryCondition & bc) override;

  virtual Real computeBoundaryRHSContribution(const LinearFVBoundaryCondition & bc) override;

  virtual void addMatrixContribution() override;

  virtual void addRightHandSideContribution() override;

protected:
  /// The functor for the scaling coefficient for the diffusion term
  const Moose::Functor<Real> & _scaling_coeff;

  /// Wall boundaries
  const std::vector<BoundaryName> & _wall_boundary_names;

  /// List for wall bounded elements
  std::unordered_set<const Elem *> _wall_bounded;
};

Description
Example input syntax
Input Parameters
Input Files
Child Objects
References