ParsedCompute

ParsedCompute object.

Overview

Evaluates a user-provided algebraic expression to compute the target buffer from named inputs and coordinates. Specify the expression via "expression" and list input buffers with "inputs". Enabling "extra_symbols" adds symbols like x, y, z.

Example Input File Syntax

[TensorComputes<<<{"href": "../../syntax/TensorComputes/index.html"}>>>]
  [Initialize<<<{"href": "../../syntax/TensorComputes/Initialize/index.html"}>>>]
    [diff]
      type = ParsedCompute<<<{"description": "ParsedCompute object.", "href": "ParsedCompute.html"}>>>
      buffer<<<{"description": "The buffer this compute is writing to"}>>> = diff
      inputs<<<{"description": "Buffer names used in the expression"}>>> = 'grad_sq c2'
      expression<<<{"description": "Parsed expression"}>>> = 'abs(grad_sq - c2)'
    []
  []
[]

Input Parameters

bufferThe buffer this compute is writing to
C++ Type:std::string
Controllable:No
Description:The buffer this compute is writing to
expressionParsed expression
C++ Type:std::string
Controllable:No
Description:Parsed expression

Required Parameters

constant_expressionsVector of values for the constants in constant_names (can be an FParser expression)
C++ Type:std::vector<std::string>
Controllable:No
Description:Vector of values for the constants in constant_names (can be an FParser expression)
constant_namesVector of constants used in the parsed function (use this for kB etc.)
C++ Type:std::vector<std::string>
Controllable:No
Description:Vector of constants used in the parsed function (use this for kB etc.)
derivativesList of inputs to take the derivative w.r.t. (or none)
C++ Type:std::vector<std::string>
Controllable:No
Description:List of inputs to take the derivative w.r.t. (or none)
enable_fpoptimizerTrueUse algebraic optimizer
Default:True
C++ Type:bool
Controllable:No
Description:Use algebraic optimizer
expandNONEExpand the tensor to full size.
Default:NONE
C++ Type:MooseEnum
Options:REAL, RECIPROCAL, NONE
Controllable:No
Description:Expand the tensor to full size.
extra_symbolsFalseProvide i (imaginary unit), kx,ky,kz (reciprocal space frequency), k2 (square of the k-vector), x,y,z (real space coordinates), time t, pi, and e.
Default:False
C++ Type:bool
Controllable:No
Description:Provide i (imaginary unit), kx,ky,kz (reciprocal space frequency), k2 (square of the k-vector), x,y,z (real space coordinates), time t, pi, and e.
inputsBuffer names used in the expression
C++ Type:std::vector<std::string>
Controllable:No
Description:Buffer names used in the expression
is_integerFalseTurn the function result into an integer tensor
Default:False
C++ Type:bool
Controllable:No
Description:Turn the function result into an integer tensor

Optional 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/cahn_hilliard/cahnhilliard4.i)
(examples/lbm/Phase-field/static_bubble_3d.i)
(examples/swift_hohenberg/rotating_grain.i)
(examples/cahn_hilliard/cahnhilliard.i)
(test/tests/real_space/laplacian3d.i)
(test/tests/solvers/etdrk4_diffusion.i)
(test/tests/typed_tensors/gradient_vector.i)
(benchmarks/01_spinodal_decomposition/1a.i)
(test/tests/tensor_compute/coupled_pf_mech.i)
(examples/cahn_hilliard/cahnhilliard3.i)
(test/tests/mechanics/mech.i)
(benchmarks/02_oswald_ripening/2a.i)
(test/tests/typed_tensors/gradient.i)
(test/tests/parsed_tensor/local_vars_derivative.i)
(benchmarks/01_spinodal_decomposition/1a_solver.i)
(test/tests/gradient/gradient.i)
(benchmarks/01_spinodal_decomposition/1a_secant.i)
(test/tests/real_space/diffusion.i)
(test/tests/tensor_compute/backandforth.i)
(test/tests/kks/KKS_no_flux_bc.i)
(test/tests/solvers/diagonal.i)
(test/tests/mechanics/mech3d.i)
(test/tests/tensor_compute/rotating_grain_secant.i)
(examples/lbm/Phase-field/droplet_impact.i)
(examples/libtorch_kks/KKS_libtorch.i)
(examples/lbm/Phase-field/spinodal_decomposition.i)
(test/tests/cahnhilliard/cahnhilliard_explicit_smooth.i)
(test/tests/cahnhilliard/cahnhilliard.i)
(test/tests/cahnhilliard/cahnhilliard_explicit.i)
(examples/swift_hohenberg/swifthohenberg.i)
(test/tests/tensor_compute/coupled_pf_mech_secant.i)
(test/tests/tensor_compute/parallel.i)
(test/tests/real_space/laplacian.i)
(examples/degeus_mechanics/mech.i)
(test/tests/solvers/nl_coupled.i)
(benchmarks/02_oswald_ripening/2a_broyden.i)
(examples/lbm/Phase-field/static_bubble.i)
(test/tests/tensor_compute/parallel_roundtrip_3d.i)
(benchmarks/02_oswald_ripening/2a_secant.i)
(test/tests/gradient/gradient_square.i)
(test/tests/solvers/coupled.i)
(test/tests/lbm/phase.i)
(test/tests/postprocessors/postprocessors.i)
(test/tests/postprocessors/interface_velocity.i)
(test/tests/lbm/stacking.i)
(examples/cahn_hilliard/cahnhilliard2.i)
(test/tests/tensor_compute/group.i)
(benchmarks/01_spinodal_decomposition/1b.i)
(test/tests/histogram/test.i)
(test/tests/tensor_compute/test.i)
(benchmarks/02_oswald_ripening/simple.i)
(examples/lbm/Rayleigh-Benard/rayleigh-benard.i)
(test/tests/lbm/phase_3D.i)
(test/tests/kks/KKS_libtorch.i)
(examples/lbm/Phase-field/layered_poiseuille.i)
(test/tests/tensor_compute/parallel_roundtrip.i)

(test/tests/gradient/gradient_square.i)

[Domain]
  dim = 3
  nx = 40
  ny = 40
  nz = 40
  xmax = ${fparse pi*2}
  ymax = ${fparse pi*4}
  zmax = ${fparse pi*6}
  mesh_mode = DUMMY
  device_names = cpu
[]

[TensorBuffers]
  [s]
  []
  [grad_sq]
  []
  [c2]
  []
  [diff]
  []
[]

[TensorComputes]
  [Initialize]
    [sin]
      type = ParsedCompute
      buffer = s
      extra_symbols = true
      expression = 'sin(x)+sin(y)+sin(z)'
    []
    [cos2]
      type = ParsedCompute
      buffer = c2
      extra_symbols = true
      expression = 'cos(x)^2+cos(y)^2+cos(z)^2'
    []
    [grad_sq]
      type = FFTGradientSquare
      buffer = grad_sq
      input = s
    []
    [diff]
      type = ParsedCompute
      buffer = diff
      inputs = 'grad_sq c2'
      expression = 'abs(grad_sq - c2)'
    []
  []
[]

[Postprocessors]
  [diff]
    type = TensorIntegralPostprocessor
    buffer = diff
  []
[]

[Problem]
  type = TensorProblem
[]

[Executioner]
  type = Transient
  num_steps = 1
[]

[Outputs]
  csv = true
[]

(examples/cahn_hilliard/cahnhilliard4.i)

[Domain]
  dim = 3
  nx = 100
  ny = 100
  nz = 100
  xmax = ${fparse pi*4}
  ymax = ${fparse pi*4}
  zmax = ${fparse pi*4}

  device_names = 'cuda'

  mesh_mode = DUMMY
  debug = true
[]

[TensorBuffers]
  [c]
  []
  [cbar]
  []
  [mu]
  []
  [mubar]
  []
  [Mbarmubar]
  []
  # constant tensors
  [Mbar]
  []
  [kappabarbar]
  []
[]

[TensorOutputs]
  [xdmf]
    type = XDMFTensorOutput
    buffer = 'c mu'
    output_mode = 'Node Cell'
    enable_hdf5 = true
  []
[]

[TensorComputes]
  [Initialize]
    [c]
      type = RandomTensor
      buffer = c
      min = 0.44
      max = 0.56
    []
    [Mbar]
      type = ReciprocalLaplacianFactor
      factor = 0.2 # Mobility
      buffer = Mbar
    []
    [kappabarbar]
      type = ReciprocalLaplacianSquareFactor
      factor = -0.001 # kappa
      buffer = kappabarbar
    []
  []

  [Solve]
    [mu]
      type = ParsedCompute
      buffer = mu
      expression = '0.1*c^2*(c-1)^2 + c*sin(x/2)*0.005'
      extra_symbols = true
      derivatives = c
      # expression = "0.4*c^3-0.6*c^2+0.2*c"
      inputs = c
    []
    [mubar]
      type = ForwardFFT
      buffer = mubar
      input = mu
    []
    [Mbarmubar]
      type = ParsedCompute
      buffer = Mbarmubar
      expression = 'Mbar*mubar'
      inputs = 'Mbar mubar'
    []
    [cbar]
      type = ForwardFFT
      buffer = cbar
      input = c
    []
  []
[]

[TensorTimeIntegrators]
  [c]
    type = FFTSemiImplicit
    buffer = c
    reciprocal_buffer = cbar
    linear_reciprocal = kappabarbar
    nonlinear_reciprocal = Mbarmubar
  []
[]

[Postprocessors]
  [min_c]
    type = TensorExtremeValuePostprocessor
    buffer = c
    value_type = MIN
    execute_on = 'TIMESTEP_END'
  []
  [max_c]
    type = TensorExtremeValuePostprocessor
    buffer = c
    value_type = MAX
    execute_on = 'TIMESTEP_END'
  []
  [C]
    type = TensorIntegralPostprocessor
    buffer = c
    execute_on = 'TIMESTEP_END'
  []
  [cavg]
    type = TensorAveragePostprocessor
    buffer = c
    execute_on = 'TIMESTEP_END'
  []
[]

[Problem]
  type = TensorProblem
  spectral_solve_substeps = 1000
[]

[Executioner]
  type = Transient
  num_steps = 50
  [TimeStepper]
    type = IterationAdaptiveDT
    growth_factor = 1.8
    dt = 0.1
  []
  dtmax = 500
[]

[Outputs]
  csv = true
  perf_graph = true
  execute_on = 'TIMESTEP_END'
[]

