TensorProblem

A normal Problem object that adds the ability to perform spectral solves.

TensorProblem is the main problem class in Marlin. It adds the concepts of a Domain, TensorBuffers, TensorComputes, and TensorOutputs.

The TensorProblem object handles the execution of TensorComputes, schedules TensorOutputs, and supports fast mapping for tensor buffers onto conforming meshes using the "map_to_aux_variable" option.

Overview

Coordinator for Marlin tensor simulations: owns the Domain, schedules TensorComputes, and manages TensorOutputs. Supports fast projection of buffers to mesh variables via "map_to_aux_variable".

Example Input File Syntax

[Problem<<<{"href": "../../syntax/Problem/index.html"}>>>]
  type = TensorProblem
[]

Input Parameters

blockANY_BLOCK_ID List of subdomains for kernel coverage and material coverage checks. Setting this parameter is equivalent to setting 'kernel_coverage_block_list' and 'material_coverage_block_list' as well as using 'ONLY_LIST' as the coverage check mode.
Default:ANY_BLOCK_ID
C++ Type:std::vector<SubdomainName>
Controllable:No
Description:List of subdomains for kernel coverage and material coverage checks. Setting this parameter is equivalent to setting 'kernel_coverage_block_list' and 'material_coverage_block_list' as well as using 'ONLY_LIST' as the coverage check mode.
linear_sys_namesThe linear system names
C++ Type:std::vector<LinearSystemName>
Controllable:No
Description:The linear system names
print_debug_outputFalseShow Tensor specific debug outputs
Default:False
C++ Type:bool
Controllable:No
Description:Show Tensor specific debug outputs
regard_general_exceptions_as_errorsFalseIf we catch an exception during residual/Jacobian evaluaton for which we don't have specific handling, immediately error instead of allowing the time step to be cut
Default:False
C++ Type:bool
Controllable:No
Description:If we catch an exception during residual/Jacobian evaluaton for which we don't have specific handling, immediately error instead of allowing the time step to be cut
scalar_constant_namesScalar constant names
C++ Type:std::vector<std::string>
Controllable:No
Description:Scalar constant names
scalar_constant_valuesScalar constant values
C++ Type:std::vector<double>
Unit:(no unit assumed)
Controllable:No
Description:Scalar constant values
solveTrueWhether or not to actually solve the Nonlinear system. This is handy in the case that all you want to do is execute AuxKernels, Transfers, etc. without actually solving anything
Default:True
C++ Type:bool
Controllable:Yes
Description:Whether or not to actually solve the Nonlinear system. This is handy in the case that all you want to do is execute AuxKernels, Transfers, etc. without actually solving anything
spectral_solve_substeps1How many substeps to divide the spectral solve for each MOOSE timestep into.
Default:1
C++ Type:unsigned int
Controllable:No
Description:How many substeps to divide the spectral solve for each MOOSE timestep into.

Optional Parameters

allow_initial_conditions_with_restartFalseTrue to allow the user to specify initial conditions when restarting. Initial conditions can override any restarted field
Default:False
C++ Type:bool
Controllable:No
Description:True to allow the user to specify initial conditions when restarting. Initial conditions can override any restarted field
force_restartFalseEXPERIMENTAL: If true, a sub_app may use a restart file instead of using of using the master backup file
Default:False
C++ Type:bool
Controllable:No
Description:EXPERIMENTAL: If true, a sub_app may use a restart file instead of using of using the master backup file
restart_file_baseFile base name used for restart (e.g. / or /LATEST to grab the latest file available)
C++ Type:FileNameNoExtension
Controllable:No
Description:File base name used for restart (e.g. / or /LATEST to grab the latest file available)

Restart Parameters

allow_invalid_solutionFalseSet to true to allow convergence even though the solution has been marked as 'invalid'
Default:False
C++ Type:bool
Controllable:No
Description:Set to true to allow convergence even though the solution has been marked as 'invalid'
immediately_print_invalid_solutionFalseWhether or not to report invalid solution warnings at the time the warning is produced instead of after the calculation
Default:False
C++ Type:bool
Controllable:No
Description:Whether or not to report invalid solution warnings at the time the warning is produced instead of after the calculation
show_invalid_solution_consoleTrueSet to true to show the invalid solution occurrence summary in console
Default:True
C++ Type:bool
Controllable:No
Description:Set to true to show the invalid solution occurrence summary in console

Solution Validity Control Parameters

boundary_restricted_elem_integrity_checkTrueSet to false to disable checking of boundary restricted elemental object variable dependencies, e.g. are the variable dependencies defined on the selected boundaries?
Default:True
C++ Type:bool
Controllable:No
Description:Set to false to disable checking of boundary restricted elemental object variable dependencies, e.g. are the variable dependencies defined on the selected boundaries?
boundary_restricted_node_integrity_checkTrueSet to false to disable checking of boundary restricted nodal object variable dependencies, e.g. are the variable dependencies defined on the selected boundaries?
Default:True
C++ Type:bool
Controllable:No
Description:Set to false to disable checking of boundary restricted nodal object variable dependencies, e.g. are the variable dependencies defined on the selected boundaries?
check_uo_aux_stateFalseTrue to turn on a check that no state presents during the evaluation of user objects and aux kernels
Default:False
C++ Type:bool
Controllable:No
Description:True to turn on a check that no state presents during the evaluation of user objects and aux kernels
error_on_jacobian_nonzero_reallocationFalseThis causes PETSc to error if it had to reallocate memory in the Jacobian matrix due to not having enough nonzeros
Default:False
C++ Type:bool
Controllable:No
Description:This causes PETSc to error if it had to reallocate memory in the Jacobian matrix due to not having enough nonzeros
fv_bcs_integrity_checkTrueSet to false to disable checking of overlapping Dirichlet and Flux BCs and/or multiple DirichletBCs per sideset
Default:True
C++ Type:bool
Controllable:No
Description:Set to false to disable checking of overlapping Dirichlet and Flux BCs and/or multiple DirichletBCs per sideset
kernel_coverage_block_listList of subdomains for kernel coverage check. The meaning of this list is controlled by the parameter 'kernel_coverage_check' (whether this is the list of subdomains to be checked, not to be checked or not taken into account).
C++ Type:std::vector<SubdomainName>
Controllable:No
Description:List of subdomains for kernel coverage check. The meaning of this list is controlled by the parameter 'kernel_coverage_check' (whether this is the list of subdomains to be checked, not to be checked or not taken into account).
kernel_coverage_checkTRUEControls, if and how a kernel subdomain coverage check is performed. With 'TRUE' or 'ON' all subdomains are checked (the default). Setting 'FALSE' or 'OFF' will disable the check for all subdomains. To exclude a predefined set of subdomains 'SKIP_LIST' is to be used, while the subdomains to skip are to be defined in the parameter 'kernel_coverage_block_list'. To limit the check to a list of subdomains, 'ONLY_LIST' is to be used (again, using the parameter 'kernel_coverage_block_list').
Default:TRUE
C++ Type:MooseEnum
Options:FALSE, TRUE, OFF, ON, SKIP_LIST, ONLY_LIST
Controllable:No
Description:Controls, if and how a kernel subdomain coverage check is performed. With 'TRUE' or 'ON' all subdomains are checked (the default). Setting 'FALSE' or 'OFF' will disable the check for all subdomains. To exclude a predefined set of subdomains 'SKIP_LIST' is to be used, while the subdomains to skip are to be defined in the parameter 'kernel_coverage_block_list'. To limit the check to a list of subdomains, 'ONLY_LIST' is to be used (again, using the parameter 'kernel_coverage_block_list').
material_coverage_block_listList of subdomains for material coverage check. The meaning of this list is controlled by the parameter 'material_coverage_check' (whether this is the list of subdomains to be checked, not to be checked or not taken into account).
C++ Type:std::vector<SubdomainName>
Controllable:No
Description:List of subdomains for material coverage check. The meaning of this list is controlled by the parameter 'material_coverage_check' (whether this is the list of subdomains to be checked, not to be checked or not taken into account).
material_coverage_checkTRUEControls, if and how a material subdomain coverage check is performed. With 'TRUE' or 'ON' all subdomains are checked (the default). Setting 'FALSE' or 'OFF' will disable the check for all subdomains. To exclude a predefined set of subdomains 'SKIP_LIST' is to be used, while the subdomains to skip are to be defined in the parameter 'material_coverage_block_list'. To limit the check to a list of subdomains, 'ONLY_LIST' is to be used (again, using the parameter 'material_coverage_block_list').
Default:TRUE
C++ Type:MooseEnum
Options:FALSE, TRUE, OFF, ON, SKIP_LIST, ONLY_LIST
Controllable:No
Description:Controls, if and how a material subdomain coverage check is performed. With 'TRUE' or 'ON' all subdomains are checked (the default). Setting 'FALSE' or 'OFF' will disable the check for all subdomains. To exclude a predefined set of subdomains 'SKIP_LIST' is to be used, while the subdomains to skip are to be defined in the parameter 'material_coverage_block_list'. To limit the check to a list of subdomains, 'ONLY_LIST' is to be used (again, using the parameter 'material_coverage_block_list').
material_dependency_checkTrueSet to false to disable material dependency check
Default:True
C++ Type:bool
Controllable:No
Description:Set to false to disable material dependency check
skip_nl_system_checkTrueTrue to skip the NonlinearSystem check for work to do (e.g. Make sure that there are variables to solve for).
Default:True
C++ Type:bool
Controllable:No
Description:True to skip the NonlinearSystem check for work to do (e.g. Make sure that there are variables to solve for).

