Solver

The solver is the part of the code that is responsible for solving a given system of non-linear equations, and thereby computing the system state in the next time step. DREAM provides two different solvers, namely a linearly implicit solver as well as a regular Newton solver. In the discussion that follows, we let \(\boldsymbol{F}(t,\boldsymbol{x(t)}) = 0\) represent the non-linear equation system that must be solved to advance the system in time.

Which solver to use is determined by calling setType() on the solver object:

import DREAM.Settings.Solver as Solver

ds = DREAMSettings()
...
ds.solver.setType(Solver.NONLINEAR)

The available solver types are

Name	Description
LINEAR_IMPLICIT	The linearly implicit solver
NONLINEAR	The Newton solver

Linearly implicit solver

The linearly implicit solver linearizes the equation system in time around the next time, \(t_{l+1}\), as

\[\boldsymbol{F}\left( t_{l+1}, \boldsymbol{x}^{(l+1)} \right) \approx \mathsf{M}\left( t_l, \boldsymbol{x}^{(l)} \right) \boldsymbol{x}^{(l+1)} + \boldsymbol{S}\left( t_l, \boldsymbol{x}^{(l)} \right),\]

where \(\boldsymbol{x}^{(l)}\equiv\boldsymbol{x}(t_l)\), \(\mathsf{M}\) is a matrix and \(\boldsymbol{S}\) is a vector. Setting both sides to zero allows us to obtain the solution in the next time step by inverting the matrix \(\mathsf{M}\):

\[\boldsymbol{x}^{(l+1)} = \mathsf{M}^{-1}\left( \boldsymbol{x}^{(l)} \right) \boldsymbol{S}\left( \boldsymbol{x}^{(l)} \right)\]

Non-linear solver

The non-linear solver linearizes the equation system around a point \(\tilde{\boldsymbol{x}}\) in the vicinity of the true solution (rather than around a particular time \(t\), as the linearly implicit solver), yielding

\[\boldsymbol{F}\left(\boldsymbol{x}\right) \approx \boldsymbol{F}\left(\tilde{\boldsymbol{x}}\right) + \mathsf{J}(\tilde{\boldsymbol{x}})\left( \boldsymbol{x} - \tilde{\boldsymbol{x}} \right)\]

Setting \(\boldsymbol{F}(\boldsymbol{x})\) and solving for \(\boldsymbol{x}\) then yields the iterative Newton’s method:

\[\boldsymbol{x}^{(l+1)}_{i+1} = \boldsymbol{x}^{(l+1)}_i - \mathsf{J}^{-1}\left( \boldsymbol{x}^{(l+1)}_i \right) \boldsymbol{F}\left( \boldsymbol{x}^{(l+1)}_i \right),\]

where \(\boldsymbol{x}^{(l+1)}_i\) denotes the \(i\) th approximation to the solution \(\boldsymbol{x}(t_{l+1})\).

Tolerance settings

As explained above, the non-linear solver uses a Newton method to iteratively solve the equation system. To determine when a solution is converged, the solver checks the obtained solution of each individual unknown quantity in every iteration and requires that

\[\left\lVert \boldsymbol{y}^{(n)}_{i+1} - \boldsymbol{y}^{(n)}_{i} \right\rVert \leq \epsilon^{(n)}_{\rm abs} + \left\lVert \boldsymbol{y}^{(n)}_{i+1}\right\rVert\epsilon^{(n)}_{\rm rel}\]

where \(\boldsymbol{y}^{(n)}_{i+1}\) denotes a subset of the elements in the full solution vector \(\boldsymbol{x}_{i+1}\), corresponding to the discretized quantity denoted with index \(n\), and \(\epsilon^{(n)}_{\rm abs}\) and \(\epsilon^{(n)}_{\rm rel}\) are the absolute and relative tolerances for the quantity respectively.

The tolerances \(\epsilon^{(n)}_{\rm abs}\) and \(\epsilon^{(n)}_{\rm rel}\) can be specified for each unknown quantity of the equation system using the method tolerance.set() of the Solver object. The method takes as arguments the name of the quantity to set tolerances for, as well as the absolute and/or relative tolerance values to apply:

ds = DREAMSettings()
...
ds.solver.tolerance.set('n_re', abstol=1e6, reltol=1e-8)

Note that it is possible to specify just one of abstol and reltol, in which case the value of the tolerance not specified remains unchanged.

All quantities have default tolerances which are set in the kernel depending on the typical scales of the quantities. These tolerances are fine in many cases but may sometimes need to be adjusted.

Warning

Specifying just one of the absolute or relative tolerance in the interface will also override any defaults set by the kernel for the other tolerance. As such, it is usually a good idea to specify both the absolute and relative tolerances when calling tolerance.set().

Monitoring the residual

As noted above, the Newton solver uses a condition on the change in the solution from iteration to iteration in order to determine if it has converged. In the vast majority of cases, this works well, but in a few outlier cases (when the jacobian is particularly ill-conditioned) we have found that the solver stops iterating even though the found solution does not satisfy the system of equations. To avoid such unphysical solutions, it is possible to let DREAM also monitor the residual and issue a warning when the solution is diverged. This functionality is enabled by default, and can be switched on/off using

ds.solver.setCheckResidual(True)   # Check the residual in each iteration
ds.solver.setCheckResidual(False)  # Do NOT check the residual in each iteration

If enabled, the solver will also output information about the convergence in the output file. This information consists of a flag for each time step and non-trivial unknown quantity in the equation system, indicating whether the residual was close to zero, as well as the maximum 2-norm of all elements in the unknown vector, normalized to the scale length used in DREAM to determine if the residual is close to zero or not.

Warning

Sadly, if the residual is found to be diverging, there is no general technique for correcting this problem. One will have to try to figure out exactly what is causing the divergence manually and try to work around the issue.

Storing iteration progress in output

In each Newton iteration, the solver will estimate how well converged the solution is, and decide to either continue iteration or stop. The progress of the solution can be printed to stdout by calling ds.solver.setVerbose(True), but the same information can also be saved to the output file. To save the iteration progress to the output file, just add the following line to your settings:

ds.solver.setSaveConvergenceInfo(True)

By default, this option is disabled.

Which solver should I use?

The difference between the two solvers is primarily that the linearly implicit solver is fast, but potentially inaccurate and prone to crashing if the time evolution is rapid or highly non-linear, while the non-linear solver is generally slower but very robust and guarantees a certain accuracy in the solution. In principle it should be possible to use either solver to solve any given problem, although certain non-linear systems have proven difficult for the linearly implicit solver.

Tip

When in doubt, we recommend using the Newton (non-linear) solver. While it is generally slower than the linearly implicit solver, we often find that its robustness makes up for this.

In practice, one should therefore consider a number of questions when deciding which solver to use:

• Is the equation system non-linear? If it is, the Newton solver will often be required to avoid numerical instabilities crashing the linearly implicit solver.

• Does the system involve rapidly time-varying parameters? If so, the linearly implicit solver may require exceedingly short time steps which could partly be avoided by the more accurate Newton solver.

Linear solvers

Both the linearly implicit and non-linear solvers involve solutions of sets of linear equations. Hence, both solvers utilize third-party linear equation solvers available through the PETSc library. Both solvers currently support the use of four different LU factorization algorithms, namely

Name	Description
LINEAR_SOLVER_LU	The LU direct solver available in PETSc.
LINEAR_SOLVER_MKL	The parallel direct solver PARDISO, part of Intel’s Math Kernel Library (MKL).
LINEAR_SOLVER_MUMPS	The MUMPS parallel sparse direct solver.
LINEAR_SOLVER_SUPERLU	The SuperLU direct LU solver.

Backup linear solver

Some linear solvers are less robust than others. The less robust solvers however usually compensate for this by being orders of magnitudes faster. To obtain the best of both worlds, DREAM provides the ability to specify a “backup” linear solver which can take over inversions if the ordinary linear solver fails for some reason. In this way, we can use a fast solver (such as MKL Pardiso) to advance for as many time steps as possible, while using a more robust solver (typically PETSc’s built-in LU solver) when the ordinary solver is no longer capable of inverting the system of equations.

The backup linear solver is specified using the setBackupSolver() as follows:

import DREAM.Settings.Solver as Solver

ds = DREAMSettings()
...
# Use PETSc's built-in LU solver if the main solver fails.
ds.solver.setBackupSolver(Solver.LINEAR_SOLVER_LU)

Note

…that the backup linear solver may not be the same as the main linear solver.

Debug settings

A number of options are available which can aid in debugging numerical issues with the solver. The options are enabled/disabled using the setDebug() method for both solvers, and they apply to one or all time steps, and one or all iterations of the non-linear Newton solver.

The following example will cause DREAM to print information about the jacobian matrix of the non-linear Newton solver in every iteration of the first time step:

ds = DREAMSettings()
...
ds.solver.setDebug(printjacobianinfo=True, timestep=1, iteration=0)

The timestep and iteration options should be non-negative integers. If the integer is 0, the requested actions will be taken in every time step and/or iteration. Otherwise, the actions are taken in the specifed time step and/or iteration.

Linearly implicit solver