(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
  []
[]

(examples/swift_hohenberg/rotating_grain.i)

#
# Solve a simple Swift-Hohenberg crystal phase field problem. The initial condition is
# a circular grain that is rotated against the surropunding matrix.
# This example demonstrates the use of the [TensorComputes/Postprocess] system to perform
# compute steps just prior to running the output objects. Here we perform a low-pass filtering
# by forward transfroming the psi amplitude field into frequency space, attenuating frequencies
# by the exponent of their wave number, and transforming back into real space. This filtering
# makes the dislocation structure in the crystal more pronounced in the visualization.
#

w=60
[Domain]
  dim = 2
  nx = 400
  ny = 400
  xmax = ${fparse pi*2*w}
  ymax = ${fparse pi*2*w}

  device_names = 'cuda'

  mesh_mode = DOMAIN
[]

[TensorBuffers]
  [psi]
    map_to_aux_variable = psi
  []
  [psibar]
  []
  [psi3]
  []
  [psi3bar]
  []
  # constant tensors
  [linear]
  []

  # output
  [filter]
    map_to_aux_variable = filter
  []
  [filterbar]
  []
[]

[AuxVariables]
  [psi]
  []
  [filter]
  []
[]

crystal = '-sin(sin(a)*y/2+cos(a)*x/2)^2*sin(sin(a+1/3*pi)*y/2+cos(a+1/3*pi)*x/2)^2*sin(sin(a-1/3*pi)*y/2+cos(a-1/3*pi)*x/2)^2'
[Functions]
  [grain1]
    type = ParsedFunction
    expression = 'a := 0; ${crystal}'
  []
  [grain2]
    type = ParsedFunction
    expression = 'a := 0.95; ${crystal}'
  []
  [domain]
    type = ParsedFunction
    expression = 'r := (x-${w}*pi)^2+(y-${w}*pi)^2; if(r<(${w}*2/3*pi)^2, grain2, grain1)'
    symbol_names = 'grain1 grain2'
    symbol_values = 'grain1 grain2'
  []
[]

[TensorComputes]
  [Initialize]
    [psi]
      type = MooseFunctionTensor
      buffer = psi
      function = domain
    []
    [linear]
      type = SwiftHohenbergLinear
      buffer = linear
      alpha = 1
      r = 0.025
    []
  []

  [Solve]
    [psi3]
      type = ParsedCompute
      buffer = psi3
      expression = "0.20*psi^2-psi^3"
      inputs = psi
    []
    [psibar]
      type = ForwardFFT
      buffer = psibar
      input = psi
    []
    [psi3bar]
      type = ForwardFFT
      buffer = psi3bar
      input = psi3
    []
  []

  [Postprocess]
    [low_pass]
      type = ParsedCompute
      buffer = filterbar
      extra_symbols = true
      expression = 'psibar * exp(-k2*10)'
      inputs = psibar
    []
    [filter]
      type = InverseFFT
      buffer = filter
      input = filterbar
    []
  []
[]

[TensorSolver]
  type = AdamsBashforthMoulton
  buffer = psi
  reciprocal_buffer = psibar
  linear_reciprocal = linear
  nonlinear_reciprocal = psi3bar
  substeps = 100
[]

[Problem]
  type = TensorProblem
[]

[Executioner]
  type = Transient
  num_steps = 120
  [TimeStepper]
    type = IterationAdaptiveDT
    growth_factor = 1.1
    dt = 5
  []
  dtmax = 500
[]

[Postprocessors]
  [min_psi]
    type = TensorExtremeValuePostprocessor
    buffer = psi
    value_type = MIN
    execute_on = 'TIMESTEP_END'
  []
  [max_psi]
    type = TensorExtremeValuePostprocessor
    buffer = psi
    value_type = MAX
    execute_on = 'TIMESTEP_END'
  []
  [Psi]
    type = TensorIntegralPostprocessor
    buffer = psi
  []
[]

[Outputs]
  exodus = true
  perf_graph = true
  execute_on = 'TIMESTEP_END'
[]

(examples/cahn_hilliard/cahnhilliard.i)

#
# Simple Cahn-Hilliard solve on a 2D grid. We create a matching (conforming)
# MOOSE mesh (with one element per FFT grid cell) and project the solution onto
# the MOOSE mesh to utilize the exodus output object.
#

[Domain]
  dim = 2
  nx = 200
  ny = 200

  xmax = ${fparse pi*8}
  ymax = ${fparse pi*8}

  # automatically create a matching mesh
  mesh_mode = DOMAIN
[]


[TensorBuffers]
  [c]
    # perform fast mapping to the matching mesh by directly writing to
    # the solution vector of the specified Auxvariable
    map_to_aux_variable = c
  []
  [cbar]
  []
  [mu]
    map_to_aux_variable = mu
  []
  [mubar]
  []
  [Mbarmubar]
  []
  # constant tensors
  [Mbar]
  []
  [kappabarbar]
  []
[]

[TensorComputes]
  [Initialize]
    [c]
      # Random initial condition around a concentration of 1/2
      type = RandomTensor
      buffer = c
      min = 0.44
      max = 0.56
    []

    # precompute fixed factors for the solve
    [Mbar]
      type = ReciprocalLaplacianFactor
      factor = 0.2 # Mobility
      buffer = Mbar
    []
    [kappabarbar]
      type = ReciprocalLaplacianSquareFactor
      factor = -0.001 # kappa
      buffer = kappabarbar
    []
    [mu_init]
      type = ConstantTensor
      buffer = mu
      real = 0
    []
  []

  [Solve]
    [cahn_hilliard]
      [mu]
        type = ParsedCompute
        buffer = mu
        expression = '0.1*c^2*(c-1)^2'
        derivatives = c
        inputs = c
      []
      [mubar]
        type = ForwardFFT
        buffer = mubar
        input = mu
      []
      [Mbarmubar]
        type = ParsedCompute
        buffer = Mbarmubar
        expression = 'Mbar*mubar'
        inputs = 'Mbar mubar'
      []
      [cbar]
        type = ForwardFFT
        buffer = cbar
        input = c
      []
    []
  []
[]

[TensorSolver]
  type = AdamsBashforthMoulton
  root_compute = cahn_hilliard
  buffer = c
  reciprocal_buffer = cbar
  linear_reciprocal = kappabarbar
  nonlinear_reciprocal = Mbarmubar
  substeps = 1000
[]

[AuxVariables]
  [mu]
    # the mu tensor  is projected onto this elemental variable
    family = MONOMIAL
    order = CONSTANT
  []
  [c]
    # the c tensor is projected onto this nodal variable
  []
[]

# a slower but more flexible alternative to `map_to_aux_variable` is running
# these `ProjectTensorAux` AuxKernels to perform the projection. This aprpoach
# also supports non-conforming meshes.
[AuxKernels]
  # [c]
  #   type = ProjectTensorAux
  #   buffer = c
  #   variable = c
  #   execute_on = final
  # []
  # [f]
  #   type = ProjectTensorAux
  #   buffer = f
  #   variable = f
  #   execute_on = TIMESTEP_END
  # []
[]

[Postprocessors]
  [min_c]
    type = ElementExtremeValue
    variable = c
    value_type = MIN
    execute_on = 'TIMESTEP_END'
  []
  [max_c]
    type = ElementExtremeValue
    variable = c
    value_type = MAX
    execute_on = 'TIMESTEP_END'
  []
  # [F]
  #   type = ElementIntegralVariablePostprocessor
  #   variable = f
  #   execute_on = 'TIMESTEP_END'
  # []
  [C]
    type = ElementIntegralVariablePostprocessor
    variable = c
    execute_on = 'TIMESTEP_END'
  []
[]

[Problem]
  type = TensorProblem
[]

[Executioner]
  type = Transient
  num_steps = 100
  [TimeStepper]
    type = IterationAdaptiveDT
    growth_factor = 1.8
    dt = 0.1
  []
  dtmax = 1000
[]

[Outputs]
  exodus = true
  csv = true
  perf_graph = true
  execute_on = 'TIMESTEP_END'
[]

(test/tests/real_space/laplacian3d.i)

[Domain]
  dim = 3
  nx = 30
  ny = 24
  nz = 17
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y Z'
  xmin = -1
  ymin = -1
  zmin = -1
  xmax = 1
  ymax = 1
  zmax = 1
[]

[TensorBuffers]
  [u]
    type = PlainTensorBuffer
  []
  [u_gold]
    type = PlainTensorBuffer
  []
  [lap]
    type = PlainTensorBuffer
  []
[]

[TensorComputes]
  [Initialize]
    [u_ic]
      type = ParsedCompute
      buffer = u
      expression = 'sin(pi*(x+0.5))*sin(pi*(y+0.5))*sin(pi*(z+0.5))'
      extra_symbols = true
    []
    [lap_gold_ic]
      type = ParsedCompute
      buffer = lap_gold
      #expression = '0'
      expression = '-3*pi*pi*sin(pi*(x+0.5))*sin(pi*(y+0.5))*sin(pi*(z+0.5))'
      extra_symbols = true
      expand = REAL
    []
  []
  [Solve]
    [lap_compute]
      type = FiniteDifferenceLaplacian
      input = u
      buffer = lap
    []
    [diff]
      type = ParsedCompute
      buffer = diff
      expression = 'abs(lap - lap_gold)'
      inputs = 'lap lap_gold'
    []
  []
[]

[Postprocessors]
  [max_error]
    type = TensorExtremeValuePostprocessor
    buffer = diff
    value_type = MAX
  []
[]

[TensorOutputs]
  # active = ''
  [lap]
    type = XDMFTensorOutput
    buffer = 'lap lap_gold diff'
    enable_hdf5 = true
  []
[]

[Executioner]
  type = Transient
  num_steps = 1
[]

[Outputs]
  [out]
    type = CSV
  []
[]

(test/tests/solvers/etdrk4_diffusion.i)

D = 0.05
k = 1.0
ss = 1
dt = '${units 10 s }'

[Domain]
  dim = 1
  nx = 64
  xmax = '${fparse 2*pi}'
  mesh_mode = DUMMY
[]

[TensorComputes]
  [Initialize]
    [u0]
      type = ParsedCompute
      buffer = u0
      extra_symbols = true
      expression = 'sin(${k}*x)'
    []
    [u]
      type = ParsedCompute
      buffer = u
      inputs = u0
      expression = 'u0'
    []

    [L]
      type = ReciprocalLaplacianFactor
      factor = ${D}
      buffer = L
    []
    [zero]
      type = ConstantReciprocalTensor
      buffer = zero
    []
  []

  [Solve]
    [u_bar]
      type = ForwardFFT
      buffer = u_bar
      input = u
    []
    [u_exact]
      type = ParsedCompute
      buffer = u_exact
      inputs = u0
      extra_symbols = true
      expression = 'u0*exp(-${D}*${k}^2*t)'
    []
    [u_diff_sq]
      type = ParsedCompute
      buffer = u_diff_sq
      inputs = 'u u_exact'
      expression = '(u - u_exact)^2'
    []
  []
[]

[TensorSolver]
  type = ETDRK4Solver
  buffer = 'u'
  reciprocal_buffer = 'u_bar'
  linear_reciprocal = 'L'
  nonlinear_reciprocal = 'zero'
  substeps = ${ss}
[]

[Problem]
  type = TensorProblem
[]

[Postprocessors]
  [mse]
    type = TensorIntegralPostprocessor
    buffer = u_diff_sq
  []
  [rmse]
    type = ParsedPostprocessor
    expression = 'sqrt(mse)'
    pp_names = 'mse'
    pp_symbols = 'mse'
  []
[]

[Executioner]
  type = Transient
  num_steps = 10
  dt = ${dt}
[]

[Outputs]
  file_base = etdrk4_diffusion_rmse
  csv = true
[]

(test/tests/typed_tensors/gradient_vector.i)

[Domain]
  dim = 3
  nx = 20
  ny = 10
  nz = 5
  mesh_mode = DUMMY
  device_names = cpu
[]

[TensorComputes]
  [Initialize]
    [c]
      type = ParsedCompute
      buffer = c
      extra_symbols = true
      expression = 'sin(x*8*pi)+cos(y*4*pi)+sin(z*2*pi)'
    []
    [grad_c]
      type = GradientVector
      buffer = grad_c
      input = c
    []
  []
[]

[Problem]
  type = TensorProblem
  print_debug_output = true
[]

[Executioner]
  type = Transient
  num_steps = 1
[]

[TensorOutputs]
  [xdmf]
    type = XDMFTensorOutput
    buffer = 'c grad_c'
    file_base = gradient
    output_mode = 'NODE NODE'
    enable_hdf5 = true
  []
[]

(benchmarks/01_spinodal_decomposition/1a.i)

[Domain]
  dim = 2
  nx = 200
  ny = 200
  xmax = 200
  ymax = 200

  device_names = 'cuda'

  mesh_mode = DOMAIN
[]


[TensorBuffers]
  [c]
    map_to_aux_variable = c
  []
  [cbar]
  []
  [mu]
    # map_to_aux_variable = mu
  []
  [mubar]
  []
  [Mbarmubar]
  []
  # constant tensors
  [Mbar]
  []
  [kappabarbar]
  []
  # postprocessing
  [F]
  []
  [Fgrad]
  []
[]

[TensorComputes]
  [Initialize]
    [c]
      type = ParsedCompute
      buffer = c
      extra_symbols = true
      expression = 'c0+epsilon*(cos(0.105*x)*cos(0.11*y)+(cos(0.13*x)*cos(0.087*y))^2+cos(0.025*x-0.15*y)*cos(0.07*x-0.02*y))'
      constant_names = 'c0 epsilon'
      constant_expressions = '0.5 0.01'
    []
    [Mbar]
      type = ReciprocalLaplacianFactor
      factor = 5 # Mobility
      buffer = Mbar
    []
    [kappabarbar]
      type = ReciprocalLaplacianSquareFactor
      factor = -10 # -kappa*M
      buffer = kappabarbar
    []
  []

  [Solve]
    [mu]
      type = ParsedCompute
      buffer = mu
      expression = 'rho_s*(c-c_alpha)^2*(c_beta-c)^2'
      constant_names =       'rho_s c_alpha c_beta'
      constant_expressions = '5     0.3     0.7'
      derivatives = c
      inputs = c
    []
    [mubar]
      type = ForwardFFT
      buffer = mubar
      input = mu
    []
    [Mbarmubar]
      type = ParsedCompute
      buffer = Mbarmubar
      expression = 'Mbar*mubar'
      inputs = 'Mbar mubar'
    []
    [cbar]
      type = ForwardFFT
      buffer = cbar
      input = c
    []
  []

  [Postprocess]
    [Fgrad]
      type = FFTGradientSquare
      buffer = Fgrad
      input = c
      factor = 1 # kappa/2
    []
    [F]
      type = ParsedCompute
      buffer = F
      expression = 'rho_s * (c-c_alpha)^2 * (c_beta-c)^2 + Fgrad'
      constant_names =       'rho_s c_alpha c_beta'
      constant_expressions = '5     0.3     0.7'
      inputs = 'c Fgrad'
    []
  []
[]

[UserObjects]
  [terminator]
    type = Terminator
    expression = change<1e-4
  []
[]

[TensorTimeIntegrators]
  [c]
    type = FFTSemiImplicit
    buffer = c
    history_size = 1
    reciprocal_buffer = cbar
    linear_reciprocal = kappabarbar
    nonlinear_reciprocal = Mbarmubar
  []
[]

[AuxVariables]
  # [mu]
  #   family = MONOMIAL
  #   order = CONSTANT
  # []
  [c]
    # family = MONOMIAL
    # order = CONSTANT
  []
[]

[Postprocessors]
  [min_c]
    type = TensorExtremeValuePostprocessor
    buffer = c
    value_type = MIN
    execute_on = 'TIMESTEP_END'
  []
  [max_c]
    type = TensorExtremeValuePostprocessor
    buffer = c
    value_type = MAX
    execute_on = 'TIMESTEP_END'
  []
  [F]
    type = TensorIntegralPostprocessor
    buffer = F
  []
  [change]
    type = TensorIntegralChangePostprocessor
    buffer = c
  []
[]

[Problem]
  type = TensorProblem
  spectral_solve_substeps = 1000
[]

[Executioner]
  type = Transient
  num_steps = 1000
  [TimeStepper]
    type = IterationAdaptiveDT
    growth_factor = 1.1
    dt = 1
  []
  dtmax = 300
[]

[Outputs]
  exodus = true
  csv = true
  perf_graph = true
  execute_on = 'TIMESTEP_END'
[]

(test/tests/tensor_compute/coupled_pf_mech.i)

[Domain]
  dim = 3
  nx = 128
  ny = 128
  nz = 128
  xmax = ${fparse pi*4}
  ymax = ${fparse pi*4}
  zmax = ${fparse pi*4}
  mesh_mode = DUMMY
[]

[TensorBuffers]
  # phase field
  [c]
  []
  [cbar]
  []
  [mu]
  []
  [mubar]
  []
  [Mbarmubar]
  []

  # mechanics
  [disp_x]
  []
  [disp_y]
  []
  [disp_z]
  []
  [mumechbar]
  []
  [mumech]
  []

  # constant tensors
  [Mbar]
  []
  [kappabarbar]
  []
[]

[TensorOutputs]
  [xdmf]
    type = XDMFTensorOutput
    buffer = 'c disp_x disp_y disp_z mu mumech'
    output_mode = 'Node Node Node Node Cell Cell'
    enable_hdf5 = true
  []
[]

[TensorComputes]
  [Initialize]
    [c]
      type = RandomTensor
      buffer = c
      min = 0.44
      max = 0.56
    []
    [disp_x]
      type = RandomTensor
      buffer = disp_x
      min = 0
      max = 0
    []
    [disp_y]
      type = RandomTensor
      buffer = disp_y
      min = 0
      max = 0
    []
    [disp_z]
      type = RandomTensor
      buffer = disp_z
      min = 0
      max = 0
    []
    [Mbar]
      type = ReciprocalLaplacianFactor
      factor = 0.2 # Mobility
      buffer = Mbar
    []
    [kappabarbar]
      type = ReciprocalLaplacianSquareFactor
      factor = -0.001 # kappa
      buffer = kappabarbar
    []
  []

  [Solve]
    [mu]
      # chemical potential (real space)
      type = ParsedCompute
      buffer = mu
      expression = '0.1*c^2*(c-1)^2' # + c*sin(x/2)*0.005'
      extra_symbols = true
      derivatives = c
      inputs = c
    []
    [mubar]
      # chemical potential (reciprocal space)
      type = ForwardFFT
      buffer = mubar
      input = mu
    []
    [mumechbar]
      # mechanical chemical potential (reciprocal space)
      type = FFTElasticChemicalPotential
      buffer = mumechbar
      cbar = cbar
      displacements = 'disp_x disp_y disp_z'
      lambda = 100
      mu = 50
      e0 = 0.02
    []
    [mumech]
      # chemical potential (reciprocal space)
      type = InverseFFT
      buffer = mumech
      input = mumechbar
    []

    [Mbarmubar]
      type = ParsedCompute
      buffer = Mbarmubar
      expression = 'Mbar*(mubar+mumechbar)'
      inputs = 'Mbar mubar mumechbar'
    []
    [cbar]
      type = ForwardFFT
      buffer = cbar
      input = c
    []

    [qsmech]
      type = FFTQuasistaticElasticity
      displacements = 'disp_x disp_y disp_z'
      cbar = cbar
      lambda = 100
      mu = 50
      e0 = 0.02
    []
  []
[]

[TensorTimeIntegrators]
  [c]
    type = FFTSemiImplicit
    buffer = c
    reciprocal_buffer = cbar
    linear_reciprocal = kappabarbar
    nonlinear_reciprocal = Mbarmubar
  []
[]

[Postprocessors]
  [min_c]
    type = TensorExtremeValuePostprocessor
    buffer = c
    value_type = MIN
    execute_on = 'TIMESTEP_END'
  []
  [max_c]
    type = TensorExtremeValuePostprocessor
    buffer = c
    value_type = MAX
    execute_on = 'TIMESTEP_END'
  []

  [min_disp_x]
    type = TensorExtremeValuePostprocessor
    buffer = disp_x
    value_type = MIN
    execute_on = 'TIMESTEP_END'
  []
  [max_disp_x]
    type = TensorExtremeValuePostprocessor
    buffer = disp_x
    value_type = MAX
    execute_on = 'TIMESTEP_END'
  []
  [min_disp_y]
    type = TensorExtremeValuePostprocessor
    buffer = disp_y
    value_type = MIN
    execute_on = 'TIMESTEP_END'
  []
  [max_disp_y]
    type = TensorExtremeValuePostprocessor
    buffer = disp_y
    value_type = MAX
    execute_on = 'TIMESTEP_END'
  []
  [min_disp_z]
    type = TensorExtremeValuePostprocessor
    buffer = disp_z
    value_type = MIN
    execute_on = 'TIMESTEP_END'
  []
  [max_disp_z]
    type = TensorExtremeValuePostprocessor
    buffer = disp_z
    value_type = MAX
    execute_on = 'TIMESTEP_END'
  []

  [C]
    type = TensorIntegralPostprocessor
    buffer = c
    execute_on = 'TIMESTEP_END'
  []
  [cavg]
    type = TensorAveragePostprocessor
    buffer = c
    execute_on = 'TIMESTEP_END'
  []
[]

[Problem]
  type = TensorProblem
  spectral_solve_substeps = 1000
  print_debug_output = true
[]

[Executioner]
  type = Transient
  num_steps = 100
  [TimeStepper]
    type = IterationAdaptiveDT
    growth_factor = 1.8
    dt = 0.1
  []
  dtmax = 1000
[]

[Outputs]
  csv = true
  perf_graph = true
  execute_on = 'TIMESTEP_END'
[]

(examples/cahn_hilliard/cahnhilliard3.i)

[Domain]
  dim = 3
  nx = 100
  ny = 100
  nz = 100
  xmax = ${fparse pi*4}
  ymax = ${fparse pi*4}
  zmax = ${fparse pi*4}

  device_names = 'cuda'

  mesh_mode = DUMMY
[]

[TensorBuffers]
  [c]
  []
  [cbar]
  []
  [mu]
  []
  [mubar]
  []
  [Mbarmubar]
  []
  # constant tensors
  [Mbar]
  []
  [kappabarbar]
  []
[]

[TensorOutputs]
  [xdmf]
    type = XDMFTensorOutput
    buffer = 'c mu'
    output_mode = 'Node Cell'
    enable_hdf5 = true
  []
[]

[TensorComputes]
  [Initialize]
    [c]
      type = RandomTensor
      buffer = c
      min = 0.44
      max = 0.56
    []
    [Mbar]
      type = ReciprocalLaplacianFactor
      factor = 0.2 # Mobility
      buffer = Mbar
    []
    [kappabarbar]
      type = ReciprocalLaplacianSquareFactor
      factor = -0.001 # kappa
      buffer = kappabarbar
    []
  []

  [Solve]
    [mu]
      type = ParsedCompute
      buffer = mu
      expression = '0.1*c^2*(c-1)^2'
      derivatives = c
      # expression = "0.4*c^3-0.6*c^2+0.2*c"
      inputs = c
    []
    [mubar]
      type = ForwardFFT
      buffer = mubar
      input = mu
    []
    [Mbarmubar]
      type = ParsedCompute
      buffer = Mbarmubar
      expression = 'Mbar*mubar'
      inputs = 'Mbar mubar'
    []
    [cbar]
      type = ForwardFFT
      buffer = cbar
      input = c
    []
  []
[]

[TensorTimeIntegrators]
  [c]
    type = FFTSemiImplicit
    buffer = c
    reciprocal_buffer = cbar
    linear_reciprocal = kappabarbar
    nonlinear_reciprocal = Mbarmubar
  []
[]

[Postprocessors]
  [min_c]
    type = TensorExtremeValuePostprocessor
    buffer = c
    value_type = MIN
    execute_on = 'TIMESTEP_END'
  []
  [max_c]
    type = TensorExtremeValuePostprocessor
    buffer = c
    value_type = MAX
    execute_on = 'TIMESTEP_END'
  []
  [C]
    type = TensorIntegralPostprocessor
    buffer = c
    execute_on = 'TIMESTEP_END'
  []
  [cavg]
    type = TensorAveragePostprocessor
    buffer = c
    execute_on = 'TIMESTEP_END'
  []
[]

[Problem]
  type = TensorProblem
  spectral_solve_substeps = 1000
[]

[Executioner]
  type = Transient
  num_steps = 20
  [TimeStepper]
    type = IterationAdaptiveDT
    growth_factor = 1.8
    dt = 0.1
  []
  dtmax = 1000
[]

[Outputs]
  csv = true
  perf_graph = true
  execute_on = 'TIMESTEP_END'
[]

(test/tests/mechanics/mech.i)

[Domain]
  dim = 2
  nx = 32
  ny = 32
  xmax = ${fparse 2*pi}
  ymax = ${fparse 2*pi}
  zmax = ${fparse 2*pi}
  mesh_mode = DUMMY
  parallel_mode = FFT_AUTO
[]

[TensorComputes]
  [Initialize]
    [phase]
      type = ParsedCompute
      expression = '(cos(x)/2+0.5)^1*(cos(y)/2+0.5)^1*(cos(z)/2+0.5)^1'
      extra_symbols = true
      buffer = phase
    []
    [K]
      type = ParsedCompute
      buffer = K
      expression = '(1-phase)*Ka + phase*Kb'
      inputs = phase
      constant_names = 'Ka Kb'
      constant_expressions = '1 10'
    []
    [mu]
      type = ParsedCompute
      buffer = mu
      expression = '(1-phase)*mua + phase*mub'
      inputs = phase
      constant_names = 'mua mub'
      constant_expressions = '0.5 5'
    []
    [Finit]
      type = RankTwoIdentity
      buffer = F
    []
  []

  [Solve]
    [hyper_elasticity]
      type = HyperElasticIsotropic
      buffer = stress
      F = Fnew
      K = K
      mu = mu
    []

    [root]
      [applied_strain]
        type = MacroscopicShearTensor
        buffer = applied_strain
      []
      [mech]
        type = FFTMechanics
        buffer = Fnew
        F = F
        K = K
        mu = mu
        l_max_its = 40
        l_tol = 1e-5
        nl_rel_tol = 2e-4
        nl_abs_tol = 2e-3
        constitutive_model = hyper_elasticity
        stress = stress
        applied_macroscopic_strain = applied_strain
      []
    []
  []

  [Postprocess]
    [displacements]
      type = ComputeDisplacements
      buffer = disp
      F = F
    []
    [vonmises]
      type = ComputeVonMisesStress
      buffer = sV
    []
  []
[]

[TensorSolver]
  # no variables are integrated by this solver (FFTMechanics performs a steady state mechanics solve)
  type = ForwardEulerSolver
  root_compute = root
  # deformation tensor is just forwarded Fnew -> F
  forward_buffer = F
  forward_buffer_new = Fnew
  substeps = 3
[]

[TensorOutputs]
  active = 'parallel' # serial

  [serial]
    type = XDMFTensorOutput
    file_base = mech_serial
    buffer = 'disp sV F phase'
    output_mode = 'OVERSIZED_NODAL CELL CELL NODE'
    enable_hdf5 = true
    execute_on = 'TIMESTEP_END'
  []
  [parallel]
    # NODE and OVERSIZED_NODAL are not yet available in parallel output
    type = XDMFTensorOutput
    file_base = mech_parallel
    buffer = 'sV F phase'
    output_mode = 'CELL CELL CELL'
    enable_hdf5 = true
    execute_on = 'TIMESTEP_END'
  []
[]


[Executioner]
  type = Transient
  num_steps = 3
  dt = 0.02
[]

(benchmarks/02_oswald_ripening/2a.i)

[Domain]
  dim = 2
  nx = 200
  ny = 200
  xmax = 200
  ymax = 200
  mesh_mode = DOMAIN
[]

fchem = 'fa:=rho^2*(c-ca)^2;
fb:=rho^2*(cb-c)^2;
h:=n1^3*(6*n1^2-15*n1+10) +
   n2^3*(6*n2^2-15*n2+10) +
   n3^3*(6*n3^2-15*n3+10) +
   n4^3*(6*n4^2-15*n4+10);
g:=n1^2*(1-n1)^2 +
   n2^2*(1-n2)^2 +
   n3^2*(1-n3)^2 +
   n4^2*(1-n4)^2 +
alpha*(
n1^2*n2^2 + n1^2*n3^2 + n1^2*n4^2 +
n2^2*n1^2 + n2^2*n3^2 + n2^2*n4^2 +
n3^2*n1^2 + n3^2*n2^2 + n3^2*n4^2 +
n4^2*n1^2 + n4^2*n2^2 + n4^2*n3^2);
(fa*(1-h) + fb*h + w*g)'

nic = 'epsilon*(cos((0.01*idx)*x-4)*cos((0.007+0.01*idx)*y)
       +cos((0.11+0.01*idx)*x)*cos((0.11+0.01*idx)*y)
       +psi*(cos((0.046+0.001*idx)*x+(0.0405+0.001*idx)*y)
       *cos((0.031+0.001*idx)*x-(0.004+0.001*idx)*y))^2)^2'

cnames = 'rho     ca  cb  alpha w L M'
cvalues = 'sqrt(2) 0.3 0.7 5     1 5 5'

[TensorBuffers]
  # variables
  [c]
    # map_to_aux_variable = c
  []
  [n1]
  []
  [n2]
  []
  [n3]
  []
  [n4]
  []

  [c_bar]
  []
  [n1_bar]
  []
  [n2_bar]
  []
  [n3_bar]
  []
  [n4_bar]
  []

  [mu_c]
    # map_to_aux_variable = mu
  []
  [mu_n1]
  []
  [mu_n2]
  []
  [mu_n3]
  []
  [mu_n4]
  []

  [mu_c_bar]
  []
  [mu_n1_bar]
  []
  [mu_n2_bar]
  []
  [mu_n3_bar]
  []
  [mu_n4_bar]
  []

  [Mbar_mu_c_bar]
  []

  # constant tensors
  [Lbar] # FFT(M*laplacian)
  []
  [MkappaL2bar] # FFT(-M*kappa*laplacian^2)
  []
  [kappaLbar] # FFT(L*kappa*laplacian)
  []

  # postprocessing
  [F]
  []
  [Fgrad_c]
  []
  [Fgrad_n1]
  []
  [Fgrad_n2]
  []
  [Fgrad_n3]
  []
  [Fgrad_n4]
  []
  [bnds]
    #map_to_aux_variable = bnds
  []
[]