Simulation Checks 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.
default_ghostingFalseWhether or not to use libMesh's default amount of algebraic and geometric ghosting
Default:False
C++ Type:bool
Controllable:No
Description:Whether or not to use libMesh's default amount of algebraic and geometric ghosting
enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Controllable:No
Description:Set the enabled status of the MooseObject.

Advanced Parameters

extra_tag_matricesExtra matrices to add to the system that can be filled by objects which compute residuals and Jacobians (Kernels, BCs, etc.) by setting tags on them. The outer index is for which nonlinear system the extra tag vectors should be added for
C++ Type:std::vector<std::vector<TagName>>
Controllable:No
Description:Extra matrices to add to the system that can be filled by objects which compute residuals and Jacobians (Kernels, BCs, etc.) by setting tags on them. The outer index is for which nonlinear system the extra tag vectors should be added for
extra_tag_solutionsExtra solution vectors to add to the system that can be used by objects for coupling variable values stored in them.
C++ Type:std::vector<TagName>
Controllable:No
Description:Extra solution vectors to add to the system that can be used by objects for coupling variable values stored in them.
extra_tag_vectorsExtra vectors to add to the system that can be filled by objects which compute residuals and Jacobians (Kernels, BCs, etc.) by setting tags on them. The outer index is for which nonlinear system the extra tag vectors should be added for
C++ Type:std::vector<std::vector<TagName>>
Controllable:No
Description:Extra vectors to add to the system that can be filled by objects which compute residuals and Jacobians (Kernels, BCs, etc.) by setting tags on them. The outer index is for which nonlinear system the extra tag vectors should be added for
not_zeroed_tag_vectorsExtra vector tags which the sytem will not zero when other vector tags are zeroed. The outer index is for which nonlinear system the extra tag vectors should be added for
C++ Type:std::vector<std::vector<TagName>>
Controllable:No
Description:Extra vector tags which the sytem will not zero when other vector tags are zeroed. The outer index is for which nonlinear system the extra tag vectors should be added for

Contribution To Tagged Field Data Parameters

identify_variable_groups_in_nlTrueWhether to identify variable groups in nonlinear systems. This affects dof ordering
Default:True
C++ Type:bool
Controllable:No
Description:Whether to identify variable groups in nonlinear systems. This affects dof ordering
ignore_zeros_in_jacobianFalseDo not explicitly store zero values in the Jacobian matrix if true
Default:False
C++ Type:bool
Controllable:No
Description:Do not explicitly store zero values in the Jacobian matrix if true
nl_sys_namesnl0 The nonlinear system names
Default:nl0
C++ Type:std::vector<NonlinearSystemName>
Controllable:No
Description:The nonlinear system names
previous_nl_solution_requiredFalseTrue to indicate that this calculation requires a solution vector for storing the previous nonlinear iteration.
Default:False
C++ Type:bool
Controllable:No
Description:True to indicate that this calculation requires a solution vector for storing the previous nonlinear iteration.
restore_original_nonzero_patternFalseWhether we should reset matrix memory for every Jacobian evaluation. This option is useful if the sparsity pattern is constantly changing and you are using hash table assembly or if you wish to continually restore the matrix to the originally preallocated sparsity pattern computed by relationship managers.
Default:False
C++ Type:bool
Controllable:No
Description:Whether we should reset matrix memory for every Jacobian evaluation. This option is useful if the sparsity pattern is constantly changing and you are using hash table assembly or if you wish to continually restore the matrix to the originally preallocated sparsity pattern computed by relationship managers.
use_hash_table_matrix_assemblyFalseWhether to assemble matrices using hash tables instead of preallocating matrix memory. This can be a good option if the sparsity pattern changes throughout the course of the simulation.
Default:False
C++ Type:bool
Controllable:No
Description:Whether to assemble matrices using hash tables instead of preallocating matrix memory. This can be a good option if the sparsity pattern changes throughout the course of the simulation.
use_nonlinearTrueDetermines whether to use a Nonlinear vs a Eigenvalue system (Automatically determined based on executioner)
Default:True
C++ Type:bool
Controllable:No
Description:Determines whether to use a Nonlinear vs a Eigenvalue system (Automatically determined based on executioner)

