DeAliasingTensor

Create a de-aliasing filter.

Creates a spectral de-aliasing filter in reciprocal space. Two methods are provided:

SHARP: a two-thirds truncation ("2/3 rule") that zeros high-frequency modes. The filter value $\sigma(\mathbf{k})$ is

$var element = document.getElementById("moose-equation-9740eebf-515e-4327-89cd-1d51613a08c9");katex.render("\\sigma(\\mathbf{k}) = \\begin{cases} 0, & \\text{if } |k_x| \\gt \\tfrac{2}{3} k_{x,\\max} \\text{ or } |k_y| \\gt \\tfrac{2}{3} k_{y,\\max} \\text{ or } |k_z| \\gt \\tfrac{2}{3} k_{z,\\max}, \\\\ 1, & \\text{otherwise.} \\end{cases}", element, {displayMode:true,throwOnError:false});$

HOULI: the Hou-Li exponential filter, which smoothly damps high-frequency content using

$var element = document.getElementById("moose-equation-54a3c95b-84f5-4df2-b384-b21dbdb9e48d");katex.render("\\sigma(\\mathbf{k}) = \\exp\\!\\left(-\\alpha\\,\\left(|k_x|^{p} + |k_y|^{p} + |k_z|^{p}\\right)\\right),", element, {displayMode:true,throwOnError:false});$

where $p$ is the exponent and $\alpha$ is the pre-factor controlling the sharpness of the damping. The resulting filter $\sigma$ is applied element-wise in reciprocal space.

Example Input File Syntax

Input Parameters

bufferThe buffer this compute is writing to
C++ Type:std::string
Controllable:No
Description:The buffer this compute is writing to
methodPrefactor
C++ Type:MooseEnum
Options:SHARP, HOULI
Controllable:No
Description:Prefactor

Required Parameters

alpha36Hou-Li filter pre-factor
Default:36
C++ Type:Real
Unit:(no unit assumed)
Controllable:No
Description:Hou-Li filter pre-factor
p16Hou-Li filter exponent
Default:16
C++ Type:Real
Unit:(no unit assumed)
Controllable:No
Description:Hou-Li filter exponent

Optional Parameters

control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Controllable:No
Description:Adds user-defined labels for accessing object parameters via control logic.
enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Controllable:No
Description:Set the enabled status of the MooseObject.

Advanced Parameters

Input Files

(test/tests/cahnhilliard/cahnhilliard_explicit_smooth.i)
(examples/libtorch_kks/KKS_libtorch.i)
(test/tests/anisotropic_grain_growth/analytic_hull.i)
(examples/phase_field_crystal/pfc_fcc_atomic_structure.i)
(examples/anisotropic_grain_growth/isotropic_five_grain.i)
(examples/anisotropic_grain_growth/five_grains_2D.i)
(examples/anisotropic_grain_growth/spherical_grain_Ni.i)
(test/tests/anisotropic_grain_growth/sigma3_circular_grain.i)
(test/tests/tensor_compute/pfc_fcc.i)
(examples/phase_field_crystal/pfc_fcc_atomic_structure_3d.i)

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

(examples/libtorch_kks/KKS_libtorch.i)

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


# Constants for Initial Conditions
r = 30
l = 4.2

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

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

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


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

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

    # automatically create a matching mesh
    mesh_mode = DUMMY
[]

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

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

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

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

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

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

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

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

(test/tests/anisotropic_grain_growth/analytic_hull.i)

interface_width = 0.8
r0 = 4

gbe_max = ${units 0.8 J/m^2}
kappa = ${fparse 0.75 * gbe_max * interface_width}
L = 1

g_eta_expr = '(eta^2*(1-eta^2)^2)'

[Domain]
  dim = 2
  nx = 32
  ny = 32
  xmax = 16
  ymax = 16
  mesh_mode = DUMMY
  floating_precision = SINGLE
[]

[TensorComputes]
  [Initialize]
    [eta]
      type = ParsedCompute
      buffer = 'eta'
      extra_symbols = 'true'
      expression = 'radius:=sqrt((x-${Domain/xmax}/2)^2+(y-${Domain/ymax}/2)^2);0.5 - 0.5*tanh(2*(radius-${r0})/${interface_width})'
    []
    [kappa_linear_term]
      type = ReciprocalLaplacianFactor
      buffer = kappa_linear_term
      factor = ${fparse ${L} * ${kappa}}
    []
    [smooth]
      type = DeAliasingTensor
      method = HOULI
      buffer = smooth
    []
  []
  [Solve]
    [grad_eta_vector]
      type = GradientVector
      buffer = grad_eta
      input = eta
    []
    [gb_energy]
      type = AnisotropicGBEnergy
      buffer = gb_energy
      gb_gradient_buffer = grad_eta
      dsigma_dn = dsigma_dn
      libtorch_model_file = 'marlin:anisotropic_gb_energy/analytic_sigma.pt'
      interface_width = ${interface_width}
      chunk_size = 256
    []
    [bulk_driving_force]
      type = ParsedCompute
      buffer = 'bulk_driving_force'
      expression = '6 * gb_energy * ${g_eta_expr} / ${interface_width}'
      inputs = 'mu eta gb_energy'
      derivatives = 'eta'
    []
    [dmu_dn]
      type = ParsedCompute
      buffer = 'dmu_dn'
      expression = '6 * dsigma_dn / ${interface_width}'
      inputs = 'dsigma_dn'
    []
    [g]
      type = ParsedCompute
      buffer = 'g'
      expression = '${g_eta_expr}'
      inputs = 'eta'
    []
    [anisotropy_torque]
      type = AnisotropyTorque
      buffer = anisotropy_torque
      dmu_dn = dmu_dn
      g = g
    []
    [total_mu_driving_force]
      type = ParsedCompute
      buffer = 'total_mu_driving_force'
      expression = '-${L}*(bulk_driving_force - anisotropy_torque)'
      inputs = 'bulk_driving_force anisotropy_torque'
    []

    [nonlinear]
      type = ForwardFFT
      buffer = 'NL'
      input = 'total_mu_driving_force'
    []
    [smooth_nonlinear]
      type = ParsedCompute
      buffer = 'smooth_NL'
      expression = 'smooth * NL'
      inputs = 'smooth NL'
    []

    [etabar]
      type = ForwardFFT
      buffer = etabar
      input = eta
    []

    [gradient_energy]
      type = GradientEnergyDensity
      buffer = gradient_energy
      gradient_buffers = 'grad_eta'
      kappa = ${kappa}
    []

    [interface_energy]
      type = ParsedCompute
      buffer = interface_energy
      expression = 'g * 6 * gb_energy / ${interface_width} + gradient_energy'
      inputs = 'g gb_energy gradient_energy'
    []
  []
[]

[TensorSolver]
  type = AdamsBashforthMoulton
  buffer = 'eta'
  reciprocal_buffer = 'etabar'
  linear_reciprocal = 'kappa_linear_term'
  nonlinear_reciprocal = 'smooth_NL'
  substeps = 1e2
  predictor_order = 1
  corrector_order = 1
  corrector_steps = 1
[]

[Postprocessors]
  [total_gb_energy]
    type = TensorIntegralPostprocessor
    buffer = interface_energy
  []
[]

[TensorOutputs]
  active = ''
  [xdmf]
    type = XDMFTensorOutput
    buffer = 'eta gb_energy'
    output_mode = 'NODE NODE'
    enable_hdf5 = true
    transpose = false
  []
[]

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

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

(examples/phase_field_crystal/pfc_fcc_atomic_structure.i)

# FCC PFC test designed to clearly show atomic-scale density fluctuations
# Using larger epsilon for stronger density modulations
# Based on PhysRevE.81.061601 two-mode FCC PFC model

# Grid parameters - fine resolution to resolve atomic peaks
N = 4  # Grid points per dimension

# PFC model parameters - LARGER epsilon for clearer atomic structure
psi_mean = -0.2         # Mean density
epsilon = 0.15          # Larger undercooling for stronger modulations
R1_param = 0.0          # Two-mode coupling (0 for max FCC stability)
Q1_param = ${fparse sqrt(4.0/3.0)}  # FCC wave number ratio

# Domain size: want several unit cells but fine enough resolution
# The (111) wavelength is 2*pi/q0 = 2*pi*sqrt(3/4) ≈ 5.44 (dimensionless)
# Let's fit about 'N' wavelengths in each direction
Lx = ${fparse N * 2 * pi / ${Q1_param}}
Ly = ${Lx}

[Domain]
    dim = 2
    nx = ${fparse ${N} * 8 }
    ny = ${fparse ${N} * 8 }
    xmax = ${Lx}
    ymax = ${Ly}
[]

[TensorComputes]
    [Initialize]
        # Random perturbations to seed crystal growth
        [psi]
            type = RandomTensor
            buffer = psi
            max = ${fparse ${psi_mean} + 0.01}
            min = ${fparse ${psi_mean} - 0.01}
            seed = 12345
        []

        # Linear operator for FCC
        [linear]
            type = FCCPFCLinear
            buffer = 'linear'
            eps = ${epsilon}
            R1 = ${R1_param}
            Q1 = ${Q1_param}
            mobility = 1.0
        []

        # Dealiasing for cubic nonlinearity
        [smooth_operator]
            type = DeAliasingTensor
            buffer = smooth_operator
            method = HOULI
        []
    []

    [Solve]
        # Nonlinear term
        [nl_div_psi_cubed]
            type = FCCPFCNonlinear
            buffer = NL
            psi = psi
            dealiasing = smooth_operator
            mobility = 1.0
        []

        # FFT for spectral solver
        [psi_hat]
            type = ForwardFFT
            buffer = psi_hat
            input = psi
        []
    []
[]

[TensorSolver]
    type = AdamsBashforthMoulton
    buffer = 'psi'
    linear_reciprocal = 'linear'
    nonlinear_reciprocal = 'NL'
    reciprocal_buffer = 'psi_hat'
    corrector_order = 1
    corrector_steps = 3
    predictor_order = 1
    substeps = 10000
[]

[Postprocessors]
    [max_psi]
        type = TensorExtremeValuePostprocessor
        buffer = psi
        value_type = MAX
        execute_on = 'initial timestep_end'
    []
    [min_psi]
        type = TensorExtremeValuePostprocessor
        buffer = psi
        value_type = MIN
        execute_on = 'initial timestep_end'
    []
    [mean_psi]
        type = TensorAveragePostprocessor
        buffer = psi
        execute_on = 'initial timestep_end'
    []
[]

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