[TensorComputes]
  [Initialize]
    [c]
      type = ParsedCompute
      buffer = c
      extra_symbols = true
      expression = 'c0+epsilon*(cos(0.105*x)*cos(0.11*y)+(cos(0.13*x)*cos(0.087*y))^2+cos(0.025*x-0.15*y)*cos(0.07*x-0.02*y))'
      constant_names = 'c0 epsilon'
      constant_expressions = '0.5 0.01'
    []
    [Lbar]
      type = ReciprocalLaplacianFactor
      # Mobility is pulled into the chemical potential below
      buffer = Lbar
    []
    [MkappaL2bar]
      type = ReciprocalLaplacianSquareFactor
      factor = -15 # -kappa_c*M
      buffer = MkappaL2bar
    []
    [kappaLbar]
      type = ReciprocalLaplacianFactor
      buffer = kappaLbar
      factor = 15 # kappa_ni*L
    []
    [n1]
      type = ParsedCompute
      buffer = n1
      expression = ${nic}
      extra_symbols = true
      constant_names = 'idx epsilon psi'
      constant_expressions = '  1     0.1 1.5'
    []
    [n2]
      type = ParsedCompute
      buffer = n2
      expression = ${nic}
      extra_symbols = true
      constant_names = 'idx epsilon psi'
      constant_expressions = '  2    0.1 1.5'
    []
    [n3]
      type = ParsedCompute
      buffer = n3
      expression = ${nic}
      extra_symbols = true
      constant_names = 'idx epsilon psi'
      constant_expressions = '  3     0.1 1.5'
    []
    [n4]
      type = ParsedCompute
      buffer = n4
      expression = ${nic}
      extra_symbols = true
      constant_names = 'idx epsilon psi'
      constant_expressions = '  4     0.1 1.5'
    []
  []

  [Solve]
    [mu_c]
      type = ParsedCompute
      buffer = mu_c
      expression = '${fchem}*M'
      constant_names = ${cnames}
      constant_expressions = ${cvalues}
      derivatives = c
      inputs = 'c n1 n2 n3 n4'
    []

    [mu_n1]
      type = ParsedCompute
      buffer = mu_n1
      expression = '${fchem}*(-L)'
      constant_names = ${cnames}
      constant_expressions = ${cvalues}
      derivatives = n1
      inputs = 'c n1 n2 n3 n4'
    []
    [mu_n2]
      type = ParsedCompute
      buffer = mu_n2
      expression = '${fchem}*(-L)'
      constant_names = ${cnames}
      constant_expressions = ${cvalues}
      derivatives = n2
      inputs = 'c n1 n2 n3 n4'
    []
    [mu_n3]
      type = ParsedCompute
      buffer = mu_n3
      expression = '${fchem}*(-L)'
      constant_names = ${cnames}
      constant_expressions = ${cvalues}
      derivatives = n3
      inputs = 'c n1 n2 n3 n4'
    []
    [mu_n4]
      type = ParsedCompute
      buffer = mu_n4
      expression = '${fchem}*(-L)'
      constant_names = ${cnames}
      constant_expressions = ${cvalues}
      derivatives = n4
      inputs = 'c n1 n2 n3 n4'
    []

    [mu_c_bar]
      type = ForwardFFT
      buffer = mu_c_bar
      input = mu_c
    []
    [mu_n1_bar]
      type = ForwardFFT
      buffer = mu_n1_bar
      input = mu_n1
    []
    [mu_n2_bar]
      type = ForwardFFT
      buffer = mu_n2_bar
      input = mu_n2
    []
    [mu_n3_bar]
      type = ForwardFFT
      buffer = mu_n3_bar
      input = mu_n3
    []
    [mu_n4_bar]
      type = ForwardFFT
      buffer = mu_n4_bar
      input = mu_n4
    []

    [Mbar_mu_c_bar]
      type = ParsedCompute
      buffer = Mbar_mu_c_bar
      expression = 'Lbar*mu_c_bar'
      inputs = 'Lbar mu_c_bar'
    []

    [c_bar]
      type = ForwardFFT
      buffer = c_bar
      input = c
    []
    [n1_bar]
      type = ForwardFFT
      buffer = n1_bar
      input = n1
    []
    [n2_bar]
      type = ForwardFFT
      buffer = n2_bar
      input = n2
    []
    [n3_bar]
      type = ForwardFFT
      buffer = n3_bar
      input = n3
    []
    [n4_bar]
      type = ForwardFFT
      buffer = n4_bar
      input = n4
    []
  []

  [Postprocess]
    [Fgrad_c]
      type = FFTGradientSquare
      buffer = Fgrad_c
      input = c
      factor = 1.5 # kappa/2
    []
    [Fgrad_n1]
      type = FFTGradientSquare
      buffer = Fgrad_n1
      input = n1
      factor = 1.5 # kappa/2
    []
    [Fgrad_n2]
      type = FFTGradientSquare
      buffer = Fgrad_n2
      input = n2
      factor = 1.5 # kappa/2
    []
    [Fgrad_n3]
      type = FFTGradientSquare
      buffer = Fgrad_n3
      input = n3
      factor = 1.5 # kappa/2
    []
    [Fgrad_n4]
      type = FFTGradientSquare
      buffer = Fgrad_n4
      input = n4
      factor = 1.5 # kappa/2
    []
    [F]
      type = ParsedCompute
      buffer = F
      expression = '${fchem} + Fgrad_c + Fgrad_n1 + Fgrad_n2 + Fgrad_n3 + Fgrad_n4'
      constant_names = ${cnames}
      constant_expressions = ${cvalues}
      inputs = 'c n1 n2 n3 n4 Fgrad_c Fgrad_n1 Fgrad_n2 Fgrad_n3 Fgrad_n4'
    []
    [bnds]
      type = ParsedCompute
      buffer = bnds
      expression = 'n1^2 + n2^2 + n3^2 + n4^2'
      inputs = 'n1 n2 n3 n4'
    []
  []
[]

[TensorSolver]
  type = AdamsBashforthMoulton
  buffer = 'c n1 n2 n3 n4'
  reciprocal_buffer = 'c_bar n1_bar n2_bar n3_bar n4_bar'
  linear_reciprocal = 'MkappaL2bar kappaLbar kappaLbar kappaLbar kappaLbar'
  nonlinear_reciprocal = 'Mbar_mu_c_bar mu_n1_bar mu_n2_bar mu_n3_bar mu_n4_bar'
  substeps = 2000
  predictor_order = 2
  corrector_order = 2
  corrector_steps = 0
[]

[AuxVariables]
  [c]
  []
  [bnds]
  []
[]

[Postprocessors]
  [min_c]
    type = TensorExtremeValuePostprocessor
    buffer = c
    value_type = MIN
    execute_on = 'TIMESTEP_END'
  []
  [max_c]
    type = TensorExtremeValuePostprocessor
    buffer = c
    value_type = MAX
    execute_on = 'TIMESTEP_END'
  []
  [F]
    type = TensorIntegralPostprocessor
    buffer = F
  []
  # [stable_dt]
  #   type = SemiImplicitCriticalTimeStep
  #   buffer = MkappaL2bar
  # []
[]

[Problem]
  type = TensorProblem
[]

[Executioner]
  type = Transient
  num_steps = 1030
  [TimeStepper]
    type = IterationAdaptiveDT
    growth_factor = 1.1
    dt = 0.001
  []
  dtmax = 10
[]

[Outputs]
  csv = true
  perf_graph = true
  execute_on = 'TIMESTEP_END'
[]

[TensorOutputs]
  [xdmf]
    type = XDMFTensorOutput
    buffer = 'c bnds'
    output_mode = 'CELL CELL'
  []
[]

(test/tests/typed_tensors/gradient.i)

[Domain]
  dim = 3
  nx = 20
  ny = 10
  nz = 5
  mesh_mode = DUMMY
  device_names = cpu
[]

[TensorComputes]
  [Initialize]
    [c]
      type = ParsedCompute
      buffer = c
      extra_symbols = true
      expression = 'sin(x*8*pi)+cos(y*4*pi)+sin(z*2*pi)'
    []
    [grad_c]
      type = NEML2GradientVector
      buffer = grad_c
      input = c
    []
  []
[]

[Problem]
  type = TensorProblem
  print_debug_output = true
[]

[Executioner]
  type = Transient
  num_steps = 1
[]

[TensorOutputs]
  [xdmf]
    type = XDMFTensorOutput
    buffer = 'c grad_c'
    file_base = gradient
    output_mode = 'NODE NODE'
    enable_hdf5 = true
  []
[]

(test/tests/parsed_tensor/local_vars_derivative.i)

#
# Test differentiation of local variables in ParsedCompute
#
# We test the expression with a buffer input: r := sqrt(a^2 + 1); r^2
# The derivative with respect to a should be: d/da(r^2) = 2*r * dr/da = 2*r * a/r = 2*a
#

[Domain]
  dim = 2
  nx = 20
  ny = 20
  xmax = 2
  ymax = 2
  mesh_mode = DUMMY
[]

[TensorBuffers]
  [a]        # Input variable
  []
  [df_da]    # Auto-differentiated derivative
  []
  [df_da_exact]  # Hand-coded exact derivative
  []
  [error]    # Absolute difference
  []
[]

[TensorComputes]
  [Solve]
    # Initialize input buffer
    [init_a]
      type = ParsedCompute
      buffer = a
      expression = 'x + 0.5*y'
      extra_symbols = true
    []

    # Auto-differentiated derivative: d/da(r^2) where r:=sqrt(a^2+1)
    [auto_derivative]
      type = ParsedCompute
      buffer = df_da
      expression = 'r:=sqrt(a^2+1); r^2'
      derivatives = 'a'
      inputs = 'a'
    []

    # Hand-coded exact derivative: d/da(r^2) = 2*a
    [exact_derivative]
      type = ParsedCompute
      buffer = df_da_exact
      expression = '2*a'
      inputs = 'a'
    []

    # Compute absolute error
    [compute_error]
      type = ParsedCompute
      buffer = error
      expression = 'abs(df_da - df_da_exact)'
      inputs = 'df_da df_da_exact'
    []
  []
[]

[Postprocessors]
  [integral_error]
    type = TensorIntegralPostprocessor
    buffer = error
    execute_on = 'INITIAL'
  []
[]

[Problem]
  type = TensorProblem
[]

[Executioner]
  type = Transient
  num_steps = 0
[]

[Outputs]
  csv = true
[]

(benchmarks/01_spinodal_decomposition/1a_solver.i)

[Domain]
  dim = 2
  nx = 200
  ny = 200
  xmax = 200
  ymax = 200

  device_names = 'cuda'

  mesh_mode = DOMAIN
[]


[TensorBuffers]
  [c]
    map_to_aux_variable = c
  []
  [cbar]
  []
  [mu]
    # map_to_aux_variable = mu
  []
  [mubar]
  []
  [Mbarmubar]
  []
  # constant tensors
  [Mbar]
  []
  [kappabarbar]
  []
  # postprocessing
  [F]
  []
  [Fgrad]
  []
[]

[TensorComputes]
  [Initialize]
    [c]
      type = ParsedCompute
      buffer = c
      extra_symbols = true
      expression = 'c0+epsilon*(cos(0.105*x)*cos(0.11*y)+(cos(0.13*x)*cos(0.087*y))^2+cos(0.025*x-0.15*y)*cos(0.07*x-0.02*y))'
      constant_names = 'c0 epsilon'
      constant_expressions = '0.5 0.01'
    []
    [Mbar]
      type = ReciprocalLaplacianFactor
      factor = 5 # Mobility
      buffer = Mbar
    []
    [kappabarbar]
      type = ReciprocalLaplacianSquareFactor
      factor = -10 # -kappa*M
      buffer = kappabarbar
    []
  []

  [Solve]
    [mu]
      type = ParsedCompute
      buffer = mu
      expression = 'rho_s*(c-c_alpha)^2*(c_beta-c)^2'
      constant_names =       'rho_s c_alpha c_beta'
      constant_expressions = '5     0.3     0.7'
      derivatives = c
      inputs = c
    []
    [mubar]
      type = ForwardFFT
      buffer = mubar
      input = mu
    []
    [Mbarmubar]
      type = ParsedCompute
      buffer = Mbarmubar
      expression = 'Mbar*mubar'
      inputs = 'Mbar mubar'
    []
    [cbar]
      type = ForwardFFT
      buffer = cbar
      input = c
    []
  []

  [Postprocess]
    [Fgrad]
      type = FFTGradientSquare
      buffer = Fgrad
      input = c
      factor = 1 # kappa/2
    []
    [F]
      type = ParsedCompute
      buffer = F
      expression = 'rho_s * (c-c_alpha)^2 * (c_beta-c)^2 + Fgrad'
      constant_names =       'rho_s c_alpha c_beta'
      constant_expressions = '5     0.3     0.7'
      inputs = 'c Fgrad'
    []
  []
[]

[UserObjects]
  [terminator]
    type = Terminator
    expression = change<1e-4
  []
[]

[TensorSolver]
  type = AdamsBashforthMoulton
  buffer = c
  substeps = 1000
  history_size = 1
  reciprocal_buffer = cbar
  linear_reciprocal = kappabarbar
  nonlinear_reciprocal = Mbarmubar
[]

[AuxVariables]
  # [mu]
  #   family = MONOMIAL
  #   order = CONSTANT
  # []
  [c]
    # family = MONOMIAL
    # order = CONSTANT
  []
[]

[Postprocessors]
  [min_c]
    type = TensorExtremeValuePostprocessor
    buffer = c
    value_type = MIN
    execute_on = 'TIMESTEP_END'
  []
  [max_c]
    type = TensorExtremeValuePostprocessor
    buffer = c
    value_type = MAX
    execute_on = 'TIMESTEP_END'
  []
  [F]
    type = TensorIntegralPostprocessor
    buffer = F
  []
  [change]
    type = TensorIntegralChangePostprocessor
    buffer = c
  []
[]

[Problem]
  type = TensorProblem
  spectral_solve_substeps = 1000
[]

[Executioner]
  type = Transient
  num_steps = 1000
  [TimeStepper]
    type = IterationAdaptiveDT
    growth_factor = 1.1
    dt = 1
  []
  dtmax = 300
[]

[Outputs]
  exodus = true
  csv = true
  perf_graph = true
  execute_on = 'TIMESTEP_END'
[]

(test/tests/gradient/gradient.i)

[Domain]
  dim = 3
  nx = 40
  ny = 40
  nz = 40
  xmax = ${fparse pi*2}
  ymax = ${fparse pi*4}
  zmax = ${fparse pi*6}
  mesh_mode = DUMMY
[]

[TensorBuffers]
  [s]
  []
  [gradx_s]
  []
  [grady_s]
  []
  [gradz_s]
  []
  [cx]
  []
  [cy]
  []
  [cz]
  []
  [diff]
  []
[]

[TensorComputes]
  [Initialize]
    [sin]
      type = ParsedCompute
      buffer = s
      extra_symbols = true
      expression = 'sin(x)+sin(y)+sin(z)'
    []
    [cosx]
      type = ParsedCompute
      buffer = cx
      extra_symbols = true
      expression = 'cos(x)'
    []
    [cosy]
      type = ParsedCompute
      buffer = cy
      extra_symbols = true
      expression = 'cos(y)'
    []
    [cosz]
      type = ParsedCompute
      buffer = cz
      extra_symbols = true
      expression = 'cos(z)'
    []
    [gradx_sin]
      type = FFTGradient
      buffer = gradx_s
      input = s
      direction = x
    []
    [grady_sin]
      type = FFTGradient
      buffer = grady_s
      input = s
      direction = y
    []
    [gradz_sin]
      type = FFTGradient
      buffer = gradz_s
      input = s
      direction = z
    []
    [diff]
      type = ParsedCompute
      buffer = diff
      inputs = 'gradx_s grady_s gradz_s cx cy  cz'
      expression = 'abs(gradx_s - cx)+abs(grady_s - cy)+abs(gradz_s - cz)'
    []
  []
[]

[Postprocessors]
  [diff]
    type = TensorIntegralPostprocessor
    buffer = diff
  []
[]

[Problem]
  type = TensorProblem
[]

[Executioner]
  type = Transient
  num_steps = 1
[]

[Outputs]
  csv = true
[]

(benchmarks/01_spinodal_decomposition/1a_secant.i)

[Domain]
  dim = 2
  nx = 200
  ny = 200
  xmax = 200
  ymax = 200
  device_names = 'cuda'
  mesh_mode = DOMAIN
[]


[TensorBuffers]
  [c]
    map_to_aux_variable = c
  []
  [cbar]
  []
  [mu]
  []
  [mubar]
  []
  [Mbarmubar]
  []
  # constant tensors
  [Mbar]
  []
  [kappabarbar]
  []
  # postprocessing
  [F]
  []
  [Fgrad]
  []
[]

[TensorComputes]
  [Initialize]
    [c]
      type = ParsedCompute
      buffer = c
      extra_symbols = true
      expression = 'c0+epsilon*(cos(0.105*x)*cos(0.11*y)+(cos(0.13*x)*cos(0.087*y))^2+cos(0.025*x-0.15*y)*cos(0.07*x-0.02*y))'
      constant_names = 'c0 epsilon'
      constant_expressions = '0.5 0.01'
    []
    [Mbar]
      type = ReciprocalLaplacianFactor
      factor = 5 # Mobility
      buffer = Mbar
    []
    [kappabarbar]
      type = ReciprocalLaplacianSquareFactor
      factor = -10 # -kappa*M
      buffer = kappabarbar
    []
  []

  [Solve]
    [mu]
      type = ParsedCompute
      buffer = mu
      expression = 'rho_s*(c-c_alpha)^2*(c_beta-c)^2'
      constant_names =       'rho_s c_alpha c_beta'
      constant_expressions = '5     0.3     0.7'
      derivatives = c
      inputs = c
    []
    [mubar]
      type = ForwardFFT
      buffer = mubar
      input = mu
    []
    [Mbarmubar]
      type = ParsedCompute
      buffer = Mbarmubar
      expression = 'Mbar*mubar'
      inputs = 'Mbar mubar'
    []
    [cbar]
      type = ForwardFFT
      buffer = cbar
      input = c
    []
  []

  [Postprocess]
    [Fgrad]
      type = FFTGradientSquare
      buffer = Fgrad
      input = c
      factor = 1 # kappa/2
    []
    [F]
      type = ParsedCompute
      buffer = F
      expression = 'rho_s * (c-c_alpha)^2 * (c_beta-c)^2 + Fgrad'
      constant_names =       'rho_s c_alpha c_beta'
      constant_expressions = '5     0.3     0.7'
      inputs = 'c Fgrad'
    []
  []
[]

[TensorSolver]
  type = SecantSolver
  substeps = 5
  max_iterations = 50
  tolerance = 1e-7
  dt_epsilon = 1e-5
  buffer = c
  reciprocal_buffer = cbar
  linear_reciprocal = kappabarbar
  nonlinear_reciprocal = Mbarmubar
[]

[AuxVariables]
  [c]
  []
[]

[Postprocessors]
  [min_c]
    type = TensorExtremeValuePostprocessor
    buffer = c
    value_type = MIN
    execute_on = 'TIMESTEP_END'
  []
  [max_c]
    type = TensorExtremeValuePostprocessor
    buffer = c
    value_type = MAX
    execute_on = 'TIMESTEP_END'
  []
  [F]
    type = TensorIntegralPostprocessor
    buffer = F
  []
[]

[Problem]
  type = TensorProblem
[]

[Executioner]
  type = Transient
  num_steps = 900
  [TimeStepper]
    type = IterationAdaptiveDT
    growth_factor = 1.01
    dt = 1
  []
  # dtmax = 200
[]

[Outputs]
  exodus = true
  csv = true
  perf_graph = true
  execute_on = 'TIMESTEP_END'
[]

(test/tests/real_space/diffusion.i)

[Domain]
  dim = 2
  nx = 80
  ny = 60
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y'
  xmin = -1
  ymin = -1
  xmax = 1
  ymax = 1
[]

[TensorComputes]
  [Initialize]
    [u_ic]
      type = ParsedCompute
      buffer = u
      expression = 'exp((-x^2-y^2)*100)'
      extra_symbols = true
    []
  []
  [Solve]
    [dt]
      type = FiniteDifferenceLaplacian
      input = u
      buffer = dudt
    []
  []
[]

[TensorSolver]
  type = RealSpaceForwardEuler
  buffer = u
  time_derivative = dudt
  substeps = 100
[]

[Postprocessors]
  [max_u]
    type = TensorExtremeValuePostprocessor
    buffer = u
    value_type = MAX
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [U]
    type = TensorIntegralPostprocessor
    buffer = u
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]

[TensorOutputs]
  # active = ''
  [u]
    type = XDMFTensorOutput
    buffer = 'u dudt'
    enable_hdf5 = true
  []
[]

[Executioner]
  type = Transient
  dt = 1e-2
  num_steps = 100
[]

[Outputs]
  [out]
    type = CSV
  []
[]

(test/tests/tensor_compute/backandforth.i)

[Domain]
  xmax = ${fparse pi*4}
  ymax = ${fparse pi*4}
  mesh_mode = DUMMY
[]

[TensorBuffers]
  [eta_gold]
  []
  [eta]
  []
  [eta_bar]
  []
  [eta2]
  []
  [zero]
  []
  [diff]
  []
[]

[TensorComputes]
  [Initialize]
    [eta_gold]
      type = ParsedCompute
      buffer = eta_gold
      expression = 'sin(x)+sin(y)+sin(z)'
      extra_symbols = true
    []
    [eta]
      type = ParsedCompute
      buffer = eta
      expression = eta_gold
      inputs = eta_gold
    []
    [eta2]
      type = ConstantTensor
      buffer = eta2
      real = 1
    []
    [zero]
      type = ConstantReciprocalTensor
      buffer = zero
      real = 0
      imaginary = 0
    []
  []

  [Solve]
    [eta_bar]
      type = ForwardFFT
      buffer = eta_bar
      input = eta
    []
    [eta_2]
      type = InverseFFT
      buffer = eta2
      input = eta_bar
    []
  []

  [Postprocess]
    [diff]
      type = ParsedCompute
      buffer = diff
      expression = 'abs(eta - eta2) + abs(eta - eta_gold)'
      inputs = 'eta eta2 eta_gold'
    []
  []
