The solver is the part of the code that is responsible for solving a given system of nonlinear 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 nonlinear equation system that must be solved to advance the system in time.
Page overview
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 


The linearly implicit solver 

The Newton solver 
The linearly implicit solver linearizes the equation system in time around the next time, \(t_{l+1}\), as
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}\):
The nonlinear 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
Setting \(\boldsymbol{F}(\boldsymbol{x})\) and solving for \(\boldsymbol{x}\) then yields the iterative Newton’s method:
where \(\boldsymbol{x}^{(l+1)}_i\) denotes the \(i\) th approximation to the solution \(\boldsymbol{x}(t_{l+1})\).
As explained above, the nonlinear 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
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=1e8)
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()
.
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 nonlinear, while the nonlinear 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 nonlinear systems have proven difficult for the linearly implicit solver.
Tip
When in doubt, we recommend using the Newton (nonlinear) 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 nonlinear? If it is, the Newton solver will often be required to avoid numerical instabilities crashing the linearly implicit solver.
Does the system involve rapidly timevarying parameters? If so, the linearly implicit solver may require exceedingly short time steps which could partly be avoided by the more accurate Newton solver.
Both the linearly implicit and nonlinear solvers involve solutions of sets of linear equations. Hence, both solvers utilize thirdparty linear equation solvers available through the PETSc library. Both solvers currently support the use of four different LU factorization algorithms, namely
Name 
Description 


The LU direct solver available in PETSc. 

The parallel direct solver PARDISO, part of Intel’s Math Kernel Library (MKL). 

The MUMPS parallel sparse direct solver. 

The SuperLU direct LU 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 builtin 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 builtin 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.
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 nonlinear Newton solver.
The following example will cause DREAM to print information about the jacobian matrix of the nonlinear 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 nonnegative 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.
The debug options available for the linearly implicit solver are
Option 
Description 


Print information about the linear operator matrix after it has been built. 

Save the linear operator matrix using the PETSc MATLAB binary viewer. 

Save the righthandside vector to a 
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.
Option 
Description 


Print information about the jacobian matrix after it has been built. 

Save the rescaled jacobian matrix/residual vector (rescaled before solution to improve condition number). 

Save the jacobian matrix using the PETSc MATLAB binary viewer. 

Save the solution vector \(\boldsymbol{x}_{i+1}^{(l+1)}\boldsymbol{x}_i^{(l+1)}\) to a 

Approximate and save the jacobian matrix numerically. The matrix is saved using the PETSc MATLAB binary viewer. 

Save the residual vector \(\boldsymbol{F}(\boldsymbol{x})\) to a 

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)
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.
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 righthand side vector to a .mat
file.
NONLINEAR 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.