LBMComputeDensity

Compute object for macroscopic density reconstruction.

This compute object sums LBM distribution functions along the last dimension to compute density or temperature or any other macroscopic scalar buffer depending on the physics of the problem.

Overview

Given a distribution tensor f(x,y,z, q), computes the macroscopic scalar quantity rho(x,y,z) = sum_q f. Values on solid cells are set to zero using the lattice mask.

Example Input File Syntax

[TensorComputes]
  [Solve]
    [rho]
      type = LBMComputeDensity
      buffer = rho
      f = f
    []
  []
[]
!listing-end

Input Parameters

bufferThe buffer this compute is writing to
C++ Type:std::string
Controllable:No
Description:The buffer this compute is writing to
fDistribution function
C++ Type:std::string
Controllable:No
Description:Distribution function

Required 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:No
Description:Set the enabled status of the MooseObject.

Advanced Parameters

Input Files

(examples/lbm/Phase-field/static_bubble_3d.i)
(test/tests/lbm/vertical_density_bcs.i)
(examples/lbm/natural_convection_3D/convection.i)
(test/tests/lbm/convective_outflow_2d_right.i)
(test/tests/lbm/nee_3d_box_first.i)
(test/tests/lbm/convective_outflow_2d_left.i)
(test/tests/lbm/mixed_bcs_d3q19_reverse.i)
(examples/lbm/rarefied_gas/channel.i)
(test/tests/lbm/convective_outflow_3d_left.i)
(test/tests/lbm/convective_outflow_3d_front.i)
(test/tests/lbm/channel3D.i)
(test/tests/lbm/convective_outflow_3d_bottom.i)
(test/tests/lbm/nee_3d_box_second.i)
(test/tests/lbm/dirichlet_wall.i)
(test/tests/lbm/mixed_bcs_d3q19.i)
(test/tests/lbm/horizontal_mixed_bcs_d2q9_reverse.i)
(examples/lbm/Rayleigh-Benard/solve.i)
(test/tests/lbm/horizontal_mixed_bcs_d2q9.i)
(test/tests/lbm/nee_2d_channel.i)
(test/tests/lbm/convective_outflow_2d_bottom.i)
(test/tests/lbm/convective_outflow_2d_top.i)
(test/tests/lbm/dirichlet_regional.i)
(examples/lbm/Phase-field/droplet_impact.i)
(examples/lbm/Phase-field/spinodal_decomposition.i)
(test/tests/lbm/channel2D.i)
(examples/lbm/Karman-vortex/cylinder.i)
(test/tests/lbm/convective_outflow_3d_right.i)
(test/tests/lbm/mixed_bcs_d3q27.i)
(test/tests/lbm/obstacle.i)
(test/tests/lbm/neumann_regional.i)
(test/tests/lbm/rarefied_gas.i)
(test/tests/lbm/nee_2d_all_walls.i)
(test/tests/lbm/isotropic_stencil_mrt.i)
(test/tests/lbm/neumann_wall.i)
(test/tests/lbm/neumann_box.i)
(examples/lbm/Phase-field/static_bubble.i)
(test/tests/lbm/phase.i)
(test/tests/lbm/vertical_velocity_bcs.i)
(test/tests/lbm/convective_outflow_3d_top.i)
(test/tests/lbm/smagorinsky_mrt.i)
(examples/lbm/Formula1-aerodynamics/f1.i)
(test/tests/lbm/mixed_bcs_d3q27_reverse.i)
(test/tests/lbm/phase_3D.i)
(test/tests/lbm/convective_outflow_3d_back.i)
(examples/lbm/Phase-field/layered_poiseuille.i)
(test/tests/lbm/dirichlet_box.i)
(examples/lbm/Pebbles/pebbles.i)

(examples/lbm/Phase-field/static_bubble_3d.i)

#
# 3D Static bubble test case
# Based on the 2D static bubble (PHYSICAL REVIEW E 97, 033309 - Section III.A)
# Small domain to verify 3D gradient/laplacian correctness
#

# Domain (small for quick testing)
Nx = 40
Ny = 40
Nz = 40

# Bubble parameters
Cx = '${Nx}/2.0'
Cy = '${Ny}/2.0'
Cz = '${Nz}/2.0'
R = 12

# Fluid properties
rho_l = 1000.0
rho_g = 1.0
nu_l = 0.1
nu_g = 0.1
sigma = 0.001

# Phase field parameters
tau_h = 0.7
D = 4

[Domain]
  dim = 3
  nx = '${Nx}'
  ny = '${Ny}'
  nz = '${Nz}'
  xmax = '${Nx}'
  ymax = '${Ny}'
  zmax = '${Nz}'
  device_names = 'cpu'
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y Z'
[]

[Stencil]
  [d3q27]
    type = LBMD3Q27
  []
[]