Nonlinear System(S) Parameters

near_null_space_dimension0The dimension of the near nullspace
Default:0
C++ Type:unsigned int
Controllable:No
Description:The dimension of the near nullspace
null_space_dimension0The dimension of the nullspace
Default:0
C++ Type:unsigned int
Controllable:No
Description:The dimension of the nullspace
transpose_null_space_dimension0The dimension of the transpose nullspace
Default:0
C++ Type:unsigned int
Controllable:No
Description:The dimension of the transpose nullspace

Null Space Removal Parameters

parallel_barrier_messagingFalseDisplays messaging from parallel barrier notifications when executing or transferring to/from Multiapps (default: false)
Default:False
C++ Type:bool
Controllable:No
Description:Displays messaging from parallel barrier notifications when executing or transferring to/from Multiapps (default: false)
verbose_multiappsFalseSet to True to enable verbose screen printing related to MultiApps
Default:False
C++ Type:bool
Controllable:No
Description:Set to True to enable verbose screen printing related to MultiApps
verbose_restoreFalseSet to True to enable verbose screen printing related to solution restoration
Default:False
C++ Type:bool
Controllable:No
Description:Set to True to enable verbose screen printing related to solution restoration
verbose_setupfalseSet to 'true' to have the problem report on any object created. Set to 'extra' to also display all parameters.
Default:false
C++ Type:MooseEnum
Options:false, true, extra
Controllable:No
Description:Set to 'true' to have the problem report on any object created. Set to 'extra' to also display all parameters.

Verbosity Parameters

Input Files

(examples/cahn_hilliard/cahnhilliard4.i)
(examples/cahn_hilliard/cahnhilliard.i)
(examples/swift_hohenberg/rotating_grain.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/tensor_ics/sineic.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/tensor_compute/backandforth.i)
(test/tests/solvers/diagonal.i)
(test/tests/tensor_compute/smooth_rectangle.i)
(test/tests/tensor_compute/rotating_grain_secant.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/problem/fftproblem.i)
(test/tests/solvers/nl_coupled.i)
(benchmarks/02_oswald_ripening/2a_broyden.i)
(benchmarks/02_oswald_ripening/2a_secant.i)
(test/tests/tensor_compute/parallel_roundtrip_3d.i)
(test/tests/gradient/gradient_square.i)
(test/tests/solvers/coupled.i)
(test/tests/neml2/scalar.i)
(test/tests/postprocessors/interface_velocity.i)
(test/tests/postprocessors/postprocessors.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)
(test/tests/tensor_compute/parallel_roundtrip.i)

Child Objects

(include/problems/LatticeBoltzmannProblem.h)

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

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

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

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

[TensorBuffers]
  [eta]
  []
[]

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

[AuxVariables]
  [eta]
  []
[]

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

[Problem]
  type = TensorProblem
[]

[Executioner]
  type = Transient
  num_steps = 2
[]

[Outputs]
  exodus = true
[]

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

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

[TensorComputes]
    [Initialize]
        [rectangle_sharp]
            type = SmoothRectangleCompute
            buffer = rectangle_sharp
            x1 = 5
            x2 = 15
            y1 = 5
            y2 = 15
            inside = -1
            outside = 3
        []
        [rectangle_cos]
            type = SmoothRectangleCompute
            buffer = rectangle_cos
            x1 = 5
            x2 = 15
            y1 = 5
            y2 = 15
            inside = -1
            outside = 3
            profile = COS
            int_width = 1
        []
        [rectangle_tanh]
            type = SmoothRectangleCompute
            buffer = rectangle_tanh
            x1 = 5
            x2 = 15
            y1 = 5
            y2 = 15
            inside = -1
            outside = 3
            profile = TANH
            int_width = 1
        []
    []
[]

[Problem]
    type = TensorProblem
[]

