LStableDirk2
Second order diagonally implicit Runge Kutta method (Dirk) with two stages.
This method can be expressed as a Runge-Kutta method with the following Butcher Tableau:
where
The stability function for this method is:
The method is L-stable:
Notes
This method is derived in detail in Alexander (1977). This method is more expensive than Crank-Nicolson, but has the advantage of being L-stable (the same type of stability as the implicit Euler method) so may be more suitable for "stiff" problems.
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.
Input Files
References
- R. Alexander.
Diagonally implicit runge-kutta methods for stiff odes.
SIAM J. Numer. Anal., 14(6):1006–1021, 1977.[BibTeX]
@article{alexander1967, author = "Alexander, R.", title = "Diagonally implicit Runge-Kutta Methods for Stiff ODEs", journal = "SIAM J. Numer. Anal.", year = "1977", pages = "1006-1021", volume = "14", number = "6" }
(test/tests/time_integrators/dirk/dirk-2d-heat-adap.i)
[Mesh]
type = GeneratedMesh
dim = 2
xmin = -1
xmax = 1
ymin = -1
ymax = 1
nx = 4
ny = 4
elem_type = QUAD4
[]
[Variables]
active = 'u'
[./u]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = ConstantIC
value = 0
[../]
[../]
[]
[Functions]
[./forcing_fn]
type = ParsedFunction
expression = 3*t*t*((x*x)+(y*y))-(4*t*t*t)
[../]
[./exact_fn]
type = ParsedFunction
expression = t*t*t*((x*x)+(y*y))
[../]
[]
[Kernels]
active = 'diff ie ffn'
[./ie]
type = TimeDerivative
variable = u
[../]
[./diff]
type = Diffusion
variable = u
[../]
[./ffn]
type = BodyForce
variable = u
function = forcing_fn
[../]
[]
[BCs]
active = 'all'
[./all]
type = FunctionDirichletBC
variable = u
boundary = '0 1 2 3'
function = exact_fn
[../]
[./left]
type = DirichletBC
variable = u
boundary = 3
value = 0
[../]
[./right]
type = DirichletBC
variable = u
boundary = 1
value = 1
[../]
[]
[Postprocessors]
[./l2_err]
type = ElementL2Error
variable = u
function = exact_fn
[../]
[]
[Executioner]
type = Transient
solve_type = 'PJFNK'
petsc_options_iname = '-pc_type -pc_hypre_type'
petsc_options_value = 'hypre boomeramg'
start_time = 0.0
num_steps = 5
dt = 0.25
[./TimeIntegrator]
type = LStableDirk2
[../]
[./Adaptivity]
refine_fraction = 0.07
coarsen_fraction = 0.
max_h_level = 4
[../]
[]
[Outputs]
execute_on = 'timestep_end'
exodus = true
[]
(test/tests/time_integrators/scalar/scalar.i)
[Mesh]
type = GeneratedMesh
dim = 2
xmin = 0
xmax = 1
ymin = 0
ymax = 1
nx = 1
ny = 1
elem_type = QUAD4
[]
[Variables]
[./n]
family = SCALAR
order = FIRST
initial_condition = 1
[../]
[]
[ScalarKernels]
[./dn]
type = ODETimeDerivative
variable = n
[../]
[./ode1]
type = ParsedODEKernel
expression = '-n'
variable = n
# implicit = false
[../]
[]
[Executioner]
type = Transient
[./TimeIntegrator]
# type = ImplicitEuler
# type = BDF2
type = CrankNicolson
# type = ImplicitMidpoint
# type = LStableDirk2
# type = LStableDirk3
# type = LStableDirk4
# type = AStableDirk4
#
# Explicit methods
# type = ExplicitEuler
# type = ExplicitMidpoint
# type = Heun
# type = Ralston
[../]
start_time = 0
end_time = 1
dt = 0.001
dtmin = 0.001 # Don't allow timestep cutting
solve_type = 'PJFNK'
nl_max_its = 2
nl_abs_tol = 1.e-12 # This is an ODE, so nl_abs_tol makes sense.
[]
[Functions]
[./exact_solution]
type = ParsedFunction
expression = exp(t)
[../]
[]
[Postprocessors]
[./error_n]
# Post processor that computes the difference between the computed
# and exact solutions. For the exact solution used here, the
# error at the final time should converge at O(dt^p), where p is
# the order of the method.
type = ScalarL2Error
variable = n
function = exact_solution
# final is not currently supported for Postprocessor execute_on...
# execute_on = 'final'
[../]
[]
[Outputs]
csv = true
[]
(test/tests/time_integrators/scalar/stiff.i)
# This is a linear model problem described in Frank et al, "Order
# results for implicit Runge-Kutta methods applied to stiff systems",
# SIAM J. Numerical Analysis, vol. 22, no. 3, 1985, pp. 515-534.
#
# Problems "PL" and "PNL" from page 527 of the paper:
# { dy1/dt = lambda*y1 + y2**p, y1(0) = -1/(lambda+p)
# { dy2/dt = -y2, y2(0) = 1
#
# The exact solution is:
# y1 = -exp(-p*t)/(lambda+p)
# y2 = exp(-t)
#
# According to the following paragraph from the reference above, the
# p=1 version of this problem should not exhibit order reductions
# regardless of stiffness, while the nonlinear version (p>=2) will
# exhibit order reductions down to the "stage order" of the method for
# lambda large, negative.
# Use Dollar Bracket Expressions (DBEs) to set the value of LAMBDA in
# a single place. You can also set this on the command line with
# e.g. LAMBDA=-4, but note that this does not seem to override the
# value set in the input file. This is a bit different from the way
# that command line values normally work...
# Note that LAMBDA == Y2_EXPONENT is not allowed!
# LAMBDA = -10
# Y2_EXPONENT = 2
[Mesh]
type = GeneratedMesh
dim = 2
xmin = 0
xmax = 1
ymin = 0
ymax = 1
nx = 1
ny = 1
elem_type = QUAD4
[]
[Variables]
[./y1]
family = SCALAR
order = FIRST
[../]
[./y2]
family = SCALAR
order = FIRST
[../]
[]
[ICs]
[./y1_init]
type = FunctionScalarIC
variable = y1
function = y1_exact
[../]
[./y2_init]
type = FunctionScalarIC
variable = y2
function = y2_exact
[../]
[]
[ScalarKernels]
[./y1_time]
type = ODETimeDerivative
variable = y1
[../]
[./y1_space]
type = ParsedODEKernel
variable = y1
expression = '-(${LAMBDA})*y1 - y2^${Y2_EXPONENT}'
coupled_variables = 'y2'
[../]
[./y2_time]
type = ODETimeDerivative
variable = y2
[../]
[./y2_space]
type = ParsedODEKernel
variable = y2
expression = 'y2'
[../]
[]
[Executioner]
type = Transient
[./TimeIntegrator]
type = LStableDirk2
[../]
start_time = 0
end_time = 1
dt = 0.125
solve_type = 'PJFNK'
nl_max_its = 6
nl_abs_tol = 1.e-13
nl_rel_tol = 1.e-32 # Force nl_abs_tol to be used.
line_search = 'none'
[]
[Functions]
[./y1_exact]
type = ParsedFunction
expression = '-exp(-${Y2_EXPONENT}*t)/(lambda+${Y2_EXPONENT})'
symbol_names = 'lambda'
symbol_values = ${LAMBDA}
[../]
[./y2_exact]
type = ParsedFunction
expression = exp(-t)
[../]
[]
[Postprocessors]
[./error_y1]
type = ScalarL2Error
variable = y1
function = y1_exact
execute_on = 'initial timestep_end'
[../]
[./error_y2]
type = ScalarL2Error
variable = y2
function = y2_exact
execute_on = 'initial timestep_end'
[../]
[./max_error_y1]
# Estimate ||e_1||_{\infty}
type = TimeExtremeValue
value_type = max
postprocessor = error_y1
execute_on = 'initial timestep_end'
[../]
[./max_error_y2]
# Estimate ||e_2||_{\infty}
type = TimeExtremeValue
value_type = max
postprocessor = error_y2
execute_on = 'initial timestep_end'
[../]
[./value_y1]
type = ScalarVariable
variable = y1
execute_on = 'initial timestep_end'
[../]
[./value_y2]
type = ScalarVariable
variable = y2
execute_on = 'initial timestep_end'
[../]
[./value_y1_abs_max]
type = TimeExtremeValue
value_type = abs_max
postprocessor = value_y1
execute_on = 'initial timestep_end'
[../]
[./value_y2_abs_max]
type = TimeExtremeValue
value_type = abs_max
postprocessor = value_y2
execute_on = 'initial timestep_end'
[../]
[]
[Outputs]
csv = true
[]
(test/tests/time_integrators/dirk/dirk-2d-heat.i)
#
# Testing a solution that is second order in space and first order in time.
#
[Mesh]
type = GeneratedMesh
dim = 2
xmin = 0
xmax = 1
ymin = 0
ymax = 1
nx = 20
ny = 20
elem_type = QUAD9
[]
[Variables]
[./u]
order = SECOND
family = LAGRANGE
[./InitialCondition]
type = FunctionIC
function = exact_fn
[../]
[../]
[]
[Functions]
[./forcing_fn]
type = ParsedFunction
expression = ((x*x)+(y*y))-(4*t)
[../]
[./exact_fn]
type = ParsedFunction
expression = t*((x*x)+(y*y))
[../]
[]
[Kernels]
[./ie]
type = TimeDerivative
variable = u
[../]
[./diff]
type = Diffusion
variable = u
[../]
[./ffn]
type = BodyForce
variable = u
function = forcing_fn
[../]
[]
[BCs]
[./all]
type = FunctionDirichletBC
variable = u
boundary = '0 1 2 3'
function = exact_fn
[../]
[]
[Postprocessors]
[./l2_err]
type = ElementL2Error
variable = u
function = exact_fn
[../]
[]
[Executioner]
type = Transient
solve_type = 'PJFNK'
petsc_options_iname = '-pc_type -pc_hypre_type'
petsc_options_value = 'hypre boomeramg'
start_time = 0.0
end_time = 1.0
dt = 1.0
nl_abs_tol=1e-13
nl_rel_tol=1e-13
[./TimeIntegrator]
type = LStableDirk2
[../]
[]
[Outputs]
execute_on = 'timestep_end'
exodus = true
[]
(test/tests/time_integrators/multi_stage_time_integrator/unconverged_1st_stage.i)
# This test is designed to check that a time step solve should stop if *any*
# time integrator solve stage fails, not just the *last* stage. If a time
# integrator does not check convergence per stage, then a time step proceeds
# past intermediate stages without checking nonlinear convergence. This test
# is designed to check that the 2nd stage is never even entered by making it
# impossible for the first stage to converge.
[Mesh]
type = GeneratedMesh
dim = 1
xmin = -1
xmax = 1
nx = 5
[]
[Functions]
[./ic]
type = ParsedFunction
expression = 0
[../]
[./forcing_fn]
type = ParsedFunction
expression = x
[../]
[./exact_fn]
type = ParsedFunction
expression = t*x
[../]
[]
[Variables]
[./u]
[../]
[]
[Kernels]
[./time]
type = TimeDerivative
variable = u
[../]
[./diff]
type = Diffusion
variable = u
[../]
[./body]
type = BodyForce
variable = u
function = forcing_fn
[../]
[]
[ICs]
[./u_ic]
type = FunctionIC
variable = u
function = ic
[../]
[]
[BCs]
[./bcs]
type = FunctionDirichletBC
variable = u
boundary = '0 1'
function = exact_fn
[../]
[]
[Executioner]
type = Transient
[./TimeIntegrator]
type = LStableDirk2
[../]
num_steps = 1
abort_on_solve_fail = true
solve_type = NEWTON
nl_max_its = 0
[]