[]

[Postprocessors]
  [norm]
    type = TensorIntegralPostprocessor
    buffer = diff
  []
[]

[TensorSolver]
  type = AdamsBashforthMoulton
  buffer = eta
  reciprocal_buffer = eta_bar
  linear_reciprocal = zero
  nonlinear_reciprocal = zero
[]

[Problem]
  type = TensorProblem
[]

[Executioner]
  type = Transient
  num_steps = 4
[]

[Outputs]
  csv = true
[]

(test/tests/kks/KKS_no_flux_bc.i)

#
# Kim-Kim-Suzuki with no-flux BC imposed using the smooth boundary method (SBM), solved on a 2D grid.
# Mask tensor 'psi' supplies the mask for the solve region to the system.
# Note: c is not directly conserved here - the masked value (psi > 0.0)*c will however be conserved.
#


# Constants for Initial Conditions
r = 30
l = 4.2

# Initial condition function for order parameter
eta_IC = '0.5*(1-tanh(2*(sqrt(x^2+y^2)-${r})/${l}))'

# Phase-field model parameters
kappa_eta = 5
rho_sq = 2
w = 1
M = 5
L = 5
c0_a = 0.3
c0_b = 0.7

# Expressions for switching function and bulk Gibbs energy
h_eta = 'eta^3*(6*eta^2-15*eta+10)'
F = '${h_eta}*(${rho_sq}*((c - (1-${h_eta})*(${c0_b} - ${c0_a}))-${c0_a})^2) + (1-${h_eta})*(${rho_sq}*((c + (${h_eta})*(${c0_b} - ${c0_a}))-${c0_b})^2 ) + ${w}*(eta^2)*(1-eta)^2'


[Domain]
  dim = 2
  nx = 20
  ny = 20

  xmin = -50
  xmax = 50
  ymin = -50
  ymax = 50

  # run on a CUDA device (adjust this to `cpu` if not available)
  device_names = 'cpu'

  # automatically create a matching mesh
  mesh_mode = DUMMY
[]

[Functions]
  [psi_func]
    type = ParsedFunction
    expression = 'if(x<x_min-${l},0,if(x>x_min+${l},1,0.5-0.5*cos(pi*(x-(x_min-${l}))/2/${l}) )) * if(x<x_max-${l},1,if(x>x_max+${l},0,0.5+0.5*cos(pi*(x-(x_max-${l}))/2/${l}) ))'
    symbol_names = 'x_min x_max y_min y_max'
    symbol_values = '30 70 0 100'
  []
[]

[TensorComputes]
  [Initialize]
    [c_IC]
      type = ParsedCompute
      buffer = c
      expression = '0.6 + (${c0_a}-0.6)*${eta_IC}'
      extra_symbols = 'true'
    []
    [eta_IC]
      type = ParsedCompute
      buffer = eta
      expression = '${eta_IC}'
      extra_symbols = 'true'
    []
    [psi_init]
       type = MooseFunctionTensor
      function = psi_func
      buffer = psi
    []
    [zero]
      type = ConstantReciprocalTensor
      buffer = zero
    []
    [M]
      type = ConstantTensor
      buffer = M
      real = ${M}
    []
    [L]
      type = ConstantTensor
      buffer = L
      real = ${L}
    []
    [L_kappa]
      type = ConstantTensor
      buffer = L_kappa
      real = ${fparse  ${L} * ${kappa_eta} }
    []
  []
  [Solve]
    [cbar]
      type = ForwardFFT
      buffer = cbar
      input = c
    []
    [etabar]
      type = ForwardFFT
      buffer = etabar
      input = eta
    []
    [mu]
      type = ParsedCompute
      buffer = 'mu'
      expression = '${F}'
      inputs = 'c eta'
      derivatives = 'c'
    []
    [div_J]
      type = ReciprocalMatDiffusion
      buffer = 'div_J'
      chemical_potential = mu
      mobility = M
      psi = psi
    []
    [domega_chem_deta]
      type = ParsedCompute
      buffer = 'domega_chem_deta'
      expression = '${F} - mu*c'
      inputs = 'mu c eta'
      derivatives = 'eta'
    []
    [AC_bulk]
      type = ReciprocalAllenCahn
      buffer = AC_bulk
      dF_chem_deta = domega_chem_deta
      L = L
      psi = psi
    []
    [kappa_grad_eta]
      type = ReciprocalMatDiffusion
      buffer = 'kappa_grad_eta'
      chemical_potential = 'eta'
      mobility = 'L_kappa'
      psi = psi
    []
    [AC_bar]
      type = ParsedCompute
      buffer = AC_bar
      expression = 'kappa_grad_eta + AC_bulk'
      inputs = 'AC_bulk kappa_grad_eta'

    []
  []
[]

[TensorSolver]
  type = AdamsBashforthMoulton
  buffer = 'c eta'
  reciprocal_buffer = 'cbar etabar'
  linear_reciprocal = 'zero zero'
  nonlinear_reciprocal = 'div_J AC_bar'
  substeps = 1e3
  predictor_order = 3
[]

[Postprocessors]
  [total_C]
    type = TensorIntegralPostprocessor
    buffer = c
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [total_eta]
    type = TensorIntegralPostprocessor
    buffer = eta
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]

[TensorOutputs]
  [xdmf]
    type = XDMFTensorOutput
    buffer = 'eta c mu psi'
    enable_hdf5 = true
    transpose = false
  []
[]

[Executioner]
  type = Transient
  dt = 0.1
  num_steps = 10
[]

[Outputs]
  csv = true
  perf_graph = true
  execute_on = 'INITIAL TIMESTEP_END'
[]

(test/tests/solvers/diagonal.i)

#
# Simple Cahn-Hilliard solve on a 2D grid.
#

[Domain]
  dim = 2
  nx = 150
  ny = 150
  xmax = '${fparse pi*2}'
  ymax = '${fparse pi*2}'
  mesh_mode = DUMMY
[]

[GlobalParams]
  constant_names = 'A B'
  constant_expressions = '1 3.5'
[]

[TensorComputes]
  [Initialize]
    [u]
      type = ParsedCompute
      buffer = u
      extra_symbols = true
      expression = 'sin(x)*sin(y)'
      expand = REAL
    []
    [v]
      type = ConstantTensor
      buffer = v
      real = 0
    []

    # precompute fixed factors for the solve
    [Du]
      type = ReciprocalLaplacianFactor
      factor = 1e-2
      buffer = Du
    []
    [Dv]
      type = ReciprocalLaplacianFactor
      factor = 1e-3
      buffer = Dv
    []
  []

  [Solve]
    [u_bar]
      type = ForwardFFT
      buffer = u_bar
      input = u
    []
    [v_bar]
      type = ForwardFFT
      buffer = v_bar
      input = v
    []

    [source_u]
      type = ParsedCompute
      buffer = source_u
      expression = 'A - (B+1)*u +u^2*v'
      inputs = 'u v'
    []
    [source_u_bar]
      type = ForwardFFT
      buffer = source_u_bar
      input = source_u
    []

    [source_v]
      type = ParsedCompute
      buffer = source_v
      expression = 'B*u - u^2*v'
      inputs = 'u v'
    []
    [source_v_bar]
      type = ForwardFFT
      buffer = source_v_bar
      input = source_v
    []
  []
[]

[TensorSolver]
  type = AdamsBashforthMoulton
  buffer = 'u v'
  reciprocal_buffer = 'u_bar v_bar'
  linear_reciprocal = 'Du Dv'
  nonlinear_reciprocal = 'source_u_bar source_v_bar'
  substeps = ${ss}
  corrector_steps = ${cs}
  predictor_order = ${order}
  corrector_order = ${order}
[]

[Problem]
  type = TensorProblem
[]

[Postprocessors]
  [u_min]
    type = TensorExtremeValuePostprocessor
    buffer = u
    value_type = MIN
  []
  [u_max]
    type = TensorExtremeValuePostprocessor
    buffer = u
    value_type = MAX
  []
  [v_min]
    type = TensorExtremeValuePostprocessor
    buffer = v
    value_type = MIN
  []
  [v_max]
    type = TensorExtremeValuePostprocessor
    buffer = v
    value_type = MAX
  []
  [U]
    type = TensorIntegralPostprocessor
    buffer = u
  []
  [V]
    type = TensorIntegralPostprocessor
    buffer = v
  []
[]

[Executioner]
  type = Transient
  num_steps = 25
  dt = 0.5
[]

[Outputs]
  file_base = diagonal_${ss}_${cs}_${order}
  csv = true
[]

(test/tests/mechanics/mech3d.i)

[Domain]
  dim = 3
  nx = 16
  ny = 16
  nz = 16
  xmax = ${fparse 2*pi}
  ymax = ${fparse 2*pi}
  zmax = ${fparse 2*pi}
  mesh_mode = DUMMY
  parallel_mode = FFT_AUTO
[]

[TensorComputes]
  [Initialize]
    [phase]
      type = ParsedCompute
      expression = '(cos(x)/2+0.5)^1*(cos(y)/2+0.5)^1*(cos(z)/2+0.5)^1'
      extra_symbols = true
      buffer = phase
    []
    [K]
      type = ParsedCompute
      buffer = K
      expression = '(1-phase)*Ka + phase*Kb'
      inputs = phase
      constant_names = 'Ka Kb'
      constant_expressions = '1 10'
    []
    [mu]
      type = ParsedCompute
      buffer = mu
      expression = '(1-phase)*mua + phase*mub'
      inputs = phase
      constant_names = 'mua mub'
      constant_expressions = '0.5 5'
    []
    [Finit]
      type = RankTwoIdentity
      buffer = F
    []
  []

  [Solve]
    [hyper_elasticity]
      type = HyperElasticIsotropic
      buffer = stress
      F = Fnew
      K = K
      mu = mu
    []

    [root]
      [applied_strain]
        type = MacroscopicShearTensor
        buffer = applied_strain
      []
      [mech]
        type = FFTMechanics
        buffer = Fnew
        F = F
        K = K
        mu = mu
        l_tol = 1e-2
        nl_rel_tol = 2e-2
        nl_abs_tol = 2e-2
        constitutive_model = hyper_elasticity
        stress = stress
        applied_macroscopic_strain = applied_strain
      []
    []
  []

  [Postprocess]
    [displacements]
      type = ComputeDisplacements
      buffer = disp
      F = F
    []
    [vonmises]
      type = ComputeVonMisesStress
      buffer = sV
    []
  []
[]

[TensorSolver]
  # no variables are integrated by this solver (FFTMechanics performs a steady state mechanics solve)
  type = ForwardEulerSolver
  root_compute = root
  # deformation tensor is just forwarded Fnew -> F
  forward_buffer = F
  forward_buffer_new = Fnew
  substeps = 10
[]

[TensorOutputs]
  active = 'parallel' # serial

  [serial]
    type = XDMFTensorOutput
    file_base = mech3d_serial
    buffer = 'disp sV F phase'
    output_mode = 'OVERSIZED_NODAL CELL CELL NODE'
    enable_hdf5 = true
    execute_on = 'TIMESTEP_END'
  []
  [parallel]
    # NODE and OVERSIZED_NODAL are not yet available in parallel output
    type = XDMFTensorOutput
    file_base = mech3d_parallel
    buffer = 'sV F phase'
    output_mode = 'CELL CELL CELL'
    enable_hdf5 = true
    execute_on = 'TIMESTEP_END'
  []
[]

[Executioner]
  type = Transient
  num_steps = 3
  dt = 0.01
[]

(test/tests/tensor_compute/rotating_grain_secant.i)

w=6

[Domain]
  dim = 2
  nx = 40
  ny = 40
  xmax = ${fparse w*pi*2}
  ymax = ${fparse w*pi*2/sin(pi/3)}
  mesh_mode = DOMAIN
[]

[AuxVariables]
  [phi]
  []
[]

[Outputs]
  exodus = false
[]

[TensorBuffers]
  [psi]
    map_to_aux_variable = phi
  []
  [psibar]
  []
  [psi3]
  []
  [psi3bar]
  []
  # constant tensors
  [linear]
  []
[]

crystal = '-(sin(sin(a)*y/2+cos(a)*x/2)^2 + sin(sin(a+1/3*pi)*y/2+cos(a+1/3*pi)*x/2)^2 + sin(sin(a-1/3*pi)*y/2+cos(a-1/3*pi)*x/2)^2 - 1.5)*0.25'
[Functions]
  [grain1]
    type = ParsedFunction
    expression = 'a := 0; ${crystal}'
  []
  [grain2]
    type = ParsedFunction
    expression = 'a := 0.95; ${crystal}'
  []
  [domain]
    type = ParsedFunction
    expression = 'r := (x-${w}*pi)^2+(y-${w}*pi)^2; if(r<(${w}*2/3*pi)^2, grain2, grain1)'
    symbol_names = 'grain1 grain2'
    symbol_values = 'grain1 grain2'
  []
[]

[TensorComputes]
  [Initialize]
    [psi]
      type = MooseFunctionTensor
      buffer = psi
      function = domain
    []
    [linear]
      type = SwiftHohenbergLinear
      buffer = linear
      alpha = 1
      r = 0.025
    []
  []

  [Solve]
    [psi3]
      type = ParsedCompute
      buffer = psi3
      expression = "0.20*psi^2-psi^3"
      inputs = psi
    []
    [psibar]
      type = ForwardFFT
      buffer = psibar
      input = psi
    []
    [psi3bar]
      type = ForwardFFT
      buffer = psi3bar
      input = psi3
    []
  []
[]

[TensorSolver]
  type = SecantSolver
  buffer = psi
  substeps = 3
  reciprocal_buffer = psibar
  linear_reciprocal = linear
  nonlinear_reciprocal = psi3bar
[]

[Problem]
  type = TensorProblem
[]

[Executioner]
  type = Transient
  num_steps = 10
  [TimeStepper]
    type = TensorSolveIterationAdaptiveDT
    dt = 1
    max_iterations = 400
    min_iterations = 100
    growth_factor = 1.4
    cutback_factor = 0.9
  []
  dtmax = 500
[]

[TensorOutputs]
  [xdmf]
    type = XDMFTensorOutput
    buffer = 'psi'
    enable_hdf5 = true
    # Do not transpose output to avoid regolding the test. In practice the default
    # of transpose = true should always be used
    transpose = false
  []
[]

(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/libtorch_kks/KKS_libtorch.i)

#
# Kim-Kim-Suzuki with Gibbs energy supplied by a torch model, solved on a 2D grid.
#


# Constants for Initial Conditions
r = 30
l = 4.2

# Initial condition function for order parameter
eta_IC = '0.5*(1-tanh(2*(sqrt(x^2+y^2)-${r})/${l}))'

# Phase-field model parameters
kappa_eta = 5
w = 1
M = 5
L = 5

# Expressions for switching function and bulk Gibbs energy
h_eta = 'eta^3*(6*eta^2-15*eta+10)'


[Domain]
    dim = 2
    nx = 100
    ny = 100

    xmin = -50
    xmax = 50
    ymin = -50
    ymax = 50

    # automatically create a matching mesh
    mesh_mode = DUMMY
[]

[TensorComputes]
    [Initialize]
        [c_IC]
            type = ParsedCompute
            buffer = c
            expression = '0.6 + (0.3-0.6)*${eta_IC}'
            extra_symbols = 'true'
        []
        [eta_IC]
            type = ParsedCompute
            buffer = eta
            expression = '${eta_IC}'
            extra_symbols = 'true'
        []
        [psi_init]
            type = ConstantTensor
            buffer = psi
            real = 1
        []

        [M]
            type = ConstantTensor
            buffer = M
            real = ${M}
        []
        [L]
            type = ConstantTensor
            buffer = L
            real = ${L}
        []
        [L_kappa]
            type = ReciprocalLaplacianFactor
            buffer = L_kappa
            factor = ${fparse  ${L} * ${kappa_eta} }
        []
        [h_eta_IC]
            type = ParsedCompute
            buffer = h_eta
            expression = '${h_eta}'
            inputs = eta
        []
        [G_func_IC]
            type = LibtorchGibbsEnergy
            buffer = 'G'
            phase_fractions = 'h_eta'
            concentrations = 'c'
            domega_detas = 'dG_dh'
            chem_pots = 'mu'
            libtorch_model_file = 'torch_NN_gibbs_model.pt'
        []
        [smooth]
            type = DeAliasingTensor
            method = HOULI
            buffer = smooth
        []
    []
    [Solve]
        [h_eta]
            type = ParsedCompute
            buffer = h_eta
            expression = '${h_eta}'
            inputs = eta
        []
        [G_func]
            type = LibtorchGibbsEnergy
            buffer = 'G'
            phase_fractions = 'h_eta'
            concentrations = 'c'
            domega_detas = 'dG_dh'
            chem_pots = 'mu'
            libtorch_model_file = 'torch_NN_gibbs_model.pt'
        []
        [dG_deta]
            type = ParsedCompute
            buffer = 'dG_deta'
            inputs = 'eta dG_dh'
            expression = 'dG_dh * ${h_eta} + ${w} * eta^2 * (1-eta^2)^2'
            derivatives = 'eta'
        []

        [etabar]
            type = ForwardFFT
            buffer = etabar
            input = eta
        []
        [AC_bulk]
            type = ReciprocalAllenCahn
            L = L
            buffer = AC_bulk
            dF_chem_deta = dG_deta
            psi = psi
        []
        [NL_eta]
            type = ParsedCompute
            buffer = NL_eta
            expression = 'AC_bulk '
            inputs = 'AC_bulk'
        []
        [cbar]
            type = ForwardFFT
            buffer = cbar
            input = c
        []
        [div_J]
            type = ReciprocalMatDiffusion
            buffer = 'div_J'
            chemical_potential = mu
            mobility = M
            psi = psi
        []
        [NL_c]
            type = ParsedCompute
            buffer = 'NL_c'
            inputs = 'div_J smooth'
            expression = 'smooth * div_J'
        []
    []
[]

[TensorSolver]
    type = AdamsBashforthMoulton
    buffer = 'c eta'
    reciprocal_buffer = 'cbar etabar'
    linear_reciprocal = '0 L_kappa'
    nonlinear_reciprocal = 'NL_c NL_eta'
    substeps = 1e3
    predictor_order = 3
    corrector_order = 1
    corrector_steps = 1
[]

[Postprocessors]
    [total_c]
        type = TensorIntegralPostprocessor
        buffer = c
    []
[]

[TensorOutputs]
    [xdmf]
        type = XDMFTensorOutput
        buffer = 'eta c mu psi dG_deta dG_dh G'
        enable_hdf5 = true
        transpose = false
    []
[]

[Executioner]
    type = Transient
    num_steps = 100
  [TimeStepper]
    type = IterationAdaptiveDT
    growth_factor = 1.25
    dt = 0.1
  []
  dtmax = 10
[]

[Outputs]
    csv = true
    perf_graph = true
    execute_on = 'INITIAL TIMESTEP_END'
[]

(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/cahnhilliard/cahnhilliard_explicit_smooth.i)

#
# Simple Cahn-Hilliard solve on a 2D grid. We create a matching (conforming)
# MOOSE mesh (with one element per FFT grid cell) and project the solution onto
# the MOOSE mesh to utilize the exodus output object.
#

[Domain]
  dim = 2
  nx = 50
  ny = 50
  xmax = 3
  ymax = 3
  mesh_mode = DOMAIN
  device_names = cpu
[]

[TensorBuffers]
  [c]
    map_to_aux_variable = c
  []
  [cbar]
  []
  [mu]
    map_to_aux_variable = mu
  []
[]

[TensorComputes]
  [Initialize]
    [c]
      # Random initial condition around a concentration of 1/2
      type = RandomTensor
      buffer = c
      min = 0.44
      max = 0.56
      seed = 0
    []
    [mu_init]
      type = ConstantTensor
      buffer = mu
    []

    # precompute fixed factors for the solve
    [Mbar]
      type = ReciprocalLaplacianFactor
      factor = 0.2 # Mobility
      buffer = Mbar
    []
    [Mkappabarbar]
      type = ReciprocalLaplacianSquareFactor
      factor = '${fparse 0.2 * 1e-4}' # M * kappa
      buffer = Mkappabarbar
    []
    [dc_dt_bar_IC]
      type = ConstantReciprocalTensor
      buffer = dc_dt_bar
    []

    [smooth]
      type = DeAliasingTensor
      buffer = smooth
    []
  []

  [Solve]
    [cahn_hilliard]
      [mu]
        type = ParsedCompute
        buffer = mu
        expression = '0.1*c^2*(c-1)^2'
        derivatives = c
        inputs = c
      []
      [mubar]
        type = ForwardFFT
        buffer = mubar
        input = mu
      []
      [dc_dt_bar]
        type = ParsedCompute
        buffer = dc_dt_bar
        expression = 'smooth * (Mbar*mubar - Mkappabarbar*cbar)'
        # expression = '(Mbar*mubar - Mkappabarbar*cbar)'
        inputs = 'Mbar mubar Mkappabarbar cbar smooth'
      []
      [cbar]
        type = ForwardFFT
        buffer = cbar
        input = c
      []
    []
  []
[]

[TensorSolver]
  type = ForwardEulerSolver
  time_derivative_reciprocal = dc_dt_bar
  root_compute = cahn_hilliard
  buffer = c
  reciprocal_buffer = cbar
  substeps = 50
[]

[AuxVariables]
  [mu]
    # the mu tensor  is projected onto this elemental variable
    family = MONOMIAL
    order = CONSTANT
  []
  [c]
    # the c tensor is projected onto this nodal variable
  []
[]

[Postprocessors]
  [C]
    type = TensorIntegralPostprocessor
    buffer = c
  []
[]

[Problem]
  type = TensorProblem
[]

[Executioner]
  type = Transient
  num_steps = 20
  dt = 0.5
[]

[Outputs]
  exodus = true
  csv = true
[]

(test/tests/cahnhilliard/cahnhilliard.i)

#
# Simple Cahn-Hilliard solve on a 2D grid. We create a matching (conforming)
# MOOSE mesh (with one element per FFT grid cell) and project the solution onto
# the MOOSE mesh to utilize the exodus output object.
#

[Domain]
  dim = 2
  nx = 20
  ny = 20
  xmax = 3
  ymax = 3
  mesh_mode = DOMAIN
[]

# In this input we fully trely on implicit TensorBuffer declaration

[TensorComputes]
  [Initialize]
    [c]
      # Random initial condition around a concentration of 1/2
      type = RandomTensor
      buffer = c
      min = 0.44
      max = 0.56
      seed = 0
    []
    [mu_init]
      type = ConstantTensor
      buffer = mu
    []

    # precompute fixed factors for the solve
    [Mbar]
      type = ReciprocalLaplacianFactor
      factor = 0.2 # Mobility
      buffer = Mbar
    []
    [kappabarbar]
      type = ReciprocalLaplacianSquareFactor
      factor = -0.001 # kappa
      buffer = kappabarbar
    []
  []

  [Solve]
    [cahn_hilliard]
      [mu]
        type = ParsedCompute
        buffer = mu
        expression = '0.1*c^2*(c-1)^2'
        derivatives = c
        inputs = c
      []
      [mubar]
        type = ForwardFFT
        buffer = mubar
        input = mu
      []
      [Mbarmubar]
        type = ParsedCompute
        buffer = Mbarmubar
        expression = 'Mbar*mubar'
        inputs = 'Mbar mubar'
      []
      [cbar]
        type = ForwardFFT
        buffer = cbar
        input = c
      []
    []
  []
[]

[TensorSolver]
  type = AdamsBashforthMoulton
  root_compute = cahn_hilliard
  buffer = c
  reciprocal_buffer = cbar
  linear_reciprocal = kappabarbar
  nonlinear_reciprocal = Mbarmubar
  substeps = 10
[]

[AuxVariables]
  [mu]
    # the mu tensor  is projected onto this elemental variable
    family = MONOMIAL
    order = CONSTANT
  []
  [c]
    # the c tensor is projected onto this nodal variable
  []
[]

[AuxKernels]
  active = ''
  [c]
    type = ProjectTensorAux
    buffer = c
    variable = c
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [mu]
    type = ProjectTensorAux
    buffer = mu
    variable = mu
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]

[Postprocessors]
  [min_c]
    type = SemiImplicitCriticalTimeStep
    buffer = kappabarbar
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [delta_int_c]
    type = TensorIntegralChangePostprocessor
    buffer = c
  []
[]

[Problem]
  type = TensorProblem
[]

[Executioner]
  type = Transient
  num_steps = 10
  dt = 1e-3
[]

[TensorOutputs]
  active = ''
  [xdmf]
    type = XDMFTensorOutput
    buffer = 'c mu'
    output_mode = 'Node Cell'
    enable_hdf5 = true
    # Do not transpose output to avoid regolding the test. In practice the default
    # of transpose = true should always be used
    transpose = false
  []

  [xdmf2]
    # second output to trigger the hdf5 thread safety error
    type = XDMFTensorOutput
    buffer = 'c'
    output_mode = 'Cell'
    enable_hdf5 = true
    # Do not transpose output to avoid regolding the test. In practice the default
    # of transpose = true should always be used
    transpose = false
  []
[]

[Outputs]
  exodus = true
  csv = true
[]

(test/tests/cahnhilliard/cahnhilliard_explicit.i)

#
# Simple Cahn-Hilliard solve on a 2D grid. We create a matching (conforming)
# MOOSE mesh (with one element per FFT grid cell) and project the solution onto
# the MOOSE mesh to utilize the exodus output object.
#

[Domain]
  dim = 2
  nx = 50
  ny = 50
  xmax = 3
  ymax = 3
  mesh_mode = DOMAIN
  device_names = cpu
[]

[TensorBuffers]
  [c]
    map_to_aux_variable = c
  []
  [cbar]
  []
  [mu]
    map_to_aux_variable = mu
  []
  [mubar]
  []
  [dc_dt_bar]
  []
  # constant tensors
  [Mbar]
  []
  [Mkappabarbar]
  []
[]

[TensorComputes]
  [Initialize]
    [c]
      # Random initial condition around a concentration of 1/2
      type = RandomTensor
      buffer = c
      min = 0.44
      max = 0.56
      seed = 0
    []
    [mu_init]
      type = ConstantTensor
      buffer = mu
    []

    # precompute fixed factors for the solve
    [Mbar]
      type = ReciprocalLaplacianFactor
      factor = 0.2 # Mobility
      buffer = Mbar
    []
    [Mkappabarbar]
      type = ReciprocalLaplacianSquareFactor
      factor = '${fparse 0.2 * 1e-4}' # M * kappa
      buffer = Mkappabarbar
    []
    [dc_dt_bar_IC]
      type = ConstantReciprocalTensor
      buffer = dc_dt_bar
    []
  []

  [Solve]
    [cahn_hilliard]
      [mu]
        type = ParsedCompute
        buffer = mu
        expression = '0.1*c^2*(c-1)^2'
        derivatives = c
        inputs = c
      []
      [mubar]
        type = ForwardFFT
        buffer = mubar
        input = mu
      []
      [dc_dt_bar]
        type = ParsedCompute
        buffer = dc_dt_bar
        expression = 'Mbar*mubar - Mkappabarbar*cbar'
        inputs = 'Mbar mubar Mkappabarbar cbar'
      []
      [cbar]
        type = ForwardFFT
        buffer = cbar
        input = c
      []
    []
  []
[]

[TensorSolver]
  type = ForwardEulerSolver
  time_derivative_reciprocal = dc_dt_bar
  root_compute = cahn_hilliard
  buffer = c
  reciprocal_buffer = cbar
  substeps = 50
[]

[AuxVariables]
  [mu]
    # the mu tensor  is projected onto this elemental variable
    family = MONOMIAL
    order = CONSTANT
  []
  [c]
    # the c tensor is projected onto this nodal variable
  []
[]

[Problem]
  type = TensorProblem
[]

[Executioner]
  type = Transient
  num_steps = 100
  dt = 1e-1
[]

[Outputs]
  exodus = true
  csv = true
[]

(examples/swift_hohenberg/swifthohenberg.i)

#
# Solve a simple Swift-Hohenberg crystal phase field problem. The initial condition is
# a random melt, which crystalizes out.
#

[Domain]
  dim = 2
  nx = 400
  ny = 400
  xmax = ${fparse pi*2*30}
  ymax = ${fparse pi*2*30}

  device_names = 'cuda'

  mesh_mode = DOMAIN
[]

[TensorBuffers]
  [psi]
    map_to_aux_variable = psi
  []
  [psibar]
  []
  [psi3]
  []
  [psi3bar]
  []
  # constant tensors
  [linear]
  []

  # output
  [filter]
    map_to_aux_variable = filter
  []
  [filterbar]
  []
[]

[AuxVariables]
  [psi]
  []
  [filter]
  []
[]

[TensorComputes]
  [Initialize]
    [psi]
      type = RandomTensor
      buffer = psi
      min = 0
      max = 0.07
    []
    [linear]
      type = SwiftHohenbergLinear
      buffer = linear
      alpha = 1
      r = 0.025
    []
  []

  [Solve]
    [psi3]
      type = ParsedCompute
      buffer = psi3
      expression = "0.20*psi^2-psi^3"
      inputs = psi
    []
    [psibar]
      type = ForwardFFT
      buffer = psibar
      input = psi
    []
    [psi3bar]
      type = ForwardFFT
      buffer = psi3bar
      input = psi3
    []
  []

  [Postprocess]
    [low_pass]
      type = ParsedCompute
      buffer = filterbar
      extra_symbols = true
      expression = 'psibar * exp(-k2*10)'
      inputs = psibar
    []
    [filter]
      type = InverseFFT
      buffer = filter
      input = filterbar
    []
  []
[]

[TensorTimeIntegrators]
  [c]
    type = FFTSemiImplicit
    buffer = psi
    reciprocal_buffer = psibar
    linear_reciprocal = linear
    nonlinear_reciprocal = psi3bar
  []
[]

[Problem]
  type = TensorProblem
  spectral_solve_substeps = 1000
[]

[Executioner]
  type = Transient
  num_steps = 300
  [TimeStepper]
    type = IterationAdaptiveDT
    growth_factor = 1.2
    dt = 10
  []
  dtmax = 1000
[]

[Postprocessors]
  [min_psi]
    type = TensorExtremeValuePostprocessor
    buffer = psi
    value_type = MIN
    execute_on = 'TIMESTEP_END'
  []
  [max_psi]
    type = TensorExtremeValuePostprocessor
    buffer = psi
    value_type = MAX
    execute_on = 'TIMESTEP_END'
  []
  [Psi]
    type = TensorIntegralPostprocessor
    buffer = psi
  []
[]

[TensorOutputs]
  [xdmf]
    type = XDMFTensorOutput
    buffer = 'psi'
    output_mode = 'Node'
    enable_hdf5 = true
  []
[]

[Outputs]
  exodus = true
  perf_graph = true
  execute_on = 'TIMESTEP_END'
[]

(test/tests/tensor_compute/coupled_pf_mech_secant.i)

[Domain]
  dim = 3
  nx = 128
  ny = 128
  nz = 128
  xmax = ${fparse pi*4}
  ymax = ${fparse pi*4}
  zmax = ${fparse pi*4}
  mesh_mode = DUMMY
[]

[TensorBuffers]
  # phase field
  [c]
  []
  [cbar]
  []
  [mu]
  []
  [mubar]
  []
  [Mbarmubar]
  []

  # mechanics
  [disp_x]
  []
  [disp_y]
  []
  [disp_z]
  []
  [mumechbar]
  []
  [mumech]
  []

  # constant tensors
  [Mbar]
  []
  [kappabarbar]
  []
[]

[TensorOutputs]
  [xdmf]
    type = XDMFTensorOutput
    buffer = 'c disp_x disp_y disp_z mu mumech'
    output_mode = 'Node Node Node Node Cell Cell'
    enable_hdf5 = true
  []
[]

[TensorComputes]
  [Initialize]
    [c]
      type = RandomTensor
      buffer = c
      min = 0.44
      max = 0.56
    []
    [mu_init]
      type = ConstantTensor
      buffer = mu
    []
    [mumech_init]
      type = ConstantTensor
      buffer = mumech
    []
    [disp_x]
      type = RandomTensor
      buffer = disp_x
      min = 0
      max = 0
    []
    [disp_y]
      type = RandomTensor
      buffer = disp_y
      min = 0
      max = 0
    []
    [disp_z]
      type = RandomTensor
      buffer = disp_z
      min = 0
      max = 0
    []
    [Mbar]
      type = ReciprocalLaplacianFactor
      factor = 0.2 # Mobility
      buffer = Mbar
    []
    [kappabarbar]
      type = ReciprocalLaplacianSquareFactor
      factor = -0.001 # kappa
      buffer = kappabarbar
    []
  []

  [Solve]
    [mu]
      # chemical potential (real space)
      type = ParsedCompute
      buffer = mu
      expression = '0.1*c^2*(c-1)^2' # + c*sin(x/2)*0.005'
      extra_symbols = true
      derivatives = c
      inputs = c
    []
    [mubar]
      # chemical potential (reciprocal space)
      type = ForwardFFT
      buffer = mubar
      input = mu
    []
    [mumechbar]
      # mechanical chemical potential (reciprocal space)
      type = FFTElasticChemicalPotential
      buffer = mumechbar
      cbar = cbar
      displacements = 'disp_x disp_y disp_z'
      lambda = 100
      mu = 50
      e0 = 0.02
    []
    [mumech]
      # chemical potential (reciprocal space)
      type = InverseFFT
      buffer = mumech
      input = mumechbar
    []

    [Mbarmubar]
      type = ParsedCompute
      buffer = Mbarmubar
      expression = 'Mbar*(mubar+mumechbar)'
      inputs = 'Mbar mubar mumechbar'
    []
    [cbar]
      type = ForwardFFT
      buffer = cbar
      input = c
    []

    [qsmech]
      type = FFTQuasistaticElasticity
      displacements = 'disp_x disp_y disp_z'
      cbar = cbar
      lambda = 100
      mu = 50
      e0 = 0.02
    []
  []
[]

[TensorSolver]
  type = SecantSolver
  substeps = 1
  max_iterations = 1000
  # damping = 0.75
  relative_tolerance = 1e-6
  absolute_tolerance = 1e-6
  buffer = c
  dt_epsilon = 1e-7
  reciprocal_buffer = cbar
  linear_reciprocal = kappabarbar
  nonlinear_reciprocal = Mbarmubar
  verbose = true
[]

[Postprocessors]
  [min_c]
    type = TensorExtremeValuePostprocessor
    buffer = c
    value_type = MIN
    execute_on = 'TIMESTEP_END'
  []
  [max_c]
    type = TensorExtremeValuePostprocessor
    buffer = c
    value_type = MAX
    execute_on = 'TIMESTEP_END'
  []

  [min_disp_x]
    type = TensorExtremeValuePostprocessor
    buffer = disp_x
    value_type = MIN
    execute_on = 'TIMESTEP_END'
  []
  [max_disp_x]
    type = TensorExtremeValuePostprocessor
    buffer = disp_x
    value_type = MAX
    execute_on = 'TIMESTEP_END'
  []
  [min_disp_y]
    type = TensorExtremeValuePostprocessor
    buffer = disp_y
    value_type = MIN
    execute_on = 'TIMESTEP_END'
  []
  [max_disp_y]
    type = TensorExtremeValuePostprocessor
    buffer = disp_y
    value_type = MAX
    execute_on = 'TIMESTEP_END'
  []
  [min_disp_z]
    type = TensorExtremeValuePostprocessor
    buffer = disp_z
    value_type = MIN
    execute_on = 'TIMESTEP_END'
  []
  [max_disp_z]
    type = TensorExtremeValuePostprocessor
    buffer = disp_z
    value_type = MAX
    execute_on = 'TIMESTEP_END'
  []

  [C]
    type = TensorIntegralPostprocessor
    buffer = c
    execute_on = 'TIMESTEP_END'
  []
  [cavg]
    type = TensorAveragePostprocessor
    buffer = c
    execute_on = 'TIMESTEP_END'
  []
[]

[Problem]
  type = TensorProblem
[]

[Executioner]
  type = Transient
  end_time = 100
  [TimeStepper]
    type = TensorSolveIterationAdaptiveDT
    dt = 0.1
    max_iterations = 500
    min_iterations = 300
    growth_factor = 1.1
    cutback_factor = 0.9
  []
[]

[Outputs]
  csv = true
  perf_graph = true
  execute_on = 'TIMESTEP_END'
[]

(test/tests/tensor_compute/parallel.i)

[Domain]
  device_names = "cuda:1 cuda:0 cpu"
  device_weights = "10 10 1"

  parallel_mode = FFT_SLAB

  dim = 2
  nx = 400
  ny = 400
  xmax = ${fparse pi*2*30}
  ymax = ${fparse pi*2*30}
[]


# [Mesh]
#   type = GeneratedMesh
#   dim =1
# []


[TensorBuffers]
  [psi]
    # map_to_aux_variable = psi
  []
  [psibar]
  []
  [psi3]
  []
  [psi3bar]
  []
  # constant tensors
  [linear]
  []

  # output
  [filter]
    # map_to_aux_variable = filter
  []
  [filterbar]
  []
[]

[AuxVariables]
  [psi]
  []
  [filter]
  []
[]

[TensorComputes]
  [Initialize]
    [psi]
      type = RandomTensor
      buffer = psi
      min = 0
      max = 0.07
    []
    [linear]
      type = SwiftHohenbergLinear
      buffer = linear
      alpha = 1
      r = 0.025
    []
  []

  [Solve]
    [psi3]
      type = ParsedCompute
      buffer = psi3
      expression = "0.20*psi^2-psi^3"
      inputs = psi
    []
    [psibar]
      type = ForwardFFT
      buffer = psibar
      input = psi
    []
    [psi3bar]
      type = ForwardFFT
      buffer = psi3bar
      input = psi3
    []
  []

  [Postprocess]
    [low_pass]
      type = ParsedCompute
      buffer = filterbar
      extra_symbols = true
      expression = 'psibar * exp(-k2*10)'
      inputs = psibar
    []
    [filter]
      type = InverseFFT
      buffer = filter
      input = filterbar
    []
  []
[]

[TensorTimeIntegrators]
  [c]
    type = FFTSemiImplicit
    buffer = psi
    reciprocal_buffer = psibar
    linear_reciprocal = linear
    nonlinear_reciprocal = psi3bar
  []
[]

[Problem]
  type = TensorProblem
  spectral_solve_substeps = 1000
[]

[Executioner]
  type = Transient
  num_steps = 300
  [TimeStepper]
    type = IterationAdaptiveDT
    growth_factor = 1.2
    dt = 10
  []
  dtmax = 1000
[]

[Postprocessors]
  [min_psi]
    type = TensorExtremeValuePostprocessor
    buffer = psi
    value_type = MIN
    execute_on = 'TIMESTEP_END'
  []
  [max_psi]
    type = TensorExtremeValuePostprocessor
    buffer = psi
    value_type = MAX
    execute_on = 'TIMESTEP_END'
  []
  [Psi]
    type = TensorIntegralPostprocessor
    buffer = psi
  []
[]

# [TensorOutputs]
#   [xdmf]
#     type = XDMFTensorOutput
#     buffer = 'psi'
#     output_mode = 'Node'
#     enable_hdf5 = true
#   []
# []

[Outputs]
  exodus = true
  perf_graph = true
  execute_on = 'TIMESTEP_END'
[]

(test/tests/real_space/laplacian.i)

[Domain]
  dim = 2
  nx = 80
  ny = 60
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y'
  xmin = -1
  ymin = -1
  xmax = 1
  ymax = 1
[]

[TensorBuffers]
  [u]
    type = PlainTensorBuffer
  []
  [u_gold]
    type = PlainTensorBuffer
  []
  [lap]
    type = PlainTensorBuffer
  []
[]

[TensorComputes]
  [Initialize]
    [u_ic]
      type = ParsedCompute
      buffer = u
      expression = 'sin(pi*(x+0.5))*sin(pi*(y+0.5))'
      extra_symbols = true
    []
    [lap_gold_ic]
      type = ParsedCompute
      buffer = lap_gold
      #expression = '0'
      expression = '-2*pi*pi*sin(pi*(x+0.5))*sin(pi*(y+0.5))'
      extra_symbols = true
      expand = REAL
    []
  []
  [Solve]
    [lap_compute]
      type = FiniteDifferenceLaplacian
      input = u
      buffer = lap
    []
    [diff]
      type = ParsedCompute
      buffer = diff
      expression = 'abs(lap - lap_gold)'
      inputs = 'lap lap_gold'
    []
  []
[]

[Postprocessors]
  [max_error]
    type = TensorExtremeValuePostprocessor
    buffer = diff
    value_type = MAX
  []
[]

[TensorOutputs]
  active = ''
  [lap]
    type = XDMFTensorOutput
    buffer = 'lap lap_gold diff'
    enable_hdf5 = true
  []
[]

[Executioner]
  type = Transient
  num_steps = 1
[]

[Outputs]
  [out]
    type = CSV
  []
[]

(examples/degeus_mechanics/mech.i)

[Domain]
  dim = 3
  nx = 32
  ny = 32
  nz = 32
  xmax = ${fparse 2*pi}
  ymax = ${fparse 2*pi}
  zmax = ${fparse 2*pi}
  mesh_mode = DUMMY
[]

