Newton solver

The Newton solver and the related parameters. This describes the reference implementation Dumux::NewtonSolver. More...

Description

The Newton solver and the related parameters. This describes the reference implementation Dumux::NewtonSolver.

Newton's method

The following describes Newton's method and the reference implementation Dumux::NewtonSolver of a multi-dimensional Newton solver to solve a non-linear system of equations in DuMu^x. For solving nonlinear equations with only one unknown, DuMu^x also provides Dumux::findScalarRootNewton and Dumux::findScalarRootBrent.

We assume that the equations (in DuMu^x usually the discretized partial differential equations) are formulated as a residual equation:

$$\begin{array}{r}F(\mathbf{\text{u}})=\mathbf{\text{r}}=\mathbf{\text{0}}\end{array}$$

with the residual operator $F:{\mathbb{R}}^n\to \mathbb{R}$, $n$ is the number of equations, $\mathbf{\text{u}}\in {\mathbb{R}}^n$ is the vector of unknowns, and $\mathbf{\text{r}}\in {\mathbb{R}}^n$ is the residual vector.

The goal of Newton's method is to find a vector $\mathbf{\text{u}}$ such that the residual equation is fulfilled (at least to a good enough approximation $F(\mathbf{\text{u}})\approx \mathbf{\text{0}}$ as we will discuss below).

Newton's method is an iterative algorithm: it generates a sequence of approximations, ${\mathbf{\text{u}}}_k$, that converges (under some circumstances) to the vector ${\mathbf{\text{u}}}_{\star }$, that satisfies $F({\mathbf{\text{u}}}_{\star })=0$. We start with an initial guess ${\mathbf{\text{u}}}_0$ and calculate the initial residual ${\mathbf{\text{r}}}_0=\mathbf{\text{F}}({\mathbf{\text{u}}}_0)$. Then, we calculate the derivative of the residual operator with respect to the unknowns, the Jacobian matrix, $J:{\mathbb{R}}^n\to {\mathbb{R}}^n$ evaluated at ${\mathbf{\text{u}}}_k$ (in the first step $k=0$):

$$\begin{array}{r}J({\mathbf{\text{u}}}_k)=\frac{\partial F(\mathbf{\text{u}})}{\partial \mathbf{\text{u}}}{|}_{\mathbf{\text{u}}={\mathbf{\text{u}}}_k}\end{array}$$

Note

Computing the residual operator $F$ and the Jacobian matrix $J$ is delegated to the assembler implementation (see e.g. Dumux::FVAssembler). Typically with DuMu^x, the Jacobian matrix is approximated by finite differences but the implementation of Newton's method in Dumux::NewtonSolver is independent of the way the Jacobian is calculated. The C++ type of the Assembler is a template parameter of Dumux::NewtonSolver and and instance of the assembler must be passed to the Newton solver in its constructor Dumux::NewtonSolver::NewtonSolver.

This allows to compute a tangent direction of residual surface in n-dimensional space as

$$\begin{array}{r}\mathrm{\Delta }\mathbf{\text{u}}={\mathbf{\text{u}}}_k-{\mathbf{\text{u}}}_{k+1}=J({\mathbf{\text{u}}}_k{)}^{-1}\cdot {\mathbf{\text{r}}}_k\end{array}$$

Note

Evaluating $\mathrm{\Delta }\mathbf{\text{u}}$ is equivalent to solving the linear system

$$\begin{array}{r}J({\mathbf{\text{u}}}_k)\mathrm{\Delta }\mathbf{\text{u}}={\mathbf{\text{r}}}_k\end{array}$$

This task is delegated to a linear solver (see e.g. Dumux::LinearSolver) which is passed the Jacobian matrix and the residual vector computed by the assembler earlier. The C++ type of the LinearSolver is a template parameter of Dumux::NewtonSolver and and instance of the linear solver must be passed to the Newton solver in its constructor Dumux::NewtonSolver::NewtonSolver.

With the tangent direction $\mathrm{\Delta }\mathbf{\text{u}}$ (also called the "update" vector), we update our current approximation ${\mathbf{\text{u}}}_k$ to get a new approximation ${\mathbf{\text{u}}}_{k+1}$:

$$\begin{array}{r}{\mathbf{\text{u}}}_{k+1}={\mathbf{\text{u}}}_k-\mathrm{\Delta }\mathbf{\text{u}}\end{array}$$

Note

The update can be customized by a user-defined update strategy or by choosing from the implemented strategies (see below for more details).

This is one Newton step (or iteration). We repeat the steps of calculating the residual and the Jacobian matrix, solve the linear system, and update the approximation until we reach a good enough solution (see below for the available termination criteria).

Note

A good initial guess is crucial to get fast convergence of Newton's method. For time-dependent problems, the initial guess ${\mathbf{\text{u}}}_0$ is usually the solution of the previous time step. The initial guess is passed to the Newton solver in the method Dumux::NewtonSolver::apply and Dumux::NewtonSolver::solve.

When choosing a zero initial guess you must be aware that the assembler (based on finite difference approximations) may have difficulties to find a good epsilon to compute the finite difference approximation. In this case, you may need to help the assembler by providing an estimate of the problem-typical magnitude of the unknowns. See Dumux::FVAssembler and the parameters Assembly.NumericDifference.PriVarMagnitude and Assembly.NumericDifference.BaseEpsilon.

The Dumux::NewtonSolver reference implementation of Newton's method has several features discussed below that can be controlled by user-set parameters (e.g. via a parameter input file or the command line, see Runtime Parameters).

Update strategy

Here we describe strategies that differ from the default update strategy, that is:

$$\begin{array}{r}{\mathbf{\text{u}}}_{k+1}={\mathbf{\text{u}}}_k-\mathrm{\Delta }\mathbf{\text{u}}\end{array}$$

The different strategies can be enabled via runtime parameters as described below.

Line search

The idea of line search is to find an optimal scaling factor $\alpha \in \mathbb{R}$ for the update vector $\mathrm{\Delta }\mathbf{\text{u}}$

$$\begin{array}{r}{\mathbf{\text{u}}}_{k+1}={\mathbf{\text{u}}}_k-\alpha \mathrm{\Delta }\mathbf{\text{u}}\end{array}$$

such that the residual norm evaluated with the new iterate ${\mathbf{\text{u}}}_{k+1}$ that is $||r_{k+1}||$ is minimized. That means using line search the following optimization problem is solved: find $\alpha$ such that

$$\begin{array}{r}\alpha =\arg \underset{\alpha }{min}||F({\mathbf{\text{u}}}_{k+1}-\alpha \mathrm{\Delta }\mathbf{\text{u}})|{|}_2^2\end{array}$$

The default implementation of Dumux::NewtonSolver uses a simple backtracking line search algorithm to approximately solve the above optimization problem. The algorithm is based on the following steps:

1. Set the initial step size $\alpha =1$.
2. Compute the residual norm $||r_{k+1}||=||F({\mathbf{\text{u}}}_{k+1}-\alpha \mathrm{\Delta }\mathbf{\text{u}})|{|}_2^2$.
3. If the residual norm is smaller than $||r_k||$ (the norm of the previous residual), $\alpha$ is accepted and the algorithm terminates. Otherwise, alpha is reduced by a factor Newton.LineSearchReductionFactor (default value $0.5$) and we continue with step 2 unless $\alpha$ is smaller or equal the Newton.LineSearchMinRelaxationFactor (default value $0.125$) in which we accept the step size, even if the residual norm increased.

Note

To customize the line search strategy, create a new class derived from Dumux::NewtonSolver and override the private virtual method Dumux::NewtonSolver::lineSearchUpdate_. Activate line search update by setting the parameter Newton.UseLineSearch = true.

Chopped update

For instance in problems where the residual has multiple roots, it can be beneficial to restrict the search space to a certain region in the parameter space. This can be achieved by a chopped update strategy, where the solution vector is projected back to a feasible region after each Newton step. A chopped update strategy can be enabled by setting Newton.EnableChop = true. As such strategies are usually problem-dependent, there is no default implementation. The parameter will merely active the hook Dumux::NewtonSolver::choppedUpdate_ which can be overridden in a derived class, see Dumux::RichardsNewtonSolver for an example implementation using this strategy.

Note

To customize the chop update strategy, create a new class derived from Dumux::NewtonSolver and override the private virtual method Dumux::NewtonSolver::choppedUpdate_. Activate the chopped update strategy by setting the parameter Newton.EnableChop = true. The line search and chopped update are mutually exclusive, i.e. if Newton.UseLineSearch is set to true, the Newton.EnableChop parameter will be ignored.

Termination criteria

Newton's method finds an approximate solution to the residual equation $F({\mathbf{\text{u}}}_k)\approx \mathbf{\text{0}}$. But when is the solution "good enough"? A good criterion makes sure that the residual is small enough but should also be easy to compute. Dumux::NewtonSolver provides several termination criteria.

The default is the relative shift criterion in addition to bounds on the minimum and maximum number of iterations. The solver never terminates if $k$ is smaller than Newton.MinSteps (default value $2$). The solver always terminates if $k$ is larger than Newton.MaxSteps (default value $18$). Dumux::NewtonSolver::apply returns a boolean value indicating whether the solver converged or not; Dumux::NewtonSolver::solve throws Dumux::NumericalProblem is the solver did not converge.

Relative shift criterion

This criterion is enabled per default. The idea is to terminate the algorithm when the solution vector does not change significantly between iterations. This criterion does not check the actual residual but is useful in practice and fast to compute. We compute the maximum difference of any degree of freedom (entry of the vector ${\mathbf{\text{u}}}_{k+1}$) between the current and the previous iteration. This uses the criterion

$$\begin{array}{r}{\text{max}}_{i=1,\dots ,n}\left(\frac{|u_{k+1,i}-u_{k,i}|}{max\{1.0,|u_{k+1,i}+u_{k,i}|\cdot 0.5\}}\right)<{ϵ}_{\mathrm{shift}}\end{array}$$

where $u_{k,i}$ is the $i$-th entry of the vector ${\mathbf{\text{u}}}_k$ and ${ϵ}_{\mathrm{shift}}$ is a user-defined parameter Newton.MaxRelativeShift (default value ${10}^{-8}$). The relative shift criterion can be disabled by setting Newton.EnableShiftCriterion = false. The relative shift criterion can be combined with a residual criterion (see below) by enabling a residual criterion without disabling the shift criterion. By default the Newton solver terminates when either of the criteria are fulfilled. If you want to require both criteria to be fulfilled, you can do so by setting Newton.SatisfyResidualAndShiftCriterion = true.

Relative residual norm criterion

This criterion is disabled by default. It can be enabled by setting Newton.EnableResidualCriterion = true. The idea is to terminate the algorithm when the norm of the residual vector is small in comparison to the initial residual norm. The relative residual norm criterion checks if

$$\begin{array}{r}\frac{||{\mathbf{\text{r}}}_{k+1}|{|}_2^2}{||{\mathbf{\text{r}}}_0|{|}_2^2}<{ϵ}_{\mathrm{r},\mathrm{rel}}^2,\end{array}$$

where $||{\mathbf{\text{r}}}_{k+1}||$ is the norm of the residual vector at iteration $k+1$ and $||{\mathbf{\text{r}}}_0||$ is the norm of the residual vector at iteration $0$. The parameter Newton.ResidualReduction (default value ${10}^{-5}$) defines the threshold ${ϵ}_{\mathrm{r},\mathrm{rel}}$.

Absolute residual norm criterion

This criterion is disabled by default. The idea is to terminate the algorithm when norm of the residual vector is small enough in comparison with a given absolute threshold. This requires knowledge about the problem because the magnitude of the residual vector has to be estimated. The criterion is enabled by setting Newton.EnableAbsoluteResidualCriterion = true and Newton.EnableResidualCriterion = true. It checks if

$$\begin{array}{r}||{\mathbf{\text{r}}}_{k+1}|{|}_2^2<{ϵ}_{\mathrm{r},\mathrm{abs}}^2,\end{array}$$

where $||{\mathbf{\text{r}}}_{k+1}||$ is the norm of the residual vector at iteration $k+1$ and ${ϵ}_{\mathrm{r},\mathrm{abs}}$ is a user-defined parameter Newton.MaxAbsoluteResidual (default value ${10}^{-5}$). The absolute residual norm criterion can be combined with a relative shift criterion by enabling both criteria.

Adaptive time stepping

The extended interface Dumux::NewtonSolver::solve with time loop in the signature provides an automated adaptive time stepping strategy. The idea is retry the entire algorithm with a lower the time step size if the Newton solver fails to converge. This is repeated until the Newton solver either converges or the maximum number of allowed retry steps is reached.

The factor by which the time step size is reduced after each non-converged try is set by the parameter Newton.RetryTimeStepReductionFactor (default value $0.5$). The maximum number of retries is set by the parameter Newton.MaxTimeStepDivisions (default value $10$). If the maximum number of retries is reached, the a Dumux::NumericalProblem is thrown.

Finally, Dumux::NewtonSolver::suggestTimeStepSize suggest a new time-step size based on the current time-step size and the convergence history of the Newton solver. In particular, the implemented heuristic determines how much the time step size is adjusted depends on how far away the Newton is from converging within $m$ iterations, where $m$ is set by the parameter Newton.TargetSteps (default value $10$).

Partial reassembly (Quasi-Newton method)

Technically, when using a finite difference approximation of the Jacobian, we use a so-called quasi-Newton method for which the Jacobian matrix is replaced by an approximation. Another strategy is useful for simulations where the Jacobian matrix is expensive to compute but not all degrees of freedom are expected to change significantly between iterations (for example if changes only happen in parts of the simulation domain).

If supported by the assembler (e.g. Dumux::FVAssembler), the Newton solver can tell the assembler to only recompute those entries of the Jacobian matrix that are associated with primary variables that have changed significantly. This strategy is enabled via the parameter Newton.EnablePartialReassembly (default value is false).

The threshold for the partial reassembly is computed can also be controlled via parameters. The threshold ${\theta }_r$ for each entry is computed as

$$\begin{array}{r}{\theta }_r=max({\theta }_{\min },min({\theta }_{\max },\omega {ϵ}_{\text{shift}}))\end{array}$$

where ${\theta }_{\min }$ is set by Newton.ReassemblyMinThreshold (default value ${10}^{-1}\mathrm{\Delta }{u}_{\mathrm{shift}}$) and ${\theta }_{\max }$ is set by Newton.ReassemblyMaxThreshold (default value ${10}^{-2}\mathrm{\Delta }{u}_{\mathrm{shift}}$), $\omega$ is a weight (default value ${10}^{-3}$), and ${ϵ}_{\text{shift}}$ is the maximum relative shift threshold parameter (Newton.MaxRelativeShift) described above.

Primary variable switching

Switching primary variables during Newton's method is a form of nonlinear preconditioning. We solve an alternative problem with different primary variables where the hope that the root is easier to find. This can be, for instance, because the finite difference approximation of the Jacobian is numerically better conditioned in the new variables. The strategy depends on the problem at hand. For an example, see for instance 2pnc for a model with activated primary variable switch.

---

🛠 Edit above doc on GitLab

Files
file	multidomain/newtonconvergencewriter.hh
	This class provides the infrastructure to write the convergence behaviour of the Newton method for multidomain simulations into a VTK file.
file	multidomain/newtonsolver.hh
file	nonlinear/newtonconvergencewriter.hh
	This class provides the infrastructure to write the convergence behaviour of the newton method into a VTK file.
file	nonlinear/newtonsolver.hh
	Reference implementation of a Newton solver.
file	staggerednewtonconvergencewriter.hh
	This class provides the infrastructure to write the convergence behaviour of the newton method for the staggered discretization scheme into a VTK file.

Classes
class	Dumux::MultiDomainNewtonConvergenceWriter< MDTraits >
	Writes the intermediate solutions for every Newton iteration. More...
class	Dumux::MultiDomainNewtonSolver< Assembler, LinearSolver, CouplingManager, Reassembler, Comm >
	Newton solver for coupled problems. More...
struct	Dumux::ConvergenceWriterInterface< SolutionVector, ResidualVector >
	A convergence writer interface Provide an interface that show the minimal requirements a convergence write passed to a newton method has to fulfill. More...
class	Dumux::NewtonConvergenceWriter< GridGeometry, SolutionVector, ResidualVector >
	Writes the intermediate solutions for every Newton iteration. More...
class	Dumux::NewtonSolver< Assembler, LinearSolver, Reassembler, Comm >
	An implementation of a Newton solver. The comprehensive documentation is in Newton solver, providing more details about the algorithm and the related parameters. More...
class	Dumux::StaggeredNewtonConvergenceWriter< GridGeometry, SolutionVector, ResidualVector >
	Writes the intermediate solutions for every Newton iteration (for staggered grid scheme) More...