FiniteDifferenceLaplacian

stencil_width

Default:3

C++ Type:MooseEnum

Options:3, 5

Controllable:No

Description:Stencil width for the Laplacian (3 => 3-point, 5 => 5-point)

(test/tests/real_space/laplacian.i)

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

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

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

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

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

[Executioner]
  type = Transient
  num_steps = 1
[]

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

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

(test/tests/real_space/semi_infinite_diffusion.i)

# Verification against Crank, "The Mathematics of Diffusion" (2nd ed.), Sec. 3.3(i),
# eq. (3.13): a semi-infinite medium (x > 0) with uniform initial concentration C0 = 0
# and constant surface concentration C1 maintained at x = 0 has the exact solution
#   C(x, t) = C1 * erfc(x / (2*sqrt(D*t)))
# (special case of (3.13) with C0 = 0, using erfc = 1 - erf).
#
# Pseudo-1D (ny = 1, Y periodic), D = 1. The domain extends 10 diffusion lengths at
# the comparison time t = 1 (erfc(10) ~ 2e-45), so the far Dirichlet u = 0 does not
# measurably affect the solution near x = 0.
[Domain]
  dim = 2
  nx = 100
  ny = 1
  parallel_mode = REAL_SPACE
  periodic_directions = 'Y'
  xmin = 0
  xmax = 10
  ymin = -0.5
  ymax = 0.5
[]

[TensorComputes]
  # no [Initialize]: buffers start at zero, the C0 = 0 initial condition
  [BCs]
    [left]
      type = TensorDirichletBC
      variable = u
      boundary = left
      value = 1
    []
    [right]
      # Far boundary, approximating the semi-infinite medium's x -> infinity
      # condition; diffusion has not reached here by the comparison time.
      type = TensorDirichletBC
      variable = u
      boundary = right
      value = 0
    []
  []
  [Solve]
    [dt]
      type = FiniteDifferenceLaplacian
      input = u
      buffer = dudt
      factor = 1
    []
    [u_exact]
      type = ParsedCompute
      buffer = u_exact
      expression = 'erfc(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 = RealSpaceForwardEuler
  buffer = u
  time_derivative = dudt
  substeps = 250
[]

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

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

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

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

# Zero-flux Neumann boundary condition: pure diffusion with zero flux on the X faces
# and periodic Y, from a non-uniform initial condition. The total integral of u must
# be conserved to near machine precision while u relaxes to its mean.
[Domain]
  dim = 2
  nx = 20
  ny = 8
  parallel_mode = REAL_SPACE
  periodic_directions = 'Y'
  xmin = 0
  xmax = 1
  ymin = -1
  ymax = 1
[]

[TensorComputes]
  [Initialize]
    [u_ic]
      type = ParsedCompute
      buffer = u
      # non-uniform in x and y, so diffusion moves mass internally
      expression = '1 + sin(2 * pi * x) + 0.3 * cos(pi * y)'
      extra_symbols = true
    []
  []
  [BCs]
    [left]
      type = TensorNeumannBC
      variable = u
      boundary = left
      flux = 0
    []
    [right]
      type = TensorNeumannBC
      variable = u
      boundary = right
      flux = 0
    []
  []
  [Solve]
    [dt]
      type = FiniteDifferenceLaplacian
      input = u
      buffer = dudt
    []
  []
[]

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

[Postprocessors]
  [mass]
    type = TensorIntegralPostprocessor
    buffer = u
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]

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

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

(test/tests/implicit_fd/crank_semi_infinite.i)

# Crank, "The Mathematics of Diffusion" (2nd ed.), Sec. 3.3(i), eq. (3.13), with the implicit
# solver: a semi-infinite medium (x > 0) with zero initial concentration and constant surface
# concentration C1 = 1 at x = 0 has the exact solution
#   C(x, t) = erfc(x / (2 sqrt(D t))).
# The same problem as real_space/semi_infinite_diffusion.i, which steps it explicitly.
#
# Pseudo-1D (ny = 1, Y periodic), D = 1. The domain extends 10 diffusion lengths 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
  nx = 100
  ny = 1
  parallel_mode = REAL_SPACE
  periodic_directions = 'Y'
  xmin = 0
  xmax = 10
  ymin = -0.5
  ymax = 0.5