[TensorComputes]
  [Initialize]
    [Finit]
      type = RankTwoIdentity
      buffer = F
    []

    [phase]
      type = PhaseMechanicsTest
      buffer = phase
    []
    [K]
      type = ParsedCompute
      buffer = K
      expression = '(1-phase)*Ka + phase*Kb'
      inputs = phase
      constant_names = 'Ka Kb'
      constant_expressions = '0.833 8.33'
    []
    [mu]
      type = ParsedCompute
      buffer = mu
      expression = '(1-phase)*mua + phase*mub'
      inputs = phase
      constant_names = 'mua mub'
      constant_expressions = '0.386 3.86'
    []
  []

  [Solve]
    [hyper_elasticity]
      type = HyperElasticIsotropic
      buffer = stress
      F = Fnew
      K = K
      mu = mu
    []

    [root]
      [applied_strain]
        type = MacroscopicShearTensor
        buffer = applied_strain
      []
      [mech]
        type = FFTMechanics
        buffer = Fnew
        F = F
        K = K
        mu = mu
        l_tol = 1e-2
        nl_rel_tol = 2e-2
        nl_abs_tol = 2e-2
        constitutive_model = hyper_elasticity
        stress = stress
        applied_macroscopic_strain = applied_strain
      []
    []
  []

  [Postprocess]
    [displacements]
      type = ComputeDisplacements
      buffer = disp
      F = F
    []
    [vonmises]
      type = ComputeVonMisesStress
      buffer = sV
    []
  []
[]

[TensorSolver]
  # no variables are integrated by this solver (FFTMechanics performs a steady state mechanics solve)
  type = ForwardEulerSolver
  root_compute = root
  # deformation tensor is just forwarded Fnew -> F
  forward_buffer = F
  forward_buffer_new = Fnew
  substeps = 10
[]

[TensorOutputs]
  [deformation_tensor]
    type = XDMFTensorOutput
    buffer = 'disp sV F'
    output_mode = 'OVERSIZED_NODAL CELL CELL'
    enable_hdf5 = true
  []
[]

[Executioner]
  type = Transient
  num_steps = 100
  dt = 0.01
[]

(test/tests/solvers/nl_coupled.i)

#
# Simple Cahn-Hilliard solve on a 2D grid.
#

[Domain]
  dim = 2
  nx = 150
  ny = 150
  xmax = '${fparse pi*2}'
  ymax = '${fparse pi*2}'
  mesh_mode = DUMMY
[]

[GlobalParams]
  constant_names = 'A B'
  constant_expressions = '1 3.5'
[]

[TensorComputes]
  [Initialize]
    [u]
      type = ParsedCompute
      buffer = u
      extra_symbols = true
      expression = 'sin(x)*sin(y)'
      expand = REAL
    []
    [v]
      type = ParsedCompute
      buffer = v
      extra_symbols = true
      expression = 'cos(x)*cos(y)'
      expand = REAL
    []
    [zero]
      type = ConstantReciprocalTensor
      buffer = zero
    []

    # precompute fixed factors for the solve
    [D1]
      type = ReciprocalLaplacianFactor
      factor = 1e-2
      buffer = D1
    []
    [D2]
      type = ReciprocalLaplacianFactor
      factor = 1e-3
      buffer = D2
    []
  []

  [Solve]
    [u_bar]
      type = ForwardFFT
      buffer = u_bar
      input = u
    []
    [v_bar]
      type = ForwardFFT
      buffer = v_bar
      input = v
    []

    [Du]
      type = ParsedCompute
      buffer = Du
      expression = 'D2*v_bar'
      inputs = 'D2 v_bar'
    []
    [Dv]
      type = ParsedCompute
      buffer = Dv
      expression = 'D2*u_bar'
      inputs = 'D2 u_bar'
    []
  []
[]

[TensorSolver]
  type = AdamsBashforthMoulton
  buffer = 'u v'
  reciprocal_buffer = 'u_bar v_bar'
  linear_reciprocal = 'D1 D1'
  nonlinear_reciprocal = 'Du Dv'
  substeps = ${ss}
  corrector_steps = ${cs}
  predictor_order = ${order}
  corrector_order = ${order}
[]

[Problem]
  type = TensorProblem
[]

[Postprocessors]
  [u_min]
    type = TensorExtremeValuePostprocessor
    buffer = u
    value_type = MIN
  []
  [u_max]
    type = TensorExtremeValuePostprocessor
    buffer = u
    value_type = MAX
  []
  [v_min]
    type = TensorExtremeValuePostprocessor
    buffer = v
    value_type = MIN
  []
  [v_max]
    type = TensorExtremeValuePostprocessor
    buffer = v
    value_type = MAX
  []
  [U]
    type = TensorIntegralPostprocessor
    buffer = u
  []
  [V]
    type = TensorIntegralPostprocessor
    buffer = v
  []
[]

[Executioner]
  type = Transient
  num_steps = 25
  dt = 10
[]

[Outputs]
  file_base = nl_coupled_${ss}_${cs}_${order}
  csv = true
[]

(benchmarks/02_oswald_ripening/2a_broyden.i)

[Domain]
  dim = 2
  nx = 200
  ny = 200
  xmax = 200
  ymax = 200

  device_names = 'cuda:0'

  mesh_mode = DOMAIN
[]

fchem='fa:=rho^2*(c-ca)^2;
fb:=rho^2*(cb-c)^2;
h:=n1^3*(6*n1^2-15*n1+10) +
   n2^3*(6*n2^2-15*n2+10) +
   n3^3*(6*n3^2-15*n3+10) +
   n4^3*(6*n4^2-15*n4+10);
g:=n1^2*(1-n1)^2 +
   n2^2*(1-n2)^2 +
   n3^2*(1-n3)^2 +
   n4^2*(1-n4)^2 +
alpha*(
n1^2*n2^2 + n1^2*n3^2 + n1^2*n4^2 +
n2^2*n1^2 + n2^2*n3^2 + n2^2*n4^2 +
n3^2*n1^2 + n3^2*n2^2 + n3^2*n4^2 +
n4^2*n1^2 + n4^2*n2^2 + n4^2*n3^2);
(fa*(1-h) + fb*h + w*g)'

nic = 'epsilon*(cos((0.01*idx)*x-4)*cos((0.007+0.01*idx)*y)
       +cos((0.11+0.01*idx)*x)*cos((0.11+0.01*idx)*y)
       +psi*(cos((0.046+0.001*idx)*x+(0.0405+0.001*idx)*y)
       *cos((0.031+0.001*idx)*x-(0.004+0.001*idx)*y))^2)^2'

cnames = 'rho     ca  cb  alpha w L M'
cvalues= 'sqrt(2) 0.3 0.7 5     1 5 5'

[TensorBuffers]
  # variables
  [c]
    map_to_aux_variable = c
  []
  [n1]
  []
  [n2]
  []
  [n3]
  []
  [n4]
  []

  [c_bar]
  []
  [n1_bar]
  []
  [n2_bar]
  []
  [n3_bar]
  []
  [n4_bar]
  []

  [mu_c]
    # map_to_aux_variable = mu
  []
  [mu_n1]
  []
  [mu_n2]
  []
  [mu_n3]
  []
  [mu_n4]
  []

  [mu_c_bar]
  []
  [mu_n1_bar]
  []
  [mu_n2_bar]
  []
  [mu_n3_bar]
  []
  [mu_n4_bar]
  []

  [Mbar_mu_c_bar]
  []

  # constant tensors
  [Lbar] # FFT(M*laplacian)
  []
  [MkappaL2bar] # FFT(-M*kappa*laplacian^2)
  []
  [kappaLbar] # FFT(L*kappa*laplacian)
  []

  # postprocessing
  [F]
  []
  [Fgrad_c]
  []
  [Fgrad_n1]
  []
  [Fgrad_n2]
  []
  [Fgrad_n3]
  []
  [Fgrad_n4]
  []
  [bnds]
    map_to_aux_variable = bnds
  []
[]

[TensorComputes]
  [Initialize]
    [c]
      type = ParsedCompute
      buffer = c
      extra_symbols = true
      expression = 'c0+epsilon*(cos(0.105*x)*cos(0.11*y)+(cos(0.13*x)*cos(0.087*y))^2+cos(0.025*x-0.15*y)*cos(0.07*x-0.02*y))'
      constant_names = 'c0 epsilon'
      constant_expressions = '0.5 0.01'
    []
    [Lbar]
      type = ReciprocalLaplacianFactor
      # Mobility is pulled into the chemical potential below
      buffer = Lbar
    []
    [MkappaL2bar]
      type = ReciprocalLaplacianSquareFactor
      factor = -15 # -kappa_c*M
      buffer = MkappaL2bar
    []
    [kappaLbar]
      type = ReciprocalLaplacianFactor
      buffer = kappaLbar
      factor = 15 # kappa_ni*L
    []
    [n1]
      type = ParsedCompute
      buffer = n1
      expression = ${nic}
      extra_symbols = true
      constant_names       = 'idx epsilon psi'
      constant_expressions = '  1     0.1 1.5'
    []
    [n2]
      type = ParsedCompute
      buffer = n2
      expression = ${nic}
      extra_symbols = true
      constant_names       = 'idx epsilon psi'
      constant_expressions = '  2    0.1 1.5'
    []
    [n3]
      type = ParsedCompute
      buffer = n3
      expression = ${nic}
      extra_symbols = true
      constant_names       = 'idx epsilon psi'
      constant_expressions = '  3     0.1 1.5'
    []
    [n4]
      type = ParsedCompute
      buffer = n4
      expression = ${nic}
      extra_symbols = true
      constant_names       = 'idx epsilon psi'
      constant_expressions = '  4     0.1 1.5'
    []
  []

  [Solve]
    [mu_c]
      type = ParsedCompute
      buffer = mu_c
      expression = '${fchem}*M'
      constant_names = ${cnames}
      constant_expressions = ${cvalues}
      derivatives = c
      inputs = 'c n1 n2 n3 n4'
    []

    [mu_n1]
      type = ParsedCompute
      buffer = mu_n1
      expression = '${fchem}*(-L)'
      constant_names = ${cnames}
      constant_expressions = ${cvalues}
      derivatives = n1
      inputs = 'c n1 n2 n3 n4'
    []
    [mu_n2]
      type = ParsedCompute
      buffer = mu_n2
      expression = '${fchem}*(-L)'
      constant_names = ${cnames}
      constant_expressions = ${cvalues}
      derivatives = n2
      inputs = 'c n1 n2 n3 n4'
    []
    [mu_n3]
      type = ParsedCompute
      buffer = mu_n3
      expression = '${fchem}*(-L)'
      constant_names = ${cnames}
      constant_expressions = ${cvalues}
      derivatives = n3
      inputs = 'c n1 n2 n3 n4'
    []
    [mu_n4]
      type = ParsedCompute
      buffer = mu_n4
      expression = '${fchem}*(-L)'
      constant_names = ${cnames}
      constant_expressions = ${cvalues}
      derivatives = n4
      inputs = 'c n1 n2 n3 n4'
    []

    [mu_c_bar]
      type = ForwardFFT
      buffer = mu_c_bar
      input = mu_c
    []
    [mu_n1_bar]
      type = ForwardFFT
      buffer = mu_n1_bar
      input = mu_n1
    []
    [mu_n2_bar]
      type = ForwardFFT
      buffer = mu_n2_bar
      input = mu_n2
    []
    [mu_n3_bar]
      type = ForwardFFT
      buffer = mu_n3_bar
      input = mu_n3
    []
    [mu_n4_bar]
      type = ForwardFFT
      buffer = mu_n4_bar
      input = mu_n4
    []

    [Mbar_mu_c_bar]
      type = ParsedCompute
      buffer = Mbar_mu_c_bar
      expression = 'Lbar*mu_c_bar'
      inputs = 'Lbar mu_c_bar'
    []

    [c_bar]
      type = ForwardFFT
      buffer = c_bar
      input = c
    []
    [n1_bar]
      type = ForwardFFT
      buffer = n1_bar
      input = n1
    []
    [n2_bar]
      type = ForwardFFT
      buffer = n2_bar
      input = n2
    []
    [n3_bar]
      type = ForwardFFT
      buffer = n3_bar
      input = n3
    []
    [n4_bar]
      type = ForwardFFT
      buffer = n4_bar
      input = n4
    []
  []

  [Postprocess]
    [Fgrad_c]
      type = FFTGradientSquare
      buffer = Fgrad_c
      input = c
      factor = 1.5 # kappa/2
    []
    [Fgrad_n1]
      type = FFTGradientSquare
      buffer = Fgrad_n1
      input = n1
      factor = 1.5 # kappa/2
    []
    [Fgrad_n2]
      type = FFTGradientSquare
      buffer = Fgrad_n2
      input = n2
      factor = 1.5 # kappa/2
    []
    [Fgrad_n3]
      type = FFTGradientSquare
      buffer = Fgrad_n3
      input = n3
      factor = 1.5 # kappa/2
    []
    [Fgrad_n4]
      type = FFTGradientSquare
      buffer = Fgrad_n4
      input = n4
      factor = 1.5 # kappa/2
    []
    [F]
      type = ParsedCompute
      buffer = F
      expression = '${fchem} + Fgrad_c + Fgrad_n1 + Fgrad_n2 + Fgrad_n3 + Fgrad_n4'
      constant_names =  ${cnames}
      constant_expressions = ${cvalues}
      inputs = 'c n1 n2 n3 n4 Fgrad_c Fgrad_n1 Fgrad_n2 Fgrad_n3 Fgrad_n4'
    []
    [bnds]
      type = ParsedCompute
      buffer = bnds
      expression = 'n1^2 + n2^2 + n3^2 + n4^2'
      inputs = 'n1 n2 n3 n4'
    []
  []
[]

[TensorSolver]
  type = BroydenSolver
  substeps = 10
  max_iterations = 1000
  damping = 0.5
  relative_tolerance = 1e-6
  absolute_tolerance = 1e-5
  buffer = 'c n1 n2 n3 n4'
  dt_epsilon = 1e-5
  reciprocal_buffer = 'c_bar n1_bar n2_bar n3_bar n4_bar'
  linear_reciprocal = 'MkappaL2bar kappaLbar kappaLbar kappaLbar kappaLbar'
  nonlinear_reciprocal = 'Mbar_mu_c_bar mu_n1_bar mu_n2_bar mu_n3_bar mu_n4_bar'
  # verbose = true
[]

[AuxVariables]
  [c]
  []
  [bnds]
  []
[]

[Postprocessors]
  [min_c]
    type = TensorExtremeValuePostprocessor
    buffer = c
    value_type = MIN
    execute_on = 'TIMESTEP_END'
  []
  [max_c]
    type = TensorExtremeValuePostprocessor
    buffer = c
    value_type = MAX
    execute_on = 'TIMESTEP_END'
  []
  [F]
    type = TensorIntegralPostprocessor
    buffer = F
  []
[]

[Problem]
  type = TensorProblem
[]

[Executioner]
  type = Transient
  # num_steps = 2
  end_time = 1e+5
  [TimeStepper]
    type = TensorSolveIterationAdaptiveDT
    dt = 1e-4
    max_iterations = 500
    min_iterations = 100
    growth_factor = 1.1
  []
[]

[Outputs]
  exodus = true
  csv = true
  perf_graph = true
  execute_on = 'TIMESTEP_END'
[]

(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/tensor_compute/parallel_roundtrip_3d.i)

[Domain]
  # Test parallel FFT round-trip with slab decomposition in 3D
  device_names = "cpu cpu cpu"
  device_weights = "1 1 1"

  dim = 3
  nx = 64
  ny = 64
  nz = 64
  xmax = ${fparse pi*4}
  ymax = ${fparse pi*4}
  zmax = ${fparse pi*4}
[]

[TensorBuffers]
  [eta_gold]
  []
  [eta]
  []
  [eta_bar]
  []
  [eta_roundtrip]
  []
  [diff]
  []
  [zero]
  []
[]

[TensorComputes]
  [Initialize]
    [eta_gold]
      type = ParsedCompute
      buffer = eta_gold
      expression = 'sin(x)+sin(y)+sin(z)+cos(2*x)*sin(3*y)*cos(z)'
      extra_symbols = true
    []
    [eta]
      type = ParsedCompute
      buffer = eta
      expression = eta_gold
      inputs = eta_gold
    []
    [zero]
      type = ConstantReciprocalTensor
      buffer = zero
      real = 0
      imaginary = 0
    []
  []

  [Solve]
    # Test: eta -> FFT -> iFFT -> eta_roundtrip
    # eta_roundtrip should equal eta (within numerical precision)
    [eta_bar]
      type = ForwardFFT
      buffer = eta_bar
      input = eta
    []
    [eta_roundtrip]
      type = InverseFFT
      buffer = eta_roundtrip
      input = eta_bar
    []
  []

  [Postprocess]
    [diff]
      type = ParsedCompute
      buffer = diff
      expression = 'abs(eta - eta_roundtrip) + abs(eta - eta_gold)'
      inputs = 'eta eta_roundtrip eta_gold'
    []
  []
[]

[Postprocessors]
  [max_error]
    type = TensorExtremeValuePostprocessor
    buffer = diff
    value_type = MAX
  []
  [l2_error]
    type = TensorIntegralPostprocessor
    buffer = diff
  []
[]

[TensorSolver]
  type = AdamsBashforthMoulton
  buffer = eta
  reciprocal_buffer = eta_bar
  linear_reciprocal = zero
  nonlinear_reciprocal = zero
[]

[TensorOutputs]
  [eta]
    type = XDMFTensorOutput
    buffer = 'eta'
    output_mode = 'CELL'
    enable_hdf5 = true
  []
[]

[Problem]
  type = TensorProblem
[]

[Executioner]
  type = Transient
  num_steps = 1
[]

[Outputs]
  csv = true
  execute_on = 'INITIAL TIMESTEP_END'
[]

(benchmarks/02_oswald_ripening/2a_secant.i)

[Domain]
  dim = 2
  nx = 120
  ny = 120
  xmax = 200
  ymax = 200

  device_names = 'cuda:1'

  mesh_mode = DOMAIN
[]

fchem='fa:=rho^2*(c-ca)^2;
fb:=rho^2*(cb-c)^2;
h:=n1^3*(6*n1^2-15*n1+10) +
   n2^3*(6*n2^2-15*n2+10) +
   n3^3*(6*n3^2-15*n3+10) +
   n4^3*(6*n4^2-15*n4+10);
g:=n1^2*(1-n1)^2 +
   n2^2*(1-n2)^2 +
   n3^2*(1-n3)^2 +
   n4^2*(1-n4)^2 +
alpha*(
n1^2*n2^2 + n1^2*n3^2 + n1^2*n4^2 +
n2^2*n1^2 + n2^2*n3^2 + n2^2*n4^2 +
n3^2*n1^2 + n3^2*n2^2 + n3^2*n4^2 +
n4^2*n1^2 + n4^2*n2^2 + n4^2*n3^2);
(fa*(1-h) + fb*h + w*g)'

nic = 'epsilon*(cos((0.01*idx)*x-4)*cos((0.007+0.01*idx)*y)
       +cos((0.11+0.01*idx)*x)*cos((0.11+0.01*idx)*y)
       +psi*(cos((0.046+0.001*idx)*x+(0.0405+0.001*idx)*y)
       *cos((0.031+0.001*idx)*x-(0.004+0.001*idx)*y))^2)^2'

cnames = 'rho     ca  cb  alpha w L M'
cvalues= 'sqrt(2) 0.3 0.7 5     1 5 5'

[TensorBuffers]
  # variables
  [c]
  []
  [n1]
  []
  [n2]
  []
  [n3]
  []
  [n4]
  []

  [c_bar]
  []
  [n1_bar]
  []
  [n2_bar]
  []
  [n3_bar]
  []
  [n4_bar]
  []

  [mu_c]
  []
  [mu_n1]
  []
  [mu_n2]
  []
  [mu_n3]
  []
  [mu_n4]
  []

  [mu_c_bar]
  []
  [mu_n1_bar]
  []
  [mu_n2_bar]
  []
  [mu_n3_bar]
  []
  [mu_n4_bar]
  []

  [Mbar_mu_c_bar]
  []

  # constant tensors
  [Lbar] # FFT(M*laplacian)
  []
  [MkappaL2bar] # FFT(-M*kappa*laplacian^2)
  []
  [kappaLbar] # FFT(L*kappa*laplacian)
  []

  # postprocessing
  [F]
  []
  [Fgrad_c]
  []
  [Fgrad_n1]
  []
  [Fgrad_n2]
  []
  [Fgrad_n3]
  []
  [Fgrad_n4]
  []
  [bnds]
  []
[]

