- bufferThe buffer this compute is writing to
C++ Type:std::string
Controllable:No
Description:The buffer this compute is writing to
- inputThe buffer being time integrated
C++ Type:std::string
Controllable:No
Description:The buffer being time integrated
TimeDerivativeOperator
Backward Euler time derivative term u/dt for implicit solves, contributing u_old/dt to the right-hand side.
Overview
The backward Euler time derivative term for implicit finite-difference solves. Its action on the solved field is , and it contributes to the right-hand side, with the current substep size. Listed together with the spatial operators of a field, an implicit solver then solves
The operator is positive definite, and its stencil is the single centre entry at every grid spacing.
Input 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/implicit_fd/be_mode.i)
- (test/tests/implicit_fd/biharmonic_mode3d.i)
- (test/tests/implicit_fd/biharmonic_mode.i)
- (test/tests/implicit_fd/bc_mms.i)
- (benchmarks/01_spinodal_decomposition/1a_implicit_fd.i)
- (test/tests/implicit_fd/crank_thin_slab.i)
- (test/tests/implicit_fd/cutback.i)
- (test/tests/implicit_fd/be_mode3d.i)
- (test/tests/implicit_fd/crank_semi_infinite.i)
(test/tests/implicit_fd/be_mode.i)
# Backward Euler diffusion u_t = D Laplacian(u) of a single periodic Fourier
# mode. On the periodic cell-centred grid the mode is an eigenvector of the
# discrete Laplacian with eigenvalue -lambda_h, so every implicit step must
# multiply it by exactly g = 1 / (1 + dt D lambda_h), up to the linear solver
# tolerance.
D = 0.5
# 5 steps of dt, then 5 steps of dt2 (dt2 != dt exercises the operator rebuild)
dt = 0.01
dt2 = 0.01
n = 32
stencil = 3
dim = 2
# Symbol of -Laplacian per axis for the wavenumber 2 pi, h = 1/n:
# 3-point: (2 - 2 cos(2 pi h)) / h^2
# 5-point: (5/2 - 8/3 cos(2 pi h) + 1/6 cos(4 pi h)) / h^2
lambda3 = '${fparse n^2 * (2 - 2 * cos(2 * pi / n))}'
lambda5 = '${fparse n^2 * (5 / 2 - 8 / 3 * cos(2 * pi / n) + 1 / 6 * cos(4 * pi / n))}'
lambda = '${fparse dim * if(stencil = 3, lambda3, lambda5)}'
[Domain]
dim = 2
nx = ${n}
ny = ${n}
xmax = 1
ymax = 1
parallel_mode = REAL_SPACE
periodic_directions = 'X Y'
[]
[TensorBuffers]
[u]
[]
[dudt_u]
[]
[lap_u]
[]
[gold]
[]
[diff]
[]
[]
[TensorComputes]
[Initialize]
[u_ic]
type = ParsedCompute
buffer = u
expression = 'sin(2*pi*x)*sin(2*pi*y)'
extra_symbols = true
[]
[]
[Solve]
[dudt]
type = TimeDerivativeOperator
buffer = dudt_u
input = u
[]
[lap]
type = FiniteDifferenceLaplacian
buffer = lap_u
input = u
factor = '${fparse -D}'
stencil_width = ${stencil}
[]
[]
[Postprocess]
[gold_compute]
type = ParsedCompute
buffer = gold
# g(dt)^(t/dt) up to t1 = 5 dt, then g(dt)^5 g(dt2)^((t - t1)/dt2)
expression = 'sin(2*pi*x)*sin(2*pi*y) * if(t < t1 + dt / 2, g1^(t / dt), g1^5 * g2^((t - t1) / dt2))'
constant_names = 'g1 g2 dt dt2 t1'
constant_expressions = '${fparse 1 / (1 + dt * D * lambda)} ${fparse 1 / (1 + dt2 * D * lambda)} ${dt} ${dt2} ${fparse 5 * dt}'
extra_symbols = true
[]
[diff_compute]
type = ParsedCompute
buffer = diff
expression = 'abs(u - gold)'
inputs = 'u gold'
[]
[]
[]
[TensorSolver]
type = ImplicitFDSolver
buffer = u
operators = 'dudt lap'
root_compute = dudt
[]
[Postprocessors]
[max_error]
type = TensorExtremeValuePostprocessor
buffer = diff
value_type = MAX
[]
[max_gold]
type = TensorExtremeValuePostprocessor
buffer = gold
value_type = MAX
[]
[]
[Problem]
type = TensorProblem
[]
[Executioner]
type = Transient
num_steps = 10
[TimeStepper]
type = TimeSequenceStepper
time_sequence = '${fparse 1 * dt} ${fparse 2 * dt} ${fparse 3 * dt} ${fparse 4 * dt} ${fparse 5 * dt}
${fparse 5 * dt + dt2} ${fparse 5 * dt + 2 * dt2} ${fparse 5 * dt + 3 * dt2}
${fparse 5 * dt + 4 * dt2} ${fparse 5 * dt + 5 * dt2}'
[]
[]
[Outputs]
csv = true
[]
(test/tests/implicit_fd/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
[]
(test/tests/implicit_fd/bc_mms.i)
# Implicit solve of u_t - D Laplacian(u) = f on the unit square with boundary faces, against the
# manufactured solution
#
# u = g + tf t q, g = cos(3 pi x / 2) cos(pi y) + x y, q = x^2 - y^2,
#
# with tf = 0 for the steady problem -D Laplacian(u) = f and tf = 1 for backward Euler.
# q is harmonic, so f = tf q - D Laplacian(g) is independent of time and the backward Euler
# step makes no time error: the error is the spatial one, provided the boundary values are
# taken at the time of the unknown, t + dt (at t it would be of order dt).
#
# Each face has a different condition and value source:
# left (x = 0) Dirichlet u = cos(pi y) - tf t y^2 function
# right (x = 1) Robin ra u + rb du/dx = c function
# or buffer (right_buffer, steady)
# bottom (y = 0) Neumann du/dn = -du/dy = -x buffer
# top (y = 1) Neumann du/dn = du/dy = x - 2 tf t function
# or buffer (top_buffer, values at t)
D = 0.5
n = 32
stencil = 3
tf = 0
ops = 'lap'
num_steps = 1
dt = 0.05
# Robin coefficients on the right face (the 5-point stencil supports ra = 1, rb = 0 only)
ra = 1
rb = 1
[Domain]
dim = 2
nx = ${n}
ny = ${n}
xmax = 1
ymax = 1
parallel_mode = REAL_SPACE
[]
[Functions]
[u_exact]
type = ParsedFunction
expression = 'cos(1.5*pi*x) * cos(pi*y) + x*y + ${tf} * t * (x^2 - y^2)'
[]
[f_func]
type = ParsedFunction
# Laplacian(g) = -13/4 pi^2 cos(3 pi x / 2) cos(pi y)
expression = '${tf} * (x^2 - y^2) + ${D} * 3.25 * pi^2 * cos(1.5*pi*x) * cos(pi*y)'
[]
[right_c]
type = ParsedFunction
# on x = 1: u = y + tf t (1 - y^2), du/dx = 3/2 pi cos(pi y) + y + 2 tf t
expression = '${ra} * (y + ${tf} * t * (1 - y^2)) + ${rb} * (1.5*pi*cos(pi*y) + y + 2 * ${tf} * t)'
[]
[bottom_flux]
type = ParsedFunction
expression = '-x'
[]
[top_flux]
type = ParsedFunction
expression = 'x - 2 * ${tf} * t'
[]
[]
[TensorBuffers]
[u]
[]
[dudt_u]
[]
[lap_u]
[]
[f]
[]
[bottom_flux]
[]
[right_c]
[]
[top_flux]
[]
[exact]
[]
[diff]
[]
[]
[TensorComputes]
[Initialize]
[u_ic]
type = MooseFunctionTensor
buffer = u
function = u_exact
[]
[f_ic]
type = MooseFunctionTensor
buffer = f
function = f_func
[]
# a face value: depends on the tangential coordinate x only
[bottom_flux]
type = MooseFunctionTensor
buffer = bottom_flux
function = bottom_flux
[]
[right_c]
type = MooseFunctionTensor
buffer = right_c
function = right_c
[]
[]
[BCs]
active = 'left right bottom top'
[left]
type = TensorDirichletBC
variable = u
boundary = left
value_function = u_exact
[]
[right]
type = TensorRobinBC
variable = u
boundary = right
a = ${ra}
b = ${rb}
c_function = right_c
[]
# the Robin c from a buffer instead (active = 'left right_buffer bottom top'), steady only
[right_buffer]
type = TensorRobinBC
variable = u
boundary = right
a = ${ra}
b = ${rb}
c_buffer = right_c
[]
[bottom]
type = TensorNeumannBC
variable = u
boundary = bottom
flux_buffer = bottom_flux
[]
[top]
type = TensorNeumannBC
variable = u
boundary = top
flux_function = top_flux
[]
# the top flux from a buffer instead (active = 'left right bottom top_buffer'); a buffer
# value is the one its compute wrote at the start of the substep, at t
[top_buffer]
type = TensorNeumannBC
variable = u
boundary = top
flux_buffer = top_flux
[]
[]
[Solve]
[top_flux]
type = MooseFunctionTensor
buffer = top_flux
function = top_flux
[]
[dudt]
type = TimeDerivativeOperator
buffer = dudt_u
input = u
[]
[lap]
type = FiniteDifferenceLaplacian
buffer = lap_u
input = u
factor = '${fparse -D}'
stencil_width = ${stencil}
[]
[]
[Postprocess]
[exact_compute]
type = MooseFunctionTensor
buffer = exact
function = u_exact
[]
[diff_compute]
type = ParsedCompute
buffer = diff
expression = 'abs(u - exact)'
inputs = 'u exact'
[]
[]
[]
[TensorSolver]
type = ImplicitFDSolver
buffer = u
operators = '${ops}'
source = f
linear_solver = TORCH_PCG
l_tol = 1e-10
[]
[Postprocessors]
[max_error]
type = TensorExtremeValuePostprocessor
buffer = diff
value_type = MAX
[]
[]
[Problem]
type = TensorProblem
[]
[Executioner]
type = Transient
num_steps = ${num_steps}
dt = ${dt}
[]
[Outputs]
csv = true
[]
(benchmarks/01_spinodal_decomposition/1a_implicit_fd.i)
#
# PFHub benchmark 1a (spinodal decomposition, square periodic domain) with the
# implicit finite-difference solver. Same problem as 1a.i:
#
# dc/dt = div(M grad(f'(c) - kappa Laplacian(c))), f = rho_s (c - c_alpha)^2 (c_beta - c)^2
#
# Backward Euler, with the biharmonic implicit and the double-well flux
# divergence explicit:
#
# (c - c^n)/dt + M kappa Laplacian(Laplacian(c)) = M Laplacian(f'(c^n))
#
# M = 5, kappa = 2, rho_s = 5, c_alpha = 0.3, c_beta = 0.7.
#
# substeps per time step, as in 1a_solver.i (each substep is one linear solve)
substeps = 1000
[Domain]
dim = 2
nx = 200
ny = 200
xmax = 200
ymax = 200
parallel_mode = REAL_SPACE
periodic_directions = 'X Y'
[]
[TensorBuffers]
[c]
[]
[dcdt]
[]
[bih_c]
[]
[mu]
[]
[rhs]
[]
# postprocessing
[lap_c]
[]
[F]
[]
[]
[TensorComputes]
[Initialize]
[c]
type = ParsedCompute
buffer = c
extra_symbols = true
expression = 'c0+epsilon*(cos(0.105*x)*cos(0.11*y)+(cos(0.13*x)*cos(0.087*y))^2+cos(0.025*x-0.15*y)*cos(0.07*x-0.02*y))'
constant_names = 'c0 epsilon'
constant_expressions = '0.5 0.01'
[]
[]
[Solve]
# implicit operators
[dcdt]
type = TimeDerivativeOperator
buffer = dcdt
input = c
[]
[bih]
type = FiniteDifferenceBiharmonic
buffer = bih_c
input = c
factor = 10 # M * kappa
[]
# explicit source M Laplacian(f'(c^n))
[explicit]
[mu]
type = ParsedCompute
buffer = mu
expression = 'rho_s*(c-c_alpha)^2*(c_beta-c)^2'
constant_names = 'rho_s c_alpha c_beta'
constant_expressions = '5 0.3 0.7'
derivatives = c
inputs = c
[]
[rhs]
type = FiniteDifferenceLaplacian
buffer = rhs
input = mu
factor = 5 # M
[]
[]
[]
[Postprocess]
# On the periodic grid sum |grad c|^2 = -sum c Laplacian(c) for the 3-point
# Laplacian, so the gradient energy (kappa/2)|grad c|^2 integrates to the same
# total as -(kappa/2) c Laplacian(c).
[lap_c]
type = FiniteDifferenceLaplacian
buffer = lap_c
input = c
[]
[F]
type = ParsedCompute
buffer = F
expression = 'rho_s * (c-c_alpha)^2 * (c_beta-c)^2 - kappa / 2 * c * lap_c'
constant_names = 'rho_s c_alpha c_beta kappa'
constant_expressions = '5 0.3 0.7 2'
inputs = 'c lap_c'
[]
[]
[]
[TensorSolver]
type = ImplicitFDSolver
buffer = c
operators = 'dcdt bih'
source = rhs
root_compute = explicit
substeps = ${substeps}
[]
[UserObjects]
[terminator]
type = Terminator
expression = change<1e-4
[]
[]
[Postprocessors]
[min_c]
type = TensorExtremeValuePostprocessor
buffer = c
value_type = MIN
execute_on = 'TIMESTEP_END'
[]
[max_c]
type = TensorExtremeValuePostprocessor
buffer = c
value_type = MAX
execute_on = 'TIMESTEP_END'
[]
[F]
type = TensorIntegralPostprocessor
buffer = F
[]
[change]
type = TensorIntegralChangePostprocessor
buffer = c
[]
# cumulative wall time of the run (s)
[wall_time]
type = PerfGraphData
section_name = Root
data_type = TOTAL
[]
[]
[Problem]
type = TensorProblem
[]
[Executioner]
type = Transient
num_steps = 1000
[TimeStepper]
type = IterationAdaptiveDT
growth_factor = 1.1
dt = 1
[]
dtmax = 300
[]
[TensorOutputs]
[xdmf]
type = XDMFTensorOutput
buffer = 'c mu'
enable_hdf5 = true
[]
[]
[Outputs]
csv = true
perf_graph = true
execute_on = 'TIMESTEP_END'
[]
(test/tests/implicit_fd/crank_thin_slab.i)
# Crank, "The Mathematics of Diffusion" (2nd ed.), Sec. 2.5, eq. (2.15), with the implicit
# solver: an instantaneous planar source C0 = 1 confined to the slab -h < x < h, h = 1,
# diffusing into an infinite medium, has the exact solution
# C(x, t) = (C0/2) (erf((h - x) / (2 sqrt(D t))) + erf((h + x) / (2 sqrt(D t)))).
# The same problem as real_space/thin_slab_diffusion.i, which steps it explicitly.
#
# Pseudo-1D (ny = 1, Y periodic), D = 1. The X faces are 10 diffusion lengths beyond the slab
# edge at the comparison time t = 1, so the far Dirichlet u = 0 does not measurably affect the
# solution.
# Backward Euler with substeps of dt / substeps: the error is first order in the substep.
substeps = 100
[Domain]
dim = 2
# dx = 0.1 puts the slab edges x = +/-1 on cell faces, so the box IC holds the slab mass
nx = 220
ny = 1
parallel_mode = REAL_SPACE
periodic_directions = 'Y'
xmin = -11
xmax = 11
ymin = -0.5
ymax = 0.5
[]
[Functions]
[slab]
type = ParsedFunction
expression = 'if(abs(x) < 1, 1, 0)'
[]
[]
[TensorBuffers]
[u]
[]
[dudt_u]
[]
[lap_u]
[]
[u_exact]
[]
[diff]
[]
[]
[TensorComputes]
[Initialize]
[u_ic]
type = MooseFunctionTensor
buffer = u
function = slab
[]
[]
[BCs]
# far boundaries, approximating the infinite medium
[left]
type = TensorDirichletBC
variable = u
boundary = left
value = 0
[]
[right]
type = TensorDirichletBC
variable = u
boundary = right
value = 0
[]
[]
[Solve]
[dudt]
type = TimeDerivativeOperator
buffer = dudt_u
input = u
[]
[lap]
type = FiniteDifferenceLaplacian
buffer = lap_u
input = u
factor = -1
[]
[]
[Postprocess]
[u_exact]
type = ParsedCompute
buffer = u_exact
expression = '0.5 * (erf((1 - x) / (2 * sqrt(t))) + erf((1 + x) / (2 * sqrt(t)))) + 0*y'
extra_symbols = true
[]
[diff]
type = ParsedCompute
buffer = diff
expression = 'abs(u - u_exact)'
inputs = 'u u_exact'
[]
[]
[]
[TensorSolver]
type = ImplicitFDSolver
buffer = u
operators = 'dudt lap'
root_compute = dudt
substeps = ${substeps}
l_tol = 1e-10
[]
[Postprocessors]
[max_error]
type = TensorExtremeValuePostprocessor
buffer = diff
value_type = MAX
[]
[]
[Problem]
type = TensorProblem
[]
[Executioner]
type = Transient
dt = 1
num_steps = 1
[]
[Outputs]
csv = true
[]
(test/tests/implicit_fd/cutback.i)
# Backward Euler diffusion of a Gaussian with a linear iteration limit that the first dt
# exceeds. The torch PCG solve needs about 44 iterations at dt = 0.008, 34 at 0.004 and 23
# at 0.002 (32 x 32, l_tol = 1e-8), so with l_max_its = 30 the step is cut twice and then
# converges.
D = 0.5
n = 32
[Domain]
dim = 2
nx = ${n}
ny = ${n}
xmax = 1
ymax = 1
parallel_mode = REAL_SPACE
periodic_directions = 'X Y'
[]
[TensorBuffers]
[u]
[]
[dudt_u]
[]
[lap_u]
[]
[]
[TensorComputes]
[Initialize]
[u_ic]
type = ParsedCompute
buffer = u
expression = 'exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.005)'
extra_symbols = true
[]
[]
[Solve]
[dudt]
type = TimeDerivativeOperator
buffer = dudt_u
input = u
[]
[lap]
type = FiniteDifferenceLaplacian
buffer = lap_u
input = u
factor = '${fparse -D}'
[]
[]
[]
[TensorSolver]
type = ImplicitFDSolver
buffer = u
operators = 'dudt lap'
linear_solver = TORCH_PCG
l_max_its = 30
[]
[Postprocessors]
[max_u]
type = TensorExtremeValuePostprocessor
buffer = u
value_type = MAX
[]
[]
[Problem]
type = TensorProblem
[]
[Executioner]
type = Transient
num_steps = 2
[TimeStepper]
type = TensorSolveIterationAdaptiveDT
dt = 0.008
# no growth or shrinking from the iteration count, only the cut on failure
min_iterations = 0
max_iterations = 1000
[]
[]
[Outputs]
csv = true
[]
(test/tests/implicit_fd/be_mode3d.i)
# Backward Euler diffusion u_t = D Laplacian(u) of a single periodic Fourier
# mode. On the periodic cell-centred grid the mode is an eigenvector of the
# discrete Laplacian with eigenvalue -lambda_h, so every implicit step must
# multiply it by exactly g = 1 / (1 + dt D lambda_h), up to the linear solver
# tolerance.
D = 0.5
# 5 steps of dt, then 5 steps of dt2 (dt2 != dt exercises the operator rebuild)
dt = 0.01
dt2 = 0.01
n = 16
stencil = 3
dim = 3
# Symbol of -Laplacian per axis for the wavenumber 2 pi, h = 1/n:
# 3-point: (2 - 2 cos(2 pi h)) / h^2
# 5-point: (5/2 - 8/3 cos(2 pi h) + 1/6 cos(4 pi h)) / h^2
lambda3 = '${fparse n^2 * (2 - 2 * cos(2 * pi / n))}'
lambda5 = '${fparse n^2 * (5 / 2 - 8 / 3 * cos(2 * pi / n) + 1 / 6 * cos(4 * pi / n))}'
lambda = '${fparse dim * if(stencil = 3, lambda3, lambda5)}'
[Domain]
dim = 3
nx = ${n}
ny = ${n}
nz = ${n}
zmax = 1
xmax = 1
ymax = 1
parallel_mode = REAL_SPACE
periodic_directions = 'X Y Z'
[]
[TensorBuffers]
[u]
[]
[dudt_u]
[]
[lap_u]
[]
[gold]
[]
[diff]
[]
[]
[TensorComputes]
[Initialize]
[u_ic]
type = ParsedCompute
buffer = u
expression = 'sin(2*pi*x)*sin(2*pi*y)*sin(2*pi*z)'
extra_symbols = true
[]
[]
[Solve]
[dudt]
type = TimeDerivativeOperator
buffer = dudt_u
input = u
[]
[lap]
type = FiniteDifferenceLaplacian
buffer = lap_u
input = u
factor = '${fparse -D}'
stencil_width = ${stencil}
[]
[]
[Postprocess]
[gold_compute]
type = ParsedCompute
buffer = gold
# g(dt)^(t/dt) up to t1 = 5 dt, then g(dt)^5 g(dt2)^((t - t1)/dt2)
expression = 'sin(2*pi*x)*sin(2*pi*y)*sin(2*pi*z) * if(t < t1 + dt / 2, g1^(t / dt), g1^5 * g2^((t - t1) / dt2))'
constant_names = 'g1 g2 dt dt2 t1'
constant_expressions = '${fparse 1 / (1 + dt * D * lambda)} ${fparse 1 / (1 + dt2 * D * lambda)} ${dt} ${dt2} ${fparse 5 * dt}'
extra_symbols = true
[]
[diff_compute]
type = ParsedCompute
buffer = diff
expression = 'abs(u - gold)'
inputs = 'u gold'
[]
[]
[]
[TensorSolver]
type = ImplicitFDSolver
buffer = u
operators = 'dudt lap'
root_compute = dudt
[]
[Postprocessors]
[max_error]
type = TensorExtremeValuePostprocessor
buffer = diff
value_type = MAX
[]
[max_gold]
type = TensorExtremeValuePostprocessor
buffer = gold
value_type = MAX
[]
[]
[Problem]
type = TensorProblem
[]
[Executioner]
type = Transient
num_steps = 10
[TimeStepper]
type = TimeSequenceStepper
time_sequence = '${fparse 1 * dt} ${fparse 2 * dt} ${fparse 3 * dt} ${fparse 4 * dt} ${fparse 5 * dt}
${fparse 5 * dt + dt2} ${fparse 5 * dt + 2 * dt2} ${fparse 5 * dt + 3 * dt2}
${fparse 5 * dt + 4 * dt2} ${fparse 5 * dt + 5 * dt2}'
[]
[]
[Outputs]
csv = true
[]
(test/tests/implicit_fd/crank_semi_infinite.i)
# Crank, "The Mathematics of Diffusion" (2nd ed.), Sec. 3.3(i), eq. (3.13), with the implicit
# solver: a semi-infinite medium (x > 0) with zero initial concentration and constant surface
# concentration C1 = 1 at x = 0 has the exact solution
# C(x, t) = erfc(x / (2 sqrt(D t))).
# The same problem as real_space/semi_infinite_diffusion.i, which steps it explicitly.
#
# Pseudo-1D (ny = 1, Y periodic), D = 1. The domain extends 10 diffusion lengths at the
# comparison time t = 1, so the far Dirichlet u = 0 does not measurably affect the solution.
# Backward Euler with substeps of dt / substeps: the error is first order in the substep.
substeps = 100
[Domain]
dim = 2
nx = 100
ny = 1
parallel_mode = REAL_SPACE
periodic_directions = 'Y'
xmin = 0
xmax = 10
ymin = -0.5
ymax = 0.5
[]
[TensorBuffers]
[u]
[]
[dudt_u]
[]
[lap_u]
[]
[u_exact]
[]
[diff]
[]
[]
[TensorComputes]
# no [Initialize]: buffers start at zero, the C0 = 0 initial condition
[BCs]
[left]
type = TensorDirichletBC
variable = u
boundary = left
value = 1
[]
[right]
# far boundary, approximating the semi-infinite medium's x -> infinity condition
type = TensorDirichletBC
variable = u
boundary = right
value = 0
[]
[]
[Solve]
[dudt]
type = TimeDerivativeOperator
buffer = dudt_u
input = u
[]
[lap]
type = FiniteDifferenceLaplacian
buffer = lap_u
input = u
factor = -1
[]
[]
[Postprocess]
[u_exact]
type = ParsedCompute
buffer = u_exact
expression = 'erfc(x / (2 * sqrt(t))) + 0*y'
extra_symbols = true
[]
[diff]
type = ParsedCompute
buffer = diff
expression = 'abs(u - u_exact)'
inputs = 'u u_exact'
[]
[]
[]
[TensorSolver]
type = ImplicitFDSolver
buffer = u
operators = 'dudt lap'
root_compute = dudt
substeps = ${substeps}
l_tol = 1e-10
[]
[Postprocessors]
[max_error]
type = TensorExtremeValuePostprocessor
buffer = diff
value_type = MAX
[]
[]
[Problem]
type = TensorProblem
[]
[Executioner]
type = Transient
dt = 1
num_steps = 1
[]
[Outputs]
csv = true
[]