[Executioner]
    type = Transient
    num_steps = 1000
    dt = 20  # Larger timestep for faster evolution
[]

[Outputs]
    perf_graph = true
    csv = true
[]

(examples/anisotropic_grain_growth/isotropic_five_grain.i)

interface_width = 1.6
side = 80

gbe_max = '${units 0.65 J/m^2}'
L = '${fparse 1.0 * 1.6 / ${interface_width} }'

# ------------------------------------------------------------------
# Isotropic reference case.
#
# The GB energy is fixed at sigma = 0.65 J/m^2 for all boundaries,
# equal to gbe_max. Consequently gamma_ij = 1.5 for all pairs (the
# symmetric-well value), there are no torque terms, the mobility is
# uniform (L_NL == L), and the mobility-split correction
# (L_NL - L)*kappa*laplacian(gr_i) vanishes identically.
# ------------------------------------------------------------------

# ------------------------------------------------------------------
# Periodic Voronoi-like IC for 5 grains.
#
# Each grain is seeded at a point (xi, yi). We use a smooth *periodic*
# distance metric based on sine functions,
#     d_i^2 = (side/pi)^2 * [ sin(pi*(x-xi)/side)^2 + sin(pi*(y-yi)/side)^2 ]
# which reduces to the Euclidean distance near the seed and wraps
# correctly across the periodic boundaries.
#
# The order parameters are then a softmax over the negative distances,
#     gr_i = exp(-a*d_i) / sum_j exp(-a*d_j),   a = 4/interface_width
# For any two competing grains this reduces exactly to the usual
# 0.5*(1 - tanh(2*(d_i - d_j)/interface_width)) profile, so the
# interfaces come out with the correct width, the partition of unity
# sum_i gr_i = 1 holds identically, and the cell structure is the
# (diffuse) Voronoi tessellation of the seeds under the periodic metric.
# ------------------------------------------------------------------
aexp = '${fparse 4.0 / ${interface_width} }'
lp = '${fparse ${side} / pi }'
kp = '${fparse pi / ${side} }'

# seed points (spread out over the periodic domain)
d0 = 'sqrt((${lp}*sin(${kp}*(x-8)))^2  + (${lp}*sin(${kp}*(y-6)))^2)'
d1 = 'sqrt((${lp}*sin(${kp}*(x-28)))^2 + (${lp}*sin(${kp}*(y-10)))^2)'
d2 = 'sqrt((${lp}*sin(${kp}*(x-20)))^2 + (${lp}*sin(${kp}*(y-24)))^2)'
d3 = 'sqrt((${lp}*sin(${kp}*(x-4)))^2  + (${lp}*sin(${kp}*(y-28)))^2)'
d4 = 'sqrt((${lp}*sin(${kp}*(x-34)))^2 + (${lp}*sin(${kp}*(y-33)))^2)'

g_gamma0 = '${fparse sqrt(2) / 3 }' # g(gamma=1.5)
f0_gamma0 = ${fparse 1/8}
kappa = '${fparse gbe_max * interface_width * sqrt(f0_gamma0) / g_gamma0 }'
mu = '${fparse gbe_max  / (g_gamma0 * interface_width * sqrt(f0_gamma0) ) }'

gamma = 1.5

f0 = '((gr0^4/4 - gr0^2/2) + (gr1^4/4 - gr1^2/2) + (gr2^4/4 - gr2^2/2) + (gr3^4/4 - gr3^2/2) + (gr4^4/4 - gr4^2/2) + ${gamma}*(gr0^2*gr1^2 + gr0^2*gr2^2 + gr0^2*gr3^2 + gr0^2*gr4^2 + gr1^2*gr2^2 + gr1^2*gr3^2 + gr1^2*gr4^2 + gr2^2*gr3^2 + gr2^2*gr4^2 + gr3^2*gr4^2) + 0.25)'

[Domain]
    dim = 2
    nx = 200
    ny = 200
    xmax = ${side}
    ymax = ${side}
    mesh_mode = DUMMY
[]

[TensorComputes]
    [Initialize]
        # unnormalized softmax weights
        [c0]
            type = ParsedCompute
            buffer = c0
            expression = 'exp(-${aexp}*${d0})'
            extra_symbols = true
        []
        [c1]
            type = ParsedCompute
            buffer = c1
            expression = 'exp(-${aexp}*${d1})'
            extra_symbols = true
        []
        [c2]
            type = ParsedCompute
            buffer = c2
            expression = 'exp(-${aexp}*${d2})'
            extra_symbols = true
        []
        [c3]
            type = ParsedCompute
            buffer = c3
            expression = 'exp(-${aexp}*${d3})'
            extra_symbols = true
        []
        [c4]
            type = ParsedCompute
            buffer = c4
            expression = 'exp(-${aexp}*${d4})'
            extra_symbols = true
        []

        # normalized order parameters (partition of unity by construction)
        [gr0]
            type = ParsedCompute
            buffer = gr0
            expression = 'c0/(c0+c1+c2+c3+c4)'
            inputs = 'c0 c1 c2 c3 c4'
        []
        [gr1]
            type = ParsedCompute
            buffer = gr1
            expression = 'c1/(c0+c1+c2+c3+c4)'
            inputs = 'c0 c1 c2 c3 c4'
        []
        [gr2]
            type = ParsedCompute
            buffer = gr2
            expression = 'c2/(c0+c1+c2+c3+c4)'
            inputs = 'c0 c1 c2 c3 c4'
        []
        [gr3]
            type = ParsedCompute
            buffer = gr3
            expression = 'c3/(c0+c1+c2+c3+c4)'
            inputs = 'c0 c1 c2 c3 c4'
        []
        [gr4]
            type = ParsedCompute
            buffer = gr4
            expression = '1 - gr0 - gr1 - gr2 - gr3'
            inputs = 'gr0 gr1 gr2 gr3'
        []

        [L_kappa_laplacian]
            type = ReciprocalLaplacianFactor
            buffer = 'L_kappa_laplacian'
            factor = '${fparse ${L} * ${kappa}}'
        []
        [smooth]
            type = DeAliasingTensor
            method = HOULI
            buffer = smooth
        []
    []
    [Solve]
        [bnds]
            type = ParsedCompute
            buffer = 'bnds'
            inputs = 'gr0 gr1 gr2 gr3 gr4'
            expression = 'gr0^2 + gr1^2 + gr2^2 + gr3^2 + gr4^2'
        []

        ## FFT
        [gr0_bar]
            type = ForwardFFT
            buffer = gr0_bar
            input = gr0
        []
        [gr1_bar]
            type = ForwardFFT
            buffer = gr1_bar
            input = gr1
        []
        [gr2_bar]
            type = ForwardFFT
            buffer = gr2_bar
            input = gr2
        []
        [gr3_bar]
            type = ForwardFFT
            buffer = gr3_bar
            input = gr3
        []
        [gr4_bar]
            type = ForwardFFT
            buffer = gr4_bar
            input = gr4
        []

        ## bulk driving forces (constant mobility L, constant gamma = ${gamma})
        [gr0_bulk_term]
            type = ParsedCompute
            buffer = 'gr0_bulk_term'
            expression = '-${L}*${mu} * (gr0^3 - gr0 + 2*${gamma}*gr0*(gr1^2 + gr2^2 + gr3^2 + gr4^2))'
            inputs = 'gr0 gr1 gr2 gr3 gr4'
        []
        [gr1_bulk_term]
            type = ParsedCompute
            buffer = 'gr1_bulk_term'
            expression = '-${L}*${mu} * (gr1^3 - gr1 + 2*${gamma}*gr1*(gr0^2 + gr2^2 + gr3^2 + gr4^2))'
            inputs = 'gr0 gr1 gr2 gr3 gr4'
        []
        [gr2_bulk_term]
            type = ParsedCompute
            buffer = 'gr2_bulk_term'
            expression = '-${L}*${mu} * (gr2^3 - gr2 + 2*${gamma}*gr2*(gr0^2 + gr1^2 + gr3^2 + gr4^2))'
            inputs = 'gr0 gr1 gr2 gr3 gr4'
        []
        [gr3_bulk_term]
            type = ParsedCompute
            buffer = 'gr3_bulk_term'
            expression = '-${L}*${mu} * (gr3^3 - gr3 + 2*${gamma}*gr3*(gr0^2 + gr1^2 + gr2^2 + gr4^2))'
            inputs = 'gr0 gr1 gr2 gr3 gr4'
        []
        [gr4_bulk_term]
            type = ParsedCompute
            buffer = 'gr4_bulk_term'
            expression = '-${L}*${mu} * (gr4^3 - gr4 + 2*${gamma}*gr4*(gr0^2 + gr1^2 + gr2^2 + gr3^2))'
            inputs = 'gr0 gr1 gr2 gr3 gr4'
        []

        [NL_gr0]
            type = ForwardFFT
            buffer = NL_gr0
            input = 'gr0_bulk_term'
        []
        [NL_gr1]
            type = ForwardFFT
            buffer = NL_gr1
            input = 'gr1_bulk_term'
        []
        [NL_gr2]
            type = ForwardFFT
            buffer = NL_gr2
            input = 'gr2_bulk_term'
        []
        [NL_gr3]
            type = ForwardFFT
            buffer = NL_gr3
            input = 'gr3_bulk_term'
        []
        [NL_gr4]
            type = ForwardFFT
            buffer = NL_gr4
            input = 'gr4_bulk_term'
        []

        [NL_gr0_smooth]
            type = ParsedCompute
            buffer = NL_gr0_smooth
            expression = 'smooth*NL_gr0'
            inputs = 'smooth NL_gr0'
        []
        [NL_gr1_smooth]
            type = ParsedCompute
            buffer = NL_gr1_smooth
            expression = 'smooth*NL_gr1'
            inputs = 'smooth NL_gr1'
        []
        [NL_gr2_smooth]
            type = ParsedCompute
            buffer = NL_gr2_smooth
            expression = 'smooth*NL_gr2'
            inputs = 'smooth NL_gr2'
        []
        [NL_gr3_smooth]
            type = ParsedCompute
            buffer = NL_gr3_smooth
            expression = 'smooth*NL_gr3'
            inputs = 'smooth NL_gr3'
        []
        [NL_gr4_smooth]
            type = ParsedCompute
            buffer = NL_gr4_smooth
            expression = 'smooth*NL_gr4'
            inputs = 'smooth NL_gr4'
        []

        # Properties for output
        [grain_ID]
            type = ParsedCompute
            buffer = 'grain_ID'
            expression = 'if(gr0>=max(gr1,max(gr2,max(gr3,gr4))), 0, if(gr1>=max(gr2,max(gr3,gr4)), 1, if(gr2>=max(gr3,gr4), 2, if(gr3>=gr4, 3, 4))))'
            inputs = 'gr0 gr1 gr2 gr3 gr4'
        []

        [grad_gr0]
            type = GradientVector
            buffer = grad_gr0
            input = gr0
        []
        [grad_gr1]
            type = GradientVector
            buffer = grad_gr1
            input = gr1
        []
        [grad_gr2]
            type = GradientVector
            buffer = grad_gr2
            input = gr2
        []
        [grad_gr3]
            type = GradientVector
            buffer = grad_gr3
            input = gr3
        []
        [grad_gr4]
            type = GradientVector
            buffer = grad_gr4
            input = gr4
        []

        [gradient_energy]
            type = GradientEnergyDensity
            buffer = gradient_energy
            gradient_buffers = 'grad_gr0 grad_gr1 grad_gr2 grad_gr3 grad_gr4'
            kappa = ${kappa}
        []
        [total_energy]
            type = ParsedCompute
            buffer = 'total_energy'
            expression = '${mu} * (${f0}) + gradient_energy'
            inputs = 'gr0 gr1 gr2 gr3 gr4 gradient_energy'
        []
    []
[]

[Postprocessors]
    [total_gb_energy]
        type = TensorIntegralPostprocessor
        buffer = total_energy
    []
[]

[TensorSolver]
    type = AdamsBashforthMoulton
    buffer = 'gr0 gr1 gr2 gr3 gr4'
    reciprocal_buffer = 'gr0_bar gr1_bar gr2_bar gr3_bar gr4_bar'
    linear_reciprocal = 'L_kappa_laplacian L_kappa_laplacian L_kappa_laplacian L_kappa_laplacian L_kappa_laplacian'
    nonlinear_reciprocal = 'NL_gr0_smooth NL_gr1_smooth NL_gr2_smooth NL_gr3_smooth NL_gr4_smooth'
    substeps = 1000
    predictor_order = 1
    corrector_order = 1
    corrector_steps = 1
[]

[TensorOutputs]
    [xdmf]
        type = XDMFTensorOutput
        buffer = 'gr0 gr1 gr2 gr3 gr4 bnds grain_ID'
        output_mode = 'NODE NODE NODE NODE NODE NODE CELL'
        enable_hdf5 = true
        transpose = true
    []
[]

[Executioner]
    type = Transient
    dt = 5
    num_steps = 50
[]

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

(examples/anisotropic_grain_growth/five_grains_2D.i)

interface_width = 1.6
side = 80

gbe_max = '${units 1.24 J/m^2}'
L = '${fparse 1.0 * 1.6 / ${interface_width} }'

# ------------------------------------------------------------------
# Periodic Voronoi-like IC for 5 grains.
#
# Each grain is seeded at a point (xi, yi). We use a smooth *periodic*
# distance metric based on sine functions,
#     d_i^2 = (side/pi)^2 * [ sin(pi*(x-xi)/side)^2 + sin(pi*(y-yi)/side)^2 ]
# which reduces to the Euclidean distance near the seed and wraps
# correctly across the periodic boundaries.
#
# The order parameters are then a softmax over the negative distances,
#     gr_i = exp(-a*d_i) / sum_j exp(-a*d_j),   a = 4/interface_width
# For any two competing grains this reduces exactly to the usual
# 0.5*(1 - tanh(2*(d_i - d_j)/interface_width)) profile, so the
# interfaces come out with the correct width, the partition of unity
# sum_i gr_i = 1 holds identically, and the cell structure is the
# (diffuse) Voronoi tessellation of the seeds under the periodic metric.
# ------------------------------------------------------------------
aexp = '${fparse 4.0 / ${interface_width} }'
lp = '${fparse ${side} / pi }'
kp = '${fparse pi / ${side} }'

# seed points (spread out over the periodic 40x40 domain)
d0 = 'sqrt((${lp}*sin(${kp}*(x-8)))^2  + (${lp}*sin(${kp}*(y-6)))^2)'
d1 = 'sqrt((${lp}*sin(${kp}*(x-28)))^2 + (${lp}*sin(${kp}*(y-10)))^2)'
d2 = 'sqrt((${lp}*sin(${kp}*(x-20)))^2 + (${lp}*sin(${kp}*(y-24)))^2)'
d3 = 'sqrt((${lp}*sin(${kp}*(x-4)))^2  + (${lp}*sin(${kp}*(y-28)))^2)'
d4 = 'sqrt((${lp}*sin(${kp}*(x-34)))^2 + (${lp}*sin(${kp}*(y-33)))^2)'

g_gamma0 = '${fparse sqrt(2) / 3 }' # g(gamma=1.5)
f0_gamma0 = ${fparse 1/8}
kappa = '${fparse gbe_max * interface_width * sqrt(f0_gamma0) / g_gamma0 }'
mu = '${fparse gbe_max  / (g_gamma0 * interface_width * sqrt(f0_gamma0) ) }'

g_to_gamma_func = '(1.0/(-3.0944 * g^8 -1.8169*g^6 + 10.323 * g^4 - 8.1819*g^2 + 2.0033))'
dgamma_dsigma = '(g * (24.7552*g^6 + 10.9014*g^4 - 41.292*g^2 + 16.3638) / (-3.0944*g^8 - 1.8169*g^6 + 10.323*g^4 - 8.1819*g^2 + 2.0033)^2 / sqrt(${kappa} * ${mu}) )'

f0 = '((gr0^4/4 - gr0^2/2) + (gr1^4/4 - gr1^2/2) + (gr2^4/4 - gr2^2/2) + (gr3^4/4 - gr3^2/2) + (gr4^4/4 - gr4^2/2) + gamma_01*gr0^2*gr1^2 + gamma_02*gr0^2*gr2^2 + gamma_03*gr0^2*gr3^2 + gamma_04*gr0^2*gr4^2 + gamma_12*gr1^2*gr2^2 + gamma_13*gr1^2*gr3^2 + gamma_14*gr1^2*gr4^2 + gamma_23*gr2^2*gr3^2 + gamma_24*gr2^2*gr4^2 + gamma_34*gr3^2*gr4^2 + 0.25)'

[Domain]
    dim = 2
    nx = 200
    ny = 200
    xmax = ${side}
    ymax = ${side}
    mesh_mode = DUMMY
[]

[TensorComputes]
    [Initialize]
        # unnormalized softmax weights
        [c0]
            type = ParsedCompute
            buffer = c0
            expression = 'exp(-${aexp}*${d0})'
            extra_symbols = true
        []
        [c1]
            type = ParsedCompute
            buffer = c1
            expression = 'exp(-${aexp}*${d1})'
            extra_symbols = true
        []
        [c2]
            type = ParsedCompute
            buffer = c2
            expression = 'exp(-${aexp}*${d2})'
            extra_symbols = true
        []
        [c3]
            type = ParsedCompute
            buffer = c3
            expression = 'exp(-${aexp}*${d3})'
            extra_symbols = true
        []
        [c4]
            type = ParsedCompute
            buffer = c4
            expression = 'exp(-${aexp}*${d4})'
            extra_symbols = true
        []

        # normalized order parameters (partition of unity by construction)
        [gr0]
            type = ParsedCompute
            buffer = gr0
            expression = 'c0/(c0+c1+c2+c3+c4)'
            inputs = 'c0 c1 c2 c3 c4'
        []
        [gr1]
            type = ParsedCompute
            buffer = gr1
            expression = 'c1/(c0+c1+c2+c3+c4)'
            inputs = 'c0 c1 c2 c3 c4'
        []
        [gr2]
            type = ParsedCompute
            buffer = gr2
            expression = 'c2/(c0+c1+c2+c3+c4)'
            inputs = 'c0 c1 c2 c3 c4'
        []
        [gr3]
            type = ParsedCompute
            buffer = gr3
            expression = 'c3/(c0+c1+c2+c3+c4)'
            inputs = 'c0 c1 c2 c3 c4'
        []
        [gr4]
            type = ParsedCompute
            buffer = gr4
            expression = '1 - gr0 - gr1 - gr2 - gr3'
            inputs = 'gr0 gr1 gr2 gr3'
        []

        [L_kappa_laplacian]
            type = ReciprocalLaplacianFactor
            buffer = 'L_kappa_laplacian'
            factor = '${fparse ${L} * ${kappa}}'
        []
        [kappa_laplacian]
            type = ReciprocalLaplacianFactor
            buffer = kappa_laplacian
            factor = '${kappa}'
        []
        [smooth]
            type = DeAliasingTensor
            method = HOULI
            buffer = smooth
        []
    []
    [Solve]
        [bnds]
            type = ParsedCompute
            buffer = 'bnds'
            inputs = 'gr0 gr1 gr2 gr3 gr4'
            expression = 'gr0^2 + gr1^2 + gr2^2 + gr3^2 + gr4^2'
        []
            ## FFT
        [gr0_bar]
            type = ForwardFFT
            buffer = gr0_bar
            input = gr0
        []
        [gr1_bar]
            type = ForwardFFT
            buffer = gr1_bar
            input = gr1
        []
        [gr2_bar]
            type = ForwardFFT
            buffer = gr2_bar
            input = gr2
        []
        [gr3_bar]
            type = ForwardFFT
            buffer = gr3_bar
            input = gr3
        []
        [gr4_bar]
            type = ForwardFFT
            buffer = gr4_bar
            input = gr4
        []

        [grad_gr0]
            type = GradientVector
            buffer = grad_gr0
            input = gr0
        []
        [grad_gr1]
            type = GradientVector
            buffer = grad_gr1
            input = gr1
        []
        [grad_gr2]
            type = GradientVector
            buffer = grad_gr2
            input = gr2
        []
        [grad_gr3]
            type = GradientVector
            buffer = grad_gr3
            input = gr3
        []
        [grad_gr4]
            type = GradientVector
            buffer = grad_gr4
            input = gr4
        []

        ## Pairwise anisotropic GB energies (10 pairs)
        # NOTE: the per-pair libtorch misorientation models (Ri_Rj.pt) will be
        # developed separately; placeholder file names below follow the
        # existing convention and location.
        [sigma_gr0_gr1]
            type = PairwiseAnisotropicGBEnergy
            buffer = 'sigma_gr0_gr1'
            dsigma_dgrad_grain1 = 'dsigma_gr0_gr1_dgrad_gr0'
            dsigma_dgrad_grain2 = 'dsigma_gr0_gr1_dgrad_gr1'
            grad_grain1_buffer = 'grad_gr0'
            grad_grain2_buffer = 'grad_gr1'
            interface_width = '${interface_width}'
            libtorch_model_file = 'marlin:anisotropic_gb_energy/Ni_R1_R2.pt'
            chunk_size = 1e5
        []
        [sigma_gr0_gr2]
            type = PairwiseAnisotropicGBEnergy
            buffer = 'sigma_gr0_gr2'
            dsigma_dgrad_grain1 = 'dsigma_gr0_gr2_dgrad_gr0'
            dsigma_dgrad_grain2 = 'dsigma_gr0_gr2_dgrad_gr2'
            grad_grain1_buffer = 'grad_gr0'
            grad_grain2_buffer = 'grad_gr2'
            interface_width = '${interface_width}'
            libtorch_model_file = 'marlin:anisotropic_gb_energy/Ni_R1_R3.pt'
            chunk_size = 1e5
        []
        [sigma_gr0_gr3]
            type = PairwiseAnisotropicGBEnergy
            buffer = 'sigma_gr0_gr3'
            dsigma_dgrad_grain1 = 'dsigma_gr0_gr3_dgrad_gr0'
            dsigma_dgrad_grain2 = 'dsigma_gr0_gr3_dgrad_gr3'
            grad_grain1_buffer = 'grad_gr0'
            grad_grain2_buffer = 'grad_gr3'
            interface_width = '${interface_width}'
            libtorch_model_file = 'marlin:anisotropic_gb_energy/Ni_R1_R4.pt'
            chunk_size = 1e5
        []
        [sigma_gr0_gr4]
            type = PairwiseAnisotropicGBEnergy
            buffer = 'sigma_gr0_gr4'
            dsigma_dgrad_grain1 = 'dsigma_gr0_gr4_dgrad_gr0'
            dsigma_dgrad_grain2 = 'dsigma_gr0_gr4_dgrad_gr4'
            grad_grain1_buffer = 'grad_gr0'
            grad_grain2_buffer = 'grad_gr4'
            interface_width = '${interface_width}'
            libtorch_model_file = 'marlin:anisotropic_gb_energy/Ni_R1_R5.pt'
            chunk_size = 1e5
        []
        [sigma_gr1_gr2]
            type = PairwiseAnisotropicGBEnergy
            buffer = 'sigma_gr1_gr2'
            dsigma_dgrad_grain1 = 'dsigma_gr1_gr2_dgrad_gr1'
            dsigma_dgrad_grain2 = 'dsigma_gr1_gr2_dgrad_gr2'
            grad_grain1_buffer = 'grad_gr1'
            grad_grain2_buffer = 'grad_gr2'
            interface_width = '${interface_width}'
            libtorch_model_file = 'marlin:anisotropic_gb_energy/Ni_R2_R3.pt'
            chunk_size = 1e5
        []
        [sigma_gr1_gr3]
            type = PairwiseAnisotropicGBEnergy
            buffer = 'sigma_gr1_gr3'
            dsigma_dgrad_grain1 = 'dsigma_gr1_gr3_dgrad_gr1'
            dsigma_dgrad_grain2 = 'dsigma_gr1_gr3_dgrad_gr3'
            grad_grain1_buffer = 'grad_gr1'
            grad_grain2_buffer = 'grad_gr3'
            interface_width = '${interface_width}'
            libtorch_model_file = 'marlin:anisotropic_gb_energy/Ni_R2_R4.pt'
            chunk_size = 1e5
        []
        [sigma_gr1_gr4]
            type = PairwiseAnisotropicGBEnergy
            buffer = 'sigma_gr1_gr4'
            dsigma_dgrad_grain1 = 'dsigma_gr1_gr4_dgrad_gr1'
            dsigma_dgrad_grain2 = 'dsigma_gr1_gr4_dgrad_gr4'
            grad_grain1_buffer = 'grad_gr1'
            grad_grain2_buffer = 'grad_gr4'
            interface_width = '${interface_width}'
            libtorch_model_file = 'marlin:anisotropic_gb_energy/Ni_R2_R5.pt'
            chunk_size = 1e5
        []
        [sigma_gr2_gr3]
            type = PairwiseAnisotropicGBEnergy
            buffer = 'sigma_gr2_gr3'
            dsigma_dgrad_grain1 = 'dsigma_gr2_gr3_dgrad_gr2'
            dsigma_dgrad_grain2 = 'dsigma_gr2_gr3_dgrad_gr3'
            grad_grain1_buffer = 'grad_gr2'
            grad_grain2_buffer = 'grad_gr3'
            interface_width = '${interface_width}'
            libtorch_model_file = 'marlin:anisotropic_gb_energy/Ni_R3_R4.pt'
            chunk_size = 1e5
        []
        [sigma_gr2_gr4]
            type = PairwiseAnisotropicGBEnergy
            buffer = 'sigma_gr2_gr4'
            dsigma_dgrad_grain1 = 'dsigma_gr2_gr4_dgrad_gr2'
            dsigma_dgrad_grain2 = 'dsigma_gr2_gr4_dgrad_gr4'
            grad_grain1_buffer = 'grad_gr2'
            grad_grain2_buffer = 'grad_gr4'
            interface_width = '${interface_width}'
            libtorch_model_file = 'marlin:anisotropic_gb_energy/Ni_R3_R5.pt'
            chunk_size = 1e5
        []
        [sigma_gr3_gr4]
            type = PairwiseAnisotropicGBEnergy
            buffer = 'sigma_gr3_gr4'
            dsigma_dgrad_grain1 = 'dsigma_gr3_gr4_dgrad_gr3'
            dsigma_dgrad_grain2 = 'dsigma_gr3_gr4_dgrad_gr4'
            grad_grain1_buffer = 'grad_gr3'
            grad_grain2_buffer = 'grad_gr4'
            interface_width = '${interface_width}'
            libtorch_model_file = 'marlin:anisotropic_gb_energy/Ni_R4_R5.pt'
            chunk_size = 1e5
        []

        ## gamma_ij(sigma_ij)
        [gamma_01]
            type = ParsedCompute
            buffer = 'gamma_01'
            expression = 'g:=sigma_gr0_gr1/sqrt(${kappa}*${mu});${g_to_gamma_func}'
            inputs = 'sigma_gr0_gr1'
        []
        [gamma_02]
            type = ParsedCompute
            buffer = 'gamma_02'
            expression = 'g:=sigma_gr0_gr2/sqrt(${kappa}*${mu});${g_to_gamma_func}'
            inputs = 'sigma_gr0_gr2'
        []
        [gamma_03]
            type = ParsedCompute
            buffer = 'gamma_03'
            expression = 'g:=sigma_gr0_gr3/sqrt(${kappa}*${mu});${g_to_gamma_func}'
            inputs = 'sigma_gr0_gr3'
        []
        [gamma_04]
            type = ParsedCompute
            buffer = 'gamma_04'
            expression = 'g:=sigma_gr0_gr4/sqrt(${kappa}*${mu});${g_to_gamma_func}'
            inputs = 'sigma_gr0_gr4'
        []
        [gamma_12]
            type = ParsedCompute
            buffer = 'gamma_12'
            expression = 'g:=sigma_gr1_gr2/sqrt(${kappa}*${mu});${g_to_gamma_func}'
            inputs = 'sigma_gr1_gr2'
        []
        [gamma_13]
            type = ParsedCompute
            buffer = 'gamma_13'
            expression = 'g:=sigma_gr1_gr3/sqrt(${kappa}*${mu});${g_to_gamma_func}'
            inputs = 'sigma_gr1_gr3'
        []
        [gamma_14]
            type = ParsedCompute
            buffer = 'gamma_14'
            expression = 'g:=sigma_gr1_gr4/sqrt(${kappa}*${mu});${g_to_gamma_func}'
            inputs = 'sigma_gr1_gr4'
        []
        [gamma_23]
            type = ParsedCompute
            buffer = 'gamma_23'
            expression = 'g:=sigma_gr2_gr3/sqrt(${kappa}*${mu});${g_to_gamma_func}'
            inputs = 'sigma_gr2_gr3'
        []
        [gamma_24]
            type = ParsedCompute
            buffer = 'gamma_24'
            expression = 'g:=sigma_gr2_gr4/sqrt(${kappa}*${mu});${g_to_gamma_func}'
            inputs = 'sigma_gr2_gr4'
        []
        [gamma_34]
            type = ParsedCompute
            buffer = 'gamma_34'
            expression = 'g:=sigma_gr3_gr4/sqrt(${kappa}*${mu});${g_to_gamma_func}'
            inputs = 'sigma_gr3_gr4'
        []

        ## torque prefactors dgamma/dsigma * gr_i^2 * gr_j^2
        [coeff_01]
            type = ParsedCompute
            buffer = 'coeff_01'
            expression = 'g:=sigma_gr0_gr1/sqrt(${kappa}*${mu});${dgamma_dsigma} * gr0^2 *gr1^2'
            inputs = 'sigma_gr0_gr1 gr0 gr1'
        []
        [coeff_02]
            type = ParsedCompute
            buffer = 'coeff_02'
            expression = 'g:=sigma_gr0_gr2/sqrt(${kappa}*${mu});${dgamma_dsigma} * gr0^2 *gr2^2'
            inputs = 'sigma_gr0_gr2 gr0 gr2'
        []
        [coeff_03]
            type = ParsedCompute
            buffer = 'coeff_03'
            expression = 'g:=sigma_gr0_gr3/sqrt(${kappa}*${mu});${dgamma_dsigma} * gr0^2 *gr3^2'
            inputs = 'sigma_gr0_gr3 gr0 gr3'
        []
        [coeff_04]
            type = ParsedCompute
            buffer = 'coeff_04'
            expression = 'g:=sigma_gr0_gr4/sqrt(${kappa}*${mu});${dgamma_dsigma} * gr0^2 *gr4^2'
            inputs = 'sigma_gr0_gr4 gr0 gr4'
        []
        [coeff_12]
            type = ParsedCompute
            buffer = 'coeff_12'
            expression = 'g:=sigma_gr1_gr2/sqrt(${kappa}*${mu});${dgamma_dsigma} * gr1^2 *gr2^2'
            inputs = 'sigma_gr1_gr2 gr1 gr2'
        []
        [coeff_13]
            type = ParsedCompute
            buffer = 'coeff_13'
            expression = 'g:=sigma_gr1_gr3/sqrt(${kappa}*${mu});${dgamma_dsigma} * gr1^2 *gr3^2'
            inputs = 'sigma_gr1_gr3 gr1 gr3'
        []
        [coeff_14]
            type = ParsedCompute
            buffer = 'coeff_14'
            expression = 'g:=sigma_gr1_gr4/sqrt(${kappa}*${mu});${dgamma_dsigma} * gr1^2 *gr4^2'
            inputs = 'sigma_gr1_gr4 gr1 gr4'
        []
        [coeff_23]
            type = ParsedCompute
            buffer = 'coeff_23'
            expression = 'g:=sigma_gr2_gr3/sqrt(${kappa}*${mu});${dgamma_dsigma} * gr2^2 *gr3^2'
            inputs = 'sigma_gr2_gr3 gr2 gr3'
        []
        [coeff_24]
            type = ParsedCompute
            buffer = 'coeff_24'
            expression = 'g:=sigma_gr2_gr4/sqrt(${kappa}*${mu});${dgamma_dsigma} * gr2^2 *gr4^2'
            inputs = 'sigma_gr2_gr4 gr2 gr4'
        []
        [coeff_34]
            type = ParsedCompute
            buffer = 'coeff_34'
            expression = 'g:=sigma_gr3_gr4/sqrt(${kappa}*${mu});${dgamma_dsigma} * gr3^2 *gr4^2'
            inputs = 'sigma_gr3_gr4 gr3 gr4'
        []

        ## anisotropy torques (2 per pair)
        [torque_01_gr0]
            type = AnisotropyTorque
            buffer = torque_01_gr0
            dmu_dn = dsigma_gr0_gr1_dgrad_gr0
            g = 'coeff_01'
        []
        [torque_01_gr1]
            type = AnisotropyTorque
            buffer = torque_01_gr1
            dmu_dn = dsigma_gr0_gr1_dgrad_gr1
            g = 'coeff_01'
        []
        [torque_02_gr0]
            type = AnisotropyTorque
            buffer = torque_02_gr0
            dmu_dn = dsigma_gr0_gr2_dgrad_gr0
            g = 'coeff_02'
        []
        [torque_02_gr2]
            type = AnisotropyTorque
            buffer = torque_02_gr2
            dmu_dn = dsigma_gr0_gr2_dgrad_gr2
            g = 'coeff_02'
        []
        [torque_03_gr0]
            type = AnisotropyTorque
            buffer = torque_03_gr0
            dmu_dn = dsigma_gr0_gr3_dgrad_gr0
            g = 'coeff_03'
        []
        [torque_03_gr3]
            type = AnisotropyTorque
            buffer = torque_03_gr3
            dmu_dn = dsigma_gr0_gr3_dgrad_gr3
            g = 'coeff_03'
        []
        [torque_04_gr0]
            type = AnisotropyTorque
            buffer = torque_04_gr0
            dmu_dn = dsigma_gr0_gr4_dgrad_gr0
            g = 'coeff_04'
        []
        [torque_04_gr4]
            type = AnisotropyTorque
            buffer = torque_04_gr4
            dmu_dn = dsigma_gr0_gr4_dgrad_gr4
            g = 'coeff_04'
        []
        [torque_12_gr1]
            type = AnisotropyTorque
            buffer = torque_12_gr1
            dmu_dn = dsigma_gr1_gr2_dgrad_gr1
            g = 'coeff_12'
        []
        [torque_12_gr2]
            type = AnisotropyTorque
            buffer = torque_12_gr2
            dmu_dn = dsigma_gr1_gr2_dgrad_gr2
            g = 'coeff_12'
        []
        [torque_13_gr1]
            type = AnisotropyTorque
            buffer = torque_13_gr1
            dmu_dn = dsigma_gr1_gr3_dgrad_gr1
            g = 'coeff_13'
        []
        [torque_13_gr3]
            type = AnisotropyTorque
            buffer = torque_13_gr3
            dmu_dn = dsigma_gr1_gr3_dgrad_gr3
            g = 'coeff_13'
        []
        [torque_14_gr1]
            type = AnisotropyTorque
            buffer = torque_14_gr1
            dmu_dn = dsigma_gr1_gr4_dgrad_gr1
            g = 'coeff_14'
        []
        [torque_14_gr4]
            type = AnisotropyTorque
            buffer = torque_14_gr4
            dmu_dn = dsigma_gr1_gr4_dgrad_gr4
            g = 'coeff_14'
        []
        [torque_23_gr2]
            type = AnisotropyTorque
            buffer = torque_23_gr2
            dmu_dn = dsigma_gr2_gr3_dgrad_gr2
            g = 'coeff_23'
        []
        [torque_23_gr3]
            type = AnisotropyTorque
            buffer = torque_23_gr3
            dmu_dn = dsigma_gr2_gr3_dgrad_gr3
            g = 'coeff_23'
        []
        [torque_24_gr2]
            type = AnisotropyTorque
            buffer = torque_24_gr2
            dmu_dn = dsigma_gr2_gr4_dgrad_gr2
            g = 'coeff_24'
        []
        [torque_24_gr4]
            type = AnisotropyTorque
            buffer = torque_24_gr4
            dmu_dn = dsigma_gr2_gr4_dgrad_gr4
            g = 'coeff_24'
        []
        [torque_34_gr3]
            type = AnisotropyTorque
            buffer = torque_34_gr3
            dmu_dn = dsigma_gr3_gr4_dgrad_gr3
            g = 'coeff_34'
        []
        [torque_34_gr4]
            type = AnisotropyTorque
            buffer = torque_34_gr4
            dmu_dn = dsigma_gr3_gr4_dgrad_gr4
            g = 'coeff_34'
        []

        ## interpolated local GB energy and nonlinear mobility
        [sigma]
            type = ParsedCompute
            buffer = 'sigma'
            expression = '(gr0^2*gr1^2*sigma_gr0_gr1 + gr0^2*gr2^2*sigma_gr0_gr2 + gr0^2*gr3^2*sigma_gr0_gr3 + gr0^2*gr4^2*sigma_gr0_gr4 + gr1^2*gr2^2*sigma_gr1_gr2 + gr1^2*gr3^2*sigma_gr1_gr3 + gr1^2*gr4^2*sigma_gr1_gr4 + gr2^2*gr3^2*sigma_gr2_gr3 + gr2^2*gr4^2*sigma_gr2_gr4 + gr3^2*gr4^2*sigma_gr3_gr4)/(gr0^2*gr1^2 + gr0^2*gr2^2 + gr0^2*gr3^2 + gr0^2*gr4^2 + gr1^2*gr2^2 + gr1^2*gr3^2 + gr1^2*gr4^2 + gr2^2*gr3^2 + gr2^2*gr4^2 + gr3^2*gr4^2 + 1e-4)'
            inputs = 'sigma_gr0_gr1 sigma_gr0_gr2 sigma_gr0_gr3 sigma_gr0_gr4 sigma_gr1_gr2 sigma_gr1_gr3 sigma_gr1_gr4 sigma_gr2_gr3 sigma_gr2_gr4 sigma_gr3_gr4 gr0 gr1 gr2 gr3 gr4'
        []
        [L_NL]
            type = ParsedCompute
            buffer = 'L_NL'
            expression = '${L} * sigma / ${gbe_max}'
            inputs = 'sigma'
        []

        ## kappa * laplacian(gr_i) in real space (for the mobility-split correction)
        [kappa_laplacian_gr0_hat]
            type = ParsedCompute
            buffer = kappa_laplacian_gr0_hat
            inputs = 'gr0_bar kappa_laplacian'
            expression = 'gr0_bar * kappa_laplacian'
        []
        [kappa_laplacian_gr1_hat]
            type = ParsedCompute
            buffer = kappa_laplacian_gr1_hat
            inputs = 'gr1_bar kappa_laplacian'
            expression = 'gr1_bar * kappa_laplacian'
        []
        [kappa_laplacian_gr2_hat]
            type = ParsedCompute
            buffer = kappa_laplacian_gr2_hat
            inputs = 'gr2_bar kappa_laplacian'
            expression = 'gr2_bar * kappa_laplacian'
        []
        [kappa_laplacian_gr3_hat]
            type = ParsedCompute
            buffer = kappa_laplacian_gr3_hat
            inputs = 'gr3_bar kappa_laplacian'
            expression = 'gr3_bar * kappa_laplacian'
        []
        [kappa_laplacian_gr4_hat]
            type = ParsedCompute
            buffer = kappa_laplacian_gr4_hat
            inputs = 'gr4_bar kappa_laplacian'
            expression = 'gr4_bar * kappa_laplacian'
        []

        [kappa_laplacian_gr0]
            type = InverseFFT
            buffer = kappa_laplacian_gr0
            input = kappa_laplacian_gr0_hat
        []
        [kappa_laplacian_gr1]
            type = InverseFFT
            buffer = kappa_laplacian_gr1
            input = kappa_laplacian_gr1_hat
        []
        [kappa_laplacian_gr2]
            type = InverseFFT
            buffer = kappa_laplacian_gr2
            input = kappa_laplacian_gr2_hat
        []
        [kappa_laplacian_gr3]
            type = InverseFFT
            buffer = kappa_laplacian_gr3
            input = kappa_laplacian_gr3_hat
        []
        [kappa_laplacian_gr4]
            type = InverseFFT
            buffer = kappa_laplacian_gr4
            input = kappa_laplacian_gr4_hat
        []

        ## bulk driving forces
        [gr0_bulk_term]
            type = ParsedCompute
            buffer = 'gr0_bulk_term'
            expression = '-L_NL*${mu} * (gr0^3 - gr0 + 2*gr0*(gamma_01*gr1^2 + gamma_02*gr2^2 + gamma_03*gr3^2 + gamma_04*gr4^2) - torque_01_gr0 - torque_02_gr0 - torque_03_gr0 - torque_04_gr0) + (L_NL - ${L}) * kappa_laplacian_gr0'
            inputs = 'gr0 gr1 gr2 gr3 gr4 gamma_01 gamma_02 gamma_03 gamma_04 torque_01_gr0 torque_02_gr0 torque_03_gr0 torque_04_gr0 L_NL kappa_laplacian_gr0'
        []
        [gr1_bulk_term]
            type = ParsedCompute
            buffer = 'gr1_bulk_term'
            expression = '-L_NL*${mu} * (gr1^3 - gr1 + 2*gr1*(gamma_01*gr0^2 + gamma_12*gr2^2 + gamma_13*gr3^2 + gamma_14*gr4^2) - torque_01_gr1 - torque_12_gr1 - torque_13_gr1 - torque_14_gr1) + (L_NL - ${L}) * kappa_laplacian_gr1'
            inputs = 'gr0 gr1 gr2 gr3 gr4 gamma_01 gamma_12 gamma_13 gamma_14 torque_01_gr1 torque_12_gr1 torque_13_gr1 torque_14_gr1 L_NL kappa_laplacian_gr1'
        []
        [gr2_bulk_term]
            type = ParsedCompute
            buffer = 'gr2_bulk_term'
            expression = '-L_NL*${mu} * (gr2^3 - gr2 + 2*gr2*(gamma_02*gr0^2 + gamma_12*gr1^2 + gamma_23*gr3^2 + gamma_24*gr4^2) - torque_02_gr2 - torque_12_gr2 - torque_23_gr2 - torque_24_gr2) + (L_NL - ${L}) * kappa_laplacian_gr2'
            inputs = 'gr0 gr1 gr2 gr3 gr4 gamma_02 gamma_12 gamma_23 gamma_24 torque_02_gr2 torque_12_gr2 torque_23_gr2 torque_24_gr2 L_NL kappa_laplacian_gr2'
        []
        [gr3_bulk_term]
            type = ParsedCompute
            buffer = 'gr3_bulk_term'
            expression = '-L_NL*${mu} * (gr3^3 - gr3 + 2*gr3*(gamma_03*gr0^2 + gamma_13*gr1^2 + gamma_23*gr2^2 + gamma_34*gr4^2) - torque_03_gr3 - torque_13_gr3 - torque_23_gr3 - torque_34_gr3) + (L_NL - ${L}) * kappa_laplacian_gr3'
            inputs = 'gr0 gr1 gr2 gr3 gr4 gamma_03 gamma_13 gamma_23 gamma_34 torque_03_gr3 torque_13_gr3 torque_23_gr3 torque_34_gr3 L_NL kappa_laplacian_gr3'
        []
        [gr4_bulk_term]
            type = ParsedCompute
            buffer = 'gr4_bulk_term'
            expression = '-L_NL*${mu} * (gr4^3 - gr4 + 2*gr4*(gamma_04*gr0^2 + gamma_14*gr1^2 + gamma_24*gr2^2 + gamma_34*gr3^2) - torque_04_gr4 - torque_14_gr4 - torque_24_gr4 - torque_34_gr4) + (L_NL - ${L}) * kappa_laplacian_gr4'
            inputs = 'gr0 gr1 gr2 gr3 gr4 gamma_04 gamma_14 gamma_24 gamma_34 torque_04_gr4 torque_14_gr4 torque_24_gr4 torque_34_gr4 L_NL kappa_laplacian_gr4'
        []

        [NL_gr0]
            type = ForwardFFT
            buffer = NL_gr0
            input = 'gr0_bulk_term'
        []
        [NL_gr1]
            type = ForwardFFT
            buffer = NL_gr1
            input = 'gr1_bulk_term'
        []
        [NL_gr2]
            type = ForwardFFT
            buffer = NL_gr2
            input = 'gr2_bulk_term'
        []
        [NL_gr3]
            type = ForwardFFT
            buffer = NL_gr3
            input = 'gr3_bulk_term'
        []
        [NL_gr4]
            type = ForwardFFT
            buffer = NL_gr4
            input = 'gr4_bulk_term'
        []

        [NL_gr0_smooth]
            type = ParsedCompute
            buffer = NL_gr0_smooth
            expression = 'smooth*NL_gr0'
            inputs = 'smooth NL_gr0'
        []
        [NL_gr1_smooth]
            type = ParsedCompute
            buffer = NL_gr1_smooth
            expression = 'smooth*NL_gr1'
            inputs = 'smooth NL_gr1'
        []
        [NL_gr2_smooth]
            type = ParsedCompute
            buffer = NL_gr2_smooth
            expression = 'smooth*NL_gr2'
            inputs = 'smooth NL_gr2'
        []
        [NL_gr3_smooth]
            type = ParsedCompute
            buffer = NL_gr3_smooth
            expression = 'smooth*NL_gr3'
            inputs = 'smooth NL_gr3'
        []
        [NL_gr4_smooth]
            type = ParsedCompute
            buffer = NL_gr4_smooth
            expression = 'smooth*NL_gr4'
            inputs = 'smooth NL_gr4'
        []

        # Properties for output
        [grain_ID]
            type = ParsedCompute
            buffer = 'grain_ID'
            expression = 'if(gr0>=max(gr1,max(gr2,max(gr3,gr4))), 0, if(gr1>=max(gr2,max(gr3,gr4)), 1, if(gr2>=max(gr3,gr4), 2, if(gr3>=gr4, 3, 4))))'
            inputs = 'gr0 gr1 gr2 gr3 gr4'
        []
        [gradient_energy]
            type = GradientEnergyDensity
            buffer = gradient_energy
            gradient_buffers = 'grad_gr0 grad_gr1 grad_gr2 grad_gr3 grad_gr4'
            kappa = ${kappa}
        []
        [total_energy]
            type = ParsedCompute
            buffer = 'total_energy'
            expression = '${mu} * (${f0}) + gradient_energy'
            inputs = 'gr0 gr1 gr2 gr3 gr4 gamma_01 gamma_02 gamma_03 gamma_04 gamma_12 gamma_13 gamma_14 gamma_23 gamma_24 gamma_34 gradient_energy'
        []
    []
[]

[Postprocessors]
    [total_gb_energy]
        type = TensorIntegralPostprocessor
        buffer = total_energy
    []
[]

[TensorSolver]
    type = AdamsBashforthMoulton
    buffer = 'gr0 gr1 gr2 gr3 gr4'
    reciprocal_buffer = 'gr0_bar gr1_bar gr2_bar gr3_bar gr4_bar'
    linear_reciprocal = 'L_kappa_laplacian L_kappa_laplacian L_kappa_laplacian L_kappa_laplacian L_kappa_laplacian'
    nonlinear_reciprocal = 'NL_gr0_smooth NL_gr1_smooth NL_gr2_smooth NL_gr3_smooth NL_gr4_smooth'
    substeps = 1000
    predictor_order = 1
    corrector_order = 1
    corrector_steps = 1
[]

[TensorOutputs]
    [xdmf]
        type = XDMFTensorOutput
        buffer = 'gr0 gr1 gr2 gr3 gr4 bnds sigma grain_ID'
        output_mode = 'NODE NODE NODE NODE NODE NODE NODE CELL'
        enable_hdf5 = true
        transpose = true
    []
[]

[Executioner]
    type = Transient
    dt = 5
    num_steps = 50
[]

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

(examples/anisotropic_grain_growth/spherical_grain_Ni.i)

interface_width = 1.6
r0 = 7.5

# Ni
gbe_max = '${units 0.8797226 J/m^2}'

L = '${fparse 1.0 * 1.6 / ${interface_width}}'

g_gamma0 = '${fparse sqrt(2) / 3}' # g(gamma=1.5)
f0_gamma0 = ${fparse 1/8}
kappa = '${fparse gbe_max * interface_width * sqrt(f0_gamma0) / g_gamma0}'
mu = '${fparse gbe_max / (g_gamma0 * interface_width * sqrt(f0_gamma0))}'

g_to_gamma_func = '(1.0/(-3.0944 * g^8 -1.8169*g^6 + 10.323 * g^4 - 8.1819*g^2 + 2.0033))'
dgamma_dsigma = '(g * (24.7552*g^6 + 10.9014*g^4 - 41.292*g^2 + 16.3638) / (-3.0944*g^8 - 1.8169*g^6 + 10.323*g^4 - 8.1819*g^2 + 2.0033)^2 / sqrt(${kappa} * ${mu}))'

f0 = '((gr0^4/4 - gr0^2/2) + (gr1^4/4 - gr1^2/2) + gamma_01*gr0^2*gr1^2 + 0.25)'

[Domain]
  dim = 3
  nx = 100
  ny = 100
  nz = 100
  xmax = 20
  ymax = 20
  zmax = 20
  mesh_mode = DUMMY
[]

[TensorComputes]
  [Initialize]
    [gr0]
      type = ParsedCompute
      buffer = 'gr0'
      extra_symbols = 'true'
      expression = 'radius:=sqrt((x-${Domain/xmax}/2)^2+(y-${Domain/ymax}/2)^2+(z-${Domain/zmax}/2)^2);0.5 - 0.5*tanh(2*(radius-${r0})/${interface_width})'
    []
    [gr1]
      type = ParsedCompute
      buffer = 'gr1'
      extra_symbols = 'true'
      expression = 'radius:=sqrt((x-${Domain/xmax}/2)^2+(y-${Domain/ymax}/2)^2+(z-${Domain/zmax}/2)^2);0.5 + 0.5*tanh(2*(radius-${r0})/${interface_width})'
    []
    [L_kappa_laplacian]
      type = ReciprocalLaplacianFactor
      buffer = L_kappa_laplacian
      factor = ${fparse ${L} * ${kappa}}
    []
    [kappa_laplacian]
      type = ReciprocalLaplacianFactor
      buffer = kappa_laplacian
      factor = '${kappa}'
    []
    [smooth]
      type = DeAliasingTensor
      method = HOULI
      buffer = smooth
    []
  []
  [Solve]
    ## FFT
    [gr0_bar]
      type = ForwardFFT
      buffer = gr0_bar
      input = gr0
    []
    [gr1_bar]
      type = ForwardFFT
      buffer = gr1_bar
      input = gr1
    []

    [gr0_bar_smooth]
      type = ParsedCompute
      buffer = gr0_bar_smooth
      expression = 'gr0_bar * smooth'
      inputs = 'gr0_bar smooth'
    []
    [gr1_bar_smooth]
      type = ParsedCompute
      buffer = gr1_bar_smooth
      expression = 'gr1_bar * smooth'
      inputs = 'gr1_bar smooth'
    []

    # Gradients
    [grad_gr0]
      type = GradientVector
      buffer = grad_gr0
      input = gr0_bar_smooth
      input_is_reciprocal = true
    []
    [grad_gr1]
      type = GradientVector
      buffer = grad_gr1
      input = gr1_bar_smooth
      input_is_reciprocal = true
    []

    # Get interface energy
    [sigma_gr0_gr1]
      type = PairwiseAnisotropicGBEnergy
      buffer = 'sigma_gr0_gr1'
      dsigma_dgrad_grain1 = 'dsigma_gr0_gr1_dgrad_gr0'
      dsigma_dgrad_grain2 = 'dsigma_gr0_gr1_dgrad_gr1'
      grad_grain1_buffer = 'grad_gr0'
      grad_grain2_buffer = 'grad_gr1'
      interface_width = '${interface_width}'
      libtorch_model_file = 'marlin:anisotropic_gb_energy/Ni_R1_R3.pt'
      chunk_size = 1e5
    []
    [gamma_01]
      type = ParsedCompute
      buffer = 'gamma_01'
      expression = 'g:=sigma_gr0_gr1/sqrt(${kappa}*${mu});${g_to_gamma_func}'
      inputs = 'sigma_gr0_gr1'
    []
    [coeff_01]
      type = ParsedCompute
      buffer = 'coeff_01'
      expression = 'g:=sigma_gr0_gr1/sqrt(${kappa}*${mu});${dgamma_dsigma} * gr0^2 *gr1^2'
      inputs = 'sigma_gr0_gr1 gr0 gr1'
    []
    [torque_01_gr0]
      type = AnisotropyTorque
      buffer = torque_01_gr0
      dmu_dn = dsigma_gr0_gr1_dgrad_gr0
      g = 'coeff_01'
    []
    [torque_01_gr1]
      type = AnisotropyTorque
      buffer = torque_01_gr1
      dmu_dn = dsigma_gr0_gr1_dgrad_gr1
      g = 'coeff_01'
    []
    [kappa_laplacian_gr0_hat]
      type = ParsedCompute
      buffer = kappa_laplacian_gr0_hat
      inputs = 'gr0_bar kappa_laplacian'
      expression = 'gr0_bar * kappa_laplacian'
    []
    [kappa_laplacian_gr1_hat]
      type = ParsedCompute
      buffer = kappa_laplacian_gr1_hat
      inputs = 'gr1_bar kappa_laplacian'
      expression = 'gr1_bar * kappa_laplacian'
    []

    [kappa_laplacian_gr0]
      type = InverseFFT
      buffer = kappa_laplacian_gr0
      input = kappa_laplacian_gr0_hat
    []
    [kappa_laplacian_gr1]
      type = InverseFFT
      buffer = kappa_laplacian_gr1
      input = kappa_laplacian_gr1_hat
    []

    [gr0_bulk_term]
      type = ParsedCompute
      buffer = 'gr0_bulk_term'
      expression = 'L_NL:=${L}*(sigma_gr0_gr1)/${gbe_max};-L_NL*${mu} * (gr0^3 - gr0 + 2*gr0*(gamma_01*gr1^2) - torque_01_gr0) + (L_NL - ${L})*kappa_laplacian_gr0'
      inputs = 'gr0 gr1 gamma_01 torque_01_gr0 sigma_gr0_gr1 kappa_laplacian_gr0'
    []
    [gr1_bulk_term]
      type = ParsedCompute
      buffer = 'gr1_bulk_term'
      expression = 'L_NL:=${L}*(sigma_gr0_gr1)/${gbe_max};-L_NL*${mu} * (gr1^3 - gr1 + 2*gr1*(gamma_01*gr0^2) - torque_01_gr1) + (L_NL - ${L})*kappa_laplacian_gr1'
      inputs = 'gr0 gr1 gamma_01 torque_01_gr1 sigma_gr0_gr1 kappa_laplacian_gr1'
    []

    ## build non-linear terms and smooth them
    [NL_gr0]
      type = ForwardFFT
      buffer = NL_gr0
      input = 'gr0_bulk_term'
    []
    [NL_gr1]
      type = ForwardFFT
      buffer = NL_gr1
      input = 'gr1_bulk_term'
    []
    [NL_gr0_smooth]
      type = ParsedCompute
      buffer = NL_gr0_smooth
      expression = 'smooth*NL_gr0'
      inputs = 'smooth NL_gr0'
    []
    [NL_gr1_smooth]
      type = ParsedCompute
      buffer = NL_gr1_smooth
      expression = 'smooth*NL_gr1'
      inputs = 'smooth NL_gr1'
    []

    # Properties for output
    [sigma]
      type = ParsedCompute
      buffer = 'sigma'
      expression = 'gr0^2*gr1^2*sigma_gr0_gr1'
      inputs = 'sigma_gr0_gr1 gr0 gr1'
    []

    [gradient_energy]
      type = GradientEnergyDensity
      buffer = gradient_energy
      gradient_buffers = 'grad_gr0 grad_gr1'
      kappa = ${kappa}
    []

    [total_energy]
      type = ParsedCompute
      buffer = 'total_energy'
      expression = '${mu} * (${f0}) + gradient_energy'
      inputs = 'gr0 gr1 gamma_01 gradient_energy'
    []
  []
[]

[Postprocessors]
  [total_gb_energy]
    type = TensorIntegralPostprocessor
    buffer = total_energy
  []
[]

[TensorSolver]
  type = AdamsBashforthMoulton
  buffer = 'gr0 gr1'
  reciprocal_buffer = 'gr0_bar gr1_bar'
  linear_reciprocal = 'L_kappa_laplacian L_kappa_laplacian'
  nonlinear_reciprocal = 'NL_gr0_smooth NL_gr1_smooth'
  substeps = 100
  predictor_order = 2
  corrector_order = 1
  corrector_steps = 2
[]

[TensorOutputs]
  [xdmf]
    type = XDMFTensorOutput
    buffer = 'gr0 gr1 sigma_gr0_gr1'
    output_mode = 'NODE NODE NODE'
    enable_hdf5 = true
    transpose = false
  []
[]

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

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

(test/tests/anisotropic_grain_growth/sigma3_circular_grain.i)

interface_width = 2.0
r0 = 3

gbe_max = '${units 0.8797226 J/m^2}'

L = '${fparse 1.0 * 1.6 / ${interface_width}}'

g_gamma0 = '${fparse sqrt(2) / 3}' # g(gamma=1.5)
f0_gamma0 = ${fparse 1/8}
kappa = '${fparse gbe_max * interface_width * sqrt(f0_gamma0) / g_gamma0}'
mu = '${fparse gbe_max / (g_gamma0 * interface_width * sqrt(f0_gamma0))}'

g_to_gamma_func = '(1.0/(-3.0944 * g^8 -1.8169*g^6 + 10.323 * g^4 - 8.1819*g^2 + 2.0033))'
dgamma_dsigma = '(g * (24.7552*g^6 + 10.9014*g^4 - 41.292*g^2 + 16.3638) / (-3.0944*g^8 - 1.8169*g^6 + 10.323*g^4 - 8.1819*g^2 + 2.0033)^2 / sqrt(${kappa} * ${mu}))'

f0 = '((gr0^4/4 - gr0^2/2) + (gr1^4/4 - gr1^2/2) + gamma_01*gr0^2*gr1^2 + 0.25)'

[Domain]
  dim = 2
  nx = 32
  ny = 32
  xmax = 8
  ymax = 8
  mesh_mode = DUMMY
  floating_precision = SINGLE
[]

[TensorComputes]
  [Initialize]
    [gr0]
      type = ParsedCompute
      buffer = 'gr0'
      extra_symbols = 'true'
      expression = 'radius:=sqrt((x-${Domain/xmax}/2)^2+(y-${Domain/ymax}/2)^2);0.5 - 0.5*tanh(2*(radius-${r0})/${interface_width})'
    []
    [gr1]
      type = ParsedCompute
      buffer = 'gr1'
      extra_symbols = 'true'
      expression = 'radius:=sqrt((x-${Domain/xmax}/2)^2+(y-${Domain/ymax}/2)^2);0.5 + 0.5*tanh(2*(radius-${r0})/${interface_width})'
    []
    [L_kappa_laplacian]
      type = ReciprocalLaplacianFactor
      buffer = L_kappa_laplacian
      factor = ${fparse ${L} * ${kappa}}
    []
    [kappa_laplacian]
      type = ReciprocalLaplacianFactor
      buffer = kappa_laplacian
      factor = '${kappa}'
    []
    [smooth]
      type = DeAliasingTensor
      method = HOULI
      buffer = smooth
    []
  []
  [Solve]
    ## FFT
    [gr0_bar]
      type = ForwardFFT
      buffer = gr0_bar
      input = gr0
    []
    [gr1_bar]
      type = ForwardFFT
      buffer = gr1_bar
      input = gr1
    []

    # Gradients
    [grad_gr0]
      type = GradientVector
      buffer = grad_gr0
      input = gr0_bar
      input_is_reciprocal = true
    []
    [grad_gr1]
      type = GradientVector
      buffer = grad_gr1
      input = gr1_bar
      input_is_reciprocal = true
    []

    # Get interface energy
    [sigma_gr0_gr1]
      type = PairwiseAnisotropicGBEnergy
      buffer = 'sigma_gr0_gr1'
      dsigma_dgrad_grain1 = 'dsigma_gr0_gr1_dgrad_gr0'
      dsigma_dgrad_grain2 = 'dsigma_gr0_gr1_dgrad_gr1'
      grad_grain1_buffer = 'grad_gr0'
      grad_grain2_buffer = 'grad_gr1'
      interface_width = '${interface_width}'
      libtorch_model_file = 'marlin:anisotropic_gb_energy/gb_energy_hull_2d.pt'
      chunk_size = 256
    []
    [gamma_01]
      type = ParsedCompute
      buffer = 'gamma_01'
      expression = 'g:=sigma_gr0_gr1/sqrt(${kappa}*${mu});${g_to_gamma_func}'
      inputs = 'sigma_gr0_gr1'
    []
    [coeff_01]
      type = ParsedCompute
      buffer = 'coeff_01'
      expression = 'g:=sigma_gr0_gr1/sqrt(${kappa}*${mu});${dgamma_dsigma} * gr0^2 *gr1^2'
      inputs = 'sigma_gr0_gr1 gr0 gr1'
    []
    [torque_01_gr0]
      type = AnisotropyTorque
      buffer = torque_01_gr0
      dmu_dn = dsigma_gr0_gr1_dgrad_gr0
      g = 'coeff_01'
    []
    [torque_01_gr1]
      type = AnisotropyTorque
      buffer = torque_01_gr1
      dmu_dn = dsigma_gr0_gr1_dgrad_gr1
      g = 'coeff_01'
    []
    [kappa_laplacian_gr0_hat]
      type = ParsedCompute
      buffer = kappa_laplacian_gr0_hat
      inputs = 'gr0_bar kappa_laplacian'
      expression = 'gr0_bar * kappa_laplacian'
    []
    [kappa_laplacian_gr1_hat]
      type = ParsedCompute
      buffer = kappa_laplacian_gr1_hat
      inputs = 'gr1_bar kappa_laplacian'
      expression = 'gr1_bar * kappa_laplacian'
    []

    [kappa_laplacian_gr0]
      type = InverseFFT
      buffer = kappa_laplacian_gr0
      input = kappa_laplacian_gr0_hat
    []
    [kappa_laplacian_gr1]
      type = InverseFFT
      buffer = kappa_laplacian_gr1
      input = kappa_laplacian_gr1_hat
    []

    [gr0_bulk_term]
      type = ParsedCompute
      buffer = 'gr0_bulk_term'
      expression = 'L_NL:=${L}*(sigma_gr0_gr1)/${gbe_max};-L_NL*${mu} * (gr0^3 - gr0 + 2*gr0*(gamma_01*gr1^2) - torque_01_gr0) + (L_NL - ${L})*kappa_laplacian_gr0'
      inputs = 'gr0 gr1 gamma_01 torque_01_gr0 sigma_gr0_gr1 kappa_laplacian_gr0'
    []
    [gr1_bulk_term]
      type = ParsedCompute
      buffer = 'gr1_bulk_term'
      expression = 'L_NL:=${L}*(sigma_gr0_gr1)/${gbe_max};-L_NL*${mu} * (gr1^3 - gr1 + 2*gr1*(gamma_01*gr0^2) - torque_01_gr1) + (L_NL - ${L})*kappa_laplacian_gr1'
      inputs = 'gr0 gr1 gamma_01 torque_01_gr1 sigma_gr0_gr1 kappa_laplacian_gr1'
    []

    ## build non-linear terms and smooth them
    [NL_gr0]
      type = ForwardFFT
      buffer = NL_gr0
      input = 'gr0_bulk_term'
    []
    [NL_gr1]
      type = ForwardFFT
      buffer = NL_gr1
      input = 'gr1_bulk_term'
    []
    [NL_gr0_smooth]
      type = ParsedCompute
      buffer = NL_gr0_smooth
      expression = 'smooth*NL_gr0'
      inputs = 'smooth NL_gr0'
    []
    [NL_gr1_smooth]
      type = ParsedCompute
      buffer = NL_gr1_smooth
      expression = 'smooth*NL_gr1'
      inputs = 'smooth NL_gr1'
    []

    # Properties for output
    [sigma]
      type = ParsedCompute
      buffer = 'sigma'
      expression = 'gr0^2*gr1^2*sigma_gr0_gr1'
      inputs = 'sigma_gr0_gr1 gr0 gr1'
    []

    [gradient_energy]
      type = GradientEnergyDensity
      buffer = gradient_energy
      gradient_buffers = 'grad_gr0 grad_gr1'
      kappa = ${kappa}
    []

    [total_energy]
      type = ParsedCompute
      buffer = 'total_energy'
      expression = '${mu} * (${f0}) + gradient_energy'
      inputs = 'gr0 gr1 gamma_01 gradient_energy'
    []
  []
[]

[Postprocessors]
  [total_gb_energy]
    type = TensorIntegralPostprocessor
    buffer = total_energy
  []
[]

[TensorSolver]
  type = AdamsBashforthMoulton
  buffer = 'gr0 gr1'
  reciprocal_buffer = 'gr0_bar gr1_bar'
  linear_reciprocal = 'L_kappa_laplacian L_kappa_laplacian'
  nonlinear_reciprocal = 'NL_gr0_smooth NL_gr1_smooth'
  substeps = 1e2
  predictor_order = 1
  corrector_order = 1
  corrector_steps = 1
[]

[TensorOutputs]
  active = ''
  [xdmf]
    type = XDMFTensorOutput
    buffer = 'gr0 gr1 sigma_gr0_gr1 sigma total_energy'
    output_mode = 'NODE NODE NODE NODE NODE'
    enable_hdf5 = true
    transpose = false
  []
[]

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

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

(test/tests/tensor_compute/pfc_fcc.i)

# FCC PFC test designed to clearly show atomic-scale density fluctuations
# Using larger epsilon for stronger density modulations
# Based on PhysRevE.81.061601 two-mode FCC PFC model

# Grid parameters - fine resolution to resolve atomic peaks
N = 1  # Grid points per dimension

# PFC model parameters - LARGER epsilon for clearer atomic structure
psi_mean = -0.2         # Mean density
epsilon = 0.15          # Larger undercooling for stronger modulations
R1_param = 0.0          # Two-mode coupling (0 for max FCC stability)
Q1_param = ${fparse sqrt(4.0/3.0)}  # FCC wave number ratio

# Domain size: want several unit cells but fine enough resolution
# The (111) wavelength is 2*pi/q0 = 2*pi*sqrt(3/4) ≈ 5.44 (dimensionless)
# Let's fit about 'N' wavelengths in each direction
Lx = ${fparse N * 2 * pi / ${Q1_param}}
Ly = ${Lx}

[Domain]
    dim = 2
    nx = ${fparse ${N} * 8 }
    ny = ${fparse ${N} * 8 }
    xmax = ${Lx}
    ymax = ${Ly}
[]

[TensorComputes]
    [Initialize]
        # Random perturbations to seed crystal growth
        [psi]
            type = RandomTensor
            buffer = psi
            max = ${fparse ${psi_mean} + 0.01}
            min = ${fparse ${psi_mean} - 0.01}
            seed = 12345
        []

        # Linear operator for FCC
        [linear]
            type = FCCPFCLinear
            buffer = 'linear'
            eps = ${epsilon}
            R1 = ${R1_param}
            Q1 = ${Q1_param}
            mobility = 1.0
        []

        # Dealiasing for cubic nonlinearity
        [smooth_operator]
            type = DeAliasingTensor
            buffer = smooth_operator
            method = HOULI
        []
    []

    [Solve]
        # Nonlinear term
        [nl_div_psi_cubed]
            type = FCCPFCNonlinear
            buffer = NL
            psi = psi
            dealiasing = smooth_operator
            mobility = 1.0
        []

        # FFT for spectral solver
        [psi_hat]
            type = ForwardFFT
            buffer = psi_hat
            input = psi
        []
    []
[]

[TensorSolver]
    type = AdamsBashforthMoulton
    buffer = 'psi'
    linear_reciprocal = 'linear'
    nonlinear_reciprocal = 'NL'
    reciprocal_buffer = 'psi_hat'
    corrector_order = 1
    corrector_steps = 3
    predictor_order = 1
    substeps = 1000
[]

[Postprocessors]
    [max_psi]
        type = TensorExtremeValuePostprocessor
        buffer = psi
        value_type = MAX
        execute_on = 'initial timestep_end'
    []
    [min_psi]
        type = TensorExtremeValuePostprocessor
        buffer = psi
        value_type = MIN
        execute_on = 'initial timestep_end'
    []
    [mean_psi]
        type = TensorAveragePostprocessor
        buffer = psi
        execute_on = 'initial timestep_end'
    []
[]

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

[Executioner]
    type = Transient
    num_steps = 100
    dt = 0.2  # Larger timestep for faster evolution
[]

[Outputs]
    perf_graph = true
    csv = true
[]

(examples/phase_field_crystal/pfc_fcc_atomic_structure_3d.i)

# 3D FCC PFC test showing atomic-scale density fluctuations
# Based on PhysRevE.81.061601 two-mode FCC PFC model

# Grid parameters - 3D resolution
N = 64  # Grid points per dimension (64^3 = 262k points)

# PFC model parameters - LARGER epsilon for clearer atomic structure
psi_mean = -0.2         # Mean density
epsilon = 0.15          # Larger undercooling for stronger modulations
R1_param = 0.0          # Two-mode coupling (0 for max FCC stability)
Q1_param = ${fparse sqrt(4.0/3.0)}  # FCC wave number ratio

# Domain size: want several unit cells
# The (111) wavelength is 2*pi/q0 = 2*pi*sqrt(3/4) ≈ 5.44 (dimensionless)
# Fit 8 wavelengths in each direction
Lx = ${fparse 8 * 2 * ${pi} / ${Q1_param}}
Ly = ${Lx}
Lz = ${Lx}

[Domain]
    dim = 3
    nx = ${N}
    ny = ${N}
    nz = ${N}
    xmax = ${Lx}
    ymax = ${Ly}
    zmax = ${Lz}
[]

[TensorComputes]
    [Initialize]
        # Random perturbations to seed crystal growth
        [psi]
            type = RandomTensor
            buffer = psi
            max = ${fparse ${psi_mean} + 0.01}
            min = ${fparse ${psi_mean} - 0.01}
            seed = 12345
        []

        # Linear operator for FCC
        [linear]
            type = FCCPFCLinear
            buffer = 'linear'
            eps = ${epsilon}
            R1 = ${R1_param}
            Q1 = ${Q1_param}
            mobility = 1.0
        []

        # Dealiasing for cubic nonlinearity
        [smooth_operator]
            type = DeAliasingTensor
            buffer = smooth_operator
            method = HOULI
        []
    []

    [Solve]
        # Nonlinear term
        [nl_div_psi_cubed]
            type = FCCPFCNonlinear
            buffer = NL
            psi = psi
            dealiasing = smooth_operator
            mobility = 1.0
        []

        # FFT for spectral solver
        [psi_hat]
            type = ForwardFFT
            buffer = psi_hat
            input = psi
        []
    []
[]

[TensorSolver]
    type = AdamsBashforthMoulton
    buffer = 'psi'
    linear_reciprocal = 'linear'
    nonlinear_reciprocal = 'NL'
    reciprocal_buffer = 'psi_hat'
    corrector_order = 1
    corrector_steps = 3
    predictor_order = 1
    substeps = 10000
[]

[Postprocessors]
    [max_psi]
        type = TensorExtremeValuePostprocessor
        buffer = psi
        value_type = MAX
        execute_on = 'initial timestep_end'
    []
    [min_psi]
        type = TensorExtremeValuePostprocessor
        buffer = psi
        value_type = MIN
        execute_on = 'initial timestep_end'
    []
    [mean_psi]
        type = TensorAveragePostprocessor
        buffer = psi
        execute_on = 'initial timestep_end'
    []
[]

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

[Executioner]
    type = Transient
    num_steps = 200  # Reduced for 3D (much more expensive)
    dt = 5
[]

[Outputs]
    perf_graph = true
    csv = true
[]

Example Input File Syntax
Input Parameters
Input Files