SMP

Single matrix preconditioner (SMP) builds a preconditioner using user defined off-diagonal parts of the Jacobian.

Overview

The Single Matrix Preconditioner (SMP) builds one matrix for preconditioning. As an example, consider the system:

var element = document.getElementById("moose-equation-9ea189da-75f9-4800-918f-0f4f6864e660");katex.render("\\begin{aligned} \\nabla \\cdot k(s,T) \\nabla T &= 0 \\\\ \\nabla \\cdot D(s,T) \\nabla s &= 0 , \\end{aligned},", element, {displayMode:true,throwOnError:false,macros:{"\\bs":"\\boldsymbol{#1}","\\normal":"\\bs{n}","\\grad":"\\bs{\\nabla}","\\divergence":"\\grad \\cdot","\\norm":"\\left\\lVert#1\\right\\rVert","\\abs":"\\left\\lvert#1\\right\\rvert","\\tr":"\\text{tr}","\\dev":"\\text{dev}","\\sym":"\\text{sym}","\\diff":"\\text{ d}#1","\\activepart":"{\\left< A \\right>}","\\inactivepart":"{\\left< I \\right>}","\\Gc":"{\\mathcal{G}_c}","\\strain":"\\bs{\\varepsilon}","\\stress":"\\bs{\\sigma}","\\macaulay":"\\left<#1\\right>","\\body":"\\Omega","\\bodyboundary":"{\\partial\\body}","\\ep":"{\\varepsilon^p}","\\ep0":"{\\varepsilon_0^p}","\\epdot":"{\\dot{\\varepsilon}}^p","\\epdot0":"{\\dot{\\varepsilon}}_0^p","\\bfa":"\\boldsymbol{a}","\\bfb":"\\boldsymbol{b}","\\bfc":"\\boldsymbol{c}","\\bfd":"\\boldsymbol{d}","\\bfe":"\\boldsymbol{e}","\\bff":"\\boldsymbol{f}","\\bfg":"\\boldsymbol{g}","\\bfh":"\\boldsymbol{h}","\\bfi":"\\boldsymbol{i}","\\bfj":"\\boldsymbol{j}","\\bfk":"\\boldsymbol{k}","\\bfl":"\\boldsymbol{l}","\\bfm":"\\boldsymbol{m}","\\bfn":"\\boldsymbol{n}","\\bfo":"\\boldsymbol{o}","\\bfp":"\\boldsymbol{p}","\\bfq":"\\boldsymbol{q}","\\bfr":"\\boldsymbol{r}","\\bfs":"\\boldsymbol{s}","\\bft":"\\boldsymbol{t}","\\bfu":"\\boldsymbol{u}","\\bfv":"\\boldsymbol{v}","\\bfw":"\\boldsymbol{w}","\\bfx":"\\boldsymbol{x}","\\bfy":"\\boldsymbol{y}","\\bfz":"\\boldsymbol{z}","\\bfA":"\\boldsymbol{A}","\\bfB":"\\boldsymbol{B}","\\bfC":"\\boldsymbol{C}","\\bfD":"\\boldsymbol{D}","\\bfE":"\\boldsymbol{E}","\\bfF":"\\boldsymbol{F}","\\bfG":"\\boldsymbol{G}","\\bfH":"\\boldsymbol{H}","\\bfI":"\\boldsymbol{I}","\\bfJ":"\\boldsymbol{J}","\\bfK":"\\boldsymbol{K}","\\bfL":"\\boldsymbol{L}","\\bfM":"\\boldsymbol{M}","\\bfN":"\\boldsymbol{N}","\\bfO":"\\boldsymbol{O}","\\bfP":"\\boldsymbol{P}","\\bfQ":"\\boldsymbol{Q}","\\bfR":"\\boldsymbol{R}","\\bfS":"\\boldsymbol{S}","\\bfT":"\\boldsymbol{T}","\\bfU":"\\boldsymbol{U}","\\bfV":"\\boldsymbol{V}","\\bfW":"\\boldsymbol{W}","\\bfX":"\\boldsymbol{X}","\\bfY":"\\boldsymbol{Y}","\\bfZ":"\\boldsymbol{Z}"}});

Users can then specify which off-diagonal blokcs of the matrix to use like


off_diag_row    = 's'
off_diag_column = 'T'

Which would produce a preconditioning matrix like this:

var element = document.getElementById("moose-equation-be928a0d-791e-4d70-a6c4-55ee549c1a24");katex.render("\\boldsymbol{M} \\equiv \\begin{bmatrix} \\left(k(s,T) \\nabla \\phi_j, \\nabla \\psi_i\\right) & \\boldsymbol{0} \\\\[3pt] \\left(\\frac{\\partial D(s,T)}{\\partial T_j} \\nabla s, \\nabla \\psi_i\\right) & \\left(D(s,T) \\nabla \\phi_j, \\nabla\\psi_i\\right) \\end{bmatrix}", element, {displayMode:true,throwOnError:false,macros:{"\\bs":"\\boldsymbol{#1}","\\normal":"\\bs{n}","\\grad":"\\bs{\\nabla}","\\divergence":"\\grad \\cdot","\\norm":"\\left\\lVert#1\\right\\rVert","\\abs":"\\left\\lvert#1\\right\\rvert","\\tr":"\\text{tr}","\\dev":"\\text{dev}","\\sym":"\\text{sym}","\\diff":"\\text{ d}#1","\\activepart":"{\\left< A \\right>}","\\inactivepart":"{\\left< I \\right>}","\\Gc":"{\\mathcal{G}_c}","\\strain":"\\bs{\\varepsilon}","\\stress":"\\bs{\\sigma}","\\macaulay":"\\left<#1\\right>","\\body":"\\Omega","\\bodyboundary":"{\\partial\\body}","\\ep":"{\\varepsilon^p}","\\ep0":"{\\varepsilon_0^p}","\\epdot":"{\\dot{\\varepsilon}}^p","\\epdot0":"{\\dot{\\varepsilon}}_0^p","\\bfa":"\\boldsymbol{a}","\\bfb":"\\boldsymbol{b}","\\bfc":"\\boldsymbol{c}","\\bfd":"\\boldsymbol{d}","\\bfe":"\\boldsymbol{e}","\\bff":"\\boldsymbol{f}","\\bfg":"\\boldsymbol{g}","\\bfh":"\\boldsymbol{h}","\\bfi":"\\boldsymbol{i}","\\bfj":"\\boldsymbol{j}","\\bfk":"\\boldsymbol{k}","\\bfl":"\\boldsymbol{l}","\\bfm":"\\boldsymbol{m}","\\bfn":"\\boldsymbol{n}","\\bfo":"\\boldsymbol{o}","\\bfp":"\\boldsymbol{p}","\\bfq":"\\boldsymbol{q}","\\bfr":"\\boldsymbol{r}","\\bfs":"\\boldsymbol{s}","\\bft":"\\boldsymbol{t}","\\bfu":"\\boldsymbol{u}","\\bfv":"\\boldsymbol{v}","\\bfw":"\\boldsymbol{w}","\\bfx":"\\boldsymbol{x}","\\bfy":"\\boldsymbol{y}","\\bfz":"\\boldsymbol{z}","\\bfA":"\\boldsymbol{A}","\\bfB":"\\boldsymbol{B}","\\bfC":"\\boldsymbol{C}","\\bfD":"\\boldsymbol{D}","\\bfE":"\\boldsymbol{E}","\\bfF":"\\boldsymbol{F}","\\bfG":"\\boldsymbol{G}","\\bfH":"\\boldsymbol{H}","\\bfI":"\\boldsymbol{I}","\\bfJ":"\\boldsymbol{J}","\\bfK":"\\boldsymbol{K}","\\bfL":"\\boldsymbol{L}","\\bfM":"\\boldsymbol{M}","\\bfN":"\\boldsymbol{N}","\\bfO":"\\boldsymbol{O}","\\bfP":"\\boldsymbol{P}","\\bfQ":"\\boldsymbol{Q}","\\bfR":"\\boldsymbol{R}","\\bfS":"\\boldsymbol{S}","\\bfT":"\\boldsymbol{T}","\\bfU":"\\boldsymbol{U}","\\bfV":"\\boldsymbol{V}","\\bfW":"\\boldsymbol{W}","\\bfX":"\\boldsymbol{X}","\\bfY":"\\boldsymbol{Y}","\\bfZ":"\\boldsymbol{Z}"}});

In order for this to work, the computeQpOffDiagJacobian() function must be provided in the kernels that computes the required partial derivatives. To use all off diagonal blocks, you can use the following input file syntax:


full = true

Example Input File Syntax

[Preconditioning]
  active = 'SMP_jfnk'

  [./SMP_jfnk]
    type = SMP

    off_diag_row = 'forced'
    off_diag_column = 'diffused'

    #Preconditioned JFNK (default)
    solve_type = 'PJFNK'

    petsc_options_iname = '-pc_type'
    petsc_options_value = 'lu'
  [../]

  [./SMP_jfnk_full]
    type = SMP

    full = true

    #Preconditioned JFNK (default)
    solve_type = 'PJFNK'

    petsc_options_iname = '-pc_type'
    petsc_options_value = 'lu'
  [../]

  [./SMP_n]
    type = SMP

    off_diag_row = 'forced'
    off_diag_column = 'diffused'

    solve_type = 'NEWTON'

    petsc_options_iname = '-pc_type'
    petsc_options_value = 'lu'
  [../]
[]

Input Parameters

coupled_groupsList multiple space separated groups of comma separated variables. Off-diagonal jacobians will be generated for all pairs within a group.
C++ Type:std::vector<NonlinearVariableName>
Controllable:No
Description:List multiple space separated groups of comma separated variables. Off-diagonal jacobians will be generated for all pairs within a group.
fullFalseSet to true if you want the full set of couplings between variables simply for convenience so you don't have to set every off_diag_row and off_diag_column combination.
Default:False
C++ Type:bool
Controllable:No
Description:Set to true if you want the full set of couplings between variables simply for convenience so you don't have to set every off_diag_row and off_diag_column combination.
ksp_normunpreconditionedSets the norm that is used for convergence testing
Default:unpreconditioned
C++ Type:MooseEnum
Options:none, preconditioned, unpreconditioned, natural, default
Controllable:No
Description:Sets the norm that is used for convergence testing
off_diag_columnThe variable names for the off-diagonal columns you want to add into the matrix; they will be associated with an off-diagonal row from the same position in off_diag_row.
C++ Type:std::vector<NonlinearVariableName>
Controllable:No
Description:The variable names for the off-diagonal columns you want to add into the matrix; they will be associated with an off-diagonal row from the same position in off_diag_row.
off_diag_rowThe variable names for the off-diagonal rows you want to add into the matrix; they will be associated with an off-diagonal column from the same position in off_diag_column.
C++ Type:std::vector<NonlinearVariableName>
Controllable:No
Description:The variable names for the off-diagonal rows you want to add into the matrix; they will be associated with an off-diagonal column from the same position in off_diag_column.
pc_sidedefaultPreconditioning side
Default:default
C++ Type:MooseEnum
Options:left, right, symmetric, default
Controllable:No
Description:Preconditioning side
trust_my_couplingFalseWhether to trust my coupling even if the framework wants to be paranoid and create a full coupling matrix, which can happen when using global AD indexing for example.
Default:False
C++ Type:bool
Controllable:No
Description:Whether to trust my coupling even if the framework wants to be paranoid and create a full coupling matrix, which can happen when using global AD indexing for example.

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

mffd_typewpSpecifies the finite differencing type for Jacobian-free solve types. Note that the default is wp (for Walker and Pernice).
Default:wp
C++ Type:MooseEnum
Options:wp, ds
Controllable:No
Description:Specifies the finite differencing type for Jacobian-free solve types. Note that the default is wp (for Walker and Pernice).
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, -ksp_snes_ew, -snes_converged_reason, -snes_ksp, -snes_ksp_ew, -snes_linesearch_monitor, -snes_mf, -snes_mf_operator, -snes_monitor, -snes_test_display, -snes_view
Controllable:No
Description:Singleton PETSc options
petsc_options_inameNames of PETSc name/value pairs
C++ Type:MultiMooseEnum
Options:-ksp_atol, -ksp_gmres_restart, -ksp_max_it, -ksp_pc_side, -ksp_rtol, -ksp_type, -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, -snes_atol, -snes_linesearch_type, -snes_ls, -snes_max_it, -snes_rtol, -snes_divergence_tolerance, -snes_type, -sub_ksp_type, -sub_pc_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"
solve_typePJFNK: Preconditioned Jacobian-Free Newton Krylov JFNK: Jacobian-Free Newton Krylov NEWTON: Full Newton Solve FD: Use finite differences to compute Jacobian LINEAR: Solving a linear problem
C++ Type:MooseEnum
Options:PJFNK, JFNK, NEWTON, FD, LINEAR
Controllable:No
Description:PJFNK: Preconditioned Jacobian-Free Newton Krylov JFNK: Jacobian-Free Newton Krylov NEWTON: Full Newton Solve FD: Use finite differences to compute Jacobian LINEAR: Solving a linear problem

Petsc Parameters

(../moose/examples/ex11_prec/smp.i)

[Mesh]
  file = square.e
[]

[Variables]
  [./diffused]
    order = FIRST
    family = LAGRANGE
  [../]

  [./forced]
    order = FIRST
    family = LAGRANGE
  [../]
[]

# The Preconditioning block
[Preconditioning]
  active = 'SMP_jfnk'

  [./SMP_jfnk]
    type = SMP

    off_diag_row    = 'forced'
    off_diag_column = 'diffused'


  #Preconditioned JFNK (default)
  solve_type = 'PJFNK'


    petsc_options_iname = '-pc_type'
    petsc_options_value = 'lu'
  [../]

  [./SMP_jfnk_full]
    type = SMP

    full = true


  #Preconditioned JFNK (default)
  solve_type = 'PJFNK'


    petsc_options_iname = '-pc_type'
    petsc_options_value = 'lu'
  [../]

  [./SMP_n]
    type = SMP

    off_diag_row    = 'forced'
    off_diag_column = 'diffused'


    solve_type = 'NEWTON'

    petsc_options_iname = '-pc_type'
    petsc_options_value = 'lu'
  [../]
[]

[Kernels]
  [./diff_diffused]
    type = Diffusion
    variable = diffused
  [../]

  [./conv_forced]
    type = CoupledForce
    variable = forced
    v = diffused
  [../]

  [./diff_forced]
    type = Diffusion
    variable = forced
  [../]
[]

[BCs]
  #Note we have active on and neglect the right_forced BC
  active = 'left_diffused right_diffused left_forced'
  [./left_diffused]
    type = DirichletBC
    variable = diffused
    boundary = 1
    value = 0
  [../]

  [./right_diffused]
    type = DirichletBC
    variable = diffused
    boundary = 2
    value = 100
  [../]

  [./left_forced]
    type = DirichletBC
    variable = forced
    boundary = 1
    value = 0
  [../]

  [./right_forced]
    type = DirichletBC
    variable = forced
    boundary = 2
    value = 0
  [../]
[]

[Executioner]
  type = Steady
[]

[Outputs]
  exodus = true
[]