MassFluxPenaltyIPHDG

Unlike the perturbation in MassFluxPenalty, this class does not introduce coupling between interior degress of freedom on adjacent elements. Instead it imposes on the momentum component equations the element face term

$(1)var element = document.getElementById("moose-equation-9f319bbd-c79f-4898-a113-ffaf52b74167");katex.render(" \\langle u - \\bar{u}, v - \\bar{v}\\rangle", 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-4f28087e-defe-4240-95c1-3d3391fe8416");katex.render("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}"}});$ is the primal finite element solution (defined on element interiors), $var element = document.getElementById("moose-equation-ccc074b1-4881-4765-84df-29b9cc135dcc");katex.render("\\bar{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}"}});$ is the facet solution, $var element = document.getElementById("moose-equation-25073471-b308-4997-b1ab-bef85b3bab35");katex.render("v", 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 the primal test functions, and $var element = document.getElementById("moose-equation-06b5b2e5-30e3-4b3c-9872-d7a217cd8864");katex.render("\\bar{v}", 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 the facet test functions. The augmented form of the IP-HDG discretized Navier-Stokes equations is algebraically equivalent to

$var element = document.getElementById("moose-equation-eaf23aeb-23d2-478c-8dc5-5a55831a8cb9");katex.render(" \\mathbb{A} = \\begin{bmatrix} A + \\gamma J & B^T \\\\ B & \\mathbf{0} \\end{bmatrix}.", 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-444c53a8-62a8-4180-bde3-e17199e29b3b");katex.render("A", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is a nonsingular $var element = document.getElementById("moose-equation-7988b5e8-a80e-4150-bb46-8d8e04e725d4");katex.render("n \\times 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}"}});$ matrix, $var element = document.getElementById("moose-equation-a490d5ec-8869-47f1-b684-49631a055143");katex.render("J \\ne 0", 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 a symmetric, singular $var element = document.getElementById("moose-equation-776c6125-d3a5-4ba0-a237-f248306ade21");katex.render("n \\times 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}"}});$ matrix, $var element = document.getElementById("moose-equation-5854962c-9437-436e-a060-0fab15626934");katex.render("B", 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 an $var element = document.getElementById("moose-equation-9569fbc5-fa07-4876-89f5-c5aef5e4bf28");katex.render("m \\times 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}"}});$ matrix ( $var element = document.getElementById("moose-equation-fa91d628-733c-4a92-95d6-ae09158fb1a8");katex.render("m < 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}"}});$ ), and $var element = document.getElementById("moose-equation-5b5e7b3a-adc1-4035-9de9-00a2b9c065c4");katex.render("\\gamma>0", 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 a scalar constant. Requirements for the singular perturbation are that $var element = document.getElementById("moose-equation-bd590220-c7e2-4a14-a46b-f580f7cec1ae");katex.render("\\mathcal{N}(J)\\subseteq \\mathcal{N}(B)", 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-598664bc-1579-454f-8b3b-219477655f60");katex.render("\\mathcal{R}(B^T) = \\mathcal{R}(J)", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ . Eq. (1), corresponding to the singular perturbation $var element = document.getElementById("moose-equation-5d3098d0-c5a7-444f-9750-66c2cfd32a3a");katex.render("J", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ , meets these requirements.

Input Parameters

componentThe velocity component this object is being applied to
C++ Type:unsigned short
Controllable:No
Description:The velocity component this object is being applied to
uThe x-velocity
C++ Type:NonlinearVariableName
Unit:(no unit assumed)
Controllable:No
Description:The x-velocity
u_faceThe x-velocity on the face
C++ Type:NonlinearVariableName
Unit:(no unit assumed)
Controllable:No
Description:The x-velocity on the face
vThe y-velocity
C++ Type:NonlinearVariableName
Unit:(no unit assumed)
Controllable:No
Description:The y-velocity
v_faceThe y-velocity on the face
C++ Type:NonlinearVariableName
Unit:(no unit assumed)
Controllable:No
Description:The y-velocity on the face
variableThe name of the variable that this residual object operates on
C++ Type:NonlinearVariableName
Unit:(no unit assumed)
Controllable:No
Description:The name of the variable that this residual object operates on

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
displacementsThe displacements
C++ Type:std::vector<VariableName>
Unit:(no unit assumed)
Controllable:No
Description:The displacements
gamma1The penalty to multiply the jump
Default:1
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:The penalty to multiply the jump
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)

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_tagsnontimeThe tag for the vectors this Kernel should fill
Default:nontime
C++ Type:MultiMooseEnum
Options:nontime, 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.
diag_save_inThe name of auxiliary variables to save this Kernel's diagonal Jacobian contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
C++ Type:std::vector<AuxVariableName>
Unit:(no unit assumed)
Controllable:No
Description:The name of auxiliary variables to save this Kernel's diagonal Jacobian contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
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
save_inThe name of auxiliary variables to save this Kernel's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
C++ Type:std::vector<AuxVariableName>
Unit:(no unit assumed)
Controllable:No
Description:The name of auxiliary variables to save this Kernel's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
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

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

Material Property Retrieval Parameters

Input Files

(modules/navier_stokes/test/tests/finite_element/ins/hdg/ip/lid-driven/lid-driven-scfsp.i)
(modules/navier_stokes/test/tests/finite_element/ins/hdg/ip/lid-driven/block-restricted-scfsp.i)

(modules/navier_stokes/test/tests/finite_element/ins/hdg/ip/lid-driven/lid-driven-scfsp.i)

final_re = 10000
starting_re = 10
rho = 1
l = 2
U = 1
n = 16
gamma = 1e4
degree = 2
alpha = '${fparse 10 * degree^2}'
num_steps = 10
step_length = '${fparse (log10(final_re) - log10(starting_re)) / (num_steps - 1)}'

[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0
    xmax = ${l}
    ymin = 0
    ymax = ${l}
    nx = ${n}
    ny = ${n}
    elem_type = TRI6
  []
[]

[Problem]
  type = NavierStokesProblem
  extra_tag_matrices = 'mass'
  mass_matrix = 'mass'
  use_pressure_mass_matrix = true
[]

[AuxVariables]
  [vel_mag]
    family = L2_HIERARCHIC
    order = SECOND
  []
[]

[AuxKernels]
  [vel_mag]
    type = VectorMagnitudeAux
    variable = vel_mag
    x = vel_x
    y = vel_y
  []
[]

[Variables]
  [vel_x]
    family = L2_HIERARCHIC
    order = SECOND
  []
  [vel_y]
    family = L2_HIERARCHIC
    order = SECOND
  []
  [pressure]
    family = L2_HIERARCHIC
    order = FIRST
  []
  [vel_bar_x]
    family = SIDE_HIERARCHIC
    order = SECOND
  []
  [vel_bar_y]
    family = SIDE_HIERARCHIC
    order = SECOND
  []
  [pressure_bar]
    family = SIDE_HIERARCHIC
    order = SECOND
  []
[]

[HDGKernels]
  [momentum_x_convection]
    type = AdvectionIPHDGKernel
    variable = vel_x
    face_variable = vel_bar_x
    velocity = 'velocity'
    coeff = ${rho}
  []
  [momentum_x_diffusion]
    type = NavierStokesStressIPHDGKernel
    variable = vel_x
    face_variable = vel_bar_x
    diffusivity = 'mu'
    alpha = ${alpha}
    pressure_variable = pressure
    pressure_face_variable = pressure_bar
    component = 0
  []
  [momentum_y_convection]
    type = AdvectionIPHDGKernel
    variable = vel_y
    face_variable = vel_bar_y
    velocity = 'velocity'
    coeff = ${rho}
  []
  [momentum_y_diffusion]
    type = NavierStokesStressIPHDGKernel
    variable = vel_y
    face_variable = vel_bar_y
    diffusivity = 'mu'
    alpha = ${alpha}
    pressure_variable = pressure
    pressure_face_variable = pressure_bar
    component = 1
  []
  [pressure_convection]
    type = AdvectionIPHDGKernel
    variable = pressure
    face_variable = pressure_bar
    velocity = 'velocity'
    coeff = '${fparse -rho}'
    self_advection = false
  []

  [u_jump]
    type = MassFluxPenaltyIPHDG
    variable = vel_x
    u = vel_x
    v = vel_y
    u_face = vel_bar_x
    v_face = vel_bar_y
    component = 0
    gamma = ${gamma}
  []
  [v_jump]
    type = MassFluxPenaltyIPHDG
    variable = vel_y
    u = vel_x
    v = vel_y
    u_face = vel_bar_x
    v_face = vel_bar_y
    component = 1
    gamma = ${gamma}
  []
  [pb_mass]
    type = MassMatrixHDG
    variable = pressure_bar
    matrix_tags = 'mass'
    density = '${fparse -1/gamma}'
  []
[]

[BCs]
  [momentum_x_diffusion_walls]
    type = NavierStokesStressIPHDGDirichletBC
    boundary = 'left bottom right'
    variable = vel_x
    face_variable = vel_bar_x
    pressure_variable = pressure
    pressure_face_variable = pressure_bar
    alpha = ${alpha}
    functor = '0'
    diffusivity = 'mu'
    component = 0
  []
  [momentum_x_diffusion_top]
    type = NavierStokesStressIPHDGDirichletBC
    boundary = 'top'
    variable = vel_x
    face_variable = vel_bar_x
    pressure_variable = pressure
    pressure_face_variable = pressure_bar
    alpha = ${alpha}
    functor = '${U}'
    diffusivity = 'mu'
    component = 0
  []
  [momentum_y_diffusion_all]
    type = NavierStokesStressIPHDGDirichletBC
    boundary = 'left bottom right top'
    variable = vel_y
    face_variable = vel_bar_y
    pressure_variable = pressure
    pressure_face_variable = pressure_bar
    alpha = ${alpha}
    functor = '0'
    diffusivity = 'mu'
    component = 1
  []

  [mass_convection]
    type = AdvectionIPHDGPrescribedFluxBC
    face_variable = pressure_bar
    variable = pressure
    velocity = 'velocity'
    coeff = '${fparse -rho}'
    self_advection = false
    boundary = 'left bottom top right'
    prescribed_normal_flux = 0
  []

  [pb_mass]
    type = MassMatrixIntegratedBC
    variable = pressure_bar
    matrix_tags = 'mass'
    boundary = 'left right bottom top'
    density = '${fparse -1/gamma}'
  []

  [u_jump_walls]
    type = MassFluxPenaltyBC
    variable = vel_x
    u = vel_x
    v = vel_y
    component = 0
    boundary = 'left right bottom'
    gamma = ${gamma}
    dirichlet_value = walls
  []
  [v_jump_walls]
    type = MassFluxPenaltyBC
    variable = vel_y
    u = vel_x
    v = vel_y
    component = 1
    boundary = 'left right bottom'
    gamma = ${gamma}
    dirichlet_value = walls
  []
  [u_jump_top]
    type = MassFluxPenaltyBC
    variable = vel_x
    u = vel_x
    v = vel_y
    component = 0
    boundary = 'top'
    gamma = ${gamma}
    dirichlet_value = top_vel
  []
  [v_jump_top]
    type = MassFluxPenaltyBC
    variable = vel_y
    u = vel_x
    v = vel_y
    component = 1
    boundary = 'top'
    gamma = ${gamma}
    dirichlet_value = top_vel
  []
[]

[Functions]
  [top_vel]
    type = ParsedVectorFunction
    expression_x = ${U}
  []
  [walls]
    type = ParsedVectorFunction
  []
  [reynolds]
    type = ParsedFunction
    expression = '10^(log10(${starting_re}) + (t - 1) * ${step_length})'
  []
[]

[Materials]
  [const]
    type = ADGenericConstantMaterial
    prop_names = 'rho'
    prop_values = '${rho}'
  []
  [vel]
    type = ADVectorFromComponentVariablesMaterial
    vector_prop_name = 'velocity'
    u = vel_x
    v = vel_y
  []
  [mu]
    type = ADParsedMaterial
    functor_names = 'reynolds'
    functor_symbols = 'reynolds'
    property_name = 'mu'
    expression = '${U} * ${l} / reynolds'
    output_properties = 'mu'
  []
[]

[Preconditioning]
  [FSP]
    type = SCFSP
    topsplit = 'up'
    [up]
      splitting = 'u p'
      splitting_type = schur
      petsc_options = '-ksp_monitor'
      petsc_options_iname = '-pc_fieldsplit_schur_fact_type  -pc_fieldsplit_schur_precondition -ksp_gmres_restart -ksp_type -ksp_pc_side -ksp_rtol -ksp_max_it -ksp_atol'
      petsc_options_value = 'full                            self                              300                fgmres    right        1e-4      30          1e-9'
    []
    [u]
      vars = 'vel_bar_x vel_bar_y'
      petsc_options = '-ksp_converged_reason'
      petsc_options_iname = '-pc_type -ksp_type -ksp_rtol -ksp_gmres_restart -ksp_pc_side -pc_factor_mat_solver_type -ksp_max_it'
      petsc_options_value = 'ilu      gmres     1e-2      300                right        strumpack                  30'
    []
    [p]
      vars = 'pressure_bar'
      petsc_options = '-ksp_converged_reason'
      petsc_options_iname = '-pc_type -ksp_type -ksp_rtol -ksp_gmres_restart -ksp_pc_side -pc_factor_mat_solver_type -ksp_max_it'
      petsc_options_value = 'ilu      gmres     1e-2      300                right        strumpack                  30'
    []
  []
[]

[Executioner]
  type = Transient
  num_steps = ${num_steps}
  petsc_options_iname = '-ksp_type'
  petsc_options_value = 'preonly'
[]

[Outputs]
  print_linear_residuals = 'false'
  [out]
    type = Exodus
    hide = 'pressure_average vel_bar_x vel_bar_y pressure_bar'
    output_material_properties = true
    execute_on = 'timestep_end'
  []
[]

[Postprocessors]
  [Re]
    type = FunctionValuePostprocessor
    function = 'reynolds'
  []
  [pressure_average]
    type = ElementAverageValue
    variable = pressure
  []
  [vel_average]
    type = ElementAverageValue
    variable = vel_mag
  []
[]

[Correctors]
  [set_pressure]
    type = NSPressurePin
    pin_type = 'average'
    variable = pressure
    pressure_average = 'pressure_average'
  []
[]

(modules/navier_stokes/test/tests/finite_element/ins/hdg/ip/lid-driven/block-restricted-scfsp.i)

final_re = 10000
starting_re = 10
rho = 1
l = 2
U = 1
n = 16
gamma = 1e4
degree = 2
alpha = '${fparse 10 * degree^2}'
num_steps = 10
step_length = '${fparse (log10(final_re) - log10(starting_re)) / (num_steps - 1)}'

[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    xmin = 0
    xmax = '${fparse 2 * l}'
    ymin = 0
    ymax = ${l}
    nx = '${fparse 2 * n}'
    ny = ${n}
    elem_type = TRI6
  []
  [remove]
    type = ParsedSubdomainMeshGenerator
    input = 'gen'
    expression = 'x > ${l}'
    block_id = 2
  []
  [redo_bottom]
    type = ParsedGenerateSideset
    input = 'remove'
    combinatorial_geometry = 'x > -1e8'
    included_subdomains = '0'
    included_boundaries = 'bottom'
    new_sideset_name = 'bottom_v2'
  []
  [redo_top]
    type = ParsedGenerateSideset
    input = 'redo_bottom'
    combinatorial_geometry = 'x > -1e8'
    included_subdomains = '0'
    included_boundaries = 'top'
    new_sideset_name = 'top_v2'
  []
  [redo_right]
    type = SideSetsBetweenSubdomainsGenerator
    input = 'redo_top'
    primary_block = '0'
    paired_block = '2'
    new_boundary = 'right_v2'
  []
[]

[Problem]
  type = NavierStokesProblem
  extra_tag_matrices = 'mass'
  mass_matrix = 'mass'
  use_pressure_mass_matrix = true
  kernel_coverage_check = false
[]

[AuxVariables]
  [vel_mag]
    family = L2_HIERARCHIC
    order = SECOND
    block = 0
  []
[]

[AuxKernels]
  [vel_mag]
    type = VectorMagnitudeAux
    variable = vel_mag
    x = vel_x
    y = vel_y
  []
[]

[Variables]
  [vel_x]
    family = L2_HIERARCHIC
    order = SECOND
    block = 0
  []
  [vel_y]
    family = L2_HIERARCHIC
    order = SECOND
    block = 0
  []
  [pressure]
    family = L2_HIERARCHIC
    order = FIRST
    block = 0
  []
  [vel_bar_x]
    family = SIDE_HIERARCHIC
    order = SECOND
    block = 0
  []
  [vel_bar_y]
    family = SIDE_HIERARCHIC
    order = SECOND
    block = 0
  []
  [pressure_bar]
    family = SIDE_HIERARCHIC
    order = SECOND
    block = 0
  []
[]

[HDGKernels]
  [momentum_x_convection]
    type = AdvectionIPHDGKernel
    variable = vel_x
    face_variable = vel_bar_x
    velocity = 'velocity'
    coeff = ${rho}
  []
  [momentum_x_diffusion]
    type = NavierStokesStressIPHDGKernel
    variable = vel_x
    face_variable = vel_bar_x
    diffusivity = 'mu'
    alpha = ${alpha}
    pressure_variable = pressure
    pressure_face_variable = pressure_bar
    component = 0
  []
  [momentum_y_convection]
    type = AdvectionIPHDGKernel
    variable = vel_y
    face_variable = vel_bar_y
    velocity = 'velocity'
    coeff = ${rho}
  []
  [momentum_y_diffusion]
    type = NavierStokesStressIPHDGKernel
    variable = vel_y
    face_variable = vel_bar_y
    diffusivity = 'mu'
    alpha = ${alpha}
    pressure_variable = pressure
    pressure_face_variable = pressure_bar
    component = 1
  []
  [pressure_convection]
    type = AdvectionIPHDGKernel
    variable = pressure
    face_variable = pressure_bar
    velocity = 'velocity'
    coeff = '${fparse -rho}'
    self_advection = false
  []

  [u_jump]
    type = MassFluxPenaltyIPHDG
    variable = vel_x
    u = vel_x
    v = vel_y
    u_face = vel_bar_x
    v_face = vel_bar_y
    component = 0
    gamma = ${gamma}
  []
  [v_jump]
    type = MassFluxPenaltyIPHDG
    variable = vel_y
    u = vel_x
    v = vel_y
    u_face = vel_bar_x
    v_face = vel_bar_y
    component = 1
    gamma = ${gamma}
  []
  [pb_mass]
    type = MassMatrixHDG
    variable = pressure_bar
    matrix_tags = 'mass'
    density = '${fparse -1/gamma}'
  []
[]

[BCs]
  [momentum_x_diffusion_walls]
    type = NavierStokesStressIPHDGDirichletBC
    boundary = 'left bottom_v2 right_v2'
    variable = vel_x
    face_variable = vel_bar_x
    pressure_variable = pressure
    pressure_face_variable = pressure_bar
    alpha = ${alpha}
    functor = '0'
    diffusivity = 'mu'
    component = 0
  []
  [momentum_x_diffusion_top]
    type = NavierStokesStressIPHDGDirichletBC
    boundary = 'top_v2'
    variable = vel_x
    face_variable = vel_bar_x
    pressure_variable = pressure
    pressure_face_variable = pressure_bar
    alpha = ${alpha}
    functor = '${U}'
    diffusivity = 'mu'
    component = 0
  []
  [momentum_y_diffusion_all]
    type = NavierStokesStressIPHDGDirichletBC
    boundary = 'left bottom_v2 right_v2 top_v2'
    variable = vel_y
    face_variable = vel_bar_y
    pressure_variable = pressure
    pressure_face_variable = pressure_bar
    alpha = ${alpha}
    functor = '0'
    diffusivity = 'mu'
    component = 1
  []

  [mass_convection]
    type = AdvectionIPHDGPrescribedFluxBC
    face_variable = pressure_bar
    variable = pressure
    velocity = 'velocity'
    coeff = '${fparse -rho}'
    self_advection = false
    boundary = 'left bottom_v2 top_v2 right_v2'
    prescribed_normal_flux = 0
  []

  [pb_mass]
    type = MassMatrixIntegratedBC
    variable = pressure_bar
    matrix_tags = 'mass'
    boundary = 'left right_v2 bottom_v2 top_v2'
    density = '${fparse -1/gamma}'
  []

  [u_jump_walls]
    type = MassFluxPenaltyBC
    variable = vel_x
    u = vel_x
    v = vel_y
    component = 0
    boundary = 'left right_v2 bottom_v2'
    gamma = ${gamma}
    dirichlet_value = walls
  []
  [v_jump_walls]
    type = MassFluxPenaltyBC
    variable = vel_y
    u = vel_x
    v = vel_y
    component = 1
    boundary = 'left right_v2 bottom_v2'
    gamma = ${gamma}
    dirichlet_value = walls
  []
  [u_jump_top]
    type = MassFluxPenaltyBC
    variable = vel_x
    u = vel_x
    v = vel_y
    component = 0
    boundary = 'top_v2'
    gamma = ${gamma}
    dirichlet_value = top_vel
  []
  [v_jump_top]
    type = MassFluxPenaltyBC
    variable = vel_y
    u = vel_x
    v = vel_y
    component = 1
    boundary = 'top_v2'
    gamma = ${gamma}
    dirichlet_value = top_vel
  []
[]

[Functions]
  [top_vel]
    type = ParsedVectorFunction
    expression_x = ${U}
  []
  [walls]
    type = ParsedVectorFunction
  []
  [reynolds]
    type = ParsedFunction
    expression = '10^(log10(${starting_re}) + (t - 1) * ${step_length})'
  []
[]

[Materials]
  [const]
    type = ADGenericConstantMaterial
    prop_names = 'rho'
    prop_values = '${rho}'
  []
  [vel]
    type = ADVectorFromComponentVariablesMaterial
    vector_prop_name = 'velocity'
    u = vel_x
    v = vel_y
    block = 0
  []
  [mu]
    type = ADParsedMaterial
    functor_names = 'reynolds'
    functor_symbols = 'reynolds'
    property_name = 'mu'
    expression = '${U} * ${l} / reynolds'
  []
[]

[Preconditioning]
  [FSP]
    type = SCFSP
    topsplit = 'up'
    [up]
      splitting = 'u p'
      splitting_type = schur
      petsc_options = '-ksp_monitor'
      petsc_options_iname = '-pc_fieldsplit_schur_fact_type  -pc_fieldsplit_schur_precondition -ksp_gmres_restart -ksp_type -ksp_pc_side -ksp_rtol -ksp_max_it -ksp_atol'
      petsc_options_value = 'full                            self                              300                fgmres    right        1e-4      30          1e-9'
    []
    [u]
      vars = 'vel_bar_x vel_bar_y'
      petsc_options = '-ksp_converged_reason'
      petsc_options_iname = '-pc_type -ksp_type -ksp_rtol -ksp_gmres_restart -ksp_pc_side -pc_factor_mat_solver_type -ksp_max_it'
      petsc_options_value = 'ilu      gmres     1e-2      300                right        strumpack                  30'
    []
    [p]
      vars = 'pressure_bar'
      petsc_options = '-ksp_converged_reason'
      petsc_options_iname = '-pc_type -ksp_type -ksp_rtol -ksp_gmres_restart -ksp_pc_side -pc_factor_mat_solver_type -ksp_max_it'
      petsc_options_value = 'ilu      gmres     1e-2      300                right        strumpack                  30'
    []
  []
[]

[Executioner]
  type = Transient
  num_steps = ${num_steps}
  petsc_options_iname = '-ksp_type'
  petsc_options_value = 'preonly'
[]

[Outputs]
  print_linear_residuals = 'false'
  csv = true
[]

[Postprocessors]
  [Re]
    type = FunctionValuePostprocessor
    function = 'reynolds'
  []
  [pressure_average]
    type = ElementAverageValue
    variable = pressure
    block = 0
  []
  [vel_average]
    type = ElementAverageValue
    variable = vel_mag
    block = 0
  []
[]

[Correctors]
  [set_pressure]
    type = NSPressurePin
    pin_type = 'average'
    variable = pressure
    pressure_average = 'pressure_average'
  []
[]

Input Parameters
Input Files