[TensorBuffers]
  # Macroscopic phase field variables
  [phi]
    type = LBMTensorBuffer
    buffer_type = ms
  []
  [grad_phi]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  [laplacian_phi]
    type = LBMTensorBuffer
    buffer_type = ms
  []
  [mu]
    type = LBMTensorBuffer
    buffer_type = ms
  []
  [forces]
    type = LBMTensorBuffer
    buffer_type = mv
  []

  # Macroscopic hydrodynamic variables
  [velocity]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  [pressure]
    type = LBMTensorBuffer
    buffer_type = ms
  []
  [rho]
    type = LBMTensorBuffer
    buffer_type = ms
  []

  # LBM phase field variables
  [h]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [h_post_collision]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [h_eq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [relaxation_tensor]
    type = LBMTensorBuffer
    buffer_type = ms
  []

  # LBM hydrodynamic variables
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [f_post_collision]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [f_eq]
    type = LBMTensorBuffer
    buffer_type = df
  []
[]

[TensorComputes/Initialize]
  [phi_init]
    type = ParsedCompute
    buffer = phi
    expression = '0.5 + 0.5 * tanh(2*(R - sqrt((x - Cx)^2 + (y - Cy)^2 + (z - Cz)^2)) / D)'
    constant_names = 'Cx Cy Cz R D'
    constant_expressions = '${Cx} ${Cy} ${Cz} ${R} ${D}'
    extra_symbols = true
  []
  [grad_phi_init]
    type = LBMIsotropicGradient
    buffer = grad_phi
    scalar_field = phi
  []
  [rho_init]
    type = ParsedCompute
    buffer = rho
    extra_symbols = true
    expression = 'phi*(rho_l - rho_g) + rho_g'
    constant_names = 'rho_l rho_g'
    constant_expressions = '${rho_l} ${rho_g}'
    inputs = phi
  []
  [pressure_init]
    type = LBMConstantTensor
    buffer = pressure
    constants = 0.3
  []
  [h_eq_init]
    type = LBMPhaseEquilibrium
    buffer = h_eq
    phi = phi
    velocity = velocity
  []
  [h_post_collision_init]
    type = LBMPhaseEquilibrium
    buffer = h_post_collision
    phi = phi
    velocity = velocity
  []
  [h_init]
    type = LBMPhaseEquilibrium
    buffer = h
    phi = phi
    velocity = velocity
  []
  [f_eq_init]
    type = LBMPressureCorrectedEquilibrium
    buffer = f_eq
    rho = rho
    velocity = velocity
    pressure = pressure
  []
  [f_post_collision_init]
    type = ParsedCompute
    buffer = f_post_collision
    expression = 'f_eq'
    inputs = f_eq
  []
  [f_init]
    type = ParsedCompute
    buffer = f
    expression = 'f_eq'
    inputs = f_eq
  []
[]

[TensorComputes/Solve]
  # Phase Field
  [compute_phi]
    type = LBMComputeDensity
    buffer = phi
    f = h
  []
  [grad_phi]
    type = LBMIsotropicGradient
    buffer = grad_phi
    scalar_field = phi
  []
  [laplacian_phi]
    type = LBMIsotropicLaplacian
    buffer = laplacian_phi
    scalar_field = phi
  []
  [potential]
    type = LBMComputeChemicalPotential
    buffer = mu
    phi = phi
    laplacian_phi = laplacian_phi
    thickness = D
    sigma = sigma
  []
  [forces]
    type = LBMComputeSurfaceForces
    buffer = forces
    chemical_potential = mu
    grad_phi = grad_phi
  []
  # Hydrodynamics
  [density]
    type = ParsedCompute
    buffer = rho
    extra_symbols = true
    expression = 'phi*(rho_l - rho_g) + rho_g'
    constant_names = 'rho_l rho_g'
    constant_expressions = '${rho_l} ${rho_g}'
    inputs = phi
  []
  [velocity]
    type = LBMComputeVelocity
    buffer = velocity
    f = f
    rho = rho
    enable_forces = true
    forces = forces
  []
  [h_eq]
    type = LBMPhaseEquilibrium
    buffer = h_eq
    phi = phi
    velocity = velocity
  []
  [phase_collision]
    type = LBMBGKCollision
    buffer = h_post_collision
    f = h
    feq = h_eq
    tau0 = tau_h
  []
  [apply_forces_phase]
    type = LBMAllenCahnSource
    buffer = h_post_collision
    phi = phi
    velocity = velocity
    grad_phi = grad_phi
    tau = tau_h
    thickness = D
  []
  [relaxation_tensor]
    type = ParsedCompute
    buffer = relaxation_tensor
    extra_symbols = true
    expression = '(phi*(nu_l - nu_g) + nu_g)/cs2+0.5'
    constant_names = 'nu_l nu_g cs2'
    constant_expressions = '${nu_l} ${nu_g} 0.3333'
    inputs = phi
  []
  [pressure]
    type = LBMPhaseFieldPressure
    buffer = pressure
    f = f
    velocity = velocity
    grad_phi = grad_phi
    rho = rho
    rho_l = '${rho_l}'
    rho_g = '${rho_g}'
  []
  [f_eq]
    type = LBMPressureCorrectedEquilibrium
    buffer = f_eq
    rho = rho
    velocity = velocity
    pressure = pressure
  []
  [collision]
    type = LBMBGKCollision
    buffer = f_post_collision
    f = f
    feq = f_eq
    tau0 = 1.0
    is_dynamic_relaxation = true
    tau_tensor = relaxation_tensor
  []
  [apply_forces_hydro]
    type = LBMForceDistribution
    buffer = f_post_collision
    grad_phi = grad_phi
    velocity = velocity
    forces = forces
    tau_tensor = relaxation_tensor
    tau = 1.0
    rho_l = '${rho_l}'
    rho_g = '${rho_g}'
    is_dynamic_relaxation = true
  []
  [residual]
    type = LBMComputeResidual
    buffer = phi
    speed = phi
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = 'h f'
  f_old = 'h_post_collision f_post_collision'
  root_compute = residual
[]

[Postprocessors]
  [velocity_min]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MIN
  []
  [velocity_max]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = rho
    value_type = MIN
  []
  [density_max]
    type = TensorExtremeValuePostprocessor
    buffer = rho
    value_type = MAX
  []
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 10
  print_debug_output = true
  scalar_constant_names = 'tau_h D sigma'
  scalar_constant_values = '${tau_h} ${D} ${sigma}'
[]

[Executioner]
  type = Transient
  num_steps = 50
[]

[TensorOutputs]
  [xdmf]
    type = XDMFTensorOutput
    buffer = 'phi rho velocity forces'
    output_mode = 'Cell Cell Cell Cell'
    enable_hdf5 = true
    # transpose = false
  []
[]

(test/tests/lbm/vertical_density_bcs.i)

[Domain]
  dim = 2
  nx = 10
  ny = 10
  mesh_mode = DUMMY
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y'
[]

[Stencil]
  [d2q9]
    type = LBMD2Q9
  []
[]

[TensorBuffers]
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [f_bounce_back]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [velocity]
    type=LBMTensorBuffer
    buffer_type = mv
  []
  [density]
    type=LBMTensorBuffer
    buffer_type = ms
  []
[]

[TensorComputes]
  [Initialize]
    [initial_density]
      type = LBMConstantTensor
      buffer = density
      constants = 1.0
    []
    [initial_velocity]
      type = LBMConstantTensor
      buffer = velocity
      constants = '0.0001 0.0005'
    []
    [initial_f]
      type = LBMEquilibrium
      buffer = f
      bulk = density
      velocity = velocity
    []
    [initial_f_bb]
      type = LBMEquilibrium
      buffer = f_bounce_back
      bulk = density
      velocity = velocity
    []
  []
  [Solve]
    [density]
      type = LBMComputeDensity
      buffer = density
      f = f
    []
    [velocity]
      type = LBMComputeVelocity
      buffer = velocity
      f = f
      rho = density
    []
  []
  [Boundary]
    [left]
      type = LBMBounceBack
      buffer = f
      f_old = f_bounce_back
      boundary = left
    []
    [right]
      type = LBMBounceBack
      buffer = f
      f_old = f_bounce_back
      boundary = right
    []
    [top]
      type = LBMFixedZerothOrderBC
      buffer = f
      f = f
      value = 1.1
      boundary = top
    []
    [bottom]
      type = LBMFixedZerothOrderBC
      buffer = f
      f = f
      value = 1.00000
      boundary = bottom
    []
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = f
  f_old = f_bounce_back
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 2
[]

[Postprocessors]
  [velocity_min]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MIN
  []
  [velocity_max]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MIN
  []
  [densty_max]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MAX
  []
[]

[Executioner]
  type = Transient
  num_steps = 5
[]

[Outputs]
  file_base = vertical_density_bcs
  csv = true
[]

(examples/lbm/natural_convection_3D/convection.i)

[Domain]
  dim = 3
  nx = 270
  ny = 270
  nz = 405
  xmax = 270
  ymax = 270
  zmax = 405
  device_names='cpu'
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y Z'
  floating_precision=SINGLE
[]

[Stencil]
  [d3q19]
    type = LBMD3Q19
  []
[]

[TensorBuffers]
  # Simulation binary media
  [binary_media]
    type = LBMTensorBuffer
    file = binary_media.h5
    buffer_type = ms
    is_integer = true
  []

  # Density distribution functions
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [feq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [fpc]
    type = LBMTensorBuffer
    buffer_type = df
  []
  # Temperature distribution functions
  [g]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [geq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [gpc]
    type = LBMTensorBuffer
    buffer_type = df
  []
  # Fluid macroscopic variables: density and velocity
  [density]
    type = LBMTensorBuffer
    buffer_type = ms
  []
  [velocity]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  [speed]
    type=LBMTensorBuffer
    buffer_type = ms
  []
  # Temperature macroscpic variables: temperature and 'velocity'
  [T]
    type = LBMTensorBuffer
    buffer_type = ms
  []
  # Forces
  [F]
    type = LBMTensorBuffer
    buffer_type = mv
  []
[]

[TensorComputes]
  [Initialize]
    [density]
      type = LBMConstantTensor
      buffer = density
      constants = 'rho0'
    []
    [velocity]
      type = LBMConstantTensor
      buffer = velocity
      constants = '0.0 0.0 0.0'
    []
    [temperature]
      type = LBMConstantTensor
      buffer = T
      constants = T_C
    []
    [equilibrium_fluid]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [equilibrium_fluid_total]
      type = LBMEquilibrium
      buffer = f
      bulk = density
      velocity = velocity
    []
    [equilibrium_fluid_pc]
      type = LBMEquilibrium
      buffer = fpc
      bulk = density
      velocity = velocity
    []
    [equilibrium_temperature]
      type = LBMEquilibrium
      buffer = geq
      bulk = T
      velocity = velocity
    []
    [equilibrium_temperature_total]
      type = LBMEquilibrium
      buffer = g
      bulk = T
      velocity = velocity
    []
    [equilibrium_temperature_pc]
      type = LBMEquilibrium
      buffer = gpc
      bulk = T
      velocity = velocity
    []
  []

  #### Compute ####
  [Solve]

    # For temperature
    [Temperature]
      type = LBMComputeDensity
      buffer = T
      f = g
    []

    # For fluid
    [Fluid_density]
      type = LBMComputeDensity
      buffer = density
      f = f
    []

    [Fluid_velocity]
      type = LBMComputeVelocity
      buffer = velocity
      f = f
      rho = density
      forces = F
      enable_forces = true
    []

    # For temperature
    [Equilibrium_temperature]
      type = LBMEquilibrium
      buffer = geq
      bulk = T
      velocity = velocity
    []

    [Collision_temperature]
      type = LBMBGKCollision
      buffer = gpc
      f = g
      feq = geq
      tau0 = tau_T
    []

    # For fluid
    [Compute_forces]
      type = LBMComputeForces
      buffer = F
      rho0 = 'rho0'
      temperature = T
      T0 = 1.00
      enable_buoyancy = true
      gravity = g
      gravity_direction = 2
    []

    [Equilibrium_fluid]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []

    [Collision_fluid]
      type = LBMBGKCollision
      buffer = fpc
      f = f
      feq = feq
      tau0 = tau_f
    []

    [Apply_forces]
      type = LBMApplyForces
      buffer = fpc
      velocity = velocity
      rho = density
      forces = F
      tau0 = tau_f
    []

    [Residual]
      type = LBMComputeResidual
      speed = T
      # TODO this buffer is redundant, but avoids 'missing parameter' error
      buffer = T
    []
  []

  #### Boundary ####
  [Boundary]
    ##### for fluid
    [wall]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = wall
    []
    ##### for temperature
    [heat_source]
      type = LBMNeumannBC
      buffer = g
      f_old = gpc
      feq=geq
      velocity = velocity
      rho = T
      gradient = 0.001
      region_id = 3
      boundary = regional
    []
    [walls]
      type = LBMDirichletBC
      buffer = g
      f_old = gpc
      feq=geq
      velocity = velocity
      rho = T
      value = 1.0
      region_id = 2
      boundary = regional
    []
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = 'f g'
  f_old = 'fpc gpc'
  root_compute = Residual
[]

[Problem]
  type = LatticeBoltzmannProblem
  binary_media = binary_media
  scalar_constant_names = 'rho0 T_C T_H tau_f tau_T g'
  scalar_constant_values = '1.0 1.0 1.05 0.55 0.55 0.01'
  substeps = 100
  print_debug_output = true
[]

[Executioner]
  type = Transient
  num_steps = 1000
[]

[TensorOutputs]
  [xdmf2]
    type = XDMFTensorOutput
    buffer = 'T density velocity'
    output_mode = 'Cell Cell Cell'
    enable_hdf5 = true
  []
[]

(test/tests/lbm/convective_outflow_2d_right.i)

[Domain]
  dim = 2
  nx = 10
  ny = 10
  mesh_mode = DUMMY
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y'
[]

[Stencil]
  [d2q9]
    type = LBMD2Q9
  []
[]

[TensorBuffers]
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [feq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [fpc]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [velocity]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  [density]
    type = LBMTensorBuffer
    buffer_type = ms
  []
[]

[TensorComputes]
  [Initialize]
    [initial_density]
      type = LBMConstantTensor
      buffer = density
      constants = 1.0
    []
    [initial_velocity]
      type = LBMConstantTensor
      buffer = velocity
      constants = '0.0 0.0'
    []
    [initial_feq]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [initial_f]
      type = LBMEquilibrium
      buffer = f
      bulk = density
      velocity = velocity
    []
    [initial_fpc]
      type = LBMEquilibrium
      buffer = fpc
      bulk = density
      velocity = velocity
    []
  []
  [Solve]
    [equilibrium]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [collision]
      type = LBMBGKCollision
      buffer = fpc
      f = f
      feq = feq
      tau0 = 1.0
    []
    [density]
      type = LBMComputeDensity
      buffer = density
      f = f
    []
    [velocity]
      type = LBMComputeVelocity
      buffer = velocity
      f = f
      rho = density
    []
  []
  [Boundary]
    [left]
      type = LBMNonEquilibriumExtrapolation
      buffer = f
      prescribe_type = velocity
      ux = 0.005
      uy = 0.0
      order = first
      boundary = left
    []
    [right]
      type = LBMConvectiveOutflow
      buffer = f
      f_old = f
      convection_velocity = auto
      boundary = right
    []
    [top]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = top
    []
    [bottom]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = bottom
    []
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = f
  f_old = fpc
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 1
[]

[Postprocessors]
  [velocity_min]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MIN
  []
  [velocity_max]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MIN
  []
  [density_max]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MAX
  []
[]

[Executioner]
  type = Transient
  num_steps = 5
[]

[Outputs]
  file_base = convective_outflow_2d_right
  csv = true
[]

(test/tests/lbm/nee_3d_box_first.i)

[Domain]
  dim = 3
  nx = 8
  ny = 8
  nz = 8
  mesh_mode = DUMMY
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y Z'
[]

[Stencil]
  [d3q19]
    type = LBMD3Q19
  []
[]

[TensorBuffers]
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [feq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [fpc]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [velocity]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  [density]
    type = LBMTensorBuffer
    buffer_type = ms
  []
[]

[TensorComputes]
  [Initialize]
    [initial_density]
      type = LBMConstantTensor
      buffer = density
      constants = 1.0
    []
    [initial_velocity]
      type = LBMConstantTensor
      buffer = velocity
      constants = '0.0 0.0 0.0'
    []
    [initial_feq]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [initial_f]
      type = LBMEquilibrium
      buffer = f
      bulk = density
      velocity = velocity
    []
    [initial_fpc]
      type = LBMEquilibrium
      buffer = fpc
      bulk = density
      velocity = velocity
    []
  []
  [Solve]
    [equilibrium]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [collision]
      type = LBMBGKCollision
      buffer = fpc
      f = f
      feq = feq
      tau0 = 1.0
    []
    [density]
      type = LBMComputeDensity
      buffer = density
      f = f
    []
    [velocity]
      type = LBMComputeVelocity
      buffer = velocity
      f = f
      rho = density
      add_body_force = true
      body_force_x = 0.0001
    []
  []
  [Boundary]
    [left]
      type = LBMNonEquilibriumExtrapolation
      buffer = f
      prescribe_type = velocity
      ux = 0.0
      uy = 0.0
      uz = 0.0
      order = first
      boundary = left
    []
    [right]
      type = LBMNonEquilibriumExtrapolation
      buffer = f
      prescribe_type = velocity
      ux = 0.0
      uy = 0.0
      uz = 0.0
      order = first
      boundary = right
    []
    [top]
      type = LBMNonEquilibriumExtrapolation
      buffer = f
      prescribe_type = velocity
      ux = 0.0
      uy = 0.0
      uz = 0.0
      order = first
      boundary = top
    []
    [bottom]
      type = LBMNonEquilibriumExtrapolation
      buffer = f
      prescribe_type = velocity
      ux = 0.0
      uy = 0.0
      uz = 0.0
      order = first
      boundary = bottom
    []
    [front]
      type = LBMNonEquilibriumExtrapolation
      buffer = f
      prescribe_type = velocity
      ux = 0.0
      uy = 0.0
      uz = 0.0
      order = first
      boundary = front
    []
    [back]
      type = LBMNonEquilibriumExtrapolation
      buffer = f
      prescribe_type = velocity
      ux = 0.0
      uy = 0.0
      uz = 0.0
      order = first
      boundary = back
    []
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = f
  f_old = fpc
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 1
[]

[Postprocessors]
  [velocity_min]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MIN
  []
  [velocity_max]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MIN
  []
  [density_max]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MAX
  []
[]

[Executioner]
  type = Transient
  num_steps = 5
[]

[Outputs]
  file_base = nee_3d_box_first
  csv = true
[]

(test/tests/lbm/convective_outflow_2d_left.i)

[Domain]
  dim = 2
  nx = 10
  ny = 10
  mesh_mode = DUMMY
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y'
[]

[Stencil]
  [d2q9]
    type = LBMD2Q9
  []
[]

[TensorBuffers]
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [feq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [fpc]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [velocity]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  [density]
    type = LBMTensorBuffer
    buffer_type = ms
  []
[]

[TensorComputes]
  [Initialize]
    [initial_density]
      type = LBMConstantTensor
      buffer = density
      constants = 1.0
    []
    [initial_velocity]
      type = LBMConstantTensor
      buffer = velocity
      constants = '0.0 0.0'
    []
    [initial_feq]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [initial_f]
      type = LBMEquilibrium
      buffer = f
      bulk = density
      velocity = velocity
    []
    [initial_fpc]
      type = LBMEquilibrium
      buffer = fpc
      bulk = density
      velocity = velocity
    []
  []
  [Solve]
    [equilibrium]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [collision]
      type = LBMBGKCollision
      buffer = fpc
      f = f
      feq = feq
      tau0 = 1.0
    []
    [density]
      type = LBMComputeDensity
      buffer = density
      f = f
    []
    [velocity]
      type = LBMComputeVelocity
      buffer = velocity
      f = f
      rho = density
    []
  []
  [Boundary]
    [right]
      type = LBMNonEquilibriumExtrapolation
      buffer = f
      prescribe_type = velocity
      ux = -0.005
      uy = 0.0
      order = first
      boundary = right
    []
    [left]
      type = LBMConvectiveOutflow
      buffer = f
      f_old = f
      convection_velocity = auto
      boundary = left
    []
    [top]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = top
    []
    [bottom]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = bottom
    []
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = f
  f_old = fpc
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 1
[]

[Postprocessors]
  [velocity_min]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MIN
  []
  [velocity_max]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MIN
  []
  [density_max]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MAX
  []
[]

[Executioner]
  type = Transient
  num_steps = 5
[]

[Outputs]
  file_base = convective_outflow_2d_left
  csv = true
[]

(test/tests/lbm/mixed_bcs_d3q19_reverse.i)

[Domain]
  dim = 3
  nx = 10
  ny = 10
  nz = 10
  mesh_mode = DUMMY
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y Z'
[]

[Stencil]
  [d3q19]
    type = LBMD3Q19
  []
[]

[TensorBuffers]
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [f_bounce_back]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [velocity]
    type=LBMTensorBuffer
    buffer_type = mv
  []
  [density]
    type=LBMTensorBuffer
    buffer_type = ms
  []
[]

[TensorComputes]
  [Initialize]
    [initial_density]
      type = LBMConstantTensor
      buffer = density
      constants = 1.0
    []
    [initial_velocity]
      type = LBMConstantTensor
      buffer = velocity
      constants = '0.0001 0.0005 0.0'
    []
    [initial_f]
      type = LBMEquilibrium
      buffer = f
      bulk = density
      velocity = velocity
    []
    [initial_f_bb]
      type = LBMEquilibrium
      buffer = f_bounce_back
      bulk = density
      velocity = velocity
    []
  []
  [Solve]
    [density]
      type = LBMComputeDensity
      buffer = density
      f = f
    []
    [velocity]
      type = LBMComputeVelocity
      buffer = velocity
      f = f
      rho = density
      add_body_force = true
      body_force_x = 0.0001
    []
  []
  [Boundary]
    [right]
      type = LBMFixedZerothOrderBC
      buffer = f
      f = f
      value = 1.11
      boundary = right
    []
    [left]
      type = LBMFixedFirstOrderBC
      buffer = f
      f = f
      value = 0.0001
      boundary = left
    []
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = f
  f_old = f_bounce_back
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 1
[]

[Postprocessors]
  [velocity_min]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MIN
  []
  [velocity_max]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MIN
  []
  [densty_max]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MAX
  []
[]

[Executioner]
  type = Transient
  num_steps = 5
[]

[Outputs]
  file_base = mixed_bcs_d3q19_reverse
  csv = true
[]

(examples/lbm/rarefied_gas/channel.i)

[Domain]
  dim = 2
  nx = 102
  ny = 102
  xmax = 102
  ymax = 102
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y'
[]

[Stencil]
  [d2q9]
    type = LBMD2Q9
  []
[]

[TensorBuffers]
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [feq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [fpc]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [velocity]
    type=LBMTensorBuffer
    buffer_type = mv
  []
  [density]
    type=LBMTensorBuffer
    buffer_type = ms
  []
  [speed]
    type=LBMTensorBuffer
    buffer_type = ms
  []
  [domain]
    type=LBMTensorBuffer
    file = binary_media.h5
    buffer_type = ms
    is_integer = true
  []
  [local_pore]
    type=LBMTensorBuffer
    file = local_pore_size.h5
    buffer_type = ms
    is_integer = false
  []
  [Kn]
    type=LBMTensorBuffer
    file = Kn.h5
    buffer_type = ms
    is_integer = false
  []
  [relaxation_matrix]
    type = LBMTensorBuffer
    buffer_type = df
  []
[]

[TensorComputes]
  [Initialize]
    [initial_density]
      type = LBMConstantTensor
      buffer = density
      constants = 0.2355545440759889
    []
    [initial_velocity]
      type = LBMConstantTensor
      buffer = velocity
      constants = '0.0 0.0'
    []
    [initial_equilibrium]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [initial_distribution]
      type = LBMEquilibrium
      buffer = f
      bulk = density
      velocity = velocity
    []
    [initial_distribution_pc]
      type = LBMEquilibrium
      buffer = fpc
      bulk = density
      velocity = velocity
    []
    [relaxation_matrix_init]
      type = LBMComputeEffectiveRelaxation
      buffer = relaxation_matrix
      local_pore_size = local_pore
      local_Knudsen_number = Kn
      mfp = 7.904614716131531e-10
      dx = 0.122e-9
      A2 = 0.82
    []
  []
  [Solve]
    [density]
      type = LBMComputeDensity
      buffer = density
      f = f
    []
    [velocity]
      type = LBMComputeVelocity
      buffer = velocity
      f = f
      rho = density
      add_body_force = true
      body_force_x = 1.0e-8
    []
    [equilibrium]
      type=LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [collision]
      type = LBMMRTCollision
      buffer = fpc
      f = f
      feq = feq
      is_dynamic_relaxation = true
      local_relaxation_matrix = relaxation_matrix
      projection = true
    []
    [speed]
      type = LBMComputeVelocityMagnitude
      buffer = speed
      velocity = velocity
    []
    [residual]
      type = LBMComputeResidual
      buffer = speed
      speed = speed
    []
  []
  [Boundary]
    [wall]
      type = LBMSpecularReflectionBoundary
      buffer = f
      f_old = fpc
      local_Knudsen_number = Kn
      boundary = wall
    []
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = f
  f_old = fpc
  root_compute = residual
[]

[Postprocessors]
  [max_u]
    type = TensorExtremeValuePostprocessor
    buffer = speed
    value_type = MAX
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 2000
  print_debug_output = true
  binary_media = domain
[]

[Executioner]
  type = Transient
  num_steps = 2
[]

[TensorOutputs]
  [xdmf2]
    type = XDMFTensorOutput
    buffer = 'density velocity speed'
    output_mode = 'Cell Cell Cell'
    enable_hdf5 = true
    transpose=false
  []
[]

(test/tests/lbm/convective_outflow_3d_left.i)

[Domain]
  dim = 3
  nx = 8
  ny = 8
  nz = 8
  mesh_mode = DUMMY
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y Z'
[]

[Stencil]
  [d3q19]
    type = LBMD3Q19
  []
[]

[TensorBuffers]
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [feq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [fpc]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [velocity]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  [density]
    type = LBMTensorBuffer
    buffer_type = ms
  []
[]

[TensorComputes]
  [Initialize]
    [initial_density]
      type = LBMConstantTensor
      buffer = density
      constants = 1.0
    []
    [initial_velocity]
      type = LBMConstantTensor
      buffer = velocity
      constants = '0.0 0.0 0.0'
    []
    [initial_feq]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [initial_f]
      type = LBMEquilibrium
      buffer = f
      bulk = density
      velocity = velocity
    []
    [initial_fpc]
      type = LBMEquilibrium
      buffer = fpc
      bulk = density
      velocity = velocity
    []
  []
  [Solve]
    [equilibrium]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [collision]
      type = LBMBGKCollision
      buffer = fpc
      f = f
      feq = feq
      tau0 = 1.0
    []
    [density]
      type = LBMComputeDensity
      buffer = density
      f = f
    []
    [velocity]
      type = LBMComputeVelocity
      buffer = velocity
      f = f
      rho = density
    []
  []
  [Boundary]
    [right]
      type = LBMNonEquilibriumExtrapolation
      buffer = f
      prescribe_type = velocity
      ux = -0.005
      uy = 0.0
      uz = 0.0
      order = first
      boundary = right
    []
    [left]
      type = LBMConvectiveOutflow
      buffer = f
      f_old = f
      convection_velocity = Uc
      boundary = left
    []
    [top]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = top
    []
    [bottom]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = bottom
    []
    [front]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = front
    []
    [back]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = back
    []
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = f
  f_old = fpc
[]

[Problem]
  type = LatticeBoltzmannProblem
  scalar_constant_names = 'Uc'
  scalar_constant_values = '0.005'
  substeps = 1
[]

[Postprocessors]
  [velocity_min]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MIN
  []
  [velocity_max]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MIN
  []
  [density_max]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MAX
  []
[]

[Executioner]
  type = Transient
  num_steps = 5
[]

[Outputs]
  file_base = convective_outflow_3d_left
  csv = true
[]

(test/tests/lbm/convective_outflow_3d_front.i)

[Domain]
  dim = 3
  nx = 8
  ny = 8
  nz = 8
  mesh_mode = DUMMY
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y Z'
[]

[Stencil]
  [d3q19]
    type = LBMD3Q19
  []
[]

[TensorBuffers]
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [feq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [fpc]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [velocity]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  [density]
    type = LBMTensorBuffer
    buffer_type = ms
  []
[]

[TensorComputes]
  [Initialize]
    [initial_density]
      type = LBMConstantTensor
      buffer = density
      constants = 1.0
    []
    [initial_velocity]
      type = LBMConstantTensor
      buffer = velocity
      constants = '0.0 0.0 0.0'
    []
    [initial_feq]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [initial_f]
      type = LBMEquilibrium
      buffer = f
      bulk = density
      velocity = velocity
    []
    [initial_fpc]
      type = LBMEquilibrium
      buffer = fpc
      bulk = density
      velocity = velocity
    []
  []
  [Solve]
    [equilibrium]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [collision]
      type = LBMBGKCollision
      buffer = fpc
      f = f
      feq = feq
      tau0 = 1.0
    []
    [density]
      type = LBMComputeDensity
      buffer = density
      f = f
    []
    [velocity]
      type = LBMComputeVelocity
      buffer = velocity
      f = f
      rho = density
    []
  []
  [Boundary]
    [back]
      type = LBMNonEquilibriumExtrapolation
      buffer = f
      prescribe_type = velocity
      ux = 0.0
      uy = 0.0
      uz = -0.005
      order = first
      boundary = back
    []
    [front]
      type = LBMConvectiveOutflow
      buffer = f
      f_old = f
      convection_velocity = Uc
      boundary = front
    []
    [left]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = left
    []
    [right]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = right
    []
    [top]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = top
    []
    [bottom]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = bottom
    []
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = f
  f_old = fpc
[]

[Problem]
  type = LatticeBoltzmannProblem
  scalar_constant_names = 'Uc'
  scalar_constant_values = '0.005'
  substeps = 1
[]

[Postprocessors]
  [velocity_min]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MIN
  []
  [velocity_max]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MIN
  []
  [density_max]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MAX
  []
[]

[Executioner]
  type = Transient
  num_steps = 5
[]

[Outputs]
  file_base = convective_outflow_3d_front
  csv = true
[]

(test/tests/lbm/channel3D.i)

[Domain]
  dim = 3
  nx = 10
  ny = 10
  nz = 10
  mesh_mode = DUMMY
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y Z'
[]

[Stencil]
  [d3q19]
    type = LBMD3Q19
  []
[]

[TensorBuffers]
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [feq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [fpc]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [velocity]
    type=LBMTensorBuffer
    buffer_type = mv
  []
  [density]
    type=LBMTensorBuffer
    buffer_type = ms
  []
  [speed]
    type=LBMTensorBuffer
    buffer_type = ms
  []
[]

[TensorComputes]
  [Initialize]
    [initial_density]
      type = LBMConstantTensor
      buffer = density
      constants = 1.0
    []
    [initial_velocity]
      type = LBMConstantTensor
      buffer = velocity
      constants = '0.0 0.0 0.0'
    []
    [initial_equilibrium]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [initial_distribution]
      type = LBMEquilibrium
      buffer = f
      bulk = density
      velocity = velocity
    []
    [initial_distribution_pc]
      type = LBMEquilibrium
      buffer = fpc
      bulk = density
      velocity = velocity
    []
  []
  [Solve]
    [equilibrium]
      type=LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [collision]
      type=LBMBGKCollision
      buffer = fpc
      f = f
      feq = feq
      tau0 = 1.0
    []
    [density]
      type = LBMComputeDensity
      buffer = density
      f = f
    []
    [velocity]
      type = LBMComputeVelocity
      buffer = velocity
      f = f
      rho = density
      add_body_force = true
      body_force_x = 0.0001
    []
    [speed]
      type = LBMComputeVelocityMagnitude
      buffer = speed
      velocity = velocity
    []
    [residual]
      type = LBMComputeResidual
      buffer = speed
      speed = speed
    []
  []
  [Boundary]
    [top]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = top
    []
    [bottom]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = bottom
    []
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = f
  f_old = fpc
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 10
[]

[Postprocessors]
  [velocity_min]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MIN
  []
  [velocity_max]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MAX
  []
  [speed_min]
    type = TensorExtremeValuePostprocessor
    buffer = speed
    value_type = MIN
  []
  [speed_max]
    type = TensorExtremeValuePostprocessor
    buffer = speed
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MIN
  []
  [densty_max]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MAX
  []
[]

[Executioner]
  type = Transient
  num_steps = 11
[]

[Outputs]
  file_base = channel3D
  csv = true
[]

(test/tests/lbm/convective_outflow_3d_bottom.i)

[Domain]
  dim = 3
  nx = 8
  ny = 8
  nz = 8
  mesh_mode = DUMMY
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y Z'
[]

[Stencil]
  [d3q19]
    type = LBMD3Q19
  []
[]

[TensorBuffers]
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [feq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [fpc]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [velocity]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  [density]
    type = LBMTensorBuffer
    buffer_type = ms
  []
[]

[TensorComputes]
  [Initialize]
    [initial_density]
      type = LBMConstantTensor
      buffer = density
      constants = 1.0
    []
    [initial_velocity]
      type = LBMConstantTensor
      buffer = velocity
      constants = '0.0 0.0 0.0'
    []
    [initial_feq]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [initial_f]
      type = LBMEquilibrium
      buffer = f
      bulk = density
      velocity = velocity
    []
    [initial_fpc]
      type = LBMEquilibrium
      buffer = fpc
      bulk = density
      velocity = velocity
    []
  []
  [Solve]
    [equilibrium]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [collision]
      type = LBMBGKCollision
      buffer = fpc
      f = f
      feq = feq
      tau0 = 1.0
    []
    [density]
      type = LBMComputeDensity
      buffer = density
      f = f
    []
    [velocity]
      type = LBMComputeVelocity
      buffer = velocity
      f = f
      rho = density
    []
  []
  [Boundary]
    [top]
      type = LBMNonEquilibriumExtrapolation
      buffer = f
      prescribe_type = velocity
      ux = 0.0
      uy = -0.005
      uz = 0.0
      order = first
      boundary = top
    []
    [bottom]
      type = LBMConvectiveOutflow
      buffer = f
      f_old = f
      convection_velocity = Uc
      boundary = bottom
    []
    [left]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = left
    []
    [right]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = right
    []
    [front]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = front
    []
    [back]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = back
    []
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = f
  f_old = fpc
[]

[Problem]
  type = LatticeBoltzmannProblem
  scalar_constant_names = 'Uc'
  scalar_constant_values = '0.005'
  substeps = 1
[]

[Postprocessors]
  [velocity_min]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MIN
  []
  [velocity_max]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MIN
  []
  [density_max]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MAX
  []
[]

[Executioner]
  type = Transient
  num_steps = 5
[]

[Outputs]
  file_base = convective_outflow_3d_bottom
  csv = true
[]

(test/tests/lbm/nee_3d_box_second.i)

[Domain]
  dim = 3
  nx = 10
  ny = 10
  nz = 10
  mesh_mode = DUMMY
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y Z'
[]

[Stencil]
  [d3q19]
    type = LBMD3Q19
  []
[]

[TensorBuffers]
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [feq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [fpc]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [velocity]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  [density]
    type = LBMTensorBuffer
    buffer_type = ms
  []
[]

[TensorComputes]
  [Initialize]
    [initial_density]
      type = LBMConstantTensor
      buffer = density
      constants = 1.0
    []
    [initial_velocity]
      type = LBMConstantTensor
      buffer = velocity
      constants = '0.0 0.0 0.0'
    []
    [initial_feq]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [initial_f]
      type = LBMEquilibrium
      buffer = f
      bulk = density
      velocity = velocity
    []
    [initial_fpc]
      type = LBMEquilibrium
      buffer = fpc
      bulk = density
      velocity = velocity
    []
  []
  [Solve]
    [equilibrium]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [collision]
      type = LBMBGKCollision
      buffer = fpc
      f = f
      feq = feq
      tau0 = 1.0
    []
    [density]
      type = LBMComputeDensity
      buffer = density
      f = f
    []
    [velocity]
      type = LBMComputeVelocity
      buffer = velocity
      f = f
      rho = density
    []
  []
  [Boundary]
    [left]
      type = LBMNonEquilibriumExtrapolation
      buffer = f
      prescribe_type = velocity
      ux = 0.005
      uy = 0.0
      uz = 0.0
      order = second
      boundary = left
    []
    [right]
      type = LBMNonEquilibriumExtrapolation
      buffer = f
      prescribe_type = density
      prescribed_rho = 1.0
      order = second
      boundary = right
    []
    [top]
      type = LBMNonEquilibriumExtrapolation
      buffer = f
      prescribe_type = velocity
      ux = 0.0
      uy = 0.0
      uz = 0.0
      order = second
      boundary = top
    []
    [bottom]
      type = LBMNonEquilibriumExtrapolation
      buffer = f
      prescribe_type = velocity
      ux = 0.0
      uy = 0.0
      uz = 0.0
      order = second
      boundary = bottom
    []
    [front]
      type = LBMNonEquilibriumExtrapolation
      buffer = f
      prescribe_type = density
      prescribed_rho = 1.0
      order = second
      boundary = front
    []
    [back]
      type = LBMNonEquilibriumExtrapolation
      buffer = f
      prescribe_type = density
      prescribed_rho = 1.0
      order = second
      boundary = back
    []
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = f
  f_old = fpc
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 1
[]

[Postprocessors]
  [velocity_min]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MIN
  []
  [velocity_max]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MIN
  []
  [density_max]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MAX
  []
[]

[Executioner]
  type = Transient
  num_steps = 5
[]

[Outputs]
  file_base = nee_3d_box_second
  csv = true
[]

(test/tests/lbm/dirichlet_wall.i)

[Domain]
  dim = 2
  nx = 11
  ny = 11
  mesh_mode=DUMMY
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y'
[]

[Stencil]
  [d2q9]
    type = LBMD2Q9
  []
[]

[TensorBuffers]
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [feq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [fpc]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [T]
    type=LBMTensorBuffer
    buffer_type = ms
  []
  [velocity]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  [binary_media]
    type = LBMTensorBuffer
    file = binary_media.h5
    is_integer = true
    buffer_type = ms
  []
[]

[TensorComputes]
  [Initialize]
    [init_T]
      type = LBMConstantTensor
      buffer = T
      constants = 1.0
    []
    [equilibrium]
      type = LBMEquilibrium
      buffer = feq
      bulk = T
      velocity = velocity
    []
    [non_equilibrium]
      type = LBMEquilibrium
      buffer = f
      bulk = T
      velocity = velocity
    []
    [post_collision_equilibrium]
      type = LBMEquilibrium
      buffer = fpc
      bulk = T
      velocity = velocity
    []
  []
  [Solve]
    [T]
      type = LBMComputeDensity
      buffer = T
      f = f
    []
  []
  [Boundary]
    [wall]
      type = LBMDirichletBC
      buffer = f
      f_old = fpc
      feq = feq
      velocity = velocity
      rho = T
      value = 1.1
      boundary = wall
    []
  []
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 1
  binary_media = binary_media
[]

[Postprocessors]
  [velocity_min]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MIN
  []
  [velocity_max]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = T
    value_type = MIN
  []
  [densty_max]
    type = TensorExtremeValuePostprocessor
    buffer = T
    value_type = MAX
  []
[]

[Executioner]
  type = Transient
  num_steps = 5
[]

[Outputs]
  file_base = dirichlet_wall
  csv = true
[]

(test/tests/lbm/mixed_bcs_d3q19.i)

[Domain]
  dim = 3
  nx = 10
  ny = 10
  nz = 10
  mesh_mode = DUMMY
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y Z'
[]

[Stencil]
  [d3q19]
    type = LBMD3Q19
  []
[]

[TensorBuffers]
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [f_bounce_back]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [velocity]
    type=LBMTensorBuffer
    buffer_type = mv
  []
  [density]
    type=LBMTensorBuffer
    buffer_type = ms
  []
[]

[TensorComputes]
  [Initialize]
    [initial_density]
      type = LBMConstantTensor
      buffer = density
      constants = 1.0
    []
    [initial_velocity]
      type = LBMConstantTensor
      buffer = velocity
      constants = '0.0001 0.0005 0.0'
    []
    [initial_f]
      type = LBMEquilibrium
      buffer = f
      bulk = density
      velocity = velocity
    []
    [initial_f_bb]
      type = LBMEquilibrium
      buffer = f_bounce_back
      bulk = density
      velocity = velocity
    []
  []
  [Solve]
    [density]
      type = LBMComputeDensity
      buffer = density
      f = f
    []
    [velocity]
      type = LBMComputeVelocity
      buffer = velocity
      f = f
      rho = density
      add_body_force = true
      body_force_x = 0.0001
    []
  []
  [Boundary]
    [left]
      type = LBMFixedZerothOrderBC
      buffer = f
      f = f
      value = 1.11
      boundary = left
    []
    [right]
      type = LBMFixedFirstOrderBC
      buffer = f
      f = f
      value = 0.0001
      boundary = right
    []
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = f
  f_old = f_bounce_back
[]

[Postprocessors]
  [velocity_min]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MIN
  []
  [velocity_max]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MIN
  []
  [densty_max]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MAX
  []
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 1
[]

[Executioner]
  type = Transient
  num_steps = 5
[]

[Outputs]
  file_base = mixed_bcs_d3q19
  csv = true
[]

(test/tests/lbm/horizontal_mixed_bcs_d2q9_reverse.i)

[Domain]
  dim = 2
  nx = 10
  ny = 10
  mesh_mode = DUMMY
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y'
[]

[Stencil]
  [d2q9]
    type = LBMD2Q9
  []
[]

[TensorBuffers]
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [f_bounce_back]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [velocity]
    type=LBMTensorBuffer
    buffer_type = mv
  []
  [density]
    type=LBMTensorBuffer
    buffer_type = ms
  []
[]

[TensorComputes]
  [Initialize]
    [initial_density]
      type = LBMConstantTensor
      buffer = density
      constants = 1.0
    []
    [initial_velocity]
      type = LBMConstantTensor
      buffer = velocity
      constants = '0.0001 0.0005'
    []
    [initial_f]
      type = LBMEquilibrium
      buffer = f
      bulk = density
      velocity = velocity
    []
    [initial_f_bb]
      type = LBMEquilibrium
      buffer = f_bounce_back
      bulk = density
      velocity = velocity
    []
  []
  [Solve]
    [density]
      type = LBMComputeDensity
      buffer = density
      f = f
    []
    [velocity]
      type = LBMComputeVelocity
      buffer = velocity
      f = f
      rho = density
      add_body_force = true
      body_force_x = 0.0001
    []
  []
  [Boundary]
    [right]
      type = LBMFixedZerothOrderBC
      buffer = f
      f = f
      value = 1.1
      boundary = right
    []
    [left]
      type = LBMMicroscopicZeroGradientBC
      buffer = f
      boundary = left
    []
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = f
  f_old = f_bounce_back
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 1
[]

[Postprocessors]
  [velocity_min]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MIN
  []
  [velocity_max]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MIN
  []
  [densty_max]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MAX
  []
[]

[Executioner]
  type = Transient
  num_steps = 5
[]

[Outputs]
  file_base = horizontal_mixed_bcs_d2q9_reverse
  csv = true
[]

(examples/lbm/Rayleigh-Benard/solve.i)

#### Compute ####
[Solve]

  # For temperature
  [Temperature]
    type = LBMComputeDensity
    buffer = T
    f = g
  []

  # For fluid
  [Fluid_density]
    type = LBMComputeDensity
    buffer = density
    f = f
  []

  [Fluid_velocity]
    type = LBMComputeVelocity
    buffer = velocity
    f = f
    rho = density
    forces = F
    enable_forces = true
  []

  # For temperature
  [Equilibrium_temperature]
    type = LBMEquilibrium
    buffer = geq
    bulk = T
    velocity = velocity
  []

  [Collision_temperature]
    type = LBMBGKCollision
    buffer = gpc
    f = g
    feq = geq
    tau0 = tau_T
  []

  # For fluid
  [Compute_forces]
    type = LBMComputeForces
    buffer = F
    rho0 = 'rho0'
    temperature = T
    T0 = T_0
    enable_buoyancy = true
    gravity = g
  []

  [Equilibrium_fluid]
    type = LBMEquilibrium
    buffer = feq
    bulk = density
    velocity = velocity
  []

  [Collision_fluid]
    type = LBMBGKCollision
    buffer = fpc
    f = f
    feq = feq
    tau0 = tau_f
  []

  [Apply_forces]
    type = LBMApplyForces
    buffer = fpc
    velocity = velocity
    rho = density
    forces = F
    tau0 = tau_f
  []

  [Residual]
    type = LBMComputeResidual
    speed = T
    # TODO this buffer is redundant, but avoids 'missing parameter' error
    buffer = T
  []
[]

(test/tests/lbm/horizontal_mixed_bcs_d2q9.i)

[Domain]
  dim = 2
  nx = 10
  ny = 10
  mesh_mode = DUMMY
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y'
[]

[Stencil]
  [d2q9]
    type = LBMD2Q9
  []
[]

[TensorBuffers]
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [f_bounce_back]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [velocity]
    type=LBMTensorBuffer
    buffer_type = mv
  []
  [density]
    type=LBMTensorBuffer
    buffer_type = ms
  []
[]

[TensorComputes]
  [Initialize]
    [initial_density]
      type = LBMConstantTensor
      buffer = density
      constants = 1.0
    []
    [initial_velocity]
      type = LBMConstantTensor
      buffer = velocity
      constants = '0.0001 0.0005'
    []
    [initial_f]
      type = LBMEquilibrium
      buffer = f
      bulk = density
      velocity = velocity
    []
    [initial_f_bb]
      type = LBMEquilibrium
      buffer = f_bounce_back
      bulk = density
      velocity = velocity
    []
  []
  [Solve]
    [density]
      type = LBMComputeDensity
      buffer = density
      f = f
    []
    [velocity]
      type = LBMComputeVelocity
      buffer = velocity
      f = f
      rho = density
      add_body_force = true
      body_force_x = 0.0001
    []
  []
  [Boundary]
    [left]
      type = LBMFixedZerothOrderBC
      buffer = f
      f = f
      value = 1.1
      boundary = left
    []
    [right]
      type = LBMFixedFirstOrderBC
      buffer = f
      f = f
      value = 0.0001
      boundary = right
    []
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = f
  f_old = f_bounce_back
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 1
[]

[Postprocessors]
  [velocity_min]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MIN
  []
  [velocity_max]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MIN
  []
  [densty_max]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MAX
  []
[]

[Executioner]
  type = Transient
  num_steps = 5
[]

[Outputs]
  file_base = horizontal_mixed_bcs_d2q9
  csv = true
[]

(test/tests/lbm/nee_2d_channel.i)

[Domain]
  dim = 2
  nx = 10
  ny = 10
  mesh_mode = DUMMY
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y'
[]

[Stencil]
  [d2q9]
    type = LBMD2Q9
  []
[]

[TensorBuffers]
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [feq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [fpc]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [velocity]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  [density]
    type = LBMTensorBuffer
    buffer_type = ms
  []
[]

[TensorComputes]
  [Initialize]
    [initial_density]
      type = LBMConstantTensor
      buffer = density
      constants = 1.0
    []
    [initial_velocity]
      type = LBMConstantTensor
      buffer = velocity
      constants = '0.0 0.0'
    []
    [initial_feq]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [initial_f]
      type = LBMEquilibrium
      buffer = f
      bulk = density
      velocity = velocity
    []
    [initial_fpc]
      type = LBMEquilibrium
      buffer = fpc
      bulk = density
      velocity = velocity
    []
  []
  [Solve]
    [equilibrium]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [collision]
      type = LBMBGKCollision
      buffer = fpc
      f = f
      feq = feq
      tau0 = 1.0
    []
    [density]
      type = LBMComputeDensity
      buffer = density
      f = f
    []
    [velocity]
      type = LBMComputeVelocity
      buffer = velocity
      f = f
      rho = density
    []
  []
  [Boundary]
    [left]
      type = LBMNonEquilibriumExtrapolation
      buffer = f
      prescribe_type = velocity
      ux = 0.005
      uy = 0.0
      order = first
      boundary = left
    []
    [right]
      type = LBMNonEquilibriumExtrapolation
      buffer = f
      prescribe_type = density
      prescribed_rho = 1.0
      order = first
      boundary = right
    []
    [top]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = top
    []
    [bottom]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = bottom
    []
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = f
  f_old = fpc
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 1
[]

[Postprocessors]
  [velocity_min]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MIN
  []
  [velocity_max]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MIN
  []
  [density_max]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MAX
  []
[]

[Executioner]
  type = Transient
  num_steps = 5
[]

[Outputs]
  file_base = nee_2d_channel
  csv = true
[]

(test/tests/lbm/convective_outflow_2d_bottom.i)

[Domain]
  dim = 2
  nx = 10
  ny = 10
  mesh_mode = DUMMY
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y'
[]

[Stencil]
  [d2q9]
    type = LBMD2Q9
  []
[]

[TensorBuffers]
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [feq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [fpc]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [velocity]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  [density]
    type = LBMTensorBuffer
    buffer_type = ms
  []
[]

[TensorComputes]
  [Initialize]
    [initial_density]
      type = LBMConstantTensor
      buffer = density
      constants = 1.0
    []
    [initial_velocity]
      type = LBMConstantTensor
      buffer = velocity
      constants = '0.0 0.0'
    []
    [initial_feq]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [initial_f]
      type = LBMEquilibrium
      buffer = f
      bulk = density
      velocity = velocity
    []
    [initial_fpc]
      type = LBMEquilibrium
      buffer = fpc
      bulk = density
      velocity = velocity
    []
  []
  [Solve]
    [equilibrium]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [collision]
      type = LBMBGKCollision
      buffer = fpc
      f = f
      feq = feq
      tau0 = 1.0
    []
    [density]
      type = LBMComputeDensity
      buffer = density
      f = f
    []
    [velocity]
      type = LBMComputeVelocity
      buffer = velocity
      f = f
      rho = density
    []
  []
  [Boundary]
    [top]
      type = LBMNonEquilibriumExtrapolation
      buffer = f
      prescribe_type = velocity
      ux = 0.0
      uy = -0.005
      order = first
      boundary = top
    []
    [bottom]
      type = LBMConvectiveOutflow
      buffer = f
      f_old = f
      convection_velocity = auto
      boundary = bottom
    []
    [left]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = left
    []
    [right]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = right
    []
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = f
  f_old = fpc
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 1
[]

[Postprocessors]
  [velocity_min]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MIN
  []
  [velocity_max]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MIN
  []
  [density_max]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MAX
  []
[]

[Executioner]
  type = Transient
  num_steps = 5
[]

[Outputs]
  file_base = convective_outflow_2d_bottom
  csv = true
[]

(test/tests/lbm/convective_outflow_2d_top.i)

[Domain]
  dim = 2
  nx = 10
  ny = 10
  mesh_mode = DUMMY
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y'
[]

[Stencil]
  [d2q9]
    type = LBMD2Q9
  []
[]

[TensorBuffers]
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [feq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [fpc]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [velocity]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  [density]
    type = LBMTensorBuffer
    buffer_type = ms
  []
[]

[TensorComputes]
  [Initialize]
    [initial_density]
      type = LBMConstantTensor
      buffer = density
      constants = 1.0
    []
    [initial_velocity]
      type = LBMConstantTensor
      buffer = velocity
      constants = '0.0 0.0'
    []
    [initial_feq]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [initial_f]
      type = LBMEquilibrium
      buffer = f
      bulk = density
      velocity = velocity
    []
    [initial_fpc]
      type = LBMEquilibrium
      buffer = fpc
      bulk = density
      velocity = velocity
    []
  []
  [Solve]
    [equilibrium]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [collision]
      type = LBMBGKCollision
      buffer = fpc
      f = f
      feq = feq
      tau0 = 1.0
    []
    [density]
      type = LBMComputeDensity
      buffer = density
      f = f
    []
    [velocity]
      type = LBMComputeVelocity
      buffer = velocity
      f = f
      rho = density
    []
  []
  [Boundary]
    [bottom]
      type = LBMNonEquilibriumExtrapolation
      buffer = f
      prescribe_type = velocity
      ux = 0.0
      uy = 0.005
      order = first
      boundary = bottom
    []
    [top]
      type = LBMConvectiveOutflow
      buffer = f
      f_old = f
      convection_velocity = auto
      boundary = top
    []
    [left]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = left
    []
    [right]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = right
    []
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = f
  f_old = fpc
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 1
[]

[Postprocessors]
  [velocity_min]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MIN
  []
  [velocity_max]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MIN
  []
  [density_max]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MAX
  []
[]

[Executioner]
  type = Transient
  num_steps = 5
[]

[Outputs]
  file_base = convective_outflow_2d_top
  csv = true
[]

(test/tests/lbm/dirichlet_regional.i)

[Domain]
  dim = 2
  nx = 11
  ny = 11
  mesh_mode=DUMMY
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y'
[]

[Stencil]
  [d2q9]
    type = LBMD2Q9
  []
[]

[TensorBuffers]
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [feq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [fpc]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [T]
    type=LBMTensorBuffer
    buffer_type = ms
  []
  [velocity]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  [binary_media]
    type = LBMTensorBuffer
    file = binary_regional.h5
    is_integer = true
    buffer_type = ms
  []
[]

[TensorComputes]
  [Initialize]
    [init_T]
      type = LBMConstantTensor
      buffer = T
      constants = 1.0
    []
    [equilibrium]
      type = LBMEquilibrium
      buffer = feq
      bulk = T
      velocity = velocity
    []
    [non_equilibrium]
      type = LBMEquilibrium
      buffer = f
      bulk = T
      velocity = velocity
    []
    [post_collision_equilibrium]
      type = LBMEquilibrium
      buffer = fpc
      bulk = T
      velocity = velocity
    []
  []
  [Solve]
    [T]
      type = LBMComputeDensity
      buffer = T
      f = f
    []
  []
  [Boundary]
    [regional]
      type = LBMDirichletBC
      buffer = f
      f_old = fpc
      feq = feq
      velocity = velocity
      rho = T
      value = 1.1
      region_id = 2
      boundary = regional
    []
  []
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 1
  binary_media = binary_media
[]

[Postprocessors]
  [velocity_min]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MIN
  []
  [velocity_max]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = T
    value_type = MIN
  []
  [densty_max]
    type = TensorExtremeValuePostprocessor
    buffer = T
    value_type = MAX
  []
[]

[Executioner]
  type = Transient
  num_steps = 5
[]

[Outputs]
  file_base = dirichlet_regional
  csv = true
[]

(examples/lbm/Phase-field/droplet_impact.i)

#
# Droplet Impact on a Thin Liquid Film
# PHYSICAL REVIEW E 97, 033309 (2018) - Section III.D
#

# Domain
Nx = 1500
Ny = 500

# Fluid properties (Re=500, We=8000, 1000:1)
rho_l = 1000.0
rho_g = 1.0
# mu_l = 20.0
# mu_g = 0.2
nu_l = 0.02
nu_g = 0.2
sigma = 0.0625

# Phase field parameters
tau_h = 0.8
D = 5

[Domain]
  dim = 2
  nx = '${Nx}'
  ny = '${Ny}'
  xmax = '${Nx}'
  ymax = '${Ny}'
  device_names = 'cuda'
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y'
[]

[Stencil]
  [d2q9]
    type = LBMD2Q9
  []
[]

[TensorBuffers]
  [phi]
    type = LBMTensorBuffer
    buffer_type = ms
    file = phi.h5
  []
  [ux]
    type = LBMTensorBuffer
    buffer_type = ms
    file = ux.h5
  []
  [uy]
    type = LBMTensorBuffer
    buffer_type = ms
    file = uy.h5
  []
  [velocity]
    type = LBMTensorBuffer
    buffer_type = mv
  []

  [grad_phi]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  [laplacian_phi]
    type = LBMTensorBuffer
    buffer_type = ms
  []
  [mu]
    type = LBMTensorBuffer
    buffer_type = ms
  []
  [forces]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  [speed]
    type = LBMTensorBuffer
    buffer_type = ms
  []
  [pressure]
    type = LBMTensorBuffer
    buffer_type = ms
  []
  [rho]
    type = LBMTensorBuffer
    buffer_type = ms
  []
  [h]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [h_post_collision]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [h_eq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [relaxation_tensor]
    type = LBMTensorBuffer
    buffer_type = ms
  []
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [f_post_collision]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [f_eq]
    type = LBMTensorBuffer
    buffer_type = df
  []
[]

[TensorComputes/Initialize]
  [u_stack]
    type = LBMStackTensors
    buffer = velocity
    inputs = 'ux uy'
  []

  [grad_phi_init]
    type = LBMIsotropicGradient
    buffer = grad_phi
    scalar_field = phi
  []
  [rho_init]
    type = ParsedCompute
    buffer = rho
    expression = 'phi * (rho_l - rho_g) + rho_g'
    constant_names = 'rho_l rho_g'
    constant_expressions = '${rho_l} ${rho_g}'
    inputs = phi
  []
  [pressure_init]
    type = LBMConstantTensor
    buffer = pressure
    constants = 0.3
  []

  # Equilibrium
  [h_eq_init]
    type = LBMPhaseEquilibrium
    buffer = h_eq
    phi = phi
    velocity = velocity
  []
  [h_init]
    type = ParsedCompute
    buffer = h
    expression = 'h_eq'
    inputs = h_eq
  []
  [h_post_collision_init]
    type = ParsedCompute
    buffer = h_post_collision
    expression = 'h_eq'
    inputs = h_eq
  []
  [f_eq_init]
    type = LBMPressureCorrectedEquilibrium
    buffer = f_eq
    rho = rho
    velocity = velocity
    pressure = pressure
  []
  [f_init]
    type = ParsedCompute
    buffer = f
    expression = 'f_eq'
    inputs = f_eq
  []
  [f_post_collision_init]
    type = ParsedCompute
    buffer = f_post_collision
    expression = 'f_eq'
    inputs = f_eq
  []
[]

[TensorComputes/Solve]
  [compute_phi]
    type = LBMComputeDensity
    buffer = phi
    f = h
  []
  [grad_phi]
    type = LBMIsotropicGradient
    buffer = grad_phi
    scalar_field = phi
  []
  [laplacian_phi]
    type = LBMIsotropicLaplacian
    buffer = laplacian_phi
    scalar_field = phi
  []
  [potential]
    type = LBMComputeChemicalPotential
    buffer = mu
    phi = phi
    laplacian_phi = laplacian_phi
    thickness = D
    sigma = sigma
  []
  [compute_forces]
    type = LBMComputeSurfaceForces
    buffer = forces
    chemical_potential = mu
    grad_phi = grad_phi
  []
  [density]
    type = ParsedCompute
    buffer = rho
    expression = 'phi * (rho_l - rho_g) + rho_g'
    constant_names = 'rho_l rho_g'
    constant_expressions = '${rho_l} ${rho_g}'
    inputs = phi
  []
  [velocity]
    type = LBMComputeVelocity
    buffer = velocity
    f = f
    rho = rho
    enable_forces = true
    forces = forces
  []
  [h_eq]
    type = LBMPhaseEquilibrium
    buffer = h_eq
    phi = phi
    velocity = velocity
  []
  [phase_collision]
    type = LBMBGKCollision
    buffer = h_post_collision
    f = h
    feq = h_eq
    tau0 = tau_h
  []
  [apply_forces_phase]
    type = LBMAllenCahnSource
    buffer = h_post_collision
    phi = phi
    velocity = velocity
    grad_phi = grad_phi
    tau = tau_h
    thickness = D
  []
  [relaxation_tensor]
    type = ParsedCompute
    buffer = relaxation_tensor
    expression = '(phi * (nu_l - nu_g) + nu_g) / cs2 + 0.5'
    constant_names = 'nu_l nu_g cs2'
    constant_expressions = '${nu_l} ${nu_g} 0.333333333333'
    inputs = 'phi'
  []
  [pressure]
    type = LBMPhaseFieldPressure
    buffer = pressure
    f = f
    velocity = velocity
    grad_phi = grad_phi
    rho = rho
    rho_l = '${rho_l}'
    rho_g = '${rho_g}'
  []
  [f_eq]
    type = LBMPressureCorrectedEquilibrium
    buffer = f_eq
    rho = rho
    velocity = velocity
    pressure = pressure
  []
  [collision]
    type = LBMBGKCollision
    buffer = f_post_collision
    f = f
    feq = f_eq
    tau0 = 1.0
    is_dynamic_relaxation = true
    tau_tensor = relaxation_tensor
  []
  [apply_forces_hydro]
    type = LBMForceDistribution
    buffer = f_post_collision
    grad_phi = grad_phi
    velocity = velocity
    forces = forces
    tau_tensor = relaxation_tensor
    tau = 1.0
    rho_l = '${rho_l}'
    rho_g = '${rho_g}'
    is_dynamic_relaxation = true
  []
  [residual]
    type = LBMComputeResidual
    buffer = phi
    speed = phi
  []
[]

[TensorComputes/Boundary]
  # Bounce back on both top and bottom
  [top_fluid]
    type = LBMBounceBack
    buffer = f
    f_old = f_post_collision
    boundary = top
  []
  [bottom_fluid]
    type = LBMBounceBack
    buffer = f
    f_old = f_post_collision
    boundary = bottom
  []
  [top_phase]
    type = LBMBounceBack
    buffer = h
    f_old = h_post_collision
    boundary = top
  []
  [bottom_phase]
    type = LBMBounceBack
    buffer = h
    f_old = h_post_collision
    boundary = bottom
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = 'h f'
  f_old = 'h_post_collision f_post_collision'
  root_compute = residual
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 200
  print_debug_output = true
  scalar_constant_names = 'tau_h D sigma'
  scalar_constant_values = '${tau_h} ${D} ${sigma}'
[]

[Executioner]
  type = Transient
  num_steps = 50
[]

[TensorOutputs]
  [xdmf]
    type = XDMFTensorOutput
    buffer = 'phi velocity rho'
    output_mode = 'Cell Cell Cell'
    enable_hdf5 = true
    transpose = false
  []
[]

(examples/lbm/Phase-field/spinodal_decomposition.i)

#
# Spinodal decomposition
# PHYSICAL REVIEW E 97, 033309 (2018) - Section III.C
#

# Domain
Nx = 200
Ny = 200

# Fluid properties
rho_l = 1000.0
rho_g = 1.0
nu_l = 0.1
nu_g = 1.0
sigma = 0.2

# Phase field parameters
tau_h = 0.67
D = 4

[Domain]
  dim = 2
  nx = '${Nx}'
  ny = '${Ny}'
  xmax = '${Nx}'
  ymax = '${Ny}'
  device_names='cpu'
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y'
[]

[Stencil]
  [d2q9]
    type = LBMD2Q9
  []
[]

[TensorBuffers]
  # Macroscopic phase field variables
  [phi]
    type = LBMTensorBuffer
    buffer_type = ms
  []
  [grad_phi]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  [laplacian_phi]
    type = LBMTensorBuffer
    buffer_type = ms
  []
  [mu]
    type = LBMTensorBuffer
    buffer_type = ms
  []
  [forces]
    type = LBMTensorBuffer
    buffer_type = mv
  []

  # Macroscopic hydrodynamic variables
  [velocity]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  [pressure]
    type = LBMTensorBuffer
    buffer_type = ms
  []
  [rho]
    type = LBMTensorBuffer
    buffer_type = ms
  []

  # LBM phase field variabels
  [h]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [h_post_collision]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [h_eq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [relaxation_tensor]
    type = LBMTensorBuffer
    buffer_type = ms
  []

  # LBM hydrodynamic variables
  [fdummy]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [f_post_collision]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [f_eq]
    type = LBMTensorBuffer
    buffer_type = df
  []
[]

[TensorComputes/Initialize]
  [phi_init]
    type = RandomTensor
    buffer = phi
    min = 0.3233
    max = 0.3433
  []
  [grad_phi_init]
    type = LBMIsotropicGradient
    buffer = grad_phi
    scalar_field = phi
  []
  [rho_init]
    type = ParsedCompute
    buffer = rho
    extra_symbols = true
    expression = 'phi*(rho_l - rho_g) + rho_g'
    constant_names = 'rho_l rho_g'
    constant_expressions = '${rho_l} ${rho_g}'
    inputs = phi
  []
  [pressure_init]
    type = LBMConstantTensor
    buffer = pressure
    constants = 0.3
  []
  # Phase field equilibrium distribution initialization
  [h_eq_init]
    type = LBMPhaseEquilibrium
    buffer = h_eq
    phi = phi
    velocity = velocity
  []
  [h_post_collision_init]
    type = LBMPhaseEquilibrium
    buffer = h_post_collision
    phi = phi
    velocity = velocity
  []
  [h_init]
    type = LBMPhaseEquilibrium
    buffer = h
    phi = phi
    velocity = velocity
  []
  # Hydrodynamic equilibrium distribution initialization
  [f_eq_init]
    type = LBMPressureCorrectedEquilibrium
    buffer = f_eq
    rho = rho
    velocity = velocity
    pressure = pressure
  []
  [f_post_collision_init]
    type = ParsedCompute
    buffer = f_post_collision
    expression = 'f_eq'
    inputs = f_eq
  []
  [f_init]
    type = ParsedCompute
    buffer = f
    expression = 'f_eq'
    inputs = f_eq
  []
[]

[TensorComputes/Solve]
  # Phase Field
  [compute_phi]
    type = LBMComputeDensity
    buffer = phi
    f = h
  []
  [grad_phi]
    type = LBMIsotropicGradient
    buffer = grad_phi
    scalar_field = phi
  []
  [laplacian_phi]
    type = LBMIsotropicLaplacian
    buffer = laplacian_phi
    scalar_field = phi
  []
  [potential]
    type = LBMComputeChemicalPotential
    buffer = mu
    phi = phi
    laplacian_phi = laplacian_phi
    thickness = D
    sigma = sigma
  []
  [forces]
    type = LBMComputeSurfaceForces
    buffer = forces
    chemical_potential = mu
    grad_phi = grad_phi
  []
  # Hydrodynamics
  [density]
    type = ParsedCompute
    buffer = rho
    extra_symbols = true
    expression = 'phi*(rho_l - rho_g) + rho_g'
    constant_names = 'rho_l rho_g'
    constant_expressions = '${rho_l} ${rho_g}'
    inputs = phi
  []
  [velocity]
    type = LBMComputeVelocity
    buffer = velocity
    f = f
    rho = rho
    enable_forces = true
    forces = forces
  []
  # Phase-field
  [h_eq]
    type = LBMPhaseEquilibrium
    buffer = h_eq
    phi = phi
    velocity = velocity
  []
  [phase_collision]
    type = LBMBGKCollision
    buffer = h_post_collision
    f = h
    feq = h_eq
    tau0 = tau_h
  []
  [apply_forces_phase]
    type = LBMAllenCahnSource
    buffer = h_post_collision
    phi = phi
    velocity = velocity
    grad_phi = grad_phi
    tau = tau_h
    thickness = D
  []
  # Hydrodynamics
  [relaxation_tensor]
    type = ParsedCompute
    buffer = relaxation_tensor
    extra_symbols = true
    expression = '(phi*(nu_l - nu_g) + nu_g)/cs2+0.5'
    constant_names = 'nu_l nu_g cs2'
    constant_expressions = '${nu_l} ${nu_g} 0.3333'
    inputs = phi
  []
  [pressure]
    type = LBMPhaseFieldPressure
    buffer = pressure
    f = f
    velocity = velocity
    grad_phi = grad_phi
    rho = rho
    rho_l = '${rho_l}'
    rho_g = '${rho_g}'
  []
  [f_eq]
    type = LBMPressureCorrectedEquilibrium
    buffer = f_eq
    rho = rho
    velocity = velocity
    pressure = pressure
  []
  [collision]
    type = LBMBGKCollision
    buffer = f_post_collision
    f = f
    feq = f_eq
    tau0 = 1.0
    is_dynamic_relaxation = true
    tau_tensor = relaxation_tensor
   []
  [apply_forces_hydro]
    type = LBMForceDistribution
    buffer = f_post_collision
    grad_phi = grad_phi
    velocity = velocity
    forces = forces
    tau_tensor = relaxation_tensor
    tau = 1.0
    rho_l = '${rho_l}'
    rho_g = '${rho_g}'
    is_dynamic_relaxation = true
  []
  [residual]
    type = LBMComputeResidual
    buffer = phi
    speed = phi
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = 'h f'
  f_old = 'h_post_collision f_post_collision'
  root_compute = residual
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 500
  print_debug_output = true
  scalar_constant_names = 'tau_h D sigma'
  scalar_constant_values = '${tau_h} ${D} ${sigma}'
  log_interval = 100
[]

[Executioner]
  type = Transient
  num_steps = 100
[]

[TensorOutputs]
  [xdmf]
    type = XDMFTensorOutput
    buffer = 'phi rho velocity'
    output_mode = 'Cell Cell Cell'
    enable_hdf5 = true
    # transpose = false
  []
[]

(test/tests/lbm/channel2D.i)

[Domain]
  dim = 2
  nx = 10
  ny = 10
  mesh_mode = DUMMY
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y'
[]

[Stencil]
  [d2q9]
    type = LBMD2Q9
  []
[]

[TensorBuffers]
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [feq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [fpc]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [velocity]
    type=LBMTensorBuffer
    buffer_type = mv
  []
  [density]
    type=LBMTensorBuffer
    buffer_type = ms
  []
  [speed]
    type=LBMTensorBuffer
    buffer_type = ms
  []
[]

[TensorComputes]
  [Initialize]
    [initial_density]
      type = LBMConstantTensor
      buffer = density
      constants = 1.0
    []
    [initial_velocity]
      type = LBMConstantTensor
      buffer = velocity
      constants = '0.0 0.0'
    []
    [initial_equilibrium]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [initial_distribution]
      type = LBMEquilibrium
      buffer = f
      bulk = density
      velocity = velocity
    []
    [initial_distribution_pc]
      type = LBMEquilibrium
      buffer = fpc
      bulk = density
      velocity = velocity
    []
  []
  [Solve]
    [equilibrium]
      type=LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [collision]
      type=LBMBGKCollision
      buffer = fpc
      f = f
      feq = feq
      tau0 = 1.0
    []
    [density]
      type = LBMComputeDensity
      buffer = density
      f = f
    []
    [velocity]
      type = LBMComputeVelocity
      buffer = velocity
      f = f
      rho = density
      add_body_force = true
      body_force_x = 0.0001
    []
    [speed]
      type = LBMComputeVelocityMagnitude
      buffer = speed
      velocity = velocity
    []
    [residual]
      type = LBMComputeResidual
      buffer = speed
      speed = speed
    []
  []
  [Boundary]
    [top]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = top
    []
    [bottom]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = bottom
    []
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = f
  f_old = fpc
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 10
[]

[Postprocessors]
  [velocity_min]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MIN
  []
  [velocity_max]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MAX
  []
  [speed_min]
    type = TensorExtremeValuePostprocessor
    buffer = speed
    value_type = MIN
  []
  [speed_max]
    type = TensorExtremeValuePostprocessor
    buffer = speed
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MIN
  []
  [densty_max]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MAX
  []
[]

[Executioner]
  type = Transient
  num_steps = 11
[]

[Outputs]
  file_base = channel2D
  csv = true
[]

(examples/lbm/Karman-vortex/cylinder.i)

[Domain]
  dim = 2
  nx = 800
  ny = 200
  xmax = 4
  ymax = 1
  device_names='cpu'
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y'
[]

[Stencil]
  [d2q9]
    type = LBMD2Q9
  []
[]

[TensorBuffers]
  [rho]
    type=LBMTensorBuffer
    buffer_type = ms
  []
  [u]
    type=LBMTensorBuffer
    buffer_type = mv
  []
  [speed]
    type=LBMTensorBuffer
    buffer_type = ms
  []
  [f]
    type=LBMTensorBuffer
    buffer_type = df
  []
  [feq]
    type=LBMTensorBuffer
    buffer_type = df
  []
  [f_post_collision]
    type=LBMTensorBuffer
    buffer_type = df
  []
  [binary_media]
    type = LBMTensorBuffer
    file = binary_media.h5
    is_integer = true
    buffer_type = ms
  []
[]

[TensorComputes]
  [Initialize]
    [rho]
      type=LBMConstantTensor
      buffer=rho
      constants = rho0
    []
    [u]
      type=LBMConstantTensor
      buffer=u
      constants = 'Ux Uy'
    []
    [speed]
      type=LBMComputeVelocityMagnitude
      buffer=speed
      velocity=u
    []
    [feq]
      type=LBMEquilibrium
      buffer=feq
      bulk=rho
      velocity=u
    []
    [f]
      type=LBMEquilibrium
      buffer=f
      bulk=rho
      velocity=u
    []
    [f_post_coll]
      type=LBMEquilibrium
      buffer=f_post_collision
      bulk=rho
      velocity=u
    []
  []
  [Solve]
    [Density]
      type = LBMComputeDensity
      buffer=rho
      f = f
    []
    [Velocity]
      type = LBMComputeVelocity
      buffer=u
      f = f
      rho = rho
    []
    [Equilibrium]
      type = LBMEquilibrium
      buffer=feq
      bulk=rho
      velocity=u
    []
    [Collision]
      type = LBMSmagorinskyCollision
      buffer = f_post_collision
      f = f
      feq = feq
      tau0 = tau
      projection=true
    []
    [Speed]
      type = LBMComputeVelocityMagnitude
      buffer = speed
      velocity = u
    []
    [Residual]
      type = LBMComputeResidual
      buffer = speed
      speed = speed
    []
  []
  [Boundary]
    [wall]
      type = LBMBounceBack
      buffer = f
      f_old = f_post_collision
      boundary = wall
    []
    [left]
      type = LBMFixedFirstOrderBC
      buffer=f
      f=f
      value='Ux'
      perturb=true
      boundary=left
    []
    [right]
      type = LBMMicroscopicZeroGradientBC
      buffer=f
      boundary=right
    []
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = f
  f_old = f_post_collision
  root_compute=residual
[]

[Postprocessors]
  [rho_avg]
    type = TensorAveragePostprocessor
    buffer = rho
    execute_on = 'TIMESTEP_END'
  []
  [speed_avg]
    type = TensorAveragePostprocessor
    buffer = speed
    execute_on = 'TIMESTEP_END'
  []
  [reynolds]
    type = ComputeReynoldsNumber
    buffer = speed
    tau = tau
    diameter = D
  []
[]

[Problem]
  type = LatticeBoltzmannProblem
  scalar_constant_names = 'rho0 Ux    Uy  tau    dx    D    Cs'
  scalar_constant_values = '1.0 0.01 0.0 0.506 0.001  40  0.15'
  substeps = 100
  print_debug_output=true
  binary_media = binary_media
[]

[Executioner]
  type = Transient
  num_steps = 10000
[]

[TensorOutputs]
  [xdmf2]
    type = XDMFTensorOutput
    buffer = 'rho u speed binary_media'
    output_mode = 'Cell Cell Cell Cell'
    enable_hdf5 = true
  []
[]

(test/tests/lbm/convective_outflow_3d_right.i)

[Domain]
  dim = 3
  nx = 8
  ny = 8
  nz = 8
  mesh_mode = DUMMY
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y Z'
[]

[Stencil]
  [d3q19]
    type = LBMD3Q19
  []
[]

[TensorBuffers]
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [feq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [fpc]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [velocity]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  [density]
    type = LBMTensorBuffer
    buffer_type = ms
  []
[]

[TensorComputes]
  [Initialize]
    [initial_density]
      type = LBMConstantTensor
      buffer = density
      constants = 1.0
    []
    [initial_velocity]
      type = LBMConstantTensor
      buffer = velocity
      constants = '0.0 0.0 0.0'
    []
    [initial_feq]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [initial_f]
      type = LBMEquilibrium
      buffer = f
      bulk = density
      velocity = velocity
    []
    [initial_fpc]
      type = LBMEquilibrium
      buffer = fpc
      bulk = density
      velocity = velocity
    []
  []
  [Solve]
    [equilibrium]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [collision]
      type = LBMBGKCollision
      buffer = fpc
      f = f
      feq = feq
      tau0 = 1.0
    []
    [density]
      type = LBMComputeDensity
      buffer = density
      f = f
    []
    [velocity]
      type = LBMComputeVelocity
      buffer = velocity
      f = f
      rho = density
    []
  []
  [Boundary]
    [left]
      type = LBMNonEquilibriumExtrapolation
      buffer = f
      prescribe_type = velocity
      ux = 0.005
      uy = 0.0
      uz = 0.0
      order = first
      boundary = left
    []
    [right]
      type = LBMConvectiveOutflow
      buffer = f
      f_old = f
      convection_velocity = Uc
      boundary = right
    []
    [top]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = top
    []
    [bottom]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = bottom
    []
    [front]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = front
    []
    [back]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = back
    []
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = f
  f_old = fpc
[]

[Problem]
  type = LatticeBoltzmannProblem
  scalar_constant_names = 'Uc'
  scalar_constant_values = '0.005'
  substeps = 1
[]

[Postprocessors]
  [velocity_min]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MIN
  []
  [velocity_max]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MIN
  []
  [density_max]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MAX
  []
[]

[Executioner]
  type = Transient
  num_steps = 5
[]

[Outputs]
  file_base = convective_outflow_3d_right
  csv = true
[]

(test/tests/lbm/mixed_bcs_d3q27.i)

[Domain]
  dim = 3
  nx = 10
  ny = 10
  nz = 10
  mesh_mode = DUMMY
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y Z'
[]

[Stencil]
  [d3q27]
    type = LBMD3Q27
  []
[]

[TensorBuffers]
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [f_bounce_back]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [velocity]
    type=LBMTensorBuffer
    buffer_type = mv
  []
  [density]
    type=LBMTensorBuffer
    buffer_type = ms
  []
[]

[TensorComputes]
  [Initialize]
    [initial_density]
      type = LBMConstantTensor
      buffer = density
      constants = 1.0
    []
    [initial_velocity]
      type = LBMConstantTensor
      buffer = velocity
      constants = '0.0001 0.0005 0.0'
    []
    [initial_f]
      type = LBMEquilibrium
      buffer = f
      bulk = density
      velocity = velocity
    []
    [initial_f_bb]
      type = LBMEquilibrium
      buffer = f_bounce_back
      bulk = density
      velocity = velocity
    []
  []
  [Solve]
    [density]
      type = LBMComputeDensity
      buffer = density
      f = f
    []
    [velocity]
      type = LBMComputeVelocity
      buffer = velocity
      f = f
      rho = density
      add_body_force = true
      body_force_x = 0.0001
    []
  []
  [Boundary]
    [left]
      type = LBMFixedZerothOrderBC
      buffer = f
      f = f
      value = 1.1
      boundary = left
    []
    [right]
      type = LBMFixedFirstOrderBC
      buffer = f
      f = f
      value = 0.0001
      boundary = right
    []
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = f
  f_old = f_bounce_back
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 1
[]

[Postprocessors]
  [velocity_min]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MIN
  []
  [velocity_max]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MIN
  []
  [densty_max]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MAX
  []
[]

[Executioner]
  type = Transient
  num_steps = 5
[]

[Outputs]
  file_base = mixed_bcs_d3q27
  csv = true
[]

(test/tests/lbm/obstacle.i)

[Domain]
  dim = 2
  nx = 11
  ny = 11
  mesh_mode=DUMMY
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y'
[]

[Stencil]
  [d2q9]
    type = LBMD2Q9
  []
[]

[TensorBuffers]
  [rho]
    type=LBMTensorBuffer
    buffer_type = ms
  []
  [u]
    type=LBMTensorBuffer
    buffer_type = mv
  []
  [speed]
    type=LBMTensorBuffer
    buffer_type = ms
  []
  [f]
    type=LBMTensorBuffer
    buffer_type = df
  []
  [feq]
    type=LBMTensorBuffer
    buffer_type = df
  []
  [f_post_collision]
    type=LBMTensorBuffer
    buffer_type = df
  []
  [binary_media]
    type = LBMTensorBuffer
    file = binary_media.h5
    is_integer = true
    buffer_type = ms
  []
[]

[TensorComputes]
  [Initialize]
    [rho]
      type=LBMConstantTensor
      buffer=rho
      constants = rho0
    []
    [u]
      type=LBMConstantTensor
      buffer=u
      constants = 'Ux Uy'
    []
    [speed]
      type=LBMComputeVelocityMagnitude
      buffer=speed
      velocity=u
    []
    [feq]
      type=LBMEquilibrium
      buffer=feq
      bulk=rho
      velocity=u
    []
    [f]
      type=LBMEquilibrium
      buffer=f
      bulk=rho
      velocity=u
    []
    [f_post_coll]
      type=LBMEquilibrium
      buffer=f_post_collision
      bulk=rho
      velocity=u
    []
  []
  [Solve]
    [Density]
      type = LBMComputeDensity
      buffer=rho
      f = f
    []
    [Velocity]
      type = LBMComputeVelocity
      buffer=u
      f = f
      rho = rho
    []
    [Equilibrium]
      type = LBMEquilibrium
      buffer=feq
      bulk=rho
      velocity=u
    []
    [Collision]
      type = LBMSmagorinskyCollision
      buffer = f_post_collision
      f = f
      feq = feq
      tau0 = tau
      projection=true
    []
    [Speed]
      type = LBMComputeVelocityMagnitude
      buffer = speed
      velocity = u
    []
    [Residual]
      type = LBMComputeResidual
      buffer = speed
      speed = speed
    []
  []
  [Boundary]
    [wall]
      type = LBMBounceBack
      buffer = f
      f_old = f_post_collision
      boundary = wall
    []
    [left]
      type = LBMFixedFirstOrderBC
      buffer=f
      f=f
      value='Ux'
      perturb=true
      boundary=left
    []
    [right]
      type = LBMMicroscopicZeroGradientBC
      buffer=f
      boundary=right
    []
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = f
  f_old = f_post_collision
[]

[Problem]
  type = LatticeBoltzmannProblem
  scalar_constant_names = 'rho0 Ux    Uy tau   dx    D    Cs'
  scalar_constant_values = '1.0 0.01 0.0 0.501 0.001 50 0.15'
  substeps = 5
  print_debug_output=true
  binary_media = binary_media
[]

[Postprocessors]
  [velocity_min]
    type = TensorExtremeValuePostprocessor
    buffer = u
    value_type = MIN
  []
  [velocity_max]
    type = TensorExtremeValuePostprocessor
    buffer = u
    value_type = MAX
  []
  [speed_min]
    type = TensorExtremeValuePostprocessor
    buffer = speed
    value_type = MIN
  []
  [speed_max]
    type = TensorExtremeValuePostprocessor
    buffer = speed
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = rho
    value_type = MIN
  []
  [densty_max]
    type = TensorExtremeValuePostprocessor
    buffer = rho
    value_type = MAX
  []
[]

[Executioner]
  type = Transient
  num_steps = 5
[]

[Outputs]
  file_base = obstacle
  csv = true
[]

(test/tests/lbm/neumann_regional.i)

[Domain]
  dim = 2
  nx = 11
  ny = 11
  mesh_mode=DUMMY
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y'
[]

[Stencil]
  [d2q9]
    type = LBMD2Q9
  []
[]

[TensorBuffers]
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [feq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [fpc]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [T]
    type=LBMTensorBuffer
    buffer_type = ms
  []
  [velocity]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  [binary_media]
    type = LBMTensorBuffer
    file = binary_regional.h5
    is_integer = true
    buffer_type = ms
  []
[]

[TensorComputes]
  [Initialize]
    [init_T]
      type = LBMConstantTensor
      buffer = T
      constants = 1.0
    []
    [equilibrium]
      type = LBMEquilibrium
      buffer = feq
      bulk = T
      velocity = velocity
    []
    [non_equilibrium]
      type = LBMEquilibrium
      buffer = f
      bulk = T
      velocity = velocity
    []
    [post_collision_equilibrium]
      type = LBMEquilibrium
      buffer = fpc
      bulk = T
      velocity = velocity
    []
  []
  [Solve]
    [T]
      type = LBMComputeDensity
      buffer = T
      f = f
    []
  []
  [Boundary]
    [regional]
      type = LBMNeumannBC
      buffer = f
      f_old = fpc
      feq = feq
      velocity = velocity
      rho = T
      gradient = 0.1
      region_id = 2
      boundary = regional
    []
  []
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 1
  binary_media = binary_media
[]

[Postprocessors]
  [velocity_min]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MIN
  []
  [velocity_max]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = T
    value_type = MIN
  []
  [densty_max]
    type = TensorExtremeValuePostprocessor
    buffer = T
    value_type = MAX
  []
[]

[Executioner]
  type = Transient
  num_steps = 5
[]

[Outputs]
  file_base = neumann_regional
  csv = true
[]

(test/tests/lbm/rarefied_gas.i)

[Domain]
  dim = 2
  nx = 48
  ny = 48
  xmax = 48
  ymax = 48
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y'
[]

[Stencil]
  [d2q9]
    type = LBMD2Q9
  []
[]

[TensorBuffers]
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [feq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [fpc]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [velocity]
    type=LBMTensorBuffer
    buffer_type = mv
  []
  [density]
    type=LBMTensorBuffer
    buffer_type = ms
  []
  [speed]
    type=LBMTensorBuffer
    buffer_type = ms
  []
  [domain]
    type=LBMTensorBuffer
    file = nanochannel.h5
    buffer_type = ms
    is_integer = true
  []
  [local_pore]
    type=LBMTensorBuffer
    file = local_pore_size.h5
    buffer_type = ms
    is_integer = false
  []
  [Kn]
    type=LBMTensorBuffer
    file = Kn.h5
    buffer_type = ms
    is_integer = false
  []
  [relaxation_matrix]
    type = LBMTensorBuffer
    buffer_type = df
  []
[]

[TensorComputes]
  [Initialize]
    [initial_density]
      type = LBMConstantTensor
      buffer = density
      constants = 0.2355545440759889
    []
    [initial_velocity]
      type = LBMConstantTensor
      buffer = velocity
      constants = '0.0 0.0'
    []
    [initial_equilibrium]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [initial_distribution]
      type = LBMEquilibrium
      buffer = f
      bulk = density
      velocity = velocity
    []
    [initial_distribution_pc]
      type = LBMEquilibrium
      buffer = fpc
      bulk = density
      velocity = velocity
    []
    [relaxation_matrix_init]
      type = LBMComputeEffectiveRelaxation
      buffer = relaxation_matrix
      local_pore_size = local_pore
      local_Knudsen_number = Kn
      mfp = 7.904614716131531e-10
      dx = 1.0e-9
      A2 = 0.82
    []
  []
  [Solve]
    [equilibrium]
      type=LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [collision]
      type = LBMMRTCollision
      buffer = fpc
      f = f
      feq = feq
      is_dynamic_relaxation = true
      local_relaxation_matrix = relaxation_matrix
      projection = true
    []
    [density]
      type = LBMComputeDensity
      buffer = density
      f = f
    []
    [velocity]
      type = LBMComputeVelocity
      buffer = velocity
      f = f
      rho = density
      add_body_force = true
      body_force_x = 1.0e-9
    []
    [speed]
      type = LBMComputeVelocityMagnitude
      buffer = speed
      velocity = velocity
    []
    [residual]
      type = LBMComputeResidual
      buffer = speed
      speed = speed
    []
  []
  [Boundary]
    [wall]
      type = LBMSpecularReflectionBoundary
      buffer = f
      f_old = fpc
      local_Knudsen_number = Kn
      boundary = wall
    []
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = f
  f_old = fpc
[]

[Postprocessors]
  [velocity_min]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MIN
  []
  [velocity_max]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MAX
  []
  [speed_min]
    type = TensorExtremeValuePostprocessor
    buffer = speed
    value_type = MIN
  []
  [speed_max]
    type = TensorExtremeValuePostprocessor
    buffer = speed
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MIN
  []
  [densty_max]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MAX
  []
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 10
  print_debug_output = true
  binary_media = domain
[]

[Executioner]
  type = Transient
  num_steps = 10
[]

[Outputs]
  file_base = rarefied_gas
  csv = true
[]

(test/tests/lbm/nee_2d_all_walls.i)

[Domain]
  dim = 2
  nx = 10
  ny = 10
  mesh_mode = DUMMY
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y'
[]

[Stencil]
  [d2q9]
    type = LBMD2Q9
  []
[]

[TensorBuffers]
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [feq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [fpc]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [velocity]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  [density]
    type = LBMTensorBuffer
    buffer_type = ms
  []
[]

[TensorComputes]
  [Initialize]
    [initial_density]
      type = LBMConstantTensor
      buffer = density
      constants = 1.0
    []
    [initial_velocity]
      type = LBMConstantTensor
      buffer = velocity
      constants = '0.0 0.0'
    []
    [initial_feq]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [initial_f]
      type = LBMEquilibrium
      buffer = f
      bulk = density
      velocity = velocity
    []
    [initial_fpc]
      type = LBMEquilibrium
      buffer = fpc
      bulk = density
      velocity = velocity
    []
  []
  [Solve]
    [equilibrium]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [collision]
      type = LBMBGKCollision
      buffer = fpc
      f = f
      feq = feq
      tau0 = 1.0
    []
    [density]
      type = LBMComputeDensity
      buffer = density
      f = f
    []
    [velocity]
      type = LBMComputeVelocity
      buffer = velocity
      f = f
      rho = density
    []
  []
  [Boundary]
    [top]
      type = LBMNonEquilibriumExtrapolation
      buffer = f
      prescribe_type = velocity
      ux = 0.005
      uy = 0.0
      order = second
      boundary = top
    []
    [bottom]
      type = LBMNonEquilibriumExtrapolation
      buffer = f
      prescribe_type = velocity
      ux = 0.0
      uy = 0.0
      order = second
      boundary = bottom
    []
    [left]
      type = LBMNonEquilibriumExtrapolation
      buffer = f
      prescribe_type = velocity
      ux = 0.0
      uy = 0.0
      order = second
      boundary = left
    []
    [right]
      type = LBMNonEquilibriumExtrapolation
      buffer = f
      prescribe_type = velocity
      ux = 0.0
      uy = 0.0
      order = second
      boundary = right
    []
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = f
  f_old = fpc
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 1
[]

[Postprocessors]
  [velocity_min]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MIN
  []
  [velocity_max]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MIN
  []
  [density_max]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MAX
  []
[]

[Executioner]
  type = Transient
  num_steps = 5
[]

[Outputs]
  file_base = nee_2d_all_walls
  csv = true
[]

(test/tests/lbm/isotropic_stencil_mrt.i)

[Domain]
  dim = 3
  nx = 10
  ny = 10
  nz = 10
  mesh_mode = DUMMY
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y Z'
[]

[Stencil]
  [d3q27]
    type = LBMD3Q27
  []
[]

[TensorBuffers]
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [feq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [fpc]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [velocity]
    type=LBMTensorBuffer
    buffer_type = mv
  []
  [density]
    type=LBMTensorBuffer
    buffer_type = ms
  []
  [speed]
    type=LBMTensorBuffer
    buffer_type = ms
  []
[]

[TensorComputes]
  [Initialize]
    [initial_density]
      type = LBMConstantTensor
      buffer = density
      constants = 1.0
    []
    [initial_velocity]
      type = LBMConstantTensor
      buffer = velocity
      constants = '0.0 0.0 0.0'
    []
    [initial_equilibrium]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [initial_distribution]
      type = LBMEquilibrium
      buffer = f
      bulk = density
      velocity = velocity
    []
    [initial_distribution_pc]
      type = LBMEquilibrium
      buffer = fpc
      bulk = density
      velocity = velocity
    []
  []
  [Solve]
    [equilibrium]
      type=LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [collision]
      type=LBMMRTCollision
      buffer = fpc
      f = f
      feq = feq
      tau0 = 1.0
    []
    [density]
      type = LBMComputeDensity
      buffer = density
      f = f
    []
    [velocity]
      type = LBMComputeVelocity
      buffer = velocity
      f = f
      rho = density
      add_body_force = true
      body_force_x = 0.0001
    []
    [speed]
      type = LBMComputeVelocityMagnitude
      buffer = speed
      velocity = velocity
    []
    [residual]
      type = LBMComputeResidual
      buffer = speed
      speed = speed
    []
  []
  [Boundary]
    [top]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = top
    []
    [bottom]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = bottom
    []
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = f
  f_old = fpc
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 11
[]

[Postprocessors]
  [velocity_min]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MIN
  []
  [velocity_max]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MAX
  []
  [speed_min]
    type = TensorExtremeValuePostprocessor
    buffer = speed
    value_type = MIN
  []
  [speed_max]
    type = TensorExtremeValuePostprocessor
    buffer = speed
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MIN
  []
  [densty_max]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MAX
  []
[]

[Executioner]
  type = Transient
  num_steps = 5
[]

[Outputs]
  file_base = isotropic_stencil_mrt
  csv = true
[]

(test/tests/lbm/neumann_wall.i)

[Domain]
  dim = 2
  nx = 11
  ny = 11
  mesh_mode=DUMMY
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y'
[]

[Stencil]
  [d2q9]
    type = LBMD2Q9
  []
[]

[TensorBuffers]
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [feq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [fpc]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [T]
    type=LBMTensorBuffer
    buffer_type = ms
  []
  [velocity]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  [binary_media]
    type = LBMTensorBuffer
    file = binary_media.h5
    is_integer = true
    buffer_type = ms
  []
[]

[TensorComputes]
  [Initialize]
    [init_T]
      type = LBMConstantTensor
      buffer = T
      constants = 1.0
    []
    [equilibrium]
      type = LBMEquilibrium
      buffer = feq
      bulk = T
      velocity = velocity
    []
    [non_equilibrium]
      type = LBMEquilibrium
      buffer = f
      bulk = T
      velocity = velocity
    []
    [post_collision_equilibrium]
      type = LBMEquilibrium
      buffer = fpc
      bulk = T
      velocity = velocity
    []
  []
  [Solve]
    [T]
      type = LBMComputeDensity
      buffer = T
      f = f
    []
  []
  [Boundary]
    [wall]
      type = LBMNeumannBC
      buffer = f
      f_old = fpc
      feq = feq
      velocity = velocity
      rho = T
      gradient = 0.1
      boundary = wall
    []
  []
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 1
  binary_media = binary_media
[]

[Postprocessors]
  [velocity_min]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MIN
  []
  [velocity_max]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = T
    value_type = MIN
  []
  [densty_max]
    type = TensorExtremeValuePostprocessor
    buffer = T
    value_type = MAX
  []
[]

[Executioner]
  type = Transient
  num_steps = 5
[]

[Outputs]
  file_base = neumann_wall
  csv = true
[]

(test/tests/lbm/neumann_box.i)

[Domain]
  dim = 3
  nx = 5
  ny = 5
  nz = 5
  mesh_mode = DUMMY
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y Z'
[]

[Stencil]
  [d3q19]
    type = LBMD3Q19
  []
[]

[TensorBuffers]
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [feq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [fpc]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [density]
    type=LBMTensorBuffer
    buffer_type = ms
  []
  [velocity]
    type = LBMTensorBuffer
    buffer_type = mv
  []
[]

[TensorComputes]
  [Initialize]
    [init_density]
      type = LBMConstantTensor
      buffer = density
      constants = 1.0
    []
    [equilibrium]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [non_equilibrium]
      type = LBMEquilibrium
      buffer = f
      bulk = density
      velocity = velocity
    []
    [post_collision_equilibrium]
      type = LBMEquilibrium
      buffer = fpc
      bulk = density
      velocity = velocity
    []
  []
  [Solve]
    [density]
      type = LBMComputeDensity
      buffer = density
      f = f
    []
  []
  [Boundary]
    [left]
      type = LBMNeumannBC
      buffer = f
      f_old = fpc
      feq = feq
      velocity = velocity
      rho = density
      gradient = 0.001
      boundary = left
    []
    [right]
      type = LBMNeumannBC
      buffer = f
      f_old = fpc
      feq = feq
      velocity = velocity
      rho = density
      gradient = 0.1
      boundary = right
    []
    [top]
      type = LBMNeumannBC
      buffer = f
      f_old = fpc
      feq = feq
      velocity = velocity
      rho = density
      gradient = 0.1
      boundary = top
    []
    [bottom]
      type = LBMNeumannBC
      buffer = f
      f_old = fpc
      feq = feq
      velocity = velocity
      rho = density
      gradient = 0.1
      boundary = bottom
    []
    [front]
      type = LBMNeumannBC
      buffer = f
      f_old = fpc
      feq = feq
      velocity = velocity
      rho = density
      gradient = 0.1
      boundary = front
    []
    [back]
      type = LBMNeumannBC
      buffer = f
      f_old = fpc
      feq = feq
      velocity = velocity
      rho = density
      gradient = 0.1
      boundary = back
    []
  []
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 1
[]

[Postprocessors]
  [velocity_min]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MIN
  []
  [velocity_max]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MIN
  []
  [densty_max]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MAX
  []
[]

[Executioner]
  type = Transient
  num_steps = 5
[]

[Outputs]
  file_base = neumann_box
  csv = true
[]

(examples/lbm/Phase-field/static_bubble.i)

#
# Static bubble test case
# PHYSICAL REVIEW E 97, 033309 (2018) - Section III.A
# Intended outcome is to have (spurious) velocity magnitude around 1.0e-10
#

# Domain
Nx = 200
Ny = 200

# Bubble parameters
Cx = '${Nx}/2.0'
Cy = '${Ny}/2.0'
R = 50

# Fluid properties
rho_l = 1000.0
rho_g = 1.0
nu_l = 0.1
nu_g = 0.1
sigma = 0.001

# Phase field parameters
tau_h = 0.7
D = 5

[Domain]
  dim = 2
  nx = '${Nx}'
  ny = '${Ny}'
  xmax = '${Nx}'
  ymax = '${Ny}'
  device_names='cpu'
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y'
[]

[Stencil]
  [d2q9]
    type = LBMD2Q9
  []
[]

[TensorBuffers]
  # Macroscopic phase field variables
  [phi]
    type = LBMTensorBuffer
    buffer_type = ms
  []
  [grad_phi]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  [laplacian_phi]
    type = LBMTensorBuffer
    buffer_type = ms
  []
  [mu]
    type = LBMTensorBuffer
    buffer_type = ms
  []
  [forces]
    type = LBMTensorBuffer
    buffer_type = mv
  []

  # Macroscopic hydrodynamic variables
  [velocity]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  [pressure]
    type = LBMTensorBuffer
    buffer_type = ms
  []
  [rho]
    type = LBMTensorBuffer
    buffer_type = ms
  []

  # LBM phase field variabels
  [h]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [h_post_collision]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [h_eq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [relaxation_tensor]
    type = LBMTensorBuffer
    buffer_type = ms
  []

  # LBM hydrodynamic variables
  [fdummy]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [f_post_collision]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [f_eq]
    type = LBMTensorBuffer
    buffer_type = df
  []
[]

[TensorComputes/Initialize]
  [phi_init]
    type = ParsedCompute
    buffer = phi
    expression = '0.5 + 0.5 * tanh(2*(R - sqrt((x - Cx)^2 + (y - Cy)^2)) / D)'
    constant_names = 'Cx Cy R D'
    constant_expressions = '${Cx} ${Cy} ${R} ${D}'
    extra_symbols = true
  []
  [grad_phi_init]
    type = LBMIsotropicGradient
    buffer = grad_phi
    scalar_field = phi
  []
  [rho_init]
    type = ParsedCompute
    buffer = rho
    extra_symbols = true
    expression = 'phi*(rho_l - rho_g) + rho_g'
    constant_names = 'rho_l rho_g'
    constant_expressions = '${rho_l} ${rho_g}'
    inputs = phi
  []
  [pressure_init]
    type = LBMConstantTensor
    buffer = pressure
    constants = 0.3
  []
  # Phase field equilibrium distribution initialization
  [h_eq_init]
    type = LBMPhaseEquilibrium
    buffer = h_eq
    phi = phi
    velocity = velocity
  []
  [h_post_collision_init]
    type = LBMPhaseEquilibrium
    buffer = h_post_collision
    phi = phi
    velocity = velocity
  []
  [h_init]
    type = LBMPhaseEquilibrium
    buffer = h
    phi = phi
    velocity = velocity
  []
  # Hydrodynamic equilibrium distribution initialization
  [f_eq_init]
    type = LBMPressureCorrectedEquilibrium
    buffer = f_eq
    rho = rho
    velocity = velocity
    pressure = pressure
  []
  [f_post_collision_init]
    type = ParsedCompute
    buffer = f_post_collision
    expression = 'f_eq'
    inputs = f_eq
  []
  [f_init]
    type = ParsedCompute
    buffer = f
    expression = 'f_eq'
    inputs = f_eq
  []
[]

[TensorComputes/Solve]
  # Phase Field
  [compute_phi]
    type = LBMComputeDensity
    buffer = phi
    f = h
  []
  [grad_phi]
    type = LBMIsotropicGradient
    buffer = grad_phi
    scalar_field = phi
  []
  [laplacian_phi]
    type = LBMIsotropicLaplacian
    buffer = laplacian_phi
    scalar_field = phi
  []
  [potential]
    type = LBMComputeChemicalPotential
    buffer = mu
    phi = phi
    laplacian_phi = laplacian_phi
    thickness = D
    sigma = sigma
  []
  [forces]
    type = LBMComputeSurfaceForces
    buffer = forces
    chemical_potential = mu
    grad_phi = grad_phi
  []
  # Hydrodynamics
  [density]
    type = ParsedCompute
    buffer = rho
    extra_symbols = true
    expression = 'phi*(rho_l - rho_g) + rho_g'
    constant_names = 'rho_l rho_g'
    constant_expressions = '${rho_l} ${rho_g}'
    inputs = phi
  []
  [velocity]
    type = LBMComputeVelocity
    buffer = velocity
    f = f
    rho = rho
    enable_forces = true
    forces = forces
  []
  # Phase-field
  [h_eq]
    type = LBMPhaseEquilibrium
    buffer = h_eq
    phi = phi
    velocity = velocity
  []
  [phase_collision]
    type = LBMBGKCollision
    buffer = h_post_collision
    f = h
    feq = h_eq
    tau0 = tau_h
  []
  [apply_forces_phase]
    type = LBMAllenCahnSource
    buffer = h_post_collision
    phi = phi
    velocity = velocity
    grad_phi = grad_phi
    tau = tau_h
    thickness = D
  []
  # Hydrodynamics
  [relaxation_tensor]
    type = ParsedCompute
    buffer = relaxation_tensor
    extra_symbols = true
    expression = '(phi*(nu_l - nu_g) + nu_g)/cs2+0.5'
    constant_names = 'nu_l nu_g cs2'
    constant_expressions = '${nu_l} ${nu_g} 0.3333'
    inputs = phi
  []
  [pressure]
    type = LBMPhaseFieldPressure
    buffer = pressure
    f = f
    velocity = velocity
    grad_phi = grad_phi
    rho = rho
    rho_l = '${rho_l}'
    rho_g = '${rho_g}'
  []
  [f_eq]
    type = LBMPressureCorrectedEquilibrium
    buffer = f_eq
    rho = rho
    velocity = velocity
    pressure = pressure
  []
  [collision]
    type = LBMBGKCollision
    buffer = f_post_collision
    f = f
    feq = f_eq
    tau0 = 1.0
    is_dynamic_relaxation = true
    tau_tensor = relaxation_tensor
   []
  [apply_forces_hydro]
    type = LBMForceDistribution
    buffer = f_post_collision
    grad_phi = grad_phi
    velocity = velocity
    forces = forces
    tau_tensor = relaxation_tensor
    tau = 1.0
    rho_l = '${rho_l}'
    rho_g = '${rho_g}'
    is_dynamic_relaxation = true
  []
  [residual]
    type = LBMComputeResidual
    buffer = phi
    speed = phi
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = 'h f'
  f_old = 'h_post_collision f_post_collision'
  root_compute = residual
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 50
  print_debug_output = true
  scalar_constant_names = 'tau_h D sigma'
  scalar_constant_values = '${tau_h} ${D} ${sigma}'
[]

[Executioner]
  type = Transient
  num_steps = 100
[]

[TensorOutputs]
  [xdmf]
    type = XDMFTensorOutput
    buffer = 'phi rho velocity forces'
    output_mode = 'Cell Cell Cell Cell'
    enable_hdf5 = true
    # transpose = false
  []
[]

(test/tests/lbm/phase.i)

# Domain
Nx = 20
Ny = 20

# Fluid properties
rho_l = 5.0
rho_g = 1.0
nu_l = 0.1
nu_g = 1.0
sigma = 0.2

# Phase field parameters
tau_h = 1.0
D = 4

[Domain]
  dim = 2
  nx = '${Nx}'
  ny = '${Ny}'
  xmax = '${Nx}'
  ymax = '${Ny}'
  device_names='cpu'
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y'
[]

[Stencil]
  [d2q9]
    type = LBMD2Q9
  []
[]

[TensorBuffers]
  # Macroscopic phase field variables
  [phi]
    type = LBMTensorBuffer
    buffer_type = ms
  []
  [grad_phi]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  [laplacian_phi]
    type = LBMTensorBuffer
    buffer_type = ms
  []
  [mu]
    type = LBMTensorBuffer
    buffer_type = ms
  []
  [forces]
    type = LBMTensorBuffer
    buffer_type = mv
  []

  # Macroscopic hydrodynamic variables
  [velocity]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  [pressure]
    type = LBMTensorBuffer
    buffer_type = ms
  []
  [rho]
    type = LBMTensorBuffer
    buffer_type = ms
  []

  # LBM phase field variabels
  [h]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [h_post_collision]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [h_eq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [relaxation_tensor]
    type = LBMTensorBuffer
    buffer_type = ms
  []

  # LBM hydrodynamic variables
  [fdummy]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [f_post_collision]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [f_eq]
    type = LBMTensorBuffer
    buffer_type = df
  []
[]

[TensorComputes/Initialize]
  [phi_init]
    type = ParsedCompute
    buffer = phi
    extra_symbols = true
    expression = '0.3333 + 0.01*sin((12.9898*x + 78.233*y)*2*pi)'
  []
  [grad_phi_init]
    type = LBMIsotropicGradient
    buffer = grad_phi
    scalar_field = phi
  []
  [rho_init]
    type = ParsedCompute
    buffer = rho
    extra_symbols = true
    expression = 'phi*(rho_l - rho_g) + rho_g'
    constant_names = 'rho_l rho_g'
    constant_expressions = '${rho_l} ${rho_g}'
    inputs = phi
  []
  [pressure_init]
    type = LBMConstantTensor
    buffer = pressure
    constants = 0.3
  []
  # Phase field equilibrium distribution initialization
  [h_eq_init]
    type = LBMPhaseEquilibrium
    buffer = h_eq
    phi = phi
    velocity = velocity
  []
  [h_post_collision_init]
    type = LBMPhaseEquilibrium
    buffer = h_post_collision
    phi = phi
    velocity = velocity
  []
  [h_init]
    type = LBMPhaseEquilibrium
    buffer = h
    phi = phi
    velocity = velocity
  []
  # Hydrodynamic equilibrium distribution initialization
  [f_eq_init]
    type = LBMPressureCorrectedEquilibrium
    buffer = f_eq
    rho = rho
    velocity = velocity
    pressure = pressure
  []
  [f_post_collision_init]
    type = ParsedCompute
    buffer = f_post_collision
    expression = 'f_eq'
    inputs = f_eq
  []
  [f_init]
    type = ParsedCompute
    buffer = f
    expression = 'f_eq'
    inputs = f_eq
  []
[]

[TensorComputes/Solve]
  # Phase Field
  [compute_phi]
    type = LBMComputeDensity
    buffer = phi
    f = h
  []
  [grad_phi]
    type = LBMIsotropicGradient
    buffer = grad_phi
    scalar_field = phi
  []
  [laplacian_phi]
    type = LBMIsotropicLaplacian
    buffer = laplacian_phi
    scalar_field = phi
  []
  [potential]
    type = LBMComputeChemicalPotential
    buffer = mu
    phi = phi
    laplacian_phi = laplacian_phi
    thickness = D
    sigma = sigma
  []
  [forces]
    type = LBMComputeSurfaceForces
    buffer = forces
    chemical_potential = mu
    grad_phi = grad_phi
  []
  # Hydrodynamics
  [density]
    type = ParsedCompute
    buffer = rho
    extra_symbols = true
    expression = 'phi*(rho_l - rho_g) + rho_g'
    constant_names = 'rho_l rho_g'
    constant_expressions = '${rho_l} ${rho_g}'
    inputs = phi
  []
  [velocity]
    type = LBMComputeVelocity
    buffer = velocity
    f = f
    rho = rho
    enable_forces = true
    forces = forces
  []
  # Phase-field
  [h_eq]
    type = LBMPhaseEquilibrium
    buffer = h_eq
    phi = phi
    velocity = velocity
  []
  [phase_collision]
    type = LBMBGKCollision
    buffer = h_post_collision
    f = h
    feq = h_eq
    tau0 = tau_h
  []
  [apply_forces_phase]
    type = LBMAllenCahnSource
    buffer = h_post_collision
    phi = phi
    velocity = velocity
    grad_phi = grad_phi
    tau = tau_h
    thickness = D
  []
  # Hydrodynamics
  [relaxation_tensor]
    type = ParsedCompute
    buffer = relaxation_tensor
    extra_symbols = true
    expression = '(phi*(nu_l - nu_g) + nu_g)/cs2+0.5'
    constant_names = 'nu_l nu_g cs2'
    constant_expressions = '${nu_l} ${nu_g} 0.3333'
    inputs = phi
  []
  [pressure]
    type = LBMPhaseFieldPressure
    buffer = pressure
    f = f
    velocity = velocity
    grad_phi = grad_phi
    rho = rho
    rho_l = '${rho_l}'
    rho_g = '${rho_g}'
  []
  [f_eq]
    type = LBMPressureCorrectedEquilibrium
    buffer = f_eq
    rho = rho
    velocity = velocity
    pressure = pressure
  []
  [collision]
    type = LBMBGKCollision
    buffer = f_post_collision
    f = f
    feq = f_eq
    tau0 = 1.0
    is_dynamic_relaxation = true
    tau_tensor = relaxation_tensor
   []
  [apply_forces_hydro]
    type = LBMForceDistribution
    buffer = f_post_collision
    grad_phi = grad_phi
    velocity = velocity
    forces = forces
    tau_tensor = relaxation_tensor
    tau = 1.0
    rho_l = '${rho_l}'
    rho_g = '${rho_g}'
    is_dynamic_relaxation = true
  []
  [residual]
    type = LBMComputeResidual
    buffer = phi
    speed = phi
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = 'h f'
  f_old = 'h_post_collision f_post_collision'
[]

[Postprocessors]
  [phi_min]
    type = TensorExtremeValuePostprocessor
    buffer = phi
    value_type = MIN
  []
  [phi_max]
    type = TensorExtremeValuePostprocessor
    buffer = phi
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = rho
    value_type = MIN
  []
  [density_max]
    type = TensorExtremeValuePostprocessor
    buffer = rho
    value_type = MAX
  []
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 5
  print_debug_output = true
  scalar_constant_names = 'tau_h D sigma'
  scalar_constant_values = '${tau_h} ${D} ${sigma}'
[]

[Executioner]
  type = Transient
  num_steps = 5
[]

[Outputs]
  file_base = phase
  csv = true
[]

(test/tests/lbm/vertical_velocity_bcs.i)

[Domain]
  dim = 2
  nx = 10
  ny = 10
  mesh_mode = DUMMY
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y'
[]

[Stencil]
  [d2q9]
    type = LBMD2Q9
  []
[]

[TensorBuffers]
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [f_bounce_back]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [velocity]
    type=LBMTensorBuffer
    buffer_type = mv
  []
  [density]
    type=LBMTensorBuffer
    buffer_type = ms
  []
[]

[TensorComputes]
  [Initialize]
    [initial_density]
      type = LBMConstantTensor
      buffer = density
      constants = 1.0
    []
    [initial_velocity]
      type = LBMConstantTensor
      buffer = velocity
      constants = '0.0001 0.0005'
    []
    [initial_f]
      type = LBMEquilibrium
      buffer = f
      bulk = density
      velocity = velocity
    []
    [initial_f_bb]
      type = LBMEquilibrium
      buffer = f_bounce_back
      bulk = density
      velocity = velocity
    []
  []
  [Solve]
    [density]
      type = LBMComputeDensity
      buffer = density
      f = f
    []
    [velocity]
      type = LBMComputeVelocity
      buffer = velocity
      f = f
      rho = density
      add_body_force = true
      body_force_x = 0.0001
    []
  []
  [Boundary]
    [left]
      type = LBMBounceBack
      buffer = f
      f_old = f_bounce_back
      boundary = left
    []
    [right]
      type = LBMBounceBack
      buffer = f
      f_old = f_bounce_back
      boundary = right
    []
    [top]
      type = LBMFixedFirstOrderBC
      buffer = f
      f = f
      value = 0.0001
      boundary = top
    []
    [bottom]
      type = LBMFixedFirstOrderBC
      buffer = f
      f = f
      value = 0.00011
      boundary = bottom
    []
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = f
  f_old = f_bounce_back
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 2
[]

[Postprocessors]
  [velocity_min]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MIN
  []
  [velocity_max]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MIN
  []
  [densty_max]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MAX
  []
[]

[Executioner]
  type = Transient
  num_steps = 5
[]

[Outputs]
  file_base = vertical_velocity_bcs
  csv = true
[]

(test/tests/lbm/convective_outflow_3d_top.i)

[Domain]
  dim = 3
  nx = 8
  ny = 8
  nz = 8
  mesh_mode = DUMMY
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y Z'
[]

[Stencil]
  [d3q19]
    type = LBMD3Q19
  []
[]

[TensorBuffers]
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [feq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [fpc]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [velocity]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  [density]
    type = LBMTensorBuffer
    buffer_type = ms
  []
[]

[TensorComputes]
  [Initialize]
    [initial_density]
      type = LBMConstantTensor
      buffer = density
      constants = 1.0
    []
    [initial_velocity]
      type = LBMConstantTensor
      buffer = velocity
      constants = '0.0 0.0 0.0'
    []
    [initial_feq]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [initial_f]
      type = LBMEquilibrium
      buffer = f
      bulk = density
      velocity = velocity
    []
    [initial_fpc]
      type = LBMEquilibrium
      buffer = fpc
      bulk = density
      velocity = velocity
    []
  []
  [Solve]
    [equilibrium]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [collision]
      type = LBMBGKCollision
      buffer = fpc
      f = f
      feq = feq
      tau0 = 1.0
    []
    [density]
      type = LBMComputeDensity
      buffer = density
      f = f
    []
    [velocity]
      type = LBMComputeVelocity
      buffer = velocity
      f = f
      rho = density
    []
  []
  [Boundary]
    [bottom]
      type = LBMNonEquilibriumExtrapolation
      buffer = f
      prescribe_type = velocity
      ux = 0.0
      uy = 0.005
      uz = 0.0
      order = first
      boundary = bottom
    []
    [top]
      type = LBMConvectiveOutflow
      buffer = f
      f_old = f
      convection_velocity = Uc
      boundary = top
    []
    [left]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = left
    []
    [right]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = right
    []
    [front]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = front
    []
    [back]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = back
    []
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = f
  f_old = fpc
[]

[Problem]
  type = LatticeBoltzmannProblem
  scalar_constant_names = 'Uc'
  scalar_constant_values = '0.005'
  substeps = 1
[]

[Postprocessors]
  [velocity_min]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MIN
  []
  [velocity_max]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MIN
  []
  [density_max]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MAX
  []
[]

[Executioner]
  type = Transient
  num_steps = 5
[]

[Outputs]
  file_base = convective_outflow_3d_top
  csv = true
[]

(test/tests/lbm/smagorinsky_mrt.i)

[Domain]
  dim = 2
  nx = 10
  ny = 10
  mesh_mode = DUMMY
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y'
[]

[Stencil]
  [d2q9]
    type = LBMD2Q9
  []
[]

[TensorBuffers]
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [feq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [fpc]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [velocity]
    type=LBMTensorBuffer
    buffer_type = mv
  []
  [density]
    type=LBMTensorBuffer
    buffer_type = ms
  []
  [speed]
    type=LBMTensorBuffer
    buffer_type = ms
  []
[]

[TensorComputes]
  [Initialize]
    [initial_density]
      type = LBMConstantTensor
      buffer = density
      constants = 1.0
    []
    [initial_velocity]
      type = LBMConstantTensor
      buffer = velocity
      constants = '0.0 0.0'
    []
    [initial_equilibrium]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [initial_distribution]
      type = LBMEquilibrium
      buffer = f
      bulk = density
      velocity = velocity
    []
    [initial_distribution_pc]
      type = LBMEquilibrium
      buffer = fpc
      bulk = density
      velocity = velocity
    []
  []
  [Solve]
    [root]
      [equilibrium]
        type=LBMEquilibrium
        buffer = feq
        bulk = density
        velocity = velocity
      []
      [collision]
        type=LBMSmagorinskyMRTCollision
        buffer = fpc
        f = f
        feq = feq
        tau0 = 1.0
      []
      [density]
        type = LBMComputeDensity
        buffer = density
        f = f
      []
      [velocity]
        type = LBMComputeVelocity
        buffer = velocity
        f = f
        rho = density
        add_body_force = true
        body_force_x = 0.0001
      []
      [speed]
        type = LBMComputeVelocityMagnitude
        buffer = speed
        velocity = velocity
      []
      [residual]
        type = LBMComputeResidual
        buffer = speed
        speed = speed
      []
    []
  []
  [Boundary]
    [top]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = top
    []
    [bottom]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = bottom
    []
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = f
  f_old = fpc
  root_compute = root
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 11
[]

[Postprocessors]
  [velocity_min]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MIN
  []
  [velocity_max]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MAX
  []
  [speed_min]
    type = TensorExtremeValuePostprocessor
    buffer = speed
    value_type = MIN
  []
  [speed_max]
    type = TensorExtremeValuePostprocessor
    buffer = speed
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MIN
  []
  [densty_max]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MAX
  []
[]

[Executioner]
  type = Transient
  num_steps = 5
[]

[Outputs]
  file_base = smagorinsky_mrt
  csv = true
[]

(examples/lbm/Formula1-aerodynamics/f1.i)

[Domain]
  dim = 3
  nx = 781
  ny = 274
  nz = 146
  xmax = 781
  ymax = 274
  zmax = 146
  device_names='cuda'
  floating_precision = 'single'
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y Z'
[]

[Stencil]
  [descriptor]
    type = LBMD3Q27
  []
[]

[TensorBuffers]
  [binary_media]
    type = LBMTensorBuffer
    buffer_type = ms
    file = 'binary_media.h5'
    is_integer = true
  []

  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [feq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [fpc]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [u]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  [speed]
    type=LBMTensorBuffer
    buffer_type = ms
  []
  [rho]
    type = LBMTensorBuffer
    buffer_type = ms
  []
  [force]
    type = LBMTensorBuffer
    buffer_type = mv
  []
[]

[TensorComputes]
  [Initialize]
    [density_initial]
      type = LBMConstantTensor
      buffer = rho
      constants = 1.0
    []

    [velocity_initial]
      type = LBMConstantTensor
      buffer = u
      constants = '0 0 0'
    []

    [equilibrium_init]
      type = LBMEquilibrium
      buffer = feq
      bulk = rho
      velocity = u
    []

    [equilibrium_f]
      type = LBMEquilibrium
      buffer = f
      bulk = rho
      velocity = u
    []

    [equilibrium_pc]
      type = LBMEquilibrium
      buffer = fpc
      bulk = rho
      velocity = u
    []
  []

  [Solve]
    [density]
      type = LBMComputeDensity
      buffer = rho
      f = f
    []

    [velocity]
      type = LBMComputeVelocity
      buffer = u
      f = f
      rho = rho
    []

    [equilibrium]
      type = LBMEquilibrium
      buffer = feq
      bulk = rho
      velocity = u
    []

    [collision]
      type = LBMSmagorinskyCollision
      buffer = fpc
      f = f
      feq = feq
      tau0 = 0.5001
      Cs = 0.15
      projection=true
    []

    [speed]
      type=LBMComputeVelocityMagnitude
      buffer=speed
      velocity=u
    []

    [residual]
      type = LBMComputeResidual
      buffer = speed
      speed = speed
    []
  []

  [Boundary]
    [wall]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = wall
    []

    [left]
      type = LBMFixedFirstOrderBC
      buffer=f
      f=f
      value=0.01
      boundary=left
      # perturb=true
    []

    [right]
      type = LBMMicroscopicZeroGradientBC
      buffer=f
      boundary = right
    []

  []
[]

[TensorSolver]
  type = LBMStream
  root_compute = residual
  buffer = f
  f_old = fpc
[]

[Postprocessors]
  [reynolds]
    type = ComputeReynoldsNumber
    buffer = speed
    tau = 0.5001
    diameter = 781
  []
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 100
  print_debug_output = true
  binary_media = binary_media
[]

[Executioner]
  type = Transient
  num_steps = 2000
[]

[TensorOutputs]
  [xdmf2]
    type = XDMFTensorOutput
    buffer = 'rho u'
    output_mode = 'Cell Cell'
    enable_hdf5 = true
  []
[]

(test/tests/lbm/mixed_bcs_d3q27_reverse.i)

[Domain]
  dim = 3
  nx = 10
  ny = 10
  nz = 10
  mesh_mode = DUMMY
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y Z'
[]

[Stencil]
  [d3q27]
    type = LBMD3Q27
  []
[]

[TensorBuffers]
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [f_bounce_back]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [velocity]
    type=LBMTensorBuffer
    buffer_type = mv
  []
  [density]
    type=LBMTensorBuffer
    buffer_type = ms
  []
[]

[TensorComputes]
  [Initialize]
    [initial_density]
      type = LBMConstantTensor
      buffer = density
      constants = 1.0
    []
    [initial_velocity]
      type = LBMConstantTensor
      buffer = velocity
      constants = '0.0001 0.0005 0.0'
    []
    [initial_f]
      type = LBMEquilibrium
      buffer = f
      bulk = density
      velocity = velocity
    []
    [initial_f_bb]
      type = LBMEquilibrium
      buffer = f_bounce_back
      bulk = density
      velocity = velocity
    []
  []
  [Solve]
    [density]
      type = LBMComputeDensity
      buffer = density
      f = f
    []
    [velocity]
      type = LBMComputeVelocity
      buffer = velocity
      f = f
      rho = density
      add_body_force = true
      body_force_x = 0.0001
    []
  []
  [Boundary]
    [right]
      type = LBMFixedZerothOrderBC
      buffer = f
      f = f
      value = 1.1
      boundary = right
    []
    [left]
      type = LBMFixedFirstOrderBC
      buffer = f
      f = f
      value = 0.0001
      boundary = left
    []
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = f
  f_old = f_bounce_back
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 1
[]

[Postprocessors]
  [velocity_min]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MIN
  []
  [velocity_max]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MIN
  []
  [densty_max]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MAX
  []
[]

[Executioner]
  type = Transient
  num_steps = 5
[]

[Outputs]
  file_base = mixed_bcs_d3q27_reverse
  csv = true
[]

(test/tests/lbm/phase_3D.i)

# Domain
Nx = 20
Ny = 20
Nz = 20

# Fluid properties
rho_l = 5.0
rho_g = 1.0
nu_l = 0.1
nu_g = 1.0
sigma = 0.2

# Phase field parameters
tau_h = 1.0
D = 4

[Domain]
  dim = 3
  nx = '${Nx}'
  ny = '${Ny}'
  nz = '${Nz}'
  xmax = '${Nx}'
  ymax = '${Ny}'
  zmax = '${Nz}'
  device_names='cpu'
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y Z'
[]

[Stencil]
  [d3q27]
    type = LBMD3Q27
  []
[]

[TensorBuffers]
  # Macroscopic phase field variables
  [phi]
    type = LBMTensorBuffer
    buffer_type = ms
  []
  [grad_phi]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  [laplacian_phi]
    type = LBMTensorBuffer
    buffer_type = ms
  []
  [mu]
    type = LBMTensorBuffer
    buffer_type = ms
  []
  [forces]
    type = LBMTensorBuffer
    buffer_type = mv
  []

  # Macroscopic hydrodynamic variables
  [velocity]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  [pressure]
    type = LBMTensorBuffer
    buffer_type = ms
  []
  [rho]
    type = LBMTensorBuffer
    buffer_type = ms
  []

  # LBM phase field variabels
  [h]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [h_post_collision]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [h_eq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [relaxation_tensor]
    type = LBMTensorBuffer
    buffer_type = ms
  []

  # LBM hydrodynamic variables
  [fdummy]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [f_post_collision]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [f_eq]
    type = LBMTensorBuffer
    buffer_type = df
  []
[]

[TensorComputes/Initialize]
  [phi_init]
    type = ParsedCompute
    buffer = phi
    extra_symbols = true
    expression = '0.3333 + 0.01*sin((12.9898*x + 78.233*y + 43.12*z)*2*pi)'
  []
  [grad_phi_init]
    type = LBMIsotropicGradient
    buffer = grad_phi
    scalar_field = phi
  []
  [rho_init]
    type = ParsedCompute
    buffer = rho
    extra_symbols = true
    expression = 'phi*(rho_l - rho_g) + rho_g'
    constant_names = 'rho_l rho_g'
    constant_expressions = '${rho_l} ${rho_g}'
    inputs = phi
  []
  [pressure_init]
    type = LBMConstantTensor
    buffer = pressure
    constants = 0.3
  []
  # Phase field equilibrium distribution initialization
  [h_eq_init]
    type = LBMPhaseEquilibrium
    buffer = h_eq
    phi = phi
    velocity = velocity
  []
  [h_post_collision_init]
    type = LBMPhaseEquilibrium
    buffer = h_post_collision
    phi = phi
    velocity = velocity
  []
  [h_init]
    type = LBMPhaseEquilibrium
    buffer = h
    phi = phi
    velocity = velocity
  []
  # Hydrodynamic equilibrium distribution initialization
  [f_eq_init]
    type = LBMPressureCorrectedEquilibrium
    buffer = f_eq
    rho = rho
    velocity = velocity
    pressure = pressure
  []
  [f_post_collision_init]
    type = ParsedCompute
    buffer = f_post_collision
    expression = 'f_eq'
    inputs = f_eq
  []
  [f_init]
    type = ParsedCompute
    buffer = f
    expression = 'f_eq'
    inputs = f_eq
  []
[]

[TensorComputes/Solve]
  # Phase Field
  [compute_phi]
    type = LBMComputeDensity
    buffer = phi
    f = h
  []
  [grad_phi]
    type = LBMIsotropicGradient
    buffer = grad_phi
    scalar_field = phi
  []
  [laplacian_phi]
    type = LBMIsotropicLaplacian
    buffer = laplacian_phi
    scalar_field = phi
  []
  [potential]
    type = LBMComputeChemicalPotential
    buffer = mu
    phi = phi
    laplacian_phi = laplacian_phi
    thickness = D
    sigma = sigma
  []
  [forces]
    type = LBMComputeSurfaceForces
    buffer = forces
    chemical_potential = mu
    grad_phi = grad_phi
  []
  # Hydrodynamics
  [density]
    type = ParsedCompute
    buffer = rho
    extra_symbols = true
    expression = 'phi*(rho_l - rho_g) + rho_g'
    constant_names = 'rho_l rho_g'
    constant_expressions = '${rho_l} ${rho_g}'
    inputs = phi
  []
  [velocity]
    type = LBMComputeVelocity
    buffer = velocity
    f = f
    rho = rho
    enable_forces = true
    forces = forces
  []
  # Phase-field
  [h_eq]
    type = LBMPhaseEquilibrium
    buffer = h_eq
    phi = phi
    velocity = velocity
  []
  [phase_collision]
    type = LBMBGKCollision
    buffer = h_post_collision
    f = h
    feq = h_eq
    tau0 = tau_h
  []
  [apply_forces_phase]
    type = LBMAllenCahnSource
    buffer = h_post_collision
    phi = phi
    velocity = velocity
    grad_phi = grad_phi
    tau = tau_h
    thickness = D
  []
  # Hydrodynamics
  [relaxation_tensor]
    type = ParsedCompute
    buffer = relaxation_tensor
    extra_symbols = true
    expression = '(phi*(nu_l - nu_g) + nu_g)/cs2+0.5'
    constant_names = 'nu_l nu_g cs2'
    constant_expressions = '${nu_l} ${nu_g} 0.3333'
    inputs = phi
  []
  [pressure]
    type = LBMPhaseFieldPressure
    buffer = pressure
    f = f
    velocity = velocity
    grad_phi = grad_phi
    rho = rho
    rho_l = '${rho_l}'
    rho_g = '${rho_g}'
  []
  [f_eq]
    type = LBMPressureCorrectedEquilibrium
    buffer = f_eq
    rho = rho
    velocity = velocity
    pressure = pressure
  []
  [collision]
    type = LBMBGKCollision
    buffer = f_post_collision
    f = f
    feq = f_eq
    tau0 = 1.0
    is_dynamic_relaxation = true
    tau_tensor = relaxation_tensor
   []
  [apply_forces_hydro]
    type = LBMForceDistribution
    buffer = f_post_collision
    grad_phi = grad_phi
    velocity = velocity
    forces = forces
    tau_tensor = relaxation_tensor
    tau = 1.0
    rho_l = '${rho_l}'
    rho_g = '${rho_g}'
    is_dynamic_relaxation = true
  []
  [residual]
    type = LBMComputeResidual
    buffer = phi
    speed = phi
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = 'h f'
  f_old = 'h_post_collision f_post_collision'
[]

[Postprocessors]
  [phi_min]
    type = TensorExtremeValuePostprocessor
    buffer = phi
    value_type = MIN
  []
  [phi_max]
    type = TensorExtremeValuePostprocessor
    buffer = phi
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = rho
    value_type = MIN
  []
  [density_max]
    type = TensorExtremeValuePostprocessor
    buffer = rho
    value_type = MAX
  []
  [velocity_min]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MIN
  []
  [velocity_max]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MAX
  []
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 5
  print_debug_output = true
  scalar_constant_names = 'tau_h D sigma'
  scalar_constant_values = '${tau_h} ${D} ${sigma}'
[]

[Executioner]
  type = Transient
  num_steps = 5
[]

[Outputs]
  file_base = phase_3D
  csv = true
[]

(test/tests/lbm/convective_outflow_3d_back.i)

[Domain]
  dim = 3
  nx = 8
  ny = 8
  nz = 8
  mesh_mode = DUMMY
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y Z'
[]

[Stencil]
  [d3q19]
    type = LBMD3Q19
  []
[]

[TensorBuffers]
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [feq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [fpc]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [velocity]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  [density]
    type = LBMTensorBuffer
    buffer_type = ms
  []
[]

[TensorComputes]
  [Initialize]
    [initial_density]
      type = LBMConstantTensor
      buffer = density
      constants = 1.0
    []
    [initial_velocity]
      type = LBMConstantTensor
      buffer = velocity
      constants = '0.0 0.0 0.0'
    []
    [initial_feq]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [initial_f]
      type = LBMEquilibrium
      buffer = f
      bulk = density
      velocity = velocity
    []
    [initial_fpc]
      type = LBMEquilibrium
      buffer = fpc
      bulk = density
      velocity = velocity
    []
  []
  [Solve]
    [equilibrium]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [collision]
      type = LBMBGKCollision
      buffer = fpc
      f = f
      feq = feq
      tau0 = 1.0
    []
    [density]
      type = LBMComputeDensity
      buffer = density
      f = f
    []
    [velocity]
      type = LBMComputeVelocity
      buffer = velocity
      f = f
      rho = density
    []
  []
  [Boundary]
    [front]
      type = LBMNonEquilibriumExtrapolation
      buffer = f
      prescribe_type = velocity
      ux = 0.0
      uy = 0.0
      uz = 0.005
      order = first
      boundary = front
    []
    [back]
      type = LBMConvectiveOutflow
      buffer = f
      f_old = f
      convection_velocity = Uc
      boundary = back
    []
    [left]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = left
    []
    [right]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = right
    []
    [top]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = top
    []
    [bottom]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = bottom
    []
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = f
  f_old = fpc
[]

[Problem]
  type = LatticeBoltzmannProblem
  scalar_constant_names = 'Uc'
  scalar_constant_values = '0.005'
  substeps = 1
[]

[Postprocessors]
  [velocity_min]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MIN
  []
  [velocity_max]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MIN
  []
  [density_max]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MAX
  []
[]

[Executioner]
  type = Transient
  num_steps = 5
[]

[Outputs]
  file_base = convective_outflow_3d_back
  csv = true
[]

(examples/lbm/Phase-field/layered_poiseuille.i)

#
# Layered Poiseuille Flow
# PHYSICAL REVIEW E 97, 033309 (2018) - Section III.B
#

# Domain
Nx = 10
Ny = 400

# Fluid properties
rho_l = 1000.0
rho_g = 1.0
# nu_l = 0.1
# nu_g = 1.0
sigma = 0.001
mu_l = 100.0 # rho_l * nu_l
mu_g = 1.0   # rho_g * nu_g

# Phase field parameters
# M = 0.1
# cs2 = 0.333333333333
tau_h = 0.8 # 0.5 + '${M}' / '${cs2}'
D = 5

# Driving force: Gx = uc * (mu_l + mu_g) / h^2
# uc = 1e-4
# h = # '${Ny} / 2'
Gx = 2.53e-07 #  '${uc} * (${mu_l}  + ${mu_g}) / (${h}^2)'

[Domain]
  dim = 2
  nx = '${Nx}'
  ny = '${Ny}'
  xmax = '${Nx}'
  ymax = '${Ny}'
  device_names = 'cpu'
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y'
[]

[Stencil]
  [d2q9]
    type = LBMD2Q9
  []
[]

[TensorBuffers]
  # Macroscopic phase field variables
  [phi]
    type = LBMTensorBuffer
    buffer_type = ms
    file = phi_init.h5
  []
  [grad_phi]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  [laplacian_phi]
    type = LBMTensorBuffer
    buffer_type = ms
  []
  [mu]
    type = LBMTensorBuffer
    buffer_type = ms
  []
  [forces]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  [body_force]
    type = LBMTensorBuffer
    buffer_type = mv
  []

  # Macroscopic hydrodynamic variables
  [velocity]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  [speed]
    type = LBMTensorBuffer
    buffer_type = ms
  []
  [pressure]
    type = LBMTensorBuffer
    buffer_type = ms
  []
  [rho]
    type = LBMTensorBuffer
    buffer_type = ms
  []

  # LBM phase field variables
  [h]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [h_post_collision]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [h_eq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [relaxation_tensor]
    type = LBMTensorBuffer
    buffer_type = ms
  []

  # LBM hydrodynamic variables
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [f_post_collision]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [f_eq]
    type = LBMTensorBuffer
    buffer_type = df
  []
[]

[TensorComputes/Initialize]
  [grad_phi_init]
    type = LBMIsotropicGradient
    buffer = grad_phi
    scalar_field = phi
  []
  [rho_init]
    type = ParsedCompute
    buffer = rho
    expression = 'phi * (rho_l - rho_g) + rho_g'
    constant_names = 'rho_l rho_g'
    constant_expressions = '${rho_l} ${rho_g}'
    inputs = phi
  []
  [pressure_init]
    type = LBMConstantTensor
    buffer = pressure
    constants = 0.3
  []
  [body_force_init]
    type = LBMConstantTensor
    buffer = body_force
    constants = '${Gx} 0.00'
  []
  # Phase field equilibrium distribution initialization
  [h_eq_init]
    type = LBMPhaseEquilibrium
    buffer = h_eq
    phi = phi
    velocity = velocity
  []
  [h_post_collision_init]
    type = LBMPhaseEquilibrium
    buffer = h_post_collision
    phi = phi
    velocity = velocity
  []
  [h_init]
    type = LBMPhaseEquilibrium
    buffer = h
    phi = phi
    velocity = velocity
  []
  # Hydrodynamic equilibrium distribution initialization
  [f_eq_init]
    type = LBMPressureCorrectedEquilibrium
    buffer = f_eq
    rho = rho
    velocity = velocity
    pressure = pressure
  []
  [f_post_collision_init]
    type = ParsedCompute
    buffer = f_post_collision
    expression = 'f_eq'
    inputs = f_eq
  []
  [f_init]
    type = ParsedCompute
    buffer = f
    expression = 'f_eq'
    inputs = f_eq
  []
[]

[TensorComputes/Solve]
  # Phase Field
  [compute_phi]
    type = LBMComputeDensity
    buffer = phi
    f = h
  []
  [grad_phi]
    type = LBMIsotropicGradient
    buffer = grad_phi
    scalar_field = phi
  []
  [laplacian_phi]
    type = LBMIsotropicLaplacian
    buffer = laplacian_phi
    scalar_field = phi
  []
  [potential]
    type = LBMComputeChemicalPotential
    buffer = mu
    phi = phi
    laplacian_phi = laplacian_phi
    thickness = D
    sigma = sigma
  []
  [compute_forces]
    type = LBMComputeSurfaceForces
    buffer = forces
    chemical_potential = mu
    grad_phi = grad_phi
  []
  [add_body_force]
    type = ParsedCompute
    buffer = forces
    expression = 'forces + body_force'
    inputs = 'forces body_force'
  []
  # Hydrodynamics
  [density]
    type = ParsedCompute
    buffer = rho
    expression = 'phi * (rho_l - rho_g) + rho_g'
    constant_names = 'rho_l rho_g'
    constant_expressions = '${rho_l} ${rho_g}'
    inputs = phi
  []
  [velocity]
    type = LBMComputeVelocity
    buffer = velocity
    f = f
    rho = rho
    enable_forces = true
    forces = forces
    # add_body_force = true
    # body_force_x = '${Gx}'
  []
  # Phase-field
  [h_eq]
    type = LBMPhaseEquilibrium
    buffer = h_eq
    phi = phi
    velocity = velocity
  []
  [phase_collision]
    type = LBMBGKCollision
    buffer = h_post_collision
    f = h
    feq = h_eq
    tau0 = tau_h
  []
  [apply_forces_phase]
    type = LBMAllenCahnSource
    buffer = h_post_collision
    phi = phi
    velocity = velocity
    grad_phi = grad_phi
    tau = tau_h
    thickness = D
  []
  # Hydrodynamics
  [relaxation_tensor]
    type = ParsedCompute
    buffer = relaxation_tensor
    # Implements Eq 26: Sharp step-function for dynamic viscosity
    expression = '(if(phi >= 0.5, mu_l, mu_g) / rho) / cs2 + 0.5'
    constant_names = 'mu_l mu_g cs2'
    constant_expressions = '${mu_l} ${mu_g} 0.333333333333'
    inputs = 'phi rho'
  []
  [pressure]
    type = LBMPhaseFieldPressure
    buffer = pressure
    f = f
    velocity = velocity
    grad_phi = grad_phi
    rho = rho
    rho_l = '${rho_l}'
    rho_g = '${rho_g}'
  []
  [f_eq]
    type = LBMPressureCorrectedEquilibrium
    buffer = f_eq
    rho = rho
    velocity = velocity
    pressure = pressure
  []
  [collision]
    type = LBMBGKCollision
    buffer = f_post_collision
    f = f
    feq = f_eq
    tau0 = 1.0
    is_dynamic_relaxation = true
    tau_tensor = relaxation_tensor
  []
  [apply_forces_hydro]
    type = LBMForceDistribution
    buffer = f_post_collision
    grad_phi = grad_phi
    velocity = velocity
    forces = forces
    tau_tensor = relaxation_tensor
    tau = 1.0
    rho_l = '${rho_l}'
    rho_g = '${rho_g}'
    is_dynamic_relaxation = true
  []
  [speed]
    type = LBMComputeVelocityMagnitude
    buffer = speed
    velocity = velocity
  []
  [residual]
    type = LBMComputeResidual
    buffer = speed
    speed = speed
  []
[]

[TensorComputes/Boundary]
  [top_fluid]
    type = LBMBounceBack
    buffer = f
    f_old = f_post_collision
    boundary = top
  []
  [bottom_fluid]
    type = LBMBounceBack
    buffer = f
    f_old = f_post_collision
    boundary = bottom
  []
  [top_phase]
    type = LBMBounceBack
    buffer = h
    f_old = h_post_collision
    boundary = top
  []
  [bottom_phase]
    type = LBMBounceBack
    buffer = h
    f_old = h_post_collision
    boundary = bottom
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = 'h f'
  f_old = 'h_post_collision f_post_collision'
  root_compute = residual
[]

[Problem]
  type = LatticeBoltzmannProblem
  # Keep this low for this setup: high substeps with top/bottom bounce-back can blow up to NaN.
  substeps = 100000
  print_debug_output = true
  scalar_constant_names = 'tau_h D sigma'
  scalar_constant_values = '${tau_h} ${D} ${sigma}'
[]

[Executioner]
  type = Transient
  num_steps = 2
[]

[TensorOutputs]
  [xdmf]
    type = XDMFTensorOutput
    buffer = 'phi velocity rho'
    output_mode = 'Cell Cell Cell'
    enable_hdf5 = true
    transpose = false
  []
[]

(test/tests/lbm/dirichlet_box.i)

[Domain]
  dim = 3
  nx = 5
  ny = 5
  nz = 5
  mesh_mode = DUMMY
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y Z'
[]

[Stencil]
  [d3q19]
    type = LBMD3Q19
  []
[]

[TensorBuffers]
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [feq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [fpc]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [density]
    type=LBMTensorBuffer
    buffer_type = ms
  []
  [velocity]
    type = LBMTensorBuffer
    buffer_type = mv
  []
[]

[TensorComputes]
  [Initialize]
    [init_density]
      type = LBMConstantTensor
      buffer = density
      constants = 1.0
    []
    [equilibrium]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []
    [non_equilibrium]
      type = LBMEquilibrium
      buffer = f
      bulk = density
      velocity = velocity
    []
    [post_collision_equilibrium]
      type = LBMEquilibrium
      buffer = fpc
      bulk = density
      velocity = velocity
    []
  []
  [Solve]
    [density]
      type = LBMComputeDensity
      buffer = density
      f = f
    []
  []
  [Boundary]
    [left]
      type = LBMDirichletBC
      buffer = f
      f_old = fpc
      feq = feq
      velocity = velocity
      rho = density
      value = 1.1
      boundary = left
    []
    [right]
      type = LBMDirichletBC
      buffer = f
      f_old = fpc
      feq = feq
      velocity = velocity
      rho = density
      value = 1.1
      boundary = right
    []
    [top]
      type = LBMDirichletBC
      buffer = f
      f_old = fpc
      feq = feq
      velocity = velocity
      rho = density
      value = 1.1
      boundary = top
    []
    [bottom]
      type = LBMDirichletBC
      buffer = f
      f_old = fpc
      feq = feq
      velocity = velocity
      rho = density
      value = 1.1
      boundary = bottom
    []
    [front]
      type = LBMDirichletBC
      buffer = f
      f_old = fpc
      feq = feq
      velocity = velocity
      rho = density
      value = 1.1
      boundary = front
    []
    [back]
      type = LBMDirichletBC
      buffer = f
      f_old = fpc
      feq = feq
      velocity = velocity
      rho = density
      value = 1.1
      boundary = back
    []
  []
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 1
[]

[Postprocessors]
  [velocity_min]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MIN
  []
  [velocity_max]
    type = TensorExtremeValuePostprocessor
    buffer = velocity
    value_type = MAX
  []
  [density_min]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MIN
  []
  [densty_max]
    type = TensorExtremeValuePostprocessor
    buffer = density
    value_type = MAX
  []
[]

[Executioner]
  type = Transient
  num_steps = 5
[]

[Outputs]
  file_base = dirichlet_box
  csv = true
[]

(examples/lbm/Pebbles/pebbles.i)

[Domain]
  dim = 3
  nx = 240
  ny = 50
  nz = 50
  xmax = 240
  ymax = 50
  zmax = 50
  device_names='cpu'
  # floating_precision = 'single'
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y'
[]

[Stencil]
  [d3q19]
    type = LBMD3Q19
  []
[]

[TensorBuffers]
  # Density distribution functions
  [f]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [feq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [fpc]
    type = LBMTensorBuffer
    buffer_type = df
  []
  # Temperature distribution functions
  [g]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [geq]
    type = LBMTensorBuffer
    buffer_type = df
  []
  [gpc]
    type = LBMTensorBuffer
    buffer_type = df
  []
  # Fluid macroscopic variables: density and velocity
  [density]
    type = LBMTensorBuffer
    buffer_type = ms
  []
  [velocity]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  [speed]
    type=LBMTensorBuffer
    buffer_type = ms
  []
  # Temperature macroscpic variables: temperature and 'velocity'
  [T]
    type = LBMTensorBuffer
    buffer_type = ms
  []
  # Forces
  [F]
    type = LBMTensorBuffer
    buffer_type = mv
  []
  # Simulation domain
  [binary_media]
    type = LBMTensorBuffer
    buffer_type = ms
    file = 'binary_media.h5'
    is_integer = true
  []
[]

[TensorComputes]

  #### Initialzie ####
  [Initialize]
    [density]
      type = LBMConstantTensor
      buffer = density
      constants = 'rho0'
    []

    [velocity]
      type = LBMConstantTensor
      buffer = velocity
      constants = '0.0 0.0'
    []

    [temperature]
      type = LBMConstantTensor
      buffer = T
      constants = T_C
    []

    [equilibrium_fluid]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []

    [equilibrium_fluid_total]
      type = LBMEquilibrium
      buffer = f
      bulk = density
      velocity = velocity
    []

    [equilibrium_fluid_pc]
      type = LBMEquilibrium
      buffer = fpc
      bulk = density
      velocity = velocity
    []

    [equilibrium_temperature]
      type = LBMEquilibrium
      buffer = geq
      bulk = T
      velocity = velocity
    []

    [equilibrium_temperature_total]
      type = LBMEquilibrium
      buffer = g
      bulk = T
      velocity = velocity
    []

    [equilibrium_temperature_pc]
      type = LBMEquilibrium
      buffer = gpc
      bulk = T
      velocity = velocity
    []

    # [t_wall_hot_init]
    #   type = LBMDirichletWallBC
    #   buffer = g
    #   f_old = gpc
    #   value = T_H
    #   velocity = velocity
    #   boundary = wall
    # []
  []

  #### Compute ####
  [Solve]

    # For temperature
    [Temperature]
      type = LBMComputeDensity
      buffer = T
      f = g
    []

    # For fluid
    [Fluid_density]
      type = LBMComputeDensity
      buffer = density
      f = f
    []

    [Fluid_velocity]
      type = LBMComputeVelocity
      buffer = velocity
      f = f
      rho = density
      forces = F
      enable_forces = true
    []

    # For temperature
    [Equilibrium_temperature]
      type = LBMEquilibrium
      buffer = geq
      bulk = T
      velocity = velocity
    []

    [Collision_temperature]
      type = LBMMRTCollision
      buffer = gpc
      f = g
      feq = geq
      tau0 = tau_T
    []

    # For fluid
    [Compute_forces]
      type = LBMComputeForces
      buffer = F
      rho0 = 'rho0'
      temperature = T
      T0 = T_H
      enable_buoyancy = true
      gravity = g
      gravity_direction=0
    []

    [Equilibrium_fluid]
      type = LBMEquilibrium
      buffer = feq
      bulk = density
      velocity = velocity
    []

    [Collision_fluid]
      type = LBMMRTCollision
      buffer = fpc
      f = f
      feq = feq
      tau0 = tau_f
    []

    [Apply_forces]
      type = LBMApplyForces
      buffer = fpc
      velocity = velocity
      rho = density
      forces = F
      tau0 = tau_f
    []

    [speed]
      type=LBMComputeVelocityMagnitude
      buffer=speed
      velocity=velocity
    []

    [residual]
      type = LBMComputeResidual
      buffer = speed
      speed = speed
    []
  []

  #### Boundary ####
  [Boundary]
    ##### for fluid
    [inlet]
      type = LBMFixedFirstOrderBC
      buffer = f
      f = f
      value = u0
      boundary = left
    []
    [outlet]
      type = LBMMicroscopicZeroGradientBC
      buffer = f
      boundary = right
    []
    [top]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = top
      # exclude_corners_x = true
    []
    [bottom]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = bottom
      # exclude_corners_x = true
    []
    [front]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = front
      # exclude_corners_x = true
    []
    [back]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = back
      # exclude_corners_x = true
    []
    [wall]
      type = LBMBounceBack
      buffer = f
      f_old = fpc
      boundary = wall
    []

    ##### for temperature
    [t_inlet]
      type = LBMFixedZerothOrderBC
      buffer = g
      f = g
      value = T_C
      boundary = left
    []
    [t_outlet]
      type = LBMMicroscopicZeroGradientBC
      buffer = g
      boundary = right
    []
    [t_top]
      type = LBMBounceBack
      buffer = g
      f_old = gpc
      boundary = top
    []
    [t_bottom]
      type = LBMBounceBack
      buffer = g
      f_old = gpc
      boundary = bottom
    []
    [t_front]
      type = LBMBounceBack
      buffer = g
      f_old = gpc
      boundary = front
    []
    [t_back]
      type = LBMBounceBack
      buffer = g
      f_old = gpc
      boundary = back
    []
    [t_wall_hot]
      type = LBMDirichletBC
      buffer = g
      f_old = gpc
      feq = geq
      velocity = velocity
      rho = T
      value = 1.2
      boundary = wall
    []
  []
[]

[TensorSolver]
  type = LBMStream
  buffer = 'f g'
  f_old = 'fpc gpc'
  root_compute = residual
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 20
  scalar_constant_names = 'rho0 T_C T_H tau_f tau_T g u0'
  scalar_constant_values = '1.0 1.0 1.2 0.9 0.9 0.001 0.001'

  is_binary_media = true
  binary_media = binary_media
[]

[Postprocessors]
  [reynolds]
    type = ComputeReynoldsNumber
    buffer = speed
    tau = tau_f
    diameter = 15
  []
[]

[Executioner]
  type = Transient
  num_steps = 2000
[]

# [TensorOutputs]
#   [xdmf2]
#     type = XDMFTensorOutput
#     buffer = 'binary_media'
#     output_mode = 'Cell'
#     enable_hdf5 = true
#   []
# []

[TensorOutputs]
  [xdmf2]
    type = XDMFTensorOutput
    buffer = 'T velocity density binary_media'
    output_mode = 'Cell Cell Cell Cell'
    enable_hdf5 = true
  []
[]

Overview
Example Input File Syntax
Input Parameters
Input Files