The debug options available for the linearly implicit solver are

Option	Description
printmatrixinfo	Print information about the linear operator matrix after it has been built.
savematrix	Save the linear operator matrix using the PETSc MATLAB binary viewer.
saverhs	Save the right-hand-side vector to a .mat file.

Example usage:

ds = DREAMSettings()
...
ds.solver.setDebug(printmatrixinfo=True, savematrix=True, saverhs=True, timestep=2)

Note

The iteration parameter of setDebug() has no effect when the linearly implicit solver is used.

Non-linear solver

Option	Description
printjacobianinfo	Print information about the jacobian matrix after it has been built.
rescaled	Save the rescaled jacobian matrix/residual vector (rescaled before solution to improve condition number).
savejacobian	Save the jacobian matrix using the PETSc MATLAB binary viewer.
savesolution	Save the solution vector \(\boldsymbol{x}_{i+1}^{(l+1)}-\boldsymbol{x}_i^{(l+1)}\) to a .mat file.
savenumericaljacobian	Approximate and save the jacobian matrix numerically. The matrix is saved using the PETSc MATLAB binary viewer.
saveresidual	Save the residual vector \(\boldsymbol{F}(\boldsymbol{x})\) to a .mat file.
savesystem	Generate a regular DREAM output file with the data in the last time step populated from the most recent Newton iteration.

Example usage:

ds = DREAMSettings()
...
ds.solver.setDebug(printjacobianinfo=True, savejacobian=True,
                   savenumericaljacobian=True, saveresidual=True,
                   savesystem=True, rescaled=True, timestep=1, iteration=4)

Class documentation

class DREAM.Settings.Solver.Solver(ttype=1, linsolv=1, maxiter=100, verbose=False, checkresidual=True)

Bases: object

__init__(ttype=1, linsolv=1, maxiter=100, verbose=False, checkresidual=True)

Constructor.

fromdict(data)

Load settings from the given dictionary.

setBackupSolver(backup)

Set the backup linear solver to use in case the main linear solver fails. Set to None to disable (default).

setCheckResidual(checkresidual)

If ‘True’, ensures that the residual is close to zero at the end of non-linear iteration in each time step.

setDebug(printmatrixinfo=False, printjacobianinfo=False, savejacobian=False, savesolution=False, savematrix=False, savenumericaljacobian=False, saverhs=False, saveresidual=False, savesystem=False, rescaled=False, timestep=0, iteration=1)

Enable output of debug information.

Parameters

• timestep (int) – Index of time step to generate debug info for. If 0, debug info is generated in every (iteration of every) time step.

• savesystem (int) – Save full equation system as a DREAMOutput file in the most recent iteration/time step.

LINEAR SOLVER :param bool printmatrixinfo: If True, calls PrintInfo() on the linear operator matrix. :param bool savematrix: If True, saves the linear operator matrix using a PETSc viewer. :param bool saverhs: If True, saves the right-hand side vector to a .mat file.

NON-LINEAR SOLVER :param bool printjacobianinfo: If True, calls PrintInfo() on the jacobian matrix. :param bool savejacobian: If True, saves the jacobian matrix using a PETSc viewer. :param bool savesolution: If True, saves the solution vector to a .mat file. :param bool savenumericaljacobian: If True, evaluates the jacobian matrix numerically and saves it using a PETSc viewer. :param bool saveresidual: If True, saves the residual vector to a .mat file. :param bool rescaled: If True, saves the rescaled versions of the jacobian/solution/residual. :param int iteration: Index of iteration to save debug info for. If 0, saves in all iterations. If timestep is 0, this parameter is always ignored.

setLinearSolver(linsolv)

Set the linear solver to use.

setMaxIterations(maxiter)

Set maximum number of allowed nonlinear iterations.

setOption(linsolv=None, maxiter=None, verbose=None, checkresidual=None, saveconvergenceinfo=None)

Sets a solver option.

setSaveConvergenceInfo(saveconvergenceinfo)

If ‘True’, saves information about the non-linear convergence to the output file (|x| and |dx| norms in each time step and iteration).

setTolerance(reltol)

Set relative tolerance for nonlinear solve.

setType(ttype)

Specifies which type of solver to use (either LINEAR_IMPLICIT or NONLINEAR).

setVerbose(verbose)

If ‘True’, generates excessive output during nonlinear solve.

todict(verify=True)

Returns a Python dictionary containing all settings of this Solver object.

verifyLinearSolverSettings()

Verifies the settings for the linear solver (which is used by both the ‘LINEAR_IMPLICIT’ and ‘NONLINEAR’ solvers).

verifySettings()

Verifies that the settings of this object are consistent.