MooseFunctionTensor

buffer

C++ Type:std::string

Controllable:No

Description:The buffer this compute is writing to

function

C++ Type:FunctionName

Unit:(no unit assumed)

Controllable:No

Description:Function to map.

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

# Verification against an independent analytic solution: the steady 1D elliptic
# Robin-BC benchmark from OpenGeoSys's documentation,
# https://www.opengeosys.org/6.5.7/docs/benchmarks/elliptic/elliptic-robin/
# ("second benchmark": Dirichlet at x=0, Robin at x=1). The governing equation is
# k*Laplacian(h) = 0 on [0,1] with k=1, boundary conditions
#   h(0) = g_D,                 dh/dn|_{x=1} = alpha*(h0 - h(1)),
# and exact linear solution
#   h(x) = [alpha*(h0 - g_D)/(1 + alpha)] * x + g_D.
# TensorRobinBC's convention is a*u + b*du/dn = c; the OGS Robin condition
# du/dn = alpha*(h0 - u) rearranges to alpha*u + du/dn = alpha*h0, i.e.
# a = alpha, b = 1, c = alpha*h0.
#
# Uses alpha = +2 rather than OGS's alpha = -2: the ghost coefficient on u0,
# -(a/2 - b/dx)/(a/2 + b/dx), has magnitude > 1 for a < 0, so explicit relaxation to
# the steady state diverges (OGS solves the steady system directly). With alpha = 2,
# h0 = 1.5 and g_D = 1 the exact solution is h(x) = x/3 + 1.
#
# Y is periodic, and the IC matches neither the steady profile nor the BCs.
[Domain]
  dim = 2
  nx = 20
  ny = 8
  parallel_mode = REAL_SPACE
  periodic_directions = 'Y'
  xmin = 0
  xmax = 1
  ymin = -1
  ymax = 1
[]

[Functions]
  [u_ic_func]
    type = ParsedFunction
    expression = '0.5'
  []
  [u_exact_func]
    type = ParsedFunction
    expression = '1 + x / 3'
  []
[]

[TensorComputes]
  [Initialize]
    [u_ic]
      type = MooseFunctionTensor
      buffer = u
      function = u_ic_func
    []
  []
  [BCs]
    [left]
      type = TensorDirichletBC
      variable = u
      boundary = left
      value = 1
    []
    [right]
      # alpha*u + du/dn = alpha*h0, with alpha = 2, h0 = 1.5.
      type = TensorRobinBC
      variable = u
      boundary = right
      a = 2
      b = 1
      c = 3
    []
  []
  [Solve]
    [dt]
      type = FiniteDifferenceLaplacian
      input = u
      buffer = dudt
    []
    [u_exact]
      type = MooseFunctionTensor
      buffer = u_exact
      function = u_exact_func
    []
    [diff]
      type = ParsedCompute
      buffer = diff
      expression = 'abs(u - u_exact)'
      inputs = 'u u_exact'
    []
  []
[]

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

[Postprocessors]
  [max_error]
    type = TensorExtremeValuePostprocessor
    buffer = diff
    value_type = MAX
    execute_on = 'TIMESTEP_END'
  []
[]

[Executioner]
  type = Transient
  # same grid and CFL margin as dirichlet_linear.i
  dt = 1
  num_steps = 4
[]

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

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

[Domain]
  dim = 2
  nx = 20
  ny = 8
  parallel_mode = REAL_SPACE
  periodic_directions = 'Y'
  xmin = 0
  xmax = 1
  ymin = -1
  ymax = 1
[]

[Functions]
  [u_func]
    type = ParsedFunction
    expression = '1 + 2 * x'
  []
[]

[TensorComputes]
  [Initialize]
    [u_ic]
      # u = 1 + 2*x satisfies the Dirichlet values below exactly, so the discrete
      # Laplacian is zero to roundoff everywhere, including the boundary-adjacent cells.
      type = MooseFunctionTensor
      buffer = u
      function = u_func
    []
  []
  [BCs]
    [left]
      type = TensorDirichletBC
      variable = u
      boundary = left
      value = 1
    []
    [right]
      type = TensorDirichletBC
      variable = u
      boundary = right
      value = 3
    []
  []
  [Solve]
    [lap_compute]
      type = FiniteDifferenceLaplacian
      input = u
      buffer = lap
    []
  []
