 |
Reference documentation for deal.II version 9.2.0
|
|
Reference documentation for deal.II version 9.2.0
|
#include <deal.II/sundials/kinsol.h>
Classes | |
| class | AdditionalData |
Public Member Functions | |
| KINSOL (const AdditionalData &data=AdditionalData(), const MPI_Comm mpi_comm=MPI_COMM_WORLD) | |
| ~KINSOL () | |
| unsigned int | solve (VectorType &initial_guess_and_solution) |
Static Public Member Functions | |
| static ::ExceptionBase & | ExcKINSOLError (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 ¤t_u, const VectorType ¤t_f)> | setup_jacobian |
| std::function< int(const VectorType &ycur, const VectorType &fcur, const VectorType &rhs, VectorType &dst)> | solve_jacobian_system |
| 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 ::ExceptionBase & | ExcFunctionNotProvided (std::string arg1) |
Private Attributes | |
| AdditionalData | data |
| void * | kinsol_mem |
| N_Vector | solution |
| N_Vector | u_scale |
| N_Vector | f_scale |
| MPI_Comm | communicator |
| GrowingVectorMemory< VectorType > | mem |
Interface to SUNDIALS non linear solver (KINSOL).
KINSOL is a solver for nonlinear algebraic systems in residual form \(F(u) = 0\) or fixed point form \(G(u) = u\). 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.
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. If a solver function is not provided, the internal dense solver of KINSOL is used. Be warned that this solver computes the Jacobian approximately, and may be efficient only for small systems.
At the highest level, KINSOL implements the following iteration scheme:
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:
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.
As a user option, KINSOL permits the application of inequality constraints, \(u_i > 0\) and \(u_i < 0\), as well as \(u_i \geq 0\) and \(u_i \leq 0\), where \(u_i\) is the \(i\)-th component of \(u\). Any such constraint, or no constraint, may be imposed on each component by providing the optional functions
KINSOL will reduce step lengths in order to ensure that no constraint is violated. Specifically, if a new Newton iterate will violate a constraint, the maximum step length along the Newton direction that will satisfy all constraints is found, and \(\delta_n\) is scaled to take a step of that length.
The basic fixed-point iteration scheme implemented in KINSOL is given by:
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:
Specifying residual() allows the user to use Newton strategies (i.e., \(F(u)=0\) will be solved), while specifying iteration_function(), fixed point iteration or Picard iteration will be used (i.e., \(G(u)=u\) will be solved).
If the use of a Newton method is desired, then the user should also supply
If the solve_jacobian_system() function is not supplied, then KINSOL will use its internal dense solver for Newton methods, with approximate Jacobian. This may be very expensive for large systems. 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:
| SUNDIALS::KINSOL< VectorType >::KINSOL | ( | const AdditionalData & | data = AdditionalData(), |
| const MPI_Comm | mpi_comm = MPI_COMM_WORLD |
||
| ) |
| SUNDIALS::KINSOL< VectorType >::~KINSOL |
| unsigned int SUNDIALS::KINSOL< VectorType >::solve | ( | VectorType & | initial_guess_and_solution | ) |
|
private |
| std::function<void(VectorType &)> SUNDIALS::KINSOL< VectorType >::reinit_vector |
| std::function<int(const VectorType &src, VectorType &dst)> 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 are specified.
This function should return:
| std::function<int(const VectorType &src, VectorType &dst)> 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 are specified.
This function should return:
| std::function<int(const VectorType ¤t_u, const VectorType ¤t_f)> 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). This data may be intended either for direct use (in a direct linear solver) or for use in a preconditioner (in a preconditioned iterative linear solver).
The 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.
| current_u | Current value of u |
| current_f | Current value of F(u) or G(u) |
This function should return:
| std::function<int(const VectorType &ycur, const VectorType &fcur, const VectorType &rhs, VectorType & dst)> SUNDIALS::KINSOL< VectorType >::solve_jacobian_system |
A function object that users may supply and that is intended to solve the Jacobian linear system. 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 amount 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 again setup_jacobian() 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() 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 src, i.e., J*dst = src. It is the users responsibility to set up proper solvers and preconditioners inside this function.
Arguments to the function are
| [in] | ycur | is the current \(y\) vector for the current KINSOL internal step |
| [in] | fcur | is 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 \(A^{-1} * src\) |
This function should return:
| std::function<VectorType &()> SUNDIALS::KINSOL< VectorType >::get_solution_scaling |
| std::function<VectorType &()> 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.
|
private |
|
private |
|
private |
|
private |
|
private |
|
private |
|
private |
1.8.17