[TensorComputes]
  [Initialize]
    [c]
      type = ParsedCompute
      buffer = c
      extra_symbols = true
      expression = 'c0+epsilon*(cos(0.105*x)*cos(0.11*y)+(cos(0.13*x)*cos(0.087*y))^2+cos(0.025*x-0.15*y)*cos(0.07*x-0.02*y))'
      constant_names = 'c0 epsilon'
      constant_expressions = '0.5 0.01'
    []
    [Lbar]
      type = ReciprocalLaplacianFactor
      # Mobility is pulled into the chemical potential below
      buffer = Lbar
    []
    [MkappaL2bar]
      type = ReciprocalLaplacianSquareFactor
      factor = -15 # -kappa_c*M
      buffer = MkappaL2bar
    []
    [kappaLbar]
      type = ReciprocalLaplacianFactor
      buffer = kappaLbar
      factor = 15 # kappa_ni*L
    []
    [n1]
      type = ParsedCompute
      buffer = n1
      expression = ${nic}
      extra_symbols = true
      constant_names       = 'idx epsilon psi'
      constant_expressions = '  1     0.1 1.5'
    []
    [n2]
      type = ParsedCompute
      buffer = n2
      expression = ${nic}
      extra_symbols = true
      constant_names       = 'idx epsilon psi'
      constant_expressions = '  2    0.1 1.5'
    []
    [n3]
      type = ParsedCompute
      buffer = n3
      expression = ${nic}
      extra_symbols = true
      constant_names       = 'idx epsilon psi'
      constant_expressions = '  3     0.1 1.5'
    []
    [n4]
      type = ParsedCompute
      buffer = n4
      expression = ${nic}
      extra_symbols = true
      constant_names       = 'idx epsilon psi'
      constant_expressions = '  4     0.1 1.5'
    []
  []

  [Solve]
    [mu_c]
      type = ParsedCompute
      buffer = mu_c
      expression = '${fchem}*M'
      constant_names = ${cnames}
      constant_expressions = ${cvalues}
      derivatives = c
      inputs = 'c n1 n2 n3 n4'
    []

    [mu_n1]
      type = ParsedCompute
      buffer = mu_n1
      expression = '${fchem}*(-L)'
      constant_names = ${cnames}
      constant_expressions = ${cvalues}
      derivatives = n1
      inputs = 'c n1 n2 n3 n4'
    []
    [mu_n2]
      type = ParsedCompute
      buffer = mu_n2
      expression = '${fchem}*(-L)'
      constant_names = ${cnames}
      constant_expressions = ${cvalues}
      derivatives = n2
      inputs = 'c n1 n2 n3 n4'
    []
    [mu_n3]
      type = ParsedCompute
      buffer = mu_n3
      expression = '${fchem}*(-L)'
      constant_names = ${cnames}
      constant_expressions = ${cvalues}
      derivatives = n3
      inputs = 'c n1 n2 n3 n4'
    []
    [mu_n4]
      type = ParsedCompute
      buffer = mu_n4
      expression = '${fchem}*(-L)'
      constant_names = ${cnames}
      constant_expressions = ${cvalues}
      derivatives = n4
      inputs = 'c n1 n2 n3 n4'
    []

    [mu_c_bar]
      type = ForwardFFT
      buffer = mu_c_bar
      input = mu_c
    []
    [mu_n1_bar]
      type = ForwardFFT
      buffer = mu_n1_bar
      input = mu_n1
    []
    [mu_n2_bar]
      type = ForwardFFT
      buffer = mu_n2_bar
      input = mu_n2
    []
    [mu_n3_bar]
      type = ForwardFFT
      buffer = mu_n3_bar
      input = mu_n3
    []
    [mu_n4_bar]
      type = ForwardFFT
      buffer = mu_n4_bar
      input = mu_n4
    []

    [Mbar_mu_c_bar]
      type = ParsedCompute
      buffer = Mbar_mu_c_bar
      expression = 'Lbar*mu_c_bar'
      inputs = 'Lbar mu_c_bar'
    []

    [c_bar]
      type = ForwardFFT
      buffer = c_bar
      input = c
    []
    [n1_bar]
      type = ForwardFFT
      buffer = n1_bar
      input = n1
    []
    [n2_bar]
      type = ForwardFFT
      buffer = n2_bar
      input = n2
    []
    [n3_bar]
      type = ForwardFFT
      buffer = n3_bar
      input = n3
    []
    [n4_bar]
      type = ForwardFFT
      buffer = n4_bar
      input = n4
    []
  []

  [Postprocess]
    [Fgrad_c]
      type = FFTGradientSquare
      buffer = Fgrad_c
      input = c
      factor = 1.5 # kappa/2
    []
    [Fgrad_n1]
      type = FFTGradientSquare
      buffer = Fgrad_n1
      input = n1
      factor = 1.5 # kappa/2
    []
    [Fgrad_n2]
      type = FFTGradientSquare
      buffer = Fgrad_n2
      input = n2
      factor = 1.5 # kappa/2
    []
    [Fgrad_n3]
      type = FFTGradientSquare
      buffer = Fgrad_n3
      input = n3
      factor = 1.5 # kappa/2
    []
    [Fgrad_n4]
      type = FFTGradientSquare
      buffer = Fgrad_n4
      input = n4
      factor = 1.5 # kappa/2
    []
    [F]
      type = ParsedCompute
      buffer = F
      expression = '${fchem} + Fgrad_c + Fgrad_n1 + Fgrad_n2 + Fgrad_n3 + Fgrad_n4'
      constant_names =  ${cnames}
      constant_expressions = ${cvalues}
      inputs = 'c n1 n2 n3 n4 Fgrad_c Fgrad_n1 Fgrad_n2 Fgrad_n3 Fgrad_n4'
    []
    [bnds]
      type = ParsedCompute
      buffer = bnds
      expression = 'n1^2 + n2^2 + n3^2 + n4^2'
      inputs = 'n1 n2 n3 n4'
    []
  []
[]

[TensorSolver]
  type = SecantSolver
  substeps = 1
  max_iterations = 1000
  # damping = 0.75
  relative_tolerance = 1e-6
  absolute_tolerance = 1e-6
  buffer = 'c n1 n2 n3 n4'
  dt_epsilon = 1e-7
  reciprocal_buffer = 'c_bar n1_bar n2_bar n3_bar n4_bar'
  linear_reciprocal = 'MkappaL2bar kappaLbar kappaLbar kappaLbar kappaLbar'
  nonlinear_reciprocal = 'Mbar_mu_c_bar mu_n1_bar mu_n2_bar mu_n3_bar mu_n4_bar'
  # verbose = true
[]

[Postprocessors]
  [min_c]
    type = TensorExtremeValuePostprocessor
    buffer = c
    value_type = MIN
    execute_on = 'TIMESTEP_END'
  []
  [max_c]
    type = TensorExtremeValuePostprocessor
    buffer = c
    value_type = MAX
    execute_on = 'TIMESTEP_END'
  []
  [F]
    type = TensorIntegralPostprocessor
    buffer = F
  []
[]

[Problem]
  type = TensorProblem
[]

[Executioner]
  type = Transient
  # num_steps = 2
  end_time = 1e+5
  [TimeStepper]
    type = TensorSolveIterationAdaptiveDT
    dt = 0.1
    max_iterations = 500
    min_iterations = 300
    growth_factor = 1.1
    cutback_factor = 0.9
  []
[]

[Outputs]
  csv = true
  perf_graph = true
  execute_on = 'TIMESTEP_END'
[]

[TensorOutputs]
  [c]
    type = XDMFTensorOutput
    buffer = 'c bnds'
    output_mode = 'NODE NODE'
    enable_hdf5 = true
  []
[]

(test/tests/gradient/gradient_square.i)

[Domain]
  dim = 3
  nx = 40
  ny = 40
  nz = 40
  xmax = ${fparse pi*2}
  ymax = ${fparse pi*4}
  zmax = ${fparse pi*6}
  mesh_mode = DUMMY
  device_names = cpu
[]

[TensorBuffers]
  [s]
  []
  [grad_sq]
  []
  [c2]
  []
  [diff]
  []
[]

[TensorComputes]
  [Initialize]
    [sin]
      type = ParsedCompute
      buffer = s
      extra_symbols = true
      expression = 'sin(x)+sin(y)+sin(z)'
    []
    [cos2]
      type = ParsedCompute
      buffer = c2
      extra_symbols = true
      expression = 'cos(x)^2+cos(y)^2+cos(z)^2'
    []
    [grad_sq]
      type = FFTGradientSquare
      buffer = grad_sq
      input = s
    []
    [diff]
      type = ParsedCompute
      buffer = diff
      inputs = 'grad_sq c2'
      expression = 'abs(grad_sq - c2)'
    []
  []
[]

[Postprocessors]
  [diff]
    type = TensorIntegralPostprocessor
    buffer = diff
  []
[]

[Problem]
  type = TensorProblem
[]

[Executioner]
  type = Transient
  num_steps = 1
[]

[Outputs]
  csv = true
[]

(test/tests/solvers/coupled.i)

#
# Simple Cahn-Hilliard solve on a 2D grid.
#

[Domain]
  dim = 2
  nx = 150
  ny = 150
  xmax = '${fparse pi*2}'
  ymax = '${fparse pi*2}'
  mesh_mode = DUMMY
[]

[GlobalParams]
  constant_names = 'A B'
  constant_expressions = '1 3.5'
[]

[TensorComputes]
  [Initialize]
    [u]
      type = ParsedCompute
      buffer = u
      extra_symbols = true
      expression = 'sin(x)*sin(y)'
      expand = REAL
    []
    [v]
      type = ParsedCompute
      buffer = v
      extra_symbols = true
      expression = 'cos(x)*cos(y)'
      expand = REAL
    []
    [zero]
      type = ConstantReciprocalTensor
      buffer = zero
    []

    # precompute fixed factors for the solve
    [D1]
      type = ReciprocalLaplacianFactor
      factor = 1e-2
      buffer = D1
    []
    [D2]
      type = ReciprocalLaplacianFactor
      factor = 1e-3
      buffer = D2
    []
  []

  [Solve]
    [u_bar]
      type = ForwardFFT
      buffer = u_bar
      input = u
    []
    [v_bar]
      type = ForwardFFT
      buffer = v_bar
      input = v
    []
  []
[]

[TensorSolver]
  type = AdamsBashforthMoultonCoupled
  buffer = 'u v'
  reciprocal_buffer = 'u_bar v_bar'
  linear_reciprocal = 'D1 D1'
  linear_offdiag_cols = '0 1'
  linear_offdiag_rows = '1 0'
  linear_offdiag = 'D2 D2'
  nonlinear_reciprocal = 'zero zero'
  substeps = ${ss}
  corrector_steps = ${cs}
  predictor_order = ${order}
  corrector_order = ${order}
[]

[Problem]
  type = TensorProblem
[]

[Postprocessors]
  [u_min]
    type = TensorExtremeValuePostprocessor
    buffer = u
    value_type = MIN
  []
  [u_max]
    type = TensorExtremeValuePostprocessor
    buffer = u
    value_type = MAX
  []
  [v_min]
    type = TensorExtremeValuePostprocessor
    buffer = v
    value_type = MIN
  []
  [v_max]
    type = TensorExtremeValuePostprocessor
    buffer = v
    value_type = MAX
  []
  [U]
    type = TensorIntegralPostprocessor
    buffer = u
  []
  [V]
    type = TensorIntegralPostprocessor
    buffer = v
  []
[]

[Executioner]
  type = Transient
  num_steps = 25
  dt = 10
[]

[Outputs]
  file_base = coupled_${ss}_${cs}_${order}
  csv = true
[]

(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/postprocessors/postprocessors.i)

[Domain]
  dim = 2
  nx = 40
  ny = 40
  xmax = 2
  ymax = 3
  mesh_mode = DUMMY
[]

[TensorBuffers]
  [c]
  []
  [c_bar]
  []
[]

[TensorComputes]
  [Initialize]
    [c]
      type = ParsedCompute
      buffer = c
      extra_symbols = true
      expression = -x+y+0.3
    []
    [c_bar]
      type = ForwardFFT
      buffer = c_bar
      input = c
    []
    [u]
      type = ConstantTensor
      buffer = u
      real = 0
    []
  []

  [Solve]
    [root]
      [test]
        type = ForwardFFT
        buffer = u_bar
        input = u
      []
    []
  []
[]

[TensorSolver]
  type = ForwardEulerSolver
  time_derivative_reciprocal = c_bar
  buffer = u
  reciprocal_buffer = u_bar
  substeps = 10
[]

[Postprocessors]
  [min_c]
    type = TensorExtremeValuePostprocessor
    buffer = c
    value_type = MIN
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [max_c]
    type = TensorExtremeValuePostprocessor
    buffer = c
    value_type = MAX
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [avg_c]
    type = TensorAveragePostprocessor
    buffer = c
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [int_c]
    type = TensorIntegralPostprocessor
    buffer = c
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [int_c_bar]
    type = ReciprocalIntegral
    buffer = c_bar
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [count]
    type = ComputeGroupExecutionCount
  []
[]

[Problem]
  type = TensorProblem
[]

[Executioner]
  type = Transient
  num_steps = 0
[]

[Outputs]
  csv = true
[]

(test/tests/postprocessors/interface_velocity.i)

[Domain]
  dim = 2
  nx = 10
  ny = 2
  xmax = '${fparse pi*4}'
  mesh_mode = DUMMY
[]

[TensorComputes]
  [Solve]
    [c]
      type = ParsedCompute
      buffer = c
      extra_symbols = true
      expression = sin(x+0.2*t)
      expand = REAL
    []
  []
[]

[Postprocessors]
  [v]
    type = TensorInterfaceVelocityPostprocessor
    buffer = c
  []
[]

[Problem]
  type = TensorProblem
[]

[Executioner]
  type = Transient
  num_steps = 10
  dt = 0.01
[]

[Outputs]
  csv = true
[]

(test/tests/lbm/stacking.i)

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

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

[TensorBuffers]
  [ux]
    type=LBMTensorBuffer
    buffer_type = ms
  []
  [uy]
    type=LBMTensorBuffer
    buffer_type = ms
  []
  [u]
    type=LBMTensorBuffer
    buffer_type = mv
  []
[]

[TensorComputes/Initialize]
  [velocity_x]
    type = ParsedCompute
    buffer = ux
    expression = '0.1*sin(x/(2*pi*4))*cos(y/(2*pi*4))'
    extra_symbols=true
  []
  [velocity_y]
    type = ParsedCompute
    buffer = uy
    expression = '-0.1*cos(x/(2*pi*4))*sin(y/(2*pi*4))'
    extra_symbols=true
  []
  [u_stack]
    type=LBMStackTensors
    buffer=u
    inputs='ux uy'
  []
[]

[Problem]
  type = LatticeBoltzmannProblem
  substeps = 1
[]

[Postprocessors]
  [velocity_x_min]
    type = TensorExtremeValuePostprocessor
    buffer = ux
    value_type = MIN
  []
  [velocity_x_max]
    type = TensorExtremeValuePostprocessor
    buffer = ux
    value_type = MAX
  []
  [velocity_y_min]
    type = TensorExtremeValuePostprocessor
    buffer = uy
    value_type = MIN
  []
  [velocity_y_max]
    type = TensorExtremeValuePostprocessor
    buffer = uy
    value_type = MAX
  []
  [u_stack_min]
    type = TensorExtremeValuePostprocessor
    buffer = u
    value_type = MIN
  []
  [u_stack_max]
    type = TensorExtremeValuePostprocessor
    buffer = u
    value_type = MAX
  []
[]

[Executioner]
  type = Transient
  num_steps = 2
[]

[Outputs]
  file_base = stacking
  csv = true
[]

(examples/cahn_hilliard/cahnhilliard2.i)

#
# The same simple Cahn-Hilliard solve as cahnhilliard.i, but on a 3D grid
# and using the faster TensorOutputs system.
#

[Domain]
  dim = 3
  nx = 200
  ny = 200
  nz = 200
  xmax = ${fparse pi*8}
  ymax = ${fparse pi*8}
  zmax = ${fparse pi*8}

  # run on a CUDA device (adjust this to `cpu` if not available)
  device_names = 'cuda'

  # create a single element dummy mesh. Output will use the custom XDMF output
  # in the `TensorOutputs` system.
  mesh_mode = DUMMY
[]

[TensorBuffers]
  [c]
  []
  [cbar]
  []
  [mu]
  []
  [mubar]
  []
  [Mbarmubar]
  []
  # constant tensors
  [Mbar]
  []
  [kappabarbar]
  []
[]

[TensorOutputs]
  # the TensorOutouts system supports asynchronous threaded output.
  # for GOU calculations a copy of the solution fields is moved to the CPU,
  # and while the output files are written the next time step is already
  # starting to compute.
  [xdmf]
    type = XDMFTensorOutput
    buffer = 'c mu'
    enable_hdf5 = true
  []
[]

[TensorComputes]
  [Initialize]
    [c]
      type = RandomTensor
      buffer = c
      min = 0.44
      max = 0.56
    []
    [Mbar]
      type = ReciprocalLaplacianFactor
      factor = 0.2 # Mobility
      buffer = Mbar
    []
    [kappabarbar]
      type = ReciprocalLaplacianSquareFactor
      factor = -0.001 # kappa
      buffer = kappabarbar
    []
  []

  [Solve]
    [mu]
      type = ParsedCompute
      buffer = mu
      expression = '0.1*c^2*(c-1)^2'
      derivatives = c
      inputs = c
    []
    [mubar]
      type = ForwardFFT
      buffer = mubar
      input = mu
    []
    [Mbarmubar]
      type = ParsedCompute
      buffer = Mbarmubar
      expression = 'Mbar*mubar'
      inputs = 'Mbar mubar'
    []
    [cbar]
      type = ForwardFFT
      buffer = cbar
      input = c
    []
  []
[]

[TensorSolver]
  type = AdamsBashforthMoulton
  buffer = c
  reciprocal_buffer = cbar
  linear_reciprocal = kappabarbar
  nonlinear_reciprocal = Mbarmubar
  substeps = 1000
[]

[Postprocessors]
  [min_c]
    type = TensorExtremeValuePostprocessor
    buffer = c
    value_type = MIN
    execute_on = 'TIMESTEP_END'
  []
  [max_c]
    type = TensorExtremeValuePostprocessor
    buffer = c
    value_type = MAX
    execute_on = 'TIMESTEP_END'
  []
  [C]
    type = TensorIntegralPostprocessor
    buffer = c
    execute_on = 'TIMESTEP_END'
  []
  [cavg]
    type = TensorAveragePostprocessor
    buffer = c
    execute_on = 'TIMESTEP_END'
  []
[]

[Problem]
  type = TensorProblem
[]

[Executioner]
  type = Transient
  num_steps = 100
  [TimeStepper]
    type = IterationAdaptiveDT
    growth_factor = 1.8
    dt = 0.1
  []
  dtmax = 1000
[]

[Outputs]
  csv = true
  perf_graph = true
  execute_on = 'TIMESTEP_END'
[]

(test/tests/tensor_compute/group.i)

[Domain]
  dim = 3
  nx = 128
  ny = 128
  nz = 128
  xmax = ${fparse pi*4}
  ymax = ${fparse pi*4}
  zmax = ${fparse pi*4}
  mesh_mode = DUMMY
[]

[TensorBuffers]
  # phase field
  [c]
  []
  [cbar]
  []
  [mu]
  []
  [mubar]
  []
  [Mbarmubar]
  []

  # mechanics
  [disp_x]
  []
  [disp_y]
  []
  [disp_z]
  []
  [mumechbar]
  []
  [mumech]
  []

  # constant tensors
  [Mbar]
  []
  [kappabarbar]
  []
[]

[TensorOutputs]
  [xdmf]
    type = XDMFTensorOutput
    buffer = 'c disp_x disp_y disp_z mu mumech'
    output_mode = 'Node Node Node Node Cell Cell'
    enable_hdf5 = true
  []
[]

[TensorComputes]
  [Initialize]
    [c]
      type = RandomTensor
      buffer = c
      min = 0.44
      max = 0.56
    []
    [disp_x]
      type = RandomTensor
      buffer = disp_x
      min = 0
      max = 0
    []
    [disp_y]
      type = RandomTensor
      buffer = disp_y
      min = 0
      max = 0
    []
    [disp_z]
      type = RandomTensor
      buffer = disp_z
      min = 0
      max = 0
    []
    [Mbar]
      type = ReciprocalLaplacianFactor
      factor = 0.2 # Mobility
      buffer = Mbar
    []
    [kappabarbar]
      type = ReciprocalLaplacianSquareFactor
      factor = -0.001 # kappa
      buffer = kappabarbar
    []
  []

  [Solve]
    [mu]
      # chemical potential (real space)
      type = ParsedCompute
      buffer = mu
      expression = '0.1*c^2*(c-1)^2' # + c*sin(x/2)*0.005'
      extra_symbols = true
      derivatives = c
      inputs = c
    []
    [mubar]
      # chemical potential (reciprocal space)
      type = ForwardFFT
      buffer = mubar
      input = mu
    []
    [mumechbar]
      # mechanical chemical potential (reciprocal space)
      type = FFTElasticChemicalPotential
      buffer = mumechbar
      cbar = cbar
      displacements = 'disp_x disp_y disp_z'
      lambda = 100
      mu = 50
      e0 = 0.02
    []
    [mumech]
      # chemical potential (reciprocal space)
      type = InverseFFT
      buffer = mumech
      input = mumechbar
    []

    [Mbarmubar]
      type = ParsedCompute
      buffer = Mbarmubar
      expression = 'Mbar*(mubar+mumechbar)'
      inputs = 'Mbar mubar mumechbar'
    []
    [cbar]
      type = ForwardFFT
      buffer = cbar
      input = c
    []

    [qsmech]
      type = FFTQuasistaticElasticity
      displacements = 'disp_x disp_y disp_z'
      cbar = cbar
      lambda = 100
      mu = 50
      e0 = 0.02
    []

    [group]
      type = ComputeGroup
      computes = 'cbar mumech mu mubar'
    []
  []
