FiniteDifferenceBiharmonic

Computes the biharmonic Laplacian(Laplacian(u)) using finite differences.

Overview

Computes factor times the biharmonic $var element = document.getElementById("moose-equation-86d0e69f-3b20-406e-a236-1a140a4cfbd6");katex.render("\\nabla^2\\nabla^2 u", element, {displayMode:false,throwOnError:false});$ of the input buffer in real space. The stencil is the discrete Laplacian of FiniteDifferenceLaplacian composed with itself ("stencil_width" selects the 3-point or 5-point Laplacian), which in 2D with the 3-point Laplacian is the 13-point stencil. It reaches twice as far as the Laplacian, so two (3-point) or four (5-point) ghost layers are exchanged before it is applied.

The operator also exposes its stencil for implicit solvers. On a periodic grid the square of the symmetric Laplacian is positive semidefinite with the constants as its nullspace, so with factor > 0 the operator can be solved implicitly with ImplicitFDSolver, for example as the $var element = document.getElementById("moose-equation-810e0705-0ee3-4ee6-96ae-bece1c4625c9");katex.render("M\\kappa\\nabla^4 c", element, {displayMode:false,throwOnError:false});$ term of the Cahn-Hilliard equation.

Example Input File Syntax

[TensorComputes<<<{"href": "../../syntax/TensorComputes/index.html"}>>>]
  [Solve<<<{"href": "../../syntax/TensorComputes/Solve/index.html"}>>>]
    [bih]
      type = FiniteDifferenceBiharmonic<<<{"description": "Computes the biharmonic Laplacian(Laplacian(u)) using finite differences.", "href": "FiniteDifferenceBiharmonic.html"}>>>
      buffer<<<{"description": "The buffer this compute is writing to"}>>> = bih_u
      input<<<{"description": "The input buffer"}>>> = u
      factor<<<{"description": "Pre-factor for the biharmonic"}>>> = ${Mkappa}
      stencil_width<<<{"description": "Stencil width of the Laplacian that is applied twice (3 => 3-point, 5 => 5-point)"}>>> = ${stencil}
    []
  []
[]

Input Parameters

bufferThe buffer this compute is writing to
C++ Type:std::string
Controllable:No
Description:The buffer this compute is writing to
inputThe input buffer
C++ Type:std::string
Controllable:No
Description:The input buffer

Required Parameters

factor1Pre-factor for the biharmonic
Default:1
C++ Type:Real
Unit:(no unit assumed)
Controllable:No
Description:Pre-factor for the biharmonic
stencil_width3Stencil width of the Laplacian that is applied twice (3 => 3-point, 5 => 5-point)
Default:3
C++ Type:MooseEnum
Options:3, 5
Controllable:No
Description:Stencil width of the Laplacian that is applied twice (3 => 3-point, 5 => 5-point)

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/real_space/deep_ghost_periodic.i)
(test/tests/implicit_fd/biharmonic_mode3d.i)
(test/tests/implicit_fd/biharmonic_mode.i)
(benchmarks/01_spinodal_decomposition/1a_implicit_fd.i)
(test/tests/real_space/deep_ghost_faces.i)

(test/tests/implicit_fd/biharmonic_mode.i)

# Backward Euler for u_t = -M kappa Laplacian(Laplacian(u)) of a single periodic
# Fourier mode. The mode is an eigenvector of the discrete Laplacian with
# eigenvalue -lambda_h, so of the discrete biharmonic with eigenvalue lambda_h^2,
# and every implicit step must multiply it by exactly
# g = 1 / (1 + dt M kappa lambda_h^2), up to the linear solver tolerance.

Mkappa = 1e-3
dt = 0.01
n = 32
stencil = 3

# Symbol of -Laplacian per axis for the wavenumber 2 pi, h = 1/n (see be_mode.i)
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 2 * 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]
  []
  [bih_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
    []
    [bih]
      type = FiniteDifferenceBiharmonic
      buffer = bih_u
      input = u
      factor = ${Mkappa}
      stencil_width = ${stencil}
    []
  []
  [Postprocess]
    [gold_compute]
      type = ParsedCompute
      buffer = gold
      expression = 'sin(2*pi*x)*sin(2*pi*y) * g^(t / dt)'
      constant_names = 'g dt'
      constant_expressions = '${fparse 1 / (1 + dt * Mkappa * lambda^2)} ${dt}'
      extra_symbols = true
    []
    [diff_compute]
      type = ParsedCompute
      buffer = diff
      expression = 'abs(u - gold)'
      inputs = 'u gold'
    []
  []
[]

[TensorSolver]
  type = ImplicitFDSolver
  buffer = u
  operators = 'dudt bih'
  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
  dt = ${dt}
[]

[Outputs]
  csv = true
[]

(test/tests/real_space/deep_ghost_periodic.i)

# The reference of deep_ghost_faces.i: the field mirrored onto the periodic domain 0 < x < 32,
# where the 5-point biharmonic reads only periodic ghosts. The postprocessors take 0 < x < 16.

u_ic = 'exp(cos(pi*x/16))*sin(2*pi*y/8)'

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

[TensorBuffers]
  [u]
  []
  [bih]
  []
  [bih_strip]
  []
  [bih_sq]
  []
[]

[TensorComputes]
  [Initialize]
    [u_ic]
      type = ParsedCompute
      buffer = u
      extra_symbols = true
      expression = ${u_ic}
    []
  []
  [Postprocess]
    [bih]
      type = FiniteDifferenceBiharmonic
      buffer = bih
      input = u
      stencil_width = 5
    []
    # 0 outside the strip, which is between the extremes of both test fields there
    [bih_strip]
      type = ParsedCompute
      buffer = bih_strip
      extra_symbols = true
      expression = 'if(x < 16, bih, 0)'
      inputs = bih
    []
    [bih_sq]
      type = ParsedCompute
      buffer = bih_sq
      expression = 'bih_strip^2'
      inputs = bih_strip
    []
  []
[]

[Postprocessors]
  [integral_bih_sq]
    type = TensorIntegralPostprocessor
    buffer = bih_sq
  []
  [max_bih]
    type = TensorExtremeValuePostprocessor
    buffer = bih_strip
    value_type = MAX
  []
  [min_bih]
    type = TensorExtremeValuePostprocessor
    buffer = bih_strip
    value_type = MIN
  []
[]

[Problem]
  type = TensorProblem
[]

[Executioner]
  type = Transient
  num_steps = 1
[]

[Outputs]
  csv = true
[]

(test/tests/implicit_fd/biharmonic_mode3d.i)

# Backward Euler for u_t = -M kappa Laplacian(Laplacian(u)) of a single periodic
# Fourier mode. The mode is an eigenvector of the discrete Laplacian with
# eigenvalue -lambda_h, so of the discrete biharmonic with eigenvalue lambda_h^2,
# and every implicit step must multiply it by exactly
# g = 1 / (1 + dt M kappa lambda_h^2), up to the linear solver tolerance.
# 3D version of biharmonic_mode.i.

Mkappa = 1e-3
dt = 0.01
n = 16
stencil = 3

# Symbol of -Laplacian per axis for the wavenumber 2 pi, h = 1/n (see be_mode.i)
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 3 * if(stencil = 3, lambda3, lambda5)}'

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

[TensorBuffers]
  [u]
  []
  [dudt_u]
  []
  [bih_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
    []
    [bih]
      type = FiniteDifferenceBiharmonic
      buffer = bih_u
      input = u
      factor = ${Mkappa}
      stencil_width = ${stencil}
    []
  []
  [Postprocess]
    [gold_compute]
      type = ParsedCompute
      buffer = gold
      expression = 'sin(2*pi*x)*sin(2*pi*y)*sin(2*pi*z) * g^(t / dt)'
      constant_names = 'g dt'
      constant_expressions = '${fparse 1 / (1 + dt * Mkappa * lambda^2)} ${dt}'
      extra_symbols = true
    []
    [diff_compute]
      type = ParsedCompute
      buffer = diff
      expression = 'abs(u - gold)'
      inputs = 'u gold'
    []
  []
[]

[TensorSolver]
  type = ImplicitFDSolver
  buffer = u
  operators = 'dudt bih'
  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
  dt = ${dt}
[]

[Outputs]
  csv = true
[]

(test/tests/implicit_fd/biharmonic_mode.i)

# Backward Euler for u_t = -M kappa Laplacian(Laplacian(u)) of a single periodic
# Fourier mode. The mode is an eigenvector of the discrete Laplacian with
# eigenvalue -lambda_h, so of the discrete biharmonic with eigenvalue lambda_h^2,
# and every implicit step must multiply it by exactly
# g = 1 / (1 + dt M kappa lambda_h^2), up to the linear solver tolerance.

Mkappa = 1e-3
dt = 0.01
n = 32
stencil = 3

# Symbol of -Laplacian per axis for the wavenumber 2 pi, h = 1/n (see be_mode.i)
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 2 * 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]
  []
  [bih_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
    []
    [bih]
      type = FiniteDifferenceBiharmonic
      buffer = bih_u
      input = u
      factor = ${Mkappa}
      stencil_width = ${stencil}
    []
  []
  [Postprocess]
    [gold_compute]
      type = ParsedCompute
      buffer = gold
      expression = 'sin(2*pi*x)*sin(2*pi*y) * g^(t / dt)'
      constant_names = 'g dt'
      constant_expressions = '${fparse 1 / (1 + dt * Mkappa * lambda^2)} ${dt}'
      extra_symbols = true
    []
    [diff_compute]
      type = ParsedCompute
      buffer = diff
      expression = 'abs(u - gold)'
      inputs = 'u gold'
    []
  []
[]

[TensorSolver]
  type = ImplicitFDSolver
  buffer = u
  operators = 'dudt bih'
  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
  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/deep_ghost_faces.i)

# A stencil reaching four cells across a face (the 5-point biharmonic) with Neumann faces in x:
# its four ghost layers are the even mirror of the owned cells about the face. The result on
# 0 < x < 16 then equals the same stencil on the periodic domain 0 < x < 32 of the mirrored field
# (deep_ghost_periodic.i), which needs no ghost fill at a face, to round-off.
#
# u = exp(cos(pi x / 16)) sin(2 pi y / 8) is even about x = 0 and x = 16. With Dirichlet faces and
# u_ic = sin(pi x / 16) exp(cos(2 pi y / 8)), odd about both, the ghosts are the odd mirror.

u_ic = 'exp(cos(pi*x/16))*sin(2*pi*y/8)'

[Domain]
  dim = 2
  nx = 16
  ny = 8
  xmax = 16
  ymax = 8
  parallel_mode = REAL_SPACE
  periodic_directions = 'Y'
[]

[TensorBuffers]
  [u]
  []
  [bih]
  []
  [bih_sq]
  []
[]

[TensorComputes]
  [Initialize]
    [u_ic]
      type = ParsedCompute
      buffer = u
      extra_symbols = true
      expression = ${u_ic}
    []
  []
  [BCs]
    [left]
      type = TensorNeumannBC
      variable = u
      boundary = left
    []
    [right]
      type = TensorNeumannBC
      variable = u
      boundary = right
    []
  []
  [Postprocess]
    [bih]
      type = FiniteDifferenceBiharmonic
      buffer = bih
      input = u
      stencil_width = 5
    []
    [bih_sq]
      type = ParsedCompute
      buffer = bih_sq
      expression = 'bih^2'
      inputs = bih
    []
  []
[]

[Postprocessors]
  [integral_bih_sq]
    type = TensorIntegralPostprocessor
    buffer = bih_sq
  []
  [max_bih]
    type = TensorExtremeValuePostprocessor
    buffer = bih
    value_type = MAX
  []
  [min_bih]
    type = TensorExtremeValuePostprocessor
    buffer = bih
    value_type = MIN
  []
[]

[Problem]
  type = TensorProblem
[]

[Executioner]
  type = Transient
  num_steps = 1
[]

[Outputs]
  csv = true
[]

Overview
Example Input File Syntax
Input Parameters
Input Files