[TensorOutputs]
    [xdmf]
        type = XDMFTensorOutput
        buffer = 'rectangle_sharp rectangle_cos rectangle_tanh'
        enable_hdf5 = true
    []
[]

[Executioner]
    type = Transient
    num_steps = 0
[]

[Outputs]
    perf_graph = true
[]

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

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

[Mesh]
  type = UniformTensorMesh
  dim = 2
  nx = 50
  ny = 50
[]

[TensorBuffers]
  [eta]
  []
  [f]
  []
[]

[Problem]
  type = TensorProblem
[]

[Executioner]
  type = Transient
  num_steps = 4
[]

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

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

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

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

[Problem]
  type = TensorProblem
[]

[TensorComputes]
  [Initialize]
    [A]
      type = ConstantTensor
      buffer = A
      real = 2

    []
    [B]
      type = ConstantTensor
      buffer = B
      real = 3
    []
    [C]
      type = NEML2TensorCompute
      neml2_input_file = neml2_input.i
      neml2_model = multiply
      marlin_inputs = 'A B'
      neml2_inputs = 'forces/A forces/B'
      neml2_outputs = 'state/C'
      marlin_outputs = 'C'
    []
  []
[]

[Postprocessors]
  [C]
    type = TensorAveragePostprocessor
    buffer = C
  []
[]

[Executioner]
  type = Transient
  num_steps = 1
[]

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

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

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

(include/problems/LatticeBoltzmannProblem.h)

/**********************************************************************/
/*                     DO NOT MODIFY THIS HEADER                      */
/*            Marlin, a Fourier spectral solver for MOOSE             */
/*                                                                    */
/*            Copyright 2024 Battelle Energy Alliance, LLC            */
/*                        ALL RIGHTS RESERVED                         */
/**********************************************************************/

#pragma once

#include "TensorProblem.h"

class LatticeBoltzmannStencilBase;

/**
 * Problem object for solving lattice Boltzmann problems
 */
class LatticeBoltzmannProblem : public TensorProblem
{
public:
  static InputParameters validParams();

  LatticeBoltzmannProblem(const InputParameters & parameters);

  void addTensorBoundaryCondition(const std::string & compute_name,
                                  const std::string & name,
                                  InputParameters & parameters);

  // setup stuff
  void init() override;

  // main loop
  void execute(const ExecFlagType & exec_type) override;

  void addStencil(const std::string & stencil_name,
                  const std::string & name,
                  InputParameters & parameters);

  const LatticeBoltzmannStencilBase & getStencil() const { return *_stencil; }

  const int & getTotalSteps() const { return _t_total; }

  const unsigned int & getGhostRadius() const { return _ghost_radius; }

  const bool & isBinaryMedia() { return _is_binary_media; }

  const torch::Tensor & getBinaryMedia() { return _binary_media; }

  const std::vector<int64_t> & getExtendedShape() { return _shape_extended; }
  const std::vector<int64_t> & getExtendedShapeQ() { return _shape_extended_to_q; }

  /// sets convergence residual
  void setSolverResidual(const Real & residual) { _convergence_residual = residual; };

  /// sets tensor to a value (normally zeros) at solid nodes
  void maskedFillSolids(torch::Tensor & t, const Real & value);

protected:
  /// LBM mesh/media
  torch::Tensor _binary_media;
  torch::Tensor _binary_media_owned;
  const bool _is_binary_media;

  ///
  std::vector<int64_t> _shape_extended;
  std::vector<int64_t> _shape_extended_to_q;

  /// LBM stencils object
  std::shared_ptr<LatticeBoltzmannStencilBase> _stencil;

  /// bc objects
  TensorComputeList _bcs;

  /// radius of ghost layers
  unsigned int _ghost_radius = 0;

  const Real _A_1 = 0.6;
  const Real _A_2 = 0.9;

  /// used to restrict construction of lbm stencils to only one
  unsigned int _stencil_counter = 0;

  /// convergence residual
  Real _convergence_residual = 1;

  /// total number of time steps taken
  int _t_total = 0;

  /// lbm substeps
  const unsigned int _lbm_substeps;

  /// log interval for substeps
  const unsigned int _log_interval;

  /// lbm convergence tolerance
  const Real _tolerance;

public:
  /// LBM constants
  const Real _cs = 1.0 / sqrt(3.0);
  const Real _cs2 = _cs * _cs;
  const Real _cs4 = _cs2 * _cs2;
};

Overview
Example Input File Syntax
Input Parameters
Input Files
Child Objects