[]

[TensorSolver]
  type = AdamsBashforthMoulton
  buffer = c
  reciprocal_buffer = cbar
  linear_reciprocal = kappabarbar
  nonlinear_reciprocal = Mbarmubar
[]

[Postprocessors]
  [min_c]
    type = TensorExtremeValuePostprocessor
    buffer = c
    value_type = MIN
    execute_on = 'TIMESTEP_END'
  []
  [max_c]
    type = TensorExtremeValuePostprocessor
    buffer = c
    value_type = MAX
    execute_on = 'TIMESTEP_END'
  []

  [min_disp_x]
    type = TensorExtremeValuePostprocessor
    buffer = disp_x
    value_type = MIN
    execute_on = 'TIMESTEP_END'
  []
  [max_disp_x]
    type = TensorExtremeValuePostprocessor
    buffer = disp_x
    value_type = MAX
    execute_on = 'TIMESTEP_END'
  []
  [min_disp_y]
    type = TensorExtremeValuePostprocessor
    buffer = disp_y
    value_type = MIN
    execute_on = 'TIMESTEP_END'
  []
  [max_disp_y]
    type = TensorExtremeValuePostprocessor
    buffer = disp_y
    value_type = MAX
    execute_on = 'TIMESTEP_END'
  []
  [min_disp_z]
    type = TensorExtremeValuePostprocessor
    buffer = disp_z
    value_type = MIN
    execute_on = 'TIMESTEP_END'
  []
  [max_disp_z]
    type = TensorExtremeValuePostprocessor
    buffer = disp_z
    value_type = MAX
    execute_on = 'TIMESTEP_END'
  []

  [C]
    type = TensorIntegralPostprocessor
    buffer = c
    execute_on = 'TIMESTEP_END'
  []
  [cavg]
    type = TensorAveragePostprocessor
    buffer = c
    execute_on = 'TIMESTEP_END'
  []
[]

[Problem]
  type = TensorProblem
  spectral_solve_substeps = 1000
  print_debug_output = true
[]

[Executioner]
  type = Transient
  num_steps = 100
  [TimeStepper]
    type = IterationAdaptiveDT
    growth_factor = 1.8
    dt = 0.1
  []
  dtmax = 1000
[]

[Outputs]
  csv = true
  perf_graph = true
  execute_on = 'TIMESTEP_END'
[]

(benchmarks/01_spinodal_decomposition/1b.i)

[Domain]
  dim = 2
  nx = 200
  ny = 200

  xmin = -10
  ymin = -10
  xmax = 210
  ymax = 210

  device_names = 'cuda'

  mesh_mode = DOMAIN
[]


[TensorBuffers]
  [domain]
    map_to_aux_variable = domain
  []
  [c]
    map_to_aux_variable = c
  []
  [cbar]
  []
  [mu]
    # map_to_aux_variable = mu
  []
  [mubar]
  []
  [Mbarmubar]
  []
  # constant tensors
  [Mbar]
  []
  [kappabarbar]
  []
  # postprocessing
  [F]
  []
  [Fgrad]
  []
[]

[TensorComputes]
  [Initialize]
    [domain]
      type = ParsedCompute
      buffer = domain
      extra_symbols = true
      expression = 'w:=3;dx:=max(tanh(-x/w),tanh((x-200)/w)); dy:=max(tanh(-y/w),tanh((y-200)/w));1-0.999999999*max(dx,dy)'
    []
    [c]
      type = ParsedCompute
      buffer = c
      extra_symbols = true
      expression = 'c0+epsilon*(cos(0.105*x)*cos(0.11*y)+(cos(0.13*x)*cos(0.087*y))^2+cos(0.025*x-0.15*y)*cos(0.07*x-0.02*y))'
      constant_names = 'c0 epsilon'
      constant_expressions = '0.5 0.01'
    []
    [Mbar]
      type = ReciprocalLaplacianFactor
      factor = 5 # Mobility
      buffer = Mbar
    []
    [kappabarbar]
      type = ReciprocalLaplacianSquareFactor
      factor = -10 # -kappa*M
      buffer = kappabarbar
    []
  []

  [Solve]
    [mu]
      type = ParsedCompute
      buffer = mu
      expression = 'rho_s*(c-c_alpha)^2*(c_beta-c)^2'
      constant_names =       'rho_s c_alpha c_beta'
      constant_expressions = '5     0.3     0.7'
      derivatives = c
      inputs = c
    []
    [mubar]
      type = ForwardFFT
      buffer = mubar
      input = mu
    []
    [Mbarmubar]
      type = ParsedCompute
      buffer = Mbarmubar
      expression = 'Mbar*mubar'
      inputs = 'Mbar mubar'
    []
    [cbar]
      type = ForwardFFT
      buffer = cbar
      input = c
    []
  []

  [Postprocess]
    [Fgrad]
      type = FFTGradientSquare
      buffer = Fgrad
      input = c
      factor = 1 # kappa/2
    []
    [F]
      type = ParsedCompute
      buffer = F
      expression = 'rho_s * (c-c_alpha)^2 * (c_beta-c)^2 + Fgrad'
      constant_names =       'rho_s c_alpha c_beta'
      constant_expressions = '5     0.3     0.7'
      inputs = 'c Fgrad'
    []
  []
[]

[UserObjects]
  [terminator]
    type = Terminator
    expression = change<1e-4
  []
[]

[TensorTimeIntegrators]
  [c]
    type = FFTSemiImplicit
    buffer = c
    history_size = 1
    reciprocal_buffer = cbar
    linear_reciprocal = kappabarbar
    nonlinear_reciprocal = Mbarmubar
  []
[]

[AuxVariables]
  # [mu]
  #   family = MONOMIAL
  #   order = CONSTANT
  # []
  [c]
    # family = MONOMIAL
    # order = CONSTANT
  []
  [domain]
  []
[]

[Postprocessors]
  [min_c]
    type = TensorExtremeValuePostprocessor
    buffer = c
    value_type = MIN
    execute_on = 'TIMESTEP_END'
  []
  [max_c]
    type = TensorExtremeValuePostprocessor
    buffer = c
    value_type = MAX
    execute_on = 'TIMESTEP_END'
  []
  [F]
    type = TensorIntegralPostprocessor
    buffer = F
  []
  [change]
    type = TensorIntegralChangePostprocessor
    buffer = c
  []
[]

[Problem]
  type = TensorProblem
  spectral_solve_substeps = 1000
[]

[Executioner]
  type = Transient
  num_steps = 1 # 1000
  [TimeStepper]
    type = IterationAdaptiveDT
    growth_factor = 1.1
    dt = 1
  []
  dtmax = 300
[]

[Outputs]
  exodus = true
  csv = true
  perf_graph = true
  execute_on = 'TIMESTEP_END'
[]

(test/tests/histogram/test.i)

[Domain]
  dim = 3
  nx = 10
  ny = 10
  nz = 10
  mesh_mode = DUMMY
[]

[TensorBuffers]
  [c]
  []
[]

[TensorComputes]
  [Initialize]
    [c]
      type = ParsedCompute
      buffer = c
      extra_symbols = true
      expression = '0.1*x^2+0.2*y^2+0.3*z^2'
    []
  []
[]

[VectorPostprocessors]
  [hist]
    type = TensorHistogram
    buffer = c
    bins = 20
    min = 0
    max = 1
    execute_on = 'TIMESTEP_END'
  []
[]

[Problem]
  type = TensorProblem
[]

[Executioner]
  type = Transient
  num_steps = 1
[]

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

(test/tests/tensor_compute/test.i)

[Mesh]
  type = UniformTensorMesh
  dim = 2
  nx = 100
  ny = 100
  nz = 100
  xmax = ${fparse pi*4}
  ymax = ${fparse pi*4}
  zmax = ${fparse pi*4}
[]

[TensorBuffers]
  [eta]
  []
  [eta_bar]
  []
  [f]
  []
  [fbar]
  []
  [kappa_k2]
  []
[]

[TensorComputes]
  [Initialize]
    [eta]
      type = ParsedTensor
      buffer = eta
      function = 'sin(x)+sin(y)+sin(z)'
    []
    [kappa_k2]
      type = ReciprocalLaplacianFactor
      factor = 0.2
      buffer = kappa_k2
    []
  []

  [Solve]
    [f]
      type = ParsedCompute
      buffer = f
      expression = '0.1*(eta+2)^2*(eta-2)^2'
      derivatives = eta
      inputs = eta
    []
    [fbar]
      type = ForwardFFT
      buffer = fbar
      input = f
    []
    [eta_bar]
      type = ForwardFFT
      buffer = eta_bar
      input = eta
    []
  []
[]

[TensorTimeIntegrators]
  [eta]
    type = FFTSemiImplicit
    buffer = eta
    reciprocal_buffer = eta_bar
    linear_reciprocal = kappa_k2
    nonlinear_reciprocal = fbar
  []
[]

[AuxVariables]
  [eta]
  []
  [f]
  []
[]

[AuxKernels]
  [eta]
    type = ProjectTensorAux
    buffer = eta
    variable = eta
    execute_on = TIMESTEP_END
  []
  [f]
    type = ProjectTensorAux
    buffer = f
    variable = f
    execute_on = TIMESTEP_END
  []
[]

[Postprocessors]
  [min_eta]
    type = ElementExtremeValue
    variable = eta
    value_type = MIN
    execute_on = 'TIMESTEP_END'
  []
  [max_eta]
    type = ElementExtremeValue
    variable = eta
    value_type = MAX
    execute_on = 'TIMESTEP_END'
  []
  [F]
    type = ElementIntegralVariablePostprocessor
    variable = f
    execute_on = 'TIMESTEP_END'
  []
  [Eta]
    type = ElementIntegralVariablePostprocessor
    variable = eta
    execute_on = 'TIMESTEP_END'
  []
[]

[Problem]
  type = TensorProblem
[]

[Executioner]
  type = Transient
  num_steps = 100
  dt = 0.1
[]

[Outputs]
  exodus = true
  csv = true
  perf_graph = true
  execute_on = 'TIMESTEP_END'
[]

(benchmarks/02_oswald_ripening/simple.i)

[Domain]
  dim = 2
  nx = 4
  ny = 4
  xmax = 1
  ymax = 1

  device_names = 'cuda:0'

  mesh_mode = DOMAIN
[]

[TensorBuffers]
  # variables
  [n1]
  []
  [n2]
  []
  [n3]
  []

  [n1_bar]
  []
  [n2_bar]
  []
  [n3_bar]
  []

  [mu_n1]
  []
  [mu_n2]
  []
  [mu_n3]
  []

  [mu_n1_bar]
  []
  [mu_n2_bar]
  []
  [mu_n3_bar]
  []

  [Lbar] # zero
  []
[]

[TensorComputes]
  [Initialize]
    [Lbar]
      type = ConstantReciprocalTensor
      buffer = Lbar
      real = 0
      imaginary = 0
    []

    [n1]
      type = ConstantTensor
      buffer = n1
      real = 1
    []
    [n2]
      type = ConstantTensor
      buffer = n2
      real = 1
    []
    [n3]
      type = ConstantTensor
      buffer = n3
      real = 1
    []

    [mu_n1]
      type = ConstantTensor
      buffer = mu_n1
      real = 1
    []
  []

  [Solve]
    [mu_n2]
      type = ParsedCompute
      buffer = mu_n2
      expression = 'n3'
      inputs = n3
    []
    [mu_n3]
      type = ParsedCompute
      buffer = mu_n3
      expression = 'n3*n3'
      inputs = n3
    []

    [mu_n1_bar]
      type = ForwardFFT
      buffer = mu_n1_bar
      input = mu_n1
    []
    [mu_n2_bar]
      type = ForwardFFT
      buffer = mu_n2_bar
      input = mu_n2
    []
    [mu_n3_bar]
      type = ForwardFFT
      buffer = mu_n3_bar
      input = mu_n3
    []

    [n1_bar]
      type = ForwardFFT
      buffer = n1_bar
      input = n1
    []
    [n2_bar]
      type = ForwardFFT
      buffer = n2_bar
      input = n2
    []
    [n3_bar]
      type = ForwardFFT
      buffer = n3_bar
      input = n3
    []
  []
[]

[TensorSolver]
  # type = SecantSolver
  type = BroydenSolver
  substeps = 1
  max_iterations = 10
  damping = 0.5
  buffer = 'n1 n2 n3'
  tolerance = 1e-5
  dt_epsilon = 1e-5
  reciprocal_buffer = 'n1_bar n2_bar n3_bar'
  linear_reciprocal = 'Lbar   Lbar   Lbar'
  nonlinear_reciprocal = 'mu_n1_bar mu_n2_bar mu_n3_bar'
  verbose = true
[]

[Postprocessors]
  [n1]
    type = TensorAveragePostprocessor
    buffer = n1
    execute_on = 'TIMESTEP_END'
  []
  [n2]
    type = TensorAveragePostprocessor
    buffer = n2
    execute_on = 'TIMESTEP_END'
  []
  [n3]
    type = TensorAveragePostprocessor
    buffer = n3
    execute_on = 'TIMESTEP_END'
  []
[]

[Problem]
  type = TensorProblem
[]

[Executioner]
  type = Transient
  num_steps = 1
  dt = 1e-2
[]

[Outputs]
  perf_graph = true
  execute_on = 'TIMESTEP_END'
[]

(examples/lbm/Rayleigh-Benard/rayleigh-benard.i)

Nx = 400
Ny = 200

TH = 1.1
TC = 1.0

frequency = '${Nx}/10.0'
amplitude = '${Ny}/100.0'

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

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

!include buffers.i

[TensorComputes]

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

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

    [temperature]
      type = ParsedCompute
      buffer = T
      expression = 'a:=abs(y - sin(x / (${frequency} * pi)) * ${amplitude}) + y - sin(x / (${frequency} * pi)) * ${amplitude};
                    b:= a / (a + 1.0e-14);
                    ${TC} * b - b *${TH} + ${TH}'
      extra_symbols = true
    []

    [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
    []
  []

  !include solve.i
  !include boundary.i
[]

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

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

[Problem]
  type = LatticeBoltzmannProblem
  scalar_constant_names = 'rho0 T_0  T_C  T_H  tau_f tau_T  g'
  scalar_constant_values = '1.0  1.05  1.0  1.1  0.7   0.7   0.0001'
  substeps = 100
  print_debug_output = true
[]

[Executioner]
  type = Transient
  num_steps = 10000
[]

(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/kks/KKS_libtorch.i)

#
# Kim-Kim-Suzuki with Gibbs energy supplied by a torch model, solved on a 2D grid.
#


# Constants for Initial Conditions
r = 30
l = 4.2

# Initial condition function for order parameter
eta_IC = '0.5*(1-tanh(2*(sqrt(x^2+y^2)-${r})/${l}))'

# Phase-field model parameters
kappa_eta = 5
w = 1
M = 5
L = 5

# Expressions for switching function and bulk Gibbs energy
h_eta = 'eta^3*(6*eta^2-15*eta+10)'


[Domain]
    dim = 2
    nx = 50
    ny = 50

    xmin = -50
    xmax = 50
    ymin = -50
    ymax = 50

    # automatically create a matching mesh
    mesh_mode = DUMMY
[]

[TensorComputes]
    [Initialize]
        [c_IC]
            type = ParsedCompute
            buffer = c
            expression = '0.7 + (0.3-0.6)*${eta_IC}'
            extra_symbols = 'true'
        []
        [eta_IC]
            type = ParsedCompute
            buffer = eta
            expression = '${eta_IC}'
            extra_symbols = 'true'
        []
        [psi_init]
            type = ConstantTensor
            buffer = psi
            real = 1
        []

        [M]
            type = ConstantTensor
            buffer = M
            real = ${M}
        []
        [L]
            type = ConstantTensor
            buffer = L
            real = ${L}
        []
        [L_kappa]
            type = ReciprocalLaplacianFactor
            buffer = L_kappa
            factor = ${fparse  ${L} * ${kappa_eta} }
        []
        [h_eta_IC]
            type = ParsedCompute
            buffer = h_eta
            expression = '${h_eta}'
            inputs = eta
        []
        [G_func_IC]
            type = LibtorchGibbsEnergy
            buffer = 'G'
            phase_fractions = 'h_eta'
            concentrations = 'c'
            domega_detas = 'dG_dh'
            chem_pots = 'mu'
            libtorch_model_file = 'marlin:libtorch_gibbs_energy/torch_NN_gibbs_model.pt'
        []
    []
    [Solve]
        [h_eta]
            type = ParsedCompute
            buffer = h_eta
            expression = '${h_eta}'
            inputs = eta
        []
        [G_func]
            type = LibtorchGibbsEnergy
            buffer = 'G'
            phase_fractions = 'h_eta'
            concentrations = 'c'
            domega_detas = 'dG_dh'
            chem_pots = 'mu'
            libtorch_model_file = 'marlin:libtorch_gibbs_energy/torch_NN_gibbs_model.pt'
        []
        [dG_deta]
            type = ParsedCompute
            buffer = 'dG_deta'
            inputs = 'eta dG_dh'
            expression = 'dG_dh * ${h_eta} + ${w} * eta^2 * (1-eta^2)^2'
            derivatives = 'eta'
        []

        [etabar]
            type = ForwardFFT
            buffer = etabar
            input = eta
        []
        [AC_bulk]
            type = ReciprocalAllenCahn
            L = L
            buffer = AC_bulk
            dF_chem_deta = dG_deta
            psi = psi
        []
        [NL_eta]
            type = ParsedCompute
            buffer = NL_eta
            expression = 'AC_bulk '
            inputs = 'AC_bulk'
        []
        [cbar]
            type = ForwardFFT
            buffer = cbar
            input = c
        []
        [div_J]
            type = ReciprocalMatDiffusion
            buffer = 'div_J'
            chemical_potential = mu
            mobility = M
            psi = psi
        []
    []
[]

[TensorSolver]
    type = AdamsBashforthMoulton
    buffer = 'c eta'
    reciprocal_buffer = 'cbar etabar'
    linear_reciprocal = '0 L_kappa'
    nonlinear_reciprocal = 'div_J NL_eta'
    substeps = 1e3
    predictor_order = 3
[]

[TensorOutputs]
    [xdmf]
        type = XDMFTensorOutput
        buffer = 'eta c mu psi dG_deta dG_dh G'
        enable_hdf5 = true
        transpose = false
    []
[]

[Executioner]
    type = Transient
    dt = 0.1
    num_steps = 10
[]

[Outputs]
    csv = true
    perf_graph = true
    execute_on = 'INITIAL TIMESTEP_END'
[]

(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/tensor_compute/parallel_roundtrip.i)

[Domain]
  # Test parallel FFT round-trip with slab decomposition
  device_names = "cpu cpu cpu"
  device_weights = "1 1 1"

  parallel_mode = FFT_SLAB

  dim = 2
  nx = 128
  ny = 128
  xmax = ${fparse pi*4}
  ymax = ${fparse pi*4}
[]

[TensorBuffers]
  [eta_gold]
  []
  [eta]
  []
  [eta_bar]
  []
  [eta_roundtrip]
  []
  [diff]
  []
  [zero]
  []
[]

[TensorComputes]
  [Initialize]
    [eta_gold]
      type = ParsedCompute
      buffer = eta_gold
      expression = 'sin(x)+sin(y)+cos(2*x)*sin(3*y)'
      extra_symbols = true
    []
    [eta]
      type = ParsedCompute
      buffer = eta
      expression = eta_gold
      inputs = eta_gold
    []
    [zero]
      type = ConstantReciprocalTensor
      buffer = zero
      real = 0
      imaginary = 0
    []
  []

  [Solve]
    # Test: eta -> FFT -> iFFT -> eta_roundtrip
    # eta_roundtrip should equal eta (within numerical precision)
    [eta_bar]
      type = ForwardFFT
      buffer = eta_bar
      input = eta
    []
    [eta_roundtrip]
      type = InverseFFT
      buffer = eta_roundtrip
      input = eta_bar
    []
  []

  [Postprocess]
    [diff]
      type = ParsedCompute
      buffer = diff
      expression = 'abs(eta - eta_roundtrip) + abs(eta - eta_gold)'
      inputs = 'eta eta_roundtrip eta_gold'
    []
  []
[]

[Postprocessors]
  [max_error]
    type = TensorExtremeValuePostprocessor
    buffer = diff
    value_type = MAX
  []
  [l2_error]
    type = TensorIntegralPostprocessor
    buffer = diff
  []
[]

[TensorSolver]
  type = AdamsBashforthMoulton
  buffer = eta
  reciprocal_buffer = eta_bar
  linear_reciprocal = zero
  nonlinear_reciprocal = zero
[]

[Problem]
  type = TensorProblem
[]

[Executioner]
  type = Transient
  num_steps = 1
[]

[Outputs]
  csv = true
  execute_on = 'INITIAL TIMESTEP_END'
[]

Overview
Example Input File Syntax
Input Parameters
Input Files