[]

[Postprocessors]
  [max_lap]
    type = TensorExtremeValuePostprocessor
    buffer = lap
    value_type = MAX
  []
  [min_lap]
    type = TensorExtremeValuePostprocessor
    buffer = lap
    value_type = MIN
  []
[]

[Executioner]
  type = Transient
  num_steps = 1
[]

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

(test/tests/implicit_fd/bc_mms.i)

# Implicit solve of u_t - D Laplacian(u) = f on the unit square with boundary faces, against the
# manufactured solution
#
#   u = g + tf t q,  g = cos(3 pi x / 2) cos(pi y) + x y,  q = x^2 - y^2,
#
# with tf = 0 for the steady problem -D Laplacian(u) = f and tf = 1 for backward Euler.
# q is harmonic, so f = tf q - D Laplacian(g) is independent of time and the backward Euler
# step makes no time error: the error is the spatial one, provided the boundary values are
# taken at the time of the unknown, t + dt (at t it would be of order dt).
#
# Each face has a different condition and value source:
#   left   (x = 0)  Dirichlet   u = cos(pi y) - tf t y^2                     function
#   right  (x = 1)  Robin       ra u + rb du/dx = c                         function
#                                                      or buffer (right_buffer, steady)
#   bottom (y = 0)  Neumann     du/dn = -du/dy = -x                          buffer
#   top    (y = 1)  Neumann     du/dn = du/dy = x - 2 tf t                   function
#                                                      or buffer (top_buffer, values at t)

D = 0.5
n = 32
stencil = 3
tf = 0
ops = 'lap'
num_steps = 1
dt = 0.05
# Robin coefficients on the right face (the 5-point stencil supports ra = 1, rb = 0 only)
ra = 1
rb = 1

[Domain]
  dim = 2
  nx = ${n}
  ny = ${n}
  xmax = 1
  ymax = 1
  parallel_mode = REAL_SPACE
[]

[Functions]
  [u_exact]
    type = ParsedFunction
    expression = 'cos(1.5*pi*x) * cos(pi*y) + x*y + ${tf} * t * (x^2 - y^2)'
  []
  [f_func]
    type = ParsedFunction
    # Laplacian(g) = -13/4 pi^2 cos(3 pi x / 2) cos(pi y)
    expression = '${tf} * (x^2 - y^2) + ${D} * 3.25 * pi^2 * cos(1.5*pi*x) * cos(pi*y)'
  []
  [right_c]
    type = ParsedFunction
    # on x = 1: u = y + tf t (1 - y^2), du/dx = 3/2 pi cos(pi y) + y + 2 tf t
    expression = '${ra} * (y + ${tf} * t * (1 - y^2)) + ${rb} * (1.5*pi*cos(pi*y) + y + 2 * ${tf} * t)'
  []
  [bottom_flux]
    type = ParsedFunction
    expression = '-x'
  []
  [top_flux]
    type = ParsedFunction
    expression = 'x - 2 * ${tf} * t'
  []
[]

[TensorBuffers]
  [u]
  []
  [dudt_u]
  []
  [lap_u]
  []
  [f]
  []
  [bottom_flux]
  []
  [right_c]
  []
  [top_flux]
  []
  [exact]
  []
  [diff]
  []
[]

[TensorComputes]
  [Initialize]
    [u_ic]
      type = MooseFunctionTensor
      buffer = u
      function = u_exact
    []
    [f_ic]
      type = MooseFunctionTensor
      buffer = f
      function = f_func
    []
    # a face value: depends on the tangential coordinate x only
    [bottom_flux]
      type = MooseFunctionTensor
      buffer = bottom_flux
      function = bottom_flux
    []
    [right_c]
      type = MooseFunctionTensor
      buffer = right_c
      function = right_c
    []
  []
  [BCs]
    active = 'left right bottom top'
    [left]
      type = TensorDirichletBC
      variable = u
      boundary = left
      value_function = u_exact
    []
    [right]
      type = TensorRobinBC
      variable = u
      boundary = right
      a = ${ra}
      b = ${rb}
      c_function = right_c
    []
    # the Robin c from a buffer instead (active = 'left right_buffer bottom top'), steady only
    [right_buffer]
      type = TensorRobinBC
      variable = u
      boundary = right
      a = ${ra}
      b = ${rb}
      c_buffer = right_c
    []
    [bottom]
      type = TensorNeumannBC
      variable = u
      boundary = bottom
      flux_buffer = bottom_flux
    []
    [top]
      type = TensorNeumannBC
      variable = u
      boundary = top
      flux_function = top_flux
    []
    # the top flux from a buffer instead (active = 'left right bottom top_buffer'); a buffer
    # value is the one its compute wrote at the start of the substep, at t
    [top_buffer]
      type = TensorNeumannBC
      variable = u
      boundary = top
      flux_buffer = top_flux
    []
  []
  [Solve]
    [top_flux]
      type = MooseFunctionTensor
      buffer = top_flux
      function = top_flux
    []
    [dudt]
      type = TimeDerivativeOperator
      buffer = dudt_u
      input = u
    []
    [lap]
      type = FiniteDifferenceLaplacian
      buffer = lap_u
      input = u
      factor = '${fparse -D}'
      stencil_width = ${stencil}
    []
  []
  [Postprocess]
    [exact_compute]
      type = MooseFunctionTensor
      buffer = exact
      function = u_exact
    []
    [diff_compute]
      type = ParsedCompute
      buffer = diff
      expression = 'abs(u - exact)'
      inputs = 'u exact'
    []
  []
[]

[TensorSolver]
  type = ImplicitFDSolver
  buffer = u
  operators = '${ops}'
  source = f
  linear_solver = TORCH_PCG
  l_tol = 1e-10
[]

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

[Problem]
  type = TensorProblem
[]

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

[Outputs]
  csv = true
[]

(test/tests/real_space/bc_varying_laplacian.i)

# Boundary condition values that vary along the face, from functions and buffers.
#
# u = (x + 1) * y is harmonic and linear in the normal direction of every face, so the ghost
# values are exact, and the discrete Laplacian must be zero everywhere, including next to the
# faces. Each face takes its value from a different source.
[Domain]
  dim = 2
  nx = 20
  ny = 10
  parallel_mode = REAL_SPACE
  xmin = 0
  xmax = 1
  ymin = 0
  ymax = 1
[]

[Functions]
  [u_func]
    type = ParsedFunction
    expression = '(x + 1) * y'
  []
  [right_value]
    # u on x = 1
    type = ParsedFunction
    expression = '2 * y'
  []
  [top_flux]
    # du/dy on y = 1
    type = ParsedFunction
    expression = 'x + 1'
  []
  [bottom_flux]
    # -du/dy on y = 0 (outward normal -y)
    type = ParsedFunction
    expression = '-(x + 1)'
  []
[]

[TensorComputes]
  [Initialize]
    [u_ic]
      type = MooseFunctionTensor
      buffer = u
      function = u_func
    []
    # buffer values are face values: these depend on the tangential coordinate only
    [bottom_flux]
      type = MooseFunctionTensor
      buffer = bottom_flux
      function = bottom_flux
    []
  []
  [BCs]
    [left]
      # evaluated at the face position x = 0, where it is y
      type = TensorDirichletBC
      variable = u
      boundary = left
      value_function = u_func
    []
    [right]
      type = TensorDirichletBC
      variable = u
      boundary = right
      value_buffer = right_value
    []
    [bottom]
      type = TensorNeumannBC
      variable = u
      boundary = bottom
      flux_buffer = bottom_flux
    []
    [top]
      type = TensorNeumannBC
      variable = u
      boundary = top
      flux_function = top_flux
    []
  []
  [Solve]
    [lap_compute]
      type = FiniteDifferenceLaplacian
      input = u
      buffer = lap
    []
    # listed after the Laplacian, whose ghost fill reads it: the dependency sort runs it first
    [right_value]
      type = MooseFunctionTensor
      buffer = right_value
      function = right_value
    []
  []
[]

[Postprocessors]
  [max_lap]
    type = TensorExtremeValuePostprocessor
    buffer = lap
    value_type = MAX
  []
  [min_lap]
    type = TensorExtremeValuePostprocessor
    buffer = lap
    value_type = MIN
  []
[]

[Executioner]
  type = Transient
  num_steps = 1
[]

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

(test/tests/implicit_fd/crank_thin_slab.i)

# Crank, "The Mathematics of Diffusion" (2nd ed.), Sec. 2.5, eq. (2.15), with the implicit
# solver: an instantaneous planar source C0 = 1 confined to the slab -h < x < h, h = 1,
# diffusing into an infinite medium, has the exact solution
#   C(x, t) = (C0/2) (erf((h - x) / (2 sqrt(D t))) + erf((h + x) / (2 sqrt(D t)))).
# The same problem as real_space/thin_slab_diffusion.i, which steps it explicitly.
#
# Pseudo-1D (ny = 1, Y periodic), D = 1. The X faces are 10 diffusion lengths beyond the slab
# edge at the comparison time t = 1, so the far Dirichlet u = 0 does not measurably affect the
# solution.
# Backward Euler with substeps of dt / substeps: the error is first order in the substep.

substeps = 100

[Domain]
  dim = 2
  # dx = 0.1 puts the slab edges x = +/-1 on cell faces, so the box IC holds the slab mass
  nx = 220
  ny = 1
  parallel_mode = REAL_SPACE
  periodic_directions = 'Y'
  xmin = -11
  xmax = 11
  ymin = -0.5
  ymax = 0.5
[]

[Functions]
  [slab]
    type = ParsedFunction
    expression = 'if(abs(x) < 1, 1, 0)'
  []
[]

[TensorBuffers]
  [u]
  []
  [dudt_u]
  []
  [lap_u]
  []
  [u_exact]
  []
  [diff]
  []
[]

[TensorComputes]
  [Initialize]
    [u_ic]
      type = MooseFunctionTensor
      buffer = u
      function = slab
    []
  []
  [BCs]
    # far boundaries, approximating the infinite medium
    [left]
      type = TensorDirichletBC
      variable = u
      boundary = left
      value = 0
    []
    [right]
      type = TensorDirichletBC
      variable = u
      boundary = right
      value = 0
    []
  []
  [Solve]
    [dudt]
      type = TimeDerivativeOperator
      buffer = dudt_u
      input = u
    []
    [lap]
      type = FiniteDifferenceLaplacian
      buffer = lap_u
      input = u
      factor = -1
    []
  []
  [Postprocess]
    [u_exact]
      type = ParsedCompute
      buffer = u_exact
      expression = '0.5 * (erf((1 - x) / (2 * sqrt(t))) + erf((1 + x) / (2 * sqrt(t)))) + 0*y'
      extra_symbols = true
    []
    [diff]
      type = ParsedCompute
      buffer = diff
      expression = 'abs(u - u_exact)'
      inputs = 'u u_exact'
    []
  []
[]

[TensorSolver]
  type = ImplicitFDSolver
  buffer = u
  operators = 'dudt lap'
  root_compute = dudt
  substeps = ${substeps}
  l_tol = 1e-10
[]

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

[Problem]
  type = TensorProblem
[]

[Executioner]
  type = Transient
  dt = 1
  num_steps = 1
[]

[Outputs]
  csv = true
[]

(test/tests/kks/KKS_no_flux_bc.i)

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


# Constants for Initial Conditions
r = 30
l = 4.2

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

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

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


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

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

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

  # automatically create a matching mesh
  mesh_mode = DUMMY
[]

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

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

    []
  []
[]

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

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

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

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

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

(test/tests/real_space/neumann_linear.i)

# Nonzero-flux Neumann boundary condition: steady 1D diffusion with u(0) = 1
# (Dirichlet, left) and du/dn = 2 (Neumann, right) has the exact solution
# u = 1 + 2*x. Y is periodic, and the IC is not the steady profile.
[Domain]
  dim = 2
  nx = 20
  ny = 8
  parallel_mode = REAL_SPACE
  periodic_directions = 'Y'
  xmin = 0
  xmax = 1
  ymin = -1
  ymax = 1
[]

[Functions]
  [u_ic_func]
    type = ParsedFunction
    expression = '0.5'
  []
  [u_exact_func]
    type = ParsedFunction
    expression = '1 + 2 * x'
  []
[]

[TensorComputes]
  [Initialize]
    [u_ic]
      type = MooseFunctionTensor
      buffer = u
      function = u_ic_func
    []
  []
  [BCs]
    [left]
      type = TensorDirichletBC
      variable = u
      boundary = left
      value = 1
    []
    [right]
      type = TensorNeumannBC
      variable = u
      boundary = right
      flux = 2
    []
  []
  [Solve]
    [dt]
      type = FiniteDifferenceLaplacian
      input = u
      buffer = dudt
    []
    [u_exact]
      type = MooseFunctionTensor
      buffer = u_exact
      function = u_exact_func
    []
    [diff]
      type = ParsedCompute
      buffer = diff
      expression = 'abs(u - u_exact)'
      inputs = 'u u_exact'
    []
  []
[]

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

[Postprocessors]
  [max_error]
    type = TensorExtremeValuePostprocessor
    buffer = diff
    value_type = MAX
    execute_on = 'TIMESTEP_END'
  []
[]

[Executioner]
  type = Transient
  # same grid and CFL margin as dirichlet_linear.i; the Neumann face makes the slowest
  # mode decay more slowly, so max_error needs 15 steps to plateau
  dt = 1
  num_steps = 15
[]

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

(test/tests/real_space/dirichlet_2d_corners.i)

# Non-periodic on both axes (Dirichlet left/right, zero-flux Neumann bottom/top),
# so the ghost fill reaches the four domain corners. exchangeGhostLayers() fills one
# axis at a time, so a corner ghost depends on the fill order.
#
# u = 1 + 2*x is uniform on every face, so it satisfies all four BCs exactly, and the
# discrete Laplacian must be zero everywhere, including the corner-adjacent cells.
[Domain]
  dim = 2
  nx = 20
  ny = 10
  parallel_mode = REAL_SPACE
  xmin = 0
  xmax = 1
  ymin = 0
  ymax = 1
[]

[Functions]
  [u_func]
    type = ParsedFunction
    expression = '1 + 2 * x'
  []
[]

[TensorComputes]
  [Initialize]
    [u_ic]
      # the steady profile implied by the four BCs below
      type = MooseFunctionTensor
      buffer = u
      function = u_func
    []
  []
  [BCs]
    [left]
      type = TensorDirichletBC
      variable = u
      boundary = left
      value = 1
    []
    [right]
      type = TensorDirichletBC
      variable = u
      boundary = right
      value = 3
    []
    [bottom]
      type = TensorNeumannBC
      variable = u
      boundary = bottom
      flux = 0
    []
    [top]
      type = TensorNeumannBC
      variable = u
      boundary = top
      flux = 0
    []
  []
  [Solve]
    [lap_compute]
      type = FiniteDifferenceLaplacian
      input = u
      buffer = lap
    []
  []
[]

[Postprocessors]
  [max_lap]
    type = TensorExtremeValuePostprocessor
    buffer = lap
    value_type = MAX
  []
  [min_lap]
    type = TensorExtremeValuePostprocessor
    buffer = lap
    value_type = MIN
  []
[]

[Executioner]
  type = Transient
  num_steps = 1
[]

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

(test/tests/real_space/dirichlet_linear.i)

[Domain]
  dim = 2
  nx = 20
  ny = 8
  parallel_mode = REAL_SPACE
  periodic_directions = 'Y'
  xmin = 0
  xmax = 1
  ymin = -1
  ymax = 1
[]

[Functions]
  [u_ic_func]
    type = ParsedFunction
    expression = '1 + 2 * x + 0.001 * y'
  []
[]

[TensorComputes]
  [Initialize]
    [u_ic]
      # Not the steady solution and not matching the BCs pointwise (it depends on y);
      # diffusion must drive u to the exact steady profile u = 1 + 2*x.
      type = MooseFunctionTensor
      buffer = u
      function = u_ic_func
    []
  []
  [BCs]
    [left]
      type = TensorDirichletBC
      variable = u
      boundary = left
      value = 1
    []
    [right]
      type = TensorDirichletBC
      variable = u
      boundary = right
      value = 3
    []
  []
  [Solve]
    [dt]
      type = FiniteDifferenceLaplacian
      input = u
      buffer = dudt
    []
  []
[]

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

[Postprocessors]
  [max_u]
    type = TensorExtremeValuePostprocessor
    buffer = u
    value_type = MAX
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [min_u]
    type = TensorExtremeValuePostprocessor
    buffer = u
    value_type = MIN
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]

[Executioner]
  type = Transient
  # dx = 0.05, dy = 0.25: explicit CFL limit 1/(2*(1/dx^2 + 1/dy^2)) ~= 1.2e-3,
  # above the substep dt/substeps = 1e-3.
  dt = 1
  num_steps = 4
[]

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