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().

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)

Bases: object

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

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).

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)

Sets a solver option.

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.