Reference documentation for deal.II version 9.4.0
SUNDIALS::KINSOL< VectorType > Class Template Reference

#include <deal.II/sundials/kinsol.h>

## Public Member Functions

KINSOL (const AdditionalData &data, const MPI_Comm &mpi_comm)

~KINSOL ()

unsigned int solve (VectorType &initial_guess_and_solution)

## Static Public Member Functions

static ::ExceptionBaseExcKINSOLError (int arg1)

## Public Attributes

std::function< void(VectorType &)> reinit_vector

std::function< int(const VectorType &src, VectorType &dst)> residual

std::function< int(const VectorType &src, VectorType &dst)> iteration_function

std::function< int(const VectorType &current_u, const VectorType &current_f)> setup_jacobian

std::function< int(const VectorType &ycur, const VectorType &fcur, const VectorType &rhs, VectorType &dst)> solve_jacobian_system

std::function< int(const VectorType &rhs, VectorType &dst, const double tolerance)> solve_with_jacobian

std::function< VectorType &()> get_solution_scaling

std::function< VectorType &()> get_function_scaling

## Private Member Functions

void set_functions_to_trigger_an_assert ()

## Static Private Member Functions

static ::ExceptionBaseExcFunctionNotProvided (std::string arg1)

## Private Attributes

MPI_Comm mpi_communicator

void * kinsol_mem

SUNContext kinsol_ctx

N_Vector solution

N_Vector u_scale

N_Vector f_scale

GrowingVectorMemory< VectorType > mem

## Detailed Description

template<typename VectorType = Vector<double>>
class SUNDIALS::KINSOL< VectorType >

Interface to SUNDIALS' nonlinear solver (KINSOL).

KINSOL is a solver for nonlinear algebraic systems in residual form $$F(u) = 0$$ or fixed point form $$G(u) = u$$, where $$u$$ is a vector which we will assume to be in $${\mathbb R}^n$$ or $${\mathbb C}^n$$, but that may also have a block structure and may be distributed in parallel computations; the functions $$F$$ and $$G$$ satisfy $$F,G:{\mathbb R}^N \to{\mathbb R}^N$$ or $$F,G:{\mathbb C}^N \to{\mathbb C}^N$$. It includes a Newton-Krylov solver as well as Picard and fixed point solvers, both of which can be accelerated with Anderson acceleration. KINSOL is based on the previous Fortran package NKSOL of Brown and Saad. An example of using KINSOL can be found in the step-77 tutorial program.

KINSOL's Newton solver employs the inexact Newton method. As this solver is intended mainly for large systems, the user is required to provide their own solver function.

At the highest level, KINSOL implements the following iteration scheme:

• set $$u_0$$ = an initial guess
• For $$n = 0, 1, 2, \ldots$$ until convergence do:
• Solve $$J(u_n)\delta_n = -F(u_n)$$
• Set $$u_{n+1} = u_n + \lambda \delta_n, 0 < \lambda \leq 1$$
• Test for convergence

Here, $$u_n$$ is the $$n$$-th iterate to $$u$$, and $$J(u) = \nabla_u F(u)$$ is the system Jacobian. At each stage in the iteration process, a scalar multiple of the step $$\delta_n$$, is added to $$u_n$$ to produce a new iterate, $$u_{n+1}$$. A test for convergence is made before the iteration continues.

Unless specified otherwise by the user, KINSOL strives to update Jacobian information as infrequently as possible to balance the high costs of matrix operations against other costs. Specifically, these updates occur when:

• the problem is initialized,
• $$\|\lambda \delta_{n-1} \|_{D_u,\infty} \geq 1.5$$ (inexact Newton only, see below for a definition of $$\| \cdot \|_{D_u,\infty}$$)
• a specified number of nonlinear iterations have passed since the last update,
• the linear solver failed recoverably with outdated Jacobian information,
• the global strategy failed with outdated Jacobian information, or
• $$\|\lambda \delta_{n} \|_{D_u,\infty} \leq$$ tolerance with outdated Jacobian information.

KINSOL allows changes to the above strategy through optional solver inputs. The user can disable the initial Jacobian information evaluation or change the default value of the number of nonlinear iterations after which a Jacobian information update is enforced.

To address the case of ill-conditioned nonlinear systems, KINSOL allows prescribing scaling factors both for the solution vector and for the residual vector. For scaling to be used, the user may supply the function get_solution_scaling(), that returns values $$D_u$$, which are diagonal elements of the scaling matrix such that $$D_u u_n$$ has all components roughly the same magnitude when $$u_n$$ is close to a solution, and get_function_scaling(), that supply values $$D_F$$, which are diagonal scaling matrix elements such that $$D_F F$$ has all components roughly the same magnitude when $$u_n$$ is not too close to a solution.

When scaling values are provided for the solution vector, these values are automatically incorporated into the calculation of the perturbations used for the default difference quotient approximations for Jacobian information if the user does not supply a Jacobian solver through the solve_jacobian_system() function.

Two methods of applying a computed step $$\delta_n$$ to the previously computed solution vector are implemented. The first and simplest is the standard Newton strategy which applies the update with a constant $$\lambda$$ always set to 1. The other method is a global strategy, which attempts to use the direction implied by $$\delta_n$$ in the most efficient way for furthering convergence of the nonlinear problem. This technique is implemented in the second strategy, called Linesearch. This option employs both the $$\alpha$$ and $$\beta$$ conditions of the Goldstein-Armijo linesearch algorithm given in J. E. Dennis and R. B. Schnabel. "Numerical Methods for Unconstrained Optimization and Nonlinear Equations." SIAM, Philadelphia, 1996., where $$\lambda$$ is chosen to guarantee a sufficient decrease in $$F$$ relative to the step length as well as a minimum step length relative to the initial rate of decrease of $$F$$. One property of the algorithm is that the full Newton step tends to be taken close to the solution.

The basic fixed-point iteration scheme implemented in KINSOL is given by:

• Set $$u_0 =$$ an initial guess
• For $$n = 0, 1, 2, \dots$$ until convergence do:
• Set $$u_{n+1} = G(u_n)$$
• Test for convergence

At each stage in the iteration process, function $$G$$ is applied to the current iterate to produce a new iterate, $$u_{n+1}$$. A test for convergence is made before the iteration continues.

For Picard iteration, as implemented in KINSOL, we consider a special form of the nonlinear function $$F$$, such that $$F(u) = Lu - N(u)$$, where $$L$$ is a constant nonsingular matrix and $$N$$ is (in general) nonlinear.

Then the fixed-point function $$G$$ is defined as $$G(u) = u - L^{-1}F(u)$$. Within each iteration, the Picard step is computed then added to $$u_n$$ to produce the new iterate. Next, the nonlinear residual function is evaluated at the new iterate, and convergence is checked. The Picard and fixed point methods can be significantly accelerated using Anderson's method.

The user has to provide the implementation of the following std::functions:

• reinit_vector; and only one of
• residual; or
• iteration_function;

Specifying residual() allows the user to use Newton and Picard strategies (i.e., $$F(u)=0$$ will be solved), while specifying iteration_function(), a fixed point iteration will be used (i.e., $$G(u)=u$$ will be solved).

If the use of a Newton or Picard method is desired, then the user should also supply

• solve_jacobian_system or solve_with_jacobian; and optionally
• setup_jacobian;

Fixed point iteration does not require the solution of any linear system.

Also the following functions could be rewritten, to provide additional scaling factors for both the solution and the residual evaluation during convergence checks:

• get_solution_scaling;
• get_function_scaling;

Definition at line 181 of file kinsol.h.

## ◆ KINSOL() [1/2]

template<typename VectorType >
 SUNDIALS::KINSOL< VectorType >::KINSOL ( const AdditionalData & data = AdditionalData() )

Constructor, with class parameters set by the AdditionalData object.

Parameters
 data KINSOL configuration data
Note
With SUNDIALS 6 and later this constructor sets up logging objects to only work on the present processor (i.e., results are only communicated over MPI_COMM_SELF).

Definition at line 298 of file kinsol.cc.

## ◆ KINSOL() [2/2]

template<typename VectorType >
 SUNDIALS::KINSOL< VectorType >::KINSOL ( const AdditionalData & data, const MPI_Comm & mpi_comm )

Constructor.

Parameters
 data KINSOL configuration data mpi_comm MPI Communicator over which logging operations are computed. Only used in SUNDIALS 6 and newer.

Definition at line 305 of file kinsol.cc.

## ◆ ~KINSOL()

template<typename VectorType >
 SUNDIALS::KINSOL< VectorType >::~KINSOL

Destructor.

Definition at line 337 of file kinsol.cc.

## ◆ solve()

template<typename VectorType >
 unsigned int SUNDIALS::KINSOL< VectorType >::solve ( VectorType & initial_guess_and_solution )

Solve the non linear system. Return the number of nonlinear steps taken to converge. KINSOL uses the content of initial_guess_and_solution as initial guess, and stores the final solution in the same vector.

Definition at line 351 of file kinsol.cc.

## ◆ set_functions_to_trigger_an_assert()

template<typename VectorType >
 void SUNDIALS::KINSOL< VectorType >::set_functions_to_trigger_an_assert
private

This function is executed at construction time to set the std::function above to trigger an assert if they are not implemented.

Definition at line 624 of file kinsol.cc.

## ◆ reinit_vector

template<typename VectorType = Vector<double>>
 std::function SUNDIALS::KINSOL< VectorType >::reinit_vector

A function object that users need to supply and that is intended to reinitize the given vector to its correct size, block structure (if block vectors are used), and MPI communicator (if the vector is distributed across multiple processors using MPI), along with any other properties necessary.

Definition at line 416 of file kinsol.h.

## ◆ residual

template<typename VectorType = Vector<double>>
 std::function SUNDIALS::KINSOL< VectorType >::residual

A function object that users should supply and that is intended to compute the residual dst = F(src). This function is only used if the SolutionStrategy::newton or SolutionStrategy::linesearch strategies were selected.

This function should return:

• 0: Success
• >0: Recoverable error (KINSOL will try to change its internal parameters and attempt a new solution step)
• <0: Unrecoverable error the computation will be aborted and an assertion will be thrown.

Definition at line 431 of file kinsol.h.

## ◆ iteration_function

template<typename VectorType = Vector<double>>
 std::function SUNDIALS::KINSOL< VectorType >::iteration_function

A function object that users should supply and that is intended to compute the iteration function $$G(u)$$ for the fixed point and Picard iteration. This function is only used if the SolutionStrategy::fixed_point or SolutionStrategy::picard strategies were selected.

This function should return:

• 0: Success
• >0: Recoverable error (KINSOL will try to change its internal parameters and attempt a new solution step)
• <0: Unrecoverable error; the computation will be aborted and an assertion will be thrown.

Definition at line 448 of file kinsol.h.

## ◆ setup_jacobian

template<typename VectorType = Vector<double>>
 std::function SUNDIALS::KINSOL< VectorType >::setup_jacobian

A function object that users may supply and that is intended to prepare the linear solver for subsequent calls to solve_jacobian_system().

The job of setup_jacobian() is to prepare the linear solver for subsequent calls to solve_jacobian_system(), in the solution of linear systems $$Ax = b$$. The exact nature of this system depends on the SolutionStrategy that has been selected.

In the cases strategy = SolutionStrategy::newton or SolutionStrategy::linesearch, $$A$$ is the Jacobian $$J = \partial F/\partial u$$. If strategy = SolutionStrategy::picard, $$A$$ is the approximate Jacobian matrix $$L$$. If strategy = SolutionStrategy::fixed_point, then linear systems do not arise, and this function is never called.

The setup_jacobian() function may call a user-supplied function, or a function within the linear solver module, to compute Jacobian-related data that is required by the linear solver. It may also preprocess that data as needed for solve_jacobian_system(), which may involve calling a generic function (such as for LU factorization) or, more generally, build preconditioners from the assembled Jacobian. In any case, the data so generated may then be used whenever a linear system is solved.

The point of this function is that setup_jacobian() function is not called at every Newton iteration, but only as frequently as the solver determines that it is appropriate to perform the setup task. In this way, Jacobian-related data generated by setup_jacobian() is expected to be used over a number of Newton iterations. KINSOL determines itself when it is beneficial to regenerate the Jacobian and associated information (such as preconditioners computed for the Jacobian), thereby saving the effort to regenerate the Jacobian matrix and a preconditioner for it whenever possible.

Parameters
 current_u Current value of $$u$$ current_f Current value of $$F(u)$$ or $$G(u)$$

This function should return:

• 0: Success
• >0: Recoverable error (KINSOL will try to change its internal parameters and attempt a new solution step)
• <0: Unrecoverable error the computation will be aborted and an assertion will be thrown.

Definition at line 496 of file kinsol.h.

## ◆ solve_jacobian_system

template<typename VectorType = Vector<double>>
 std::function SUNDIALS::KINSOL< VectorType >::solve_jacobian_system
Deprecated:
Versions of SUNDIALS after 4.0 no longer provide all of the information necessary for this callback (see below). Use the solve_with_jacobian callback described below.

A function object that users may supply and that is intended to solve a linear system with the Jacobian matrix. This function will be called by KINSOL (possibly several times) after setup_jacobian() has been called at least once. KINSOL tries to do its best to call setup_jacobian() the minimum number of times. If convergence can be achieved without updating the Jacobian, then KINSOL does not call setup_jacobian() again. If, on the contrary, internal KINSOL convergence tests fail, then KINSOL calls setup_jacobian() again with updated vectors and coefficients so that successive calls to solve_jacobian_systems() lead to better convergence in the Newton process.

If you do not specify a solve_jacobian_system or solve_with_jacobian function, then only a fixed point iteration strategy can be used. Notice that this may not converge, or may converge very slowly.

A call to this function should store in dst the result of $$J^{-1}$$ applied to rhs, i.e., J*dst = rhs. It is the user's responsibility to set up proper solvers and preconditioners inside this function (or in the setup_jacobian callback above).

Arguments to the function are:

Parameters
 [in] ycur The current $$y$$ vector for the current KINSOL internal step. In the documentation above, this $$y$$ vector is generally denoted by $$u$$. [in] fcur The current value of the implicit right-hand side at ycur, $$f_I (t_n, ypred)$$. [in] rhs The system right hand side to solve for [out] dst The solution of $$J^{-1} * src$$

This function should return:

• 0: Success
• >0: Recoverable error (KINSOL will try to change its internal parameters and attempt a new solution step)
• <0: Unrecoverable error the computation will be aborted and an assertion will be thrown.
Warning
Starting with SUNDIALS 4.1, SUNDIALS no longer provides the ycur and fcur variables – only rhs is provided and dst needs to be returned. The first two arguments will therefore be empty vectors in that case. In practice, that means that one can no longer compute a Jacobian matrix for the current iterate within this function. Rather, this has to happen inside the setup_jacobian function above that receives this information. If it is important that the Jacobian corresponds to the current iterate (rather than a re-used Jacobian matrix that had been computed in a previous iteration and that therefore corresponds to a previous iterate), then you will also have to set the AdditionalData::maximum_newton_step variable to one, indicating that the Jacobian should be re-computed in every iteration.

Definition at line 560 of file kinsol.h.

## ◆ solve_with_jacobian

template<typename VectorType = Vector<double>>
 std::function< int(const VectorType &rhs, VectorType &dst, const double tolerance)> SUNDIALS::KINSOL< VectorType >::solve_with_jacobian

A function object that users may supply and that is intended to solve a linear system with the Jacobian matrix. This function will be called by KINSOL (possibly several times) after setup_jacobian() has been called at least once. KINSOL tries to do its best to call setup_jacobian() the minimum number of times. If convergence can be achieved without updating the Jacobian, then KINSOL does not call setup_jacobian() again. If, on the contrary, internal KINSOL convergence tests fail, then KINSOL calls setup_jacobian() again with updated vectors and coefficients so that successive calls to solve_jacobian_systems() lead to better convergence in the Newton process.

If you do not specify a solve_with_jacobian function, then only a fixed point iteration strategy can be used. Notice that this may not converge, or may converge very slowly.

A call to this function should store in dst the result of $$J^{-1}$$ applied to rhs, i.e., J*dst = rhs. It is the user's responsibility to set up proper solvers and preconditioners inside this function (or in the setup_jacobian callback above). The function attached to this callback is also provided with a tolerance to the linear solver, indicating that it is not necessary to solve the linear system with the Jacobian matrix exactly, but only to a tolerance that KINSOL will adapt over time.

Arguments to the function are:

Parameters
 [in] rhs The system right hand side to solve for. [out] dst The solution of $$J^{-1} * src$$. [in] tolerance The tolerance with which to solve the linear system of equations.

This function should return:

• 0: Success
• >0: Recoverable error (KINSOL will try to change its internal parameters and attempt a new solution step)
• <0: Unrecoverable error the computation will be aborted and an assertion will be thrown.

Definition at line 603 of file kinsol.h.

## ◆ get_solution_scaling

template<typename VectorType = Vector<double>>
 std::function SUNDIALS::KINSOL< VectorType >::get_solution_scaling

A function object that users may supply and that is intended to return a vector whose components are the weights used by KINSOL to compute the vector norm of the solution. The implementation of this function is optional, and it is used only if implemented.

The intent for this scaling factor is for problems in which the different components of a solution have vastly different numerical magnitudes – typically because they have different physical units and represent different things. For example, if one were to solve a nonlinear Stokes problem, the solution vector has components that correspond to velocities and other components that correspond to pressures. These have different physical units and depending on which units one chooses, they may have roughly comparable numerical sizes or maybe they don't. To give just one example, in simulations of flow in the Earth's interior, one has velocities on the order of maybe ten centimeters per year, and pressures up to around 100 GPa. If one expresses this in SI units, this corresponds to velocities of around $$0.000,000,003=3 \times 10^{-9}$$ m/s, and pressures around $$10^9 \text{kg}/\text{m}/\text{s}^2$$, i.e., vastly different. In such cases, computing the $$l_2$$ norm of a solution-type vector (e.g., the difference between the previous and the current solution) makes no sense because the norm will either be dominated by the velocity components or the pressure components. The scaling vector this function returns is intended to provide each component of the solution with a scaling factor that is generally chosen as as the inverse of a "typical velocity" or "typical pressure" so that upon multiplication of a vector component by the corresponding scaling vector component, one obtains a number that is of order of magnitude of one (i.e., a reasonably small multiple of one times the typical velocity/pressure). The KINSOL manual states this as follows: "The user should supply values \_form#2588, which are diagonal elements of the scaling matrix such that \_form#2605 has all components roughly the same magnitude when \_form#262 is close to a solution".

If no function is provided to a KINSOL object, then this is interpreted as implicitly saying that all of these scaling factors should be considered as one.

Definition at line 643 of file kinsol.h.

## ◆ get_function_scaling

template<typename VectorType = Vector<double>>
 std::function SUNDIALS::KINSOL< VectorType >::get_function_scaling

A function object that users may supply and that is intended to return a vector whose components are the weights used by KINSOL to compute the vector norm of the function evaluation away from the solution. The implementation of this function is optional, and it is used only if implemented.

The point of this function and the scaling vector it returns is similar to the one discussed above for get_solution_scaling, except that it is for a vector that scales the components of the function $$F(U)$$, rather than the components of $$U$$, when computing norms. As above, if no function is provided, then this is equivalent to using a scaling vector whose components are all equal to one.

Definition at line 659 of file kinsol.h.

## ◆ data

template<typename VectorType = Vector<double>>
private

KINSOL configuration data.

Definition at line 692 of file kinsol.h.

## ◆ mpi_communicator

template<typename VectorType = Vector<double>>
 MPI_Comm SUNDIALS::KINSOL< VectorType >::mpi_communicator
private

The MPI communicator to be used by this solver, if any.

Definition at line 697 of file kinsol.h.

## ◆ kinsol_mem

template<typename VectorType = Vector<double>>
 void* SUNDIALS::KINSOL< VectorType >::kinsol_mem
private

KINSOL memory object.

Definition at line 702 of file kinsol.h.

## ◆ kinsol_ctx

template<typename VectorType = Vector<double>>
 SUNContext SUNDIALS::KINSOL< VectorType >::kinsol_ctx
private

A context object associated with the KINSOL solver.

Definition at line 708 of file kinsol.h.

## ◆ solution

template<typename VectorType = Vector<double>>
 N_Vector SUNDIALS::KINSOL< VectorType >::solution
private

KINSOL solution vector.

Definition at line 714 of file kinsol.h.

## ◆ u_scale

template<typename VectorType = Vector<double>>
 N_Vector SUNDIALS::KINSOL< VectorType >::u_scale
private

KINSOL solution scale.

Definition at line 719 of file kinsol.h.

## ◆ f_scale

template<typename VectorType = Vector<double>>
 N_Vector SUNDIALS::KINSOL< VectorType >::f_scale
private

KINSOL f scale.

Definition at line 724 of file kinsol.h.

## ◆ mem

template<typename VectorType = Vector<double>>
 GrowingVectorMemory SUNDIALS::KINSOL< VectorType >::mem
private

Memory pool of vectors.

Definition at line 729 of file kinsol.h.

The documentation for this class was generated from the following files: