ImplicitFDSolver

Implicit solver for fields whose linear terms are finite-difference operators, solving sum_i A_i u = sum_i r_i + f for each field with a Krylov method.

Overview

Solves each field in "buffer" implicitly from

where the are the finite-difference operators listed in "operators" (as written, with no sign change), the their right-hand side contributions, and the optional explicit "source". Backward Euler diffusion is written with a TimeDerivativeOperator (action , contribution ) and a FiniteDifferenceLaplacian with factor = -D.

The operator is applied matrix-free with each operator's own kernel. The system must be symmetric positive definite, or positive semidefinite with the constants as its nullspace (a steady periodic Laplacian). In the singular case the right-hand side must have zero mean up to round-off, and the zero-mean solution is returned.

The solve stops when the unpreconditioned residual satisfies "l_tol" , starting from the previous solution. When a solve does not converge within "l_max_its" iterations, a TensorSolveIterationAdaptiveDT time stepper retries the time step with a smaller time step (with one substep per time step); with any other time stepper it is an error. After the solve, the operators' output buffers hold the operators applied to the new solution.

The "root_compute" is evaluated before every solve to update the explicit terms, such as the source. The field operators in it are skipped, since the solve applies them itself.

The PETSC linear solver (the default except on MPS) is PETSc's conjugate gradient method preconditioned by geometric multigrid (KSPCG + PCMG), in double precision on CPU or CUDA. The fine level applies the operators matrix-free; the coarse levels are the operators rediscretised on grids coarsened by two per axis and assembled, connected by cell-centred multilinear interpolation. Every level but the coarsest, which is solved by LU, is smoothed by two sweeps of the smoother. The default, CHEBYSHEV, is PCMG's default Chebyshev smoother, preconditioned by SOR on the coarse levels on CPU and by Jacobi on the matrix-free fine level and on every level on CUDA (PETSc's SOR runs on the host). Chebyshev estimates each level's eigenvalue bounds with a few Krylov steps, again on every level whenever the operators change (e.g. with the time step). smoother = DAMPED_JACOBI needs no estimate: it is damped Jacobi (KSPRICHARDSON + PCJACOBI) with the weight , where is the Gershgorin bound on the Jacobi scaled spectrum from the level's stencil coefficients , the maximum over the rows including those next to a boundary face. The weight is recomputed when the operators change (e.g. with the time step), unless -<buffer>_mg_levels_ksp_richardson_scale is set. It takes more iterations on the biharmonic, and their number grows with the grid size. The divergence test is off by default (-<buffer>_ksp_divtol sets one), as in TORCH_PCG: the solve is warm started from the old solution and its residual compared to , which scales as , so at large time steps the first residuals exceed PETSc's default divergence tolerance on a converging solve. PETSc options apply with the field name as prefix: petsc_options_iname = '-pc_type' with petsc_options_value = 'jacobi' sets -u_pc_type jacobi for field u. The block options apply to every field, and command line options (e.g. -u_mg_levels_ksp_max_it 3) override them.

The TORCH_PCG linear solver is a Jacobi preconditioned conjugate gradient method in torch. It runs on every device, in single precision on devices without double precision support.

Boundary conditions

On a non-periodic axis, each face of has a TensorDirichletBC, TensorNeumannBC or TensorRobinBC (see the BCs system for face values and their timing). With boundary values the operator action is affine, , and the solver solves , with the action on , recomputed every solve. A Dirichlet face, or a Robin face with , makes the system definite; with Neumann faces only it keeps the constants as its nullspace, and the right-hand side, including the fluxes, must then have zero mean up to round-off. Source and fluxes that balance in the continuous problem generally leave a discrete mean of order , which is an error: the data must balance on the grid. Supported are the time derivative and the Laplacian, with Dirichlet, Neumann and Robin faces for the 3-point stencil and Dirichlet and Neumann faces for the 5-point stencil, whose second ghost layer makes a Robin operator nonsymmetric. A Robin face with is an error: its ghost term turns negative once , after which the operator may be indefinite, depending on the rest of the problem, so whether it could be solved would depend on the grid. The biharmonic on a non-periodic domain is an error too.

On the coarse multigrid levels, a stencil entry that reaches a ghost cell is folded onto the mirrored interior cell, times the ghost factor of the homogeneous condition with for ghost layer and level spacing ( for Dirichlet, for Neumann). The interpolation folds the same way. This affects only the preconditioner.

Example Input File Syntax

[TensorSolver<<<{"href": "../../syntax/TensorSolver/index.html"}>>>]
  type = ImplicitFDSolver
  buffer = u
  operators = 'dudt lap'
  root_compute = dudt
[]
(test/tests/implicit_fd/be_mode.i)

With boundary conditions whose values vary along the faces:

[TensorComputes<<<{"href": "../../syntax/TensorComputes/index.html"}>>>]
  [BCs<<<{"href": "../../syntax/TensorComputes/BCs/index.html"}>>>]
    active<<<{"description": "If specified only the blocks named will be visited and made active"}>>> = 'left right bottom top'
    [left]
      type = TensorDirichletBC<<<{"description": "Dirichlet boundary condition (u = value) on a non-periodic domain face.", "href": "../tensor_computes/TensorDirichletBC.html"}>>>
      variable<<<{"description": "The field this boundary condition applies to."}>>> = u
      boundary<<<{"description": "Which face of the domain this boundary condition applies to."}>>> = left
      value_function<<<{"description": "Function giving the field value along the face, evaluated at the face position at the solve time."}>>> = u_exact
    []
    [right]
      type = TensorRobinBC<<<{"description": "General Robin boundary condition a*u + b*du/dn = c on a non-periodic domain face.", "href": "../tensor_computes/TensorRobinBC.html"}>>>
      variable<<<{"description": "The field this boundary condition applies to."}>>> = u
      boundary<<<{"description": "Which face of the domain this boundary condition applies to."}>>> = right
      a<<<{"description": "Robin coefficient on u (a*u + b*du/dn = c)."}>>> = ${ra}
      b<<<{"description": "Robin coefficient on du/dn (a*u + b*du/dn = c)."}>>> = ${rb}
      c_function<<<{"description": "Function giving c along the face, evaluated at the face position at the solve time."}>>> = right_c
    []
    # the Robin c from a buffer instead (active = 'left right_buffer bottom top'), steady only
    [right_buffer]
      type = TensorRobinBC<<<{"description": "General Robin boundary condition a*u + b*du/dn = c on a non-periodic domain face.", "href": "../tensor_computes/TensorRobinBC.html"}>>>
      variable<<<{"description": "The field this boundary condition applies to."}>>> = u
      boundary<<<{"description": "Which face of the domain this boundary condition applies to."}>>> = right
      a<<<{"description": "Robin coefficient on u (a*u + b*du/dn = c)."}>>> = ${ra}
      b<<<{"description": "Robin coefficient on du/dn (a*u + b*du/dn = c)."}>>> = ${rb}
      c_buffer<<<{"description": "Buffer whose first owned cell layer next to the face holds c along the face. The compute that writes it must evaluate c at the face position, not at the cell centres."}>>> = right_c
    []
    [bottom]
      type = TensorNeumannBC<<<{"description": "Neumann boundary condition (du/dn = flux) on a non-periodic domain face.", "href": "../tensor_computes/TensorNeumannBC.html"}>>>
      variable<<<{"description": "The field this boundary condition applies to."}>>> = u
      boundary<<<{"description": "Which face of the domain this boundary condition applies to."}>>> = bottom
      flux_buffer<<<{"description": "Buffer whose first owned cell layer next to the face holds the outward normal flux along the face. The compute that writes it must evaluate it at the face position, not at the cell centres."}>>> = bottom_flux
    []
    [top]
      type = TensorNeumannBC<<<{"description": "Neumann boundary condition (du/dn = flux) on a non-periodic domain face.", "href": "../tensor_computes/TensorNeumannBC.html"}>>>
      variable<<<{"description": "The field this boundary condition applies to."}>>> = u
      boundary<<<{"description": "Which face of the domain this boundary condition applies to."}>>> = top
      flux_function<<<{"description": "Function giving the outward normal flux along the face, evaluated at the face position at the solve time."}>>> = 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<<<{"description": "Neumann boundary condition (du/dn = flux) on a non-periodic domain face.", "href": "../tensor_computes/TensorNeumannBC.html"}>>>
      variable<<<{"description": "The field this boundary condition applies to."}>>> = u
      boundary<<<{"description": "Which face of the domain this boundary condition applies to."}>>> = top
      flux_buffer<<<{"description": "Buffer whose first owned cell layer next to the face holds the outward normal flux along the face. The compute that writes it must evaluate it at the face position, not at the cell centres."}>>> = top_flux
    []
  []
[]
(test/tests/implicit_fd/bc_mms.i)

Input Parameters

  • bufferThe solved fields

    C++ Type:std::vector<std::string>

    Controllable:No

    Description:The solved fields

  • operatorsFor each field (separated by ';'), the LinearFDOperator computes whose sum is the system operator

    C++ Type:std::vector<std::vector<std::string>>

    Controllable:No

    Description:For each field (separated by ';'), the LinearFDOperator computes whose sum is the system operator

Required Parameters

  • forward_bufferThese buffers are updated with the corresponding buffers from forward_buffer_old. No integration is performed. Buffer forwarding is used only to resolve cyclic dependencies.

    C++ Type:std::vector<std::string>

    Controllable:No

    Description:These buffers are updated with the corresponding buffers from forward_buffer_old. No integration is performed. Buffer forwarding is used only to resolve cyclic dependencies.

  • forward_buffer_newNew values to update `forward_buffer` with.

    C++ Type:std::vector<std::string>

    Controllable:No

    Description:New values to update `forward_buffer` with.

  • l_max_its10000Maximum number of linear iterations

    Default:10000

    C++ Type:unsigned int

    Range:l_max_its > 0

    Controllable:No

    Description:Maximum number of linear iterations

  • l_tol1e-08Relative tolerance on the unpreconditioned residual norm

    Default:1e-08

    C++ Type:Real

    Unit:(no unit assumed)

    Range:l_tol > 0

    Controllable:No

    Description:Relative tolerance on the unpreconditioned residual norm

  • linear_solverAUTOLinear solver backend. AUTO selects TORCH_PCG on MPS and PETSC on all other devices.

    Default:AUTO

    C++ Type:MooseEnum

    Options:AUTO, PETSC, TORCH_PCG

    Controllable:No

    Description:Linear solver backend. AUTO selects TORCH_PCG on MPS and PETSC on all other devices.

  • root_computePrimary compute object that updates the buffers. This is usually a ComputeGroup object. A ComputeGroup encompassing all computes will be generated automatically if the user does not provide this parameter.

    C++ Type:std::string

    Controllable:No

    Description:Primary compute object that updates the buffers. This is usually a ComputeGroup object. A ComputeGroup encompassing all computes will be generated automatically if the user does not provide this parameter.

  • smootherCHEBYSHEVMultigrid level smoother of the PETSC linear solver. CHEBYSHEV is PCMG's default, with SOR on the coarse levels on CPU and Jacobi on the fine level and on CUDA; DAMPED_JACOBI is Jacobi weighted by the Gershgorin bound of each level operator.

    Default:CHEBYSHEV

    C++ Type:MooseEnum

    Options:CHEBYSHEV, DAMPED_JACOBI

    Controllable:No

    Description:Multigrid level smoother of the PETSC linear solver. CHEBYSHEV is PCMG's default, with SOR on the coarse levels on CPU and Jacobi on the fine level and on CUDA; DAMPED_JACOBI is Jacobi weighted by the Gershgorin bound of each level operator.

  • sourceExplicit source term f for each field (empty for no sources)

    C++ Type:std::vector<std::string>

    Controllable:No

    Description:Explicit source term f for each field (empty for no sources)

  • substeps1Solver substeps per time step.

    Default:1

    C++ Type:unsigned int

    Controllable:No

    Description:Solver substeps per time step.

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

  • petsc_optionsSingleton PETSc options

    C++ Type:MultiMooseEnum

    Options:-dm_moose_print_embedding, -dm_view, -KSP_CONVERGED_REASON, -KSP_GMRES_MODIFIEDGRAMSCHMIDT, -KSP_MONITOR, -KSP_MONITOR_SNES_LG, -SNES_KSP_EW, -SNES_CONVERGED_REASON, -SNES_KSP, -SNES_LINESEARCH_MONITOR, -SNES_MF, -SNES_MF_OPERATOR, -SNES_MONITOR, -SNES_TEST_DISPLAY, -SNES_VIEW, -SNES_MONITOR_CANCEL

    Controllable:No

    Description:Singleton PETSc options

  • petsc_options_inameNames of PETSc name/value pairs

    C++ Type:MultiMooseEnum

    Options:-mat_fd_coloring_err, -mat_fd_type, -mat_mffd_type, -pc_asm_overlap, -pc_factor_levels, -pc_factor_mat_ordering_type, -pc_hypre_boomeramg_grid_sweeps_all, -pc_hypre_boomeramg_max_iter, -pc_hypre_boomeramg_strong_threshold, -pc_hypre_type, -pc_type, -sub_pc_type, -KSP_ATOL, -KSP_GMRES_RESTART, -KSP_MAX_IT, -KSP_PC_SIDE, -KSP_RTOL, -KSP_TYPE, -SUB_KSP_TYPE, -SNES_ATOL, -SNES_LINESEARCH_TYPE, -SNES_LS, -SNES_MAX_IT, -SNES_RTOL, -SNES_DIVERGENCE_TOLERANCE, -SNES_TYPE

    Controllable:No

    Description:Names of PETSc name/value pairs

  • petsc_options_valueValues of PETSc name/value pairs (must correspond with "petsc_options_iname"

    C++ Type:std::vector<std::string>

    Controllable:No

    Description:Values of PETSc name/value pairs (must correspond with "petsc_options_iname"

Petsc Parameters

Input Files