[]

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

[TensorComputes]
  # no [Initialize]: buffers start at zero, the C0 = 0 initial condition
  [BCs]
    [left]
      type = TensorDirichletBC
      variable = u
      boundary = left
      value = 1
    []
    [right]
      # far boundary, approximating the semi-infinite medium's x -> infinity condition
      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 = 'erfc(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/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
  []
[]

(test/tests/real_space/diffusion.i)

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

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

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

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

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

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

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

(test/tests/real_space/thin_slab_diffusion.i)

# Verification against Crank, "The Mathematics of Diffusion" (2nd ed.), Sec. 2.5,
# eq. (2.15): an instantaneous planar source of concentration C0 confined initially
# to a thin slab -h < x < h, 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))] }
#
# Pseudo-1D (ny = 1, Y periodic), D = 1, h = 1. The X faces are 10 diffusion lengths
# beyond the slab edge at the comparison time t = 1, so the far-field Dirichlet u = 0
# does not measurably affect the solution near the slab.
[Domain]
  dim = 2
  # dx = 0.1 puts the slab edges x = +/-1 on cell faces, so the box IC holds exactly
  # the slab mass 2*h*C0 (nx = 200 leaves a ~1% deficit that does not refine away)
  nx = 220
  ny = 1
  parallel_mode = REAL_SPACE
  periodic_directions = 'Y'
  xmin = -11
  xmax = 11
  ymin = -0.5
  ymax = 0.5
[]

[TensorComputes]
  [Initialize]
    [u_ic]
      type = ParsedCompute
      buffer = u
      # C0 = 1 inside the slab |x| < h = 1, zero outside
      expression = 'if(abs(x) < 1, 1, 0) + 0*y'
      extra_symbols = true
    []
  []
  [BCs]
    [left]
      type = TensorDirichletBC
      variable = u
      boundary = left
      value = 0
    []
    [right]
      type = TensorDirichletBC
      variable = u
      boundary = right
      value = 0
    []
  []
  [Solve]
    [dt]
      type = FiniteDifferenceLaplacian
      input = u
      buffer = dudt
      factor = 1
    []
    [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 = RealSpaceForwardEuler
  buffer = u
  time_derivative = dudt
  # dx = 0.1; explicit-diffusion CFL limit dt_cfl = dx^2/(2*D) = 5e-3,
  # requiring >= dt/dt_cfl = 200 substeps for dt = 1; 250 gives a safety margin.
  substeps = 250
[]

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

[Executioner]
  type = Transient
  dt = 1
  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/implicit_fd/cutback.i)

# Backward Euler diffusion of a Gaussian with a linear iteration limit that the first dt
# exceeds. The torch PCG solve needs about 44 iterations at dt = 0.008, 34 at 0.004 and 23
# at 0.002 (32 x 32, l_tol = 1e-8), so with l_max_its = 30 the step is cut twice and then
# converges.

D = 0.5
n = 32

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

[TensorBuffers]
  [u]
  []
  [dudt_u]
  []
  [lap_u]
  []
[]

[TensorComputes]
  [Initialize]
    [u_ic]
      type = ParsedCompute
      buffer = u
      expression = 'exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.005)'
      extra_symbols = true
    []
  []
  [Solve]
    [dudt]
      type = TimeDerivativeOperator
      buffer = dudt_u
      input = u
    []
    [lap]
      type = FiniteDifferenceLaplacian
      buffer = lap_u
      input = u
      factor = '${fparse -D}'
    []
  []
[]

[TensorSolver]
  type = ImplicitFDSolver
  buffer = u
  operators = 'dudt lap'
  linear_solver = TORCH_PCG
  l_max_its = 30
[]

[Postprocessors]
  [max_u]
    type = TensorExtremeValuePostprocessor
    buffer = u
    value_type = MAX
  []
[]

[Problem]
  type = TensorProblem
[]

[Executioner]
  type = Transient
  num_steps = 2
  [TimeStepper]
    type = TensorSolveIterationAdaptiveDT
    dt = 0.008
    # no growth or shrinking from the iteration count, only the cut on failure
    min_iterations = 0
    max_iterations = 1000
  []
[]

[Outputs]
  csv = true
[]

(test/tests/implicit_fd/poisson_mms.i)

# Steady periodic Poisson problem -D Laplacian(u) = f with the manufactured
# solution u = sin(2 pi x) cos(4 pi y) + cos(6 pi x) / 2 (zero mean). The
# operator is singular (constants), and the solver returns the zero-mean
# solution, so the error against u converges at the stencil order (2 for the
# 3-point, 4 for the 5-point stencil).

D = 0.5
n = 32

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

[TensorBuffers]
  [u]
  []
  [lap_u]
  []
  [f]
  []
  [exact]
  []
  [diff]
  []
[]

[TensorComputes]
  [Initialize]
    [u_ic]
      type = ParsedCompute
      buffer = u
      # initial guess (any field; the solve returns the zero-mean solution)
      expression = 'sin(2*pi*x) * sin(2*pi*y)'
      extra_symbols = true
    []
    [f_ic]
      type = ParsedCompute
      buffer = f
      # -D Laplacian(u): the two modes have |k|^2 = 20 pi^2 and 36 pi^2
      expression = 'D * (20 * pi^2 * sin(2*pi*x) * cos(4*pi*y) + 18 * pi^2 * cos(6*pi*x))'
      constant_names = 'D'
      constant_expressions = '${D}'
      extra_symbols = true
    []
  []
  [Solve]
    [lap]
      type = FiniteDifferenceLaplacian
      buffer = lap_u
      input = u
      factor = '${fparse -D}'
    []
  []
  [Postprocess]
    [exact_compute]
      type = ParsedCompute
      buffer = exact
      expression = 'sin(2*pi*x) * cos(4*pi*y) + cos(6*pi*x) / 2'
      extra_symbols = true
    []
    [diff_compute]
      type = ParsedCompute
      buffer = diff
      expression = 'abs(u - exact)'
      inputs = 'u exact'
    []
  []
[]

[TensorSolver]
  type = ImplicitFDSolver
  buffer = u
  operators = 'lap'
  source = f
[]

[Postprocessors]
  [max_error]
    type = TensorExtremeValuePostprocessor
    buffer = diff
    value_type = MAX
  []
  [mean_u]
    type = TensorAveragePostprocessor
    buffer = u
  []
[]

[Problem]
  type = TensorProblem
[]

[Executioner]
  type = Transient
  num_steps = 1
[]

[Outputs]
  csv = true
[]

(test/tests/real_space/laplacian.i)

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

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

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

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

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

[Executioner]
  type = Transient
  num_steps = 1
[]

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

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

# Backward Euler diffusion u_t = D Laplacian(u) of a single periodic Fourier
# mode. On the periodic cell-centred grid the mode is an eigenvector of the
# discrete Laplacian with eigenvalue -lambda_h, so every implicit step must
# multiply it by exactly g = 1 / (1 + dt D lambda_h), up to the linear solver
# tolerance.

D = 0.5
# 5 steps of dt, then 5 steps of dt2 (dt2 != dt exercises the operator rebuild)
dt = 0.01
dt2 = 0.01
n = 16
stencil = 3
dim = 3

# Symbol of -Laplacian per axis for the wavenumber 2 pi, h = 1/n:
#   3-point: (2 - 2 cos(2 pi h)) / h^2
#   5-point: (5/2 - 8/3 cos(2 pi h) + 1/6 cos(4 pi h)) / h^2
lambda3 = '${fparse n^2 * (2 - 2 * cos(2 * pi / n))}'
lambda5 = '${fparse n^2 * (5 / 2 - 8 / 3 * cos(2 * pi / n) + 1 / 6 * cos(4 * pi / n))}'
lambda = '${fparse dim * if(stencil = 3, lambda3, lambda5)}'

[Domain]
  dim = 3
  nx = ${n}
  ny = ${n}
  nz = ${n}
  zmax = 1
  xmax = 1
  ymax = 1
  parallel_mode = REAL_SPACE
  periodic_directions = 'X Y Z'
[]

[TensorBuffers]
  [u]
  []
  [dudt_u]
  []
  [lap_u]
  []
  [gold]
  []
  [diff]
  []
[]

[TensorComputes]
  [Initialize]
    [u_ic]
      type = ParsedCompute
      buffer = u
      expression = 'sin(2*pi*x)*sin(2*pi*y)*sin(2*pi*z)'
      extra_symbols = true
    []
  []
  [Solve]
    [dudt]
      type = TimeDerivativeOperator
      buffer = dudt_u
      input = u
    []
    [lap]
      type = FiniteDifferenceLaplacian
      buffer = lap_u
      input = u
      factor = '${fparse -D}'
      stencil_width = ${stencil}
    []
  []
  [Postprocess]
    [gold_compute]
      type = ParsedCompute
      buffer = gold
      # g(dt)^(t/dt) up to t1 = 5 dt, then g(dt)^5 g(dt2)^((t - t1)/dt2)
      expression = 'sin(2*pi*x)*sin(2*pi*y)*sin(2*pi*z) * if(t < t1 + dt / 2, g1^(t / dt), g1^5 * g2^((t - t1) / dt2))'
      constant_names = 'g1 g2 dt dt2 t1'
      constant_expressions = '${fparse 1 / (1 + dt * D * lambda)} ${fparse 1 / (1 + dt2 * D * lambda)} ${dt} ${dt2} ${fparse 5 * dt}'
      extra_symbols = true
    []
    [diff_compute]
      type = ParsedCompute
      buffer = diff
      expression = 'abs(u - gold)'
      inputs = 'u gold'
    []
  []
[]

[TensorSolver]
  type = ImplicitFDSolver
  buffer = u
  operators = 'dudt lap'
  root_compute = dudt
[]

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

[Problem]
  type = TensorProblem
[]

[Executioner]
  type = Transient
  num_steps = 10
  [TimeStepper]
    type = TimeSequenceStepper
    time_sequence = '${fparse 1 * dt} ${fparse 2 * dt} ${fparse 3 * dt} ${fparse 4 * dt} ${fparse 5 * dt}
                     ${fparse 5 * dt + dt2} ${fparse 5 * dt + 2 * dt2} ${fparse 5 * dt + 3 * dt2}
                     ${fparse 5 * dt + 4 * dt2} ${fparse 5 * dt + 5 * dt2}'
  []
[]

[Outputs]
  csv = true
[]

(test/tests/implicit_fd/be_mode.i)

# Backward Euler diffusion u_t = D Laplacian(u) of a single periodic Fourier
# mode. On the periodic cell-centred grid the mode is an eigenvector of the
# discrete Laplacian with eigenvalue -lambda_h, so every implicit step must
# multiply it by exactly g = 1 / (1 + dt D lambda_h), up to the linear solver
# tolerance.

D = 0.5
# 5 steps of dt, then 5 steps of dt2 (dt2 != dt exercises the operator rebuild)
dt = 0.01
dt2 = 0.01
n = 32
stencil = 3
dim = 2

# Symbol of -Laplacian per axis for the wavenumber 2 pi, h = 1/n:
#   3-point: (2 - 2 cos(2 pi h)) / h^2
#   5-point: (5/2 - 8/3 cos(2 pi h) + 1/6 cos(4 pi h)) / h^2
lambda3 = '${fparse n^2 * (2 - 2 * cos(2 * pi / n))}'
lambda5 = '${fparse n^2 * (5 / 2 - 8 / 3 * cos(2 * pi / n) + 1 / 6 * cos(4 * pi / n))}'
lambda = '${fparse dim * if(stencil = 3, lambda3, lambda5)}'

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

[TensorBuffers]
  [u]
  []
  [dudt_u]
  []
  [lap_u]
  []
  [gold]
  []
  [diff]
  []
[]

[TensorComputes]
  [Initialize]
    [u_ic]
      type = ParsedCompute
      buffer = u
      expression = 'sin(2*pi*x)*sin(2*pi*y)'
      extra_symbols = true
    []
  []
  [Solve]
    [dudt]
      type = TimeDerivativeOperator
      buffer = dudt_u
      input = u
    []
    [lap]
      type = FiniteDifferenceLaplacian
      buffer = lap_u
      input = u
      factor = '${fparse -D}'
      stencil_width = ${stencil}
    []
  []
  [Postprocess]
    [gold_compute]
      type = ParsedCompute
      buffer = gold
      # g(dt)^(t/dt) up to t1 = 5 dt, then g(dt)^5 g(dt2)^((t - t1)/dt2)
      expression = 'sin(2*pi*x)*sin(2*pi*y) * if(t < t1 + dt / 2, g1^(t / dt), g1^5 * g2^((t - t1) / dt2))'
      constant_names = 'g1 g2 dt dt2 t1'
      constant_expressions = '${fparse 1 / (1 + dt * D * lambda)} ${fparse 1 / (1 + dt2 * D * lambda)} ${dt} ${dt2} ${fparse 5 * dt}'
      extra_symbols = true
    []
    [diff_compute]
      type = ParsedCompute
      buffer = diff
      expression = 'abs(u - gold)'
      inputs = 'u gold'
    []
  []
[]

[TensorSolver]
  type = ImplicitFDSolver
  buffer = u
  operators = 'dudt lap'
  root_compute = dudt
[]

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

[Problem]
  type = TensorProblem
[]

[Executioner]
  type = Transient
  num_steps = 10
  [TimeStepper]
    type = TimeSequenceStepper
    time_sequence = '${fparse 1 * dt} ${fparse 2 * dt} ${fparse 3 * dt} ${fparse 4 * dt} ${fparse 5 * dt}
                     ${fparse 5 * dt + dt2} ${fparse 5 * dt + 2 * dt2} ${fparse 5 * dt + 3 * dt2}
                     ${fparse 5 * dt + 4 * dt2} ${fparse 5 * dt + 5 * dt2}'
  []
[]

[Outputs]
  csv = true
[]

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

(benchmarks/01_spinodal_decomposition/1a_implicit_fd.i)

#
# PFHub benchmark 1a (spinodal decomposition, square periodic domain) with the
# implicit finite-difference solver. Same problem as 1a.i:
#
#   dc/dt = div(M grad(f'(c) - kappa Laplacian(c))),  f = rho_s (c - c_alpha)^2 (c_beta - c)^2
#
# Backward Euler, with the biharmonic implicit and the double-well flux
# divergence explicit:
#
#   (c - c^n)/dt + M kappa Laplacian(Laplacian(c)) = M Laplacian(f'(c^n))
#
# M = 5, kappa = 2, rho_s = 5, c_alpha = 0.3, c_beta = 0.7.
#

# substeps per time step, as in 1a_solver.i (each substep is one linear solve)
substeps = 1000

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

[TensorBuffers]
  [c]
  []
  [dcdt]
  []
  [bih_c]
  []
  [mu]
  []
  [rhs]
  []
  # postprocessing
  [lap_c]
  []
  [F]
  []
[]

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

  [Solve]
    # implicit operators
    [dcdt]
      type = TimeDerivativeOperator
      buffer = dcdt
      input = c
    []
    [bih]
      type = FiniteDifferenceBiharmonic
      buffer = bih_c
      input = c
      factor = 10 # M * kappa
    []

    # explicit source M Laplacian(f'(c^n))
    [explicit]
      [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
      []
      [rhs]
        type = FiniteDifferenceLaplacian
        buffer = rhs
        input = mu
        factor = 5 # M
      []
    []
  []

  [Postprocess]
    # On the periodic grid sum |grad c|^2 = -sum c Laplacian(c) for the 3-point
    # Laplacian, so the gradient energy (kappa/2)|grad c|^2 integrates to the same
    # total as -(kappa/2) c Laplacian(c).
    [lap_c]
      type = FiniteDifferenceLaplacian
      buffer = lap_c
      input = c
    []
    [F]
      type = ParsedCompute
      buffer = F
      expression = 'rho_s * (c-c_alpha)^2 * (c_beta-c)^2 - kappa / 2 * c * lap_c'
      constant_names = 'rho_s c_alpha c_beta kappa'
      constant_expressions = '5     0.3     0.7    2'
      inputs = 'c lap_c'
    []
  []
[]

[TensorSolver]
  type = ImplicitFDSolver
  buffer = c
  operators = 'dcdt bih'
  source = rhs
  root_compute = explicit
  substeps = ${substeps}
[]

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

[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
  []
  # cumulative wall time of the run (s)
  [wall_time]
    type = PerfGraphData
    section_name = Root
    data_type = TOTAL
  []
[]

[Problem]
  type = TensorProblem
[]

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

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

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

(test/tests/real_space/laplacian3d.i)

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

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

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

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

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

[Executioner]
  type = Transient
  num_steps = 1
[]

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

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