(sect:linearsolver)=
# Linear Solvers

Many solution methods rely on reliable and efficient linear solvers.
However, there are features in models that may either impede convergence
or degrade accuracy. The Helmholtz linear solver is discussed separately
in Section {numref}`%s <sect:helmholtzlinearsolver>`. In this section,
common issues are tabulated and an example with before and after
configurations is reviewed.

1.  Some problems occur only for models with lots of constraint
    equations, due to large surfaces that are tied together (e.g. one
    large sideset constrained to another with many nodes). A way to
    confirm that this is the issue is the check whether the problem is
    mitigated if tied contact over large surfaces is turned off.

2.  Decreasing the time step (e.g. halving) can mitigate convergence
    issues.

3.  Suppose there are accuracy issues. Note that the tolerance on the
    residual is always larger than the uncertainty in the solution
    vector. A linear system has a condition number, which is always
    greater than $1$. The uncertainty in the solution vector is the
    product of the condition number and the tolerance on the residual.

4.  There are alternative to `GDSW`. Sierra/SD provides serial sparse
    linear solvers, `sparsepak` for symmetric positive definite systems,
    and `SuperLU` for other systems. In addition, `Pardiso` is a
    general-purpose sparse solver that is available on Intel platforms.
    These solvers are at least as robust as the iterative methods. It
    can be enlightening to try to use the appropriate serial sparse
    linear solver as problem size permits.

Consider, for example, the following user provided configuration of the
GDSW linear solver.

:::{input}
GDSW
  prt_summary = 3
  solver_tol = 1.0e-5
  max_iter = 5000
  orthog = 200
  overlap = 1
  diag_scaling = diagonal
  scale_option = 1
END
:::

The options are generally intuitive. If the solver converges, and
accuracy issues arise, then trying a smaller `solver_tol`, and a larger
`max_iter` is recommended. If the solver diverges, then trying a larger
`solver_tol` or a larger `max_iter` is recommended. A larger `orthog` is
also recommended. However, there are memory usage limitations. If there
is an immediate error that could be related to running out of memory,
then try a smaller value of `orthog` or use more processors. See the
discussion of reducing memory usage in the training documents for
details.

There is a hidden constraint on these options. With some Krylov methods,
e.g. the default of `krylov_method` = 1 (GMRES), it turns out that
`orthog` $\ge$ `max_iter`. For this reason, when divergence is a
problem, users often switch to `gmresClassic`, which allows `orthog` $<$
`max_iter`.

In this example, `overlap = 1` is a small value for overlap. If you are
running out of memory with a higher value, then this might be a great
idea. If the linear solver is diverging, you might try a larger value
(the default is 2).

The `diag_scaling = diagonal` option can be used either to find a
convergent solver, or to find a more accurate solver. On the other hand,
there are cases in which selecting the option decreases accurate.

In this case study, the user ultimately changed the GDSW configuration
to the following to address convergence issues.

:::{input}
GDSW
    solver_tol = 1e-12
    overlap = 2
    num_vectors_keep = 0
    orthog = 4000
    max_iter = 4000
    krylov_method = gmresClassic
END
:::

The option `num_vectors_keep` can only be used with the classic version
of GMRES (krylov_method gmresClassic). The parameter orthog controls how
many search direction are stored. We store search directions to make the
linear solver faster. More is generally better. The point to understand
is which search directions are stored. In this example, the first $4000$
search directions are stored. On later solves, the first
` num_vectors_keep` are saved and recycled. The default value of
` num_vectors_keep` is `orthog`/2. In this case the solution has changed
significantly and you don't want to use any of the old search
directions. ` num_vectors_keep = 0` tells GDSW to start afresh and
remove all search directions every time the maximum is reached. Thus,
the benefits of recycling are still retained, but the entire search
space is periodically purged of older search directions.

:::{toctree}
:hidden:
:maxdepth: 2

Accuracy
Helmholtz
:::
