From 16ef332b660fd180e79ab9d759b145cf1050ffb1 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Thu, 28 Sep 2017 16:48:04 +0200 Subject: [PATCH] KINSOL support. --- doc/news/changes/major/20170928LucaHeltai | 7 + include/deal.II/sundials/kinsol.h | 664 ++++++++++++++++++++++ source/sundials/CMakeLists.txt | 1 + source/sundials/kinsol.cc | 336 +++++++++++ tests/sundials/kinsol_01.cc | 82 +++ tests/sundials/kinsol_01.output | 4 + tests/sundials/kinsol_01.prm | 16 + 7 files changed, 1110 insertions(+) create mode 100644 doc/news/changes/major/20170928LucaHeltai create mode 100644 include/deal.II/sundials/kinsol.h create mode 100644 source/sundials/kinsol.cc create mode 100644 tests/sundials/kinsol_01.cc create mode 100644 tests/sundials/kinsol_01.output create mode 100644 tests/sundials/kinsol_01.prm diff --git a/doc/news/changes/major/20170928LucaHeltai b/doc/news/changes/major/20170928LucaHeltai new file mode 100644 index 0000000000..f37bcb293d --- /dev/null +++ b/doc/news/changes/major/20170928LucaHeltai @@ -0,0 +1,7 @@ +New: Added support for the KINSOL solver of the SUNDIALS +library. KINSOL is a solver for nonlinear algebraic systems. +It includes a Newton-Krylov solver as well as Picard and +fixed point solvers, both of which can be accelerated with +Anderson acceleration. +
+(Luca Heltai, 2017/09/28) diff --git a/include/deal.II/sundials/kinsol.h b/include/deal.II/sundials/kinsol.h new file mode 100644 index 0000000000..1e2b7a3a83 --- /dev/null +++ b/include/deal.II/sundials/kinsol.h @@ -0,0 +1,664 @@ +//----------------------------------------------------------- +// +// Copyright (C) 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +//--------------------------------------------------------------- + +#ifndef dealii_sundials_kinsol_h +#define dealii_sundials_kinsol_h + +#include +#ifdef DEAL_II_WITH_SUNDIALS + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include + + +DEAL_II_NAMESPACE_OPEN + +// Shorthand notation for KINSOL error codes. +#define AssertKINSOL(code) Assert(code >= 0, ExcKINSOLError(code)) + +namespace SUNDIALS +{ + /** + * Interface to SUNDIALS non linear solver (KINSOL). + * + * KINSOL is a solver for nonlinear algebraic systems. 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 its + * 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: + * - set u0 = 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 \detla_n, 0 < \lambda \leq 1$ + * - Test for convergence + * + * Here, $u_n$ is the $n$-th iterate to $u$, and $J(u) = \partial_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 un 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) + * - 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_residual_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 + * - get_lower_than_zero_constrained_entries() + * - get_greater_than_zero_constrained_entries() + * - get_lower_equal_than_zero_constrained_entries() + * - get_greater_or_equal_than_zero_constrained_entries() + * + * 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: + * - 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 strategies (i.e., + * $F(u)=0$ will be solved), while specifying iteration_function(), fixed + * point iteration or Pircard 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 + * - solve_jacobian_system; + * and optionally + * - setup_jacobian; + * + * 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: + * - get_solution_scaling; + * - get_function_scaling; + * + * @author Luca Heltai, 2017. + */ + template > + class KINSOL + { + public: + + /** + * Additional parameters that can be passed to the KINSOL class. + */ + class AdditionalData + { + public: + /** + * KINSOL solution strategy. KINSOL includes a Newton-Krylov solver (both + * local and global) as well as Picard and fixed point solvers. + */ + enum SolutionStrategy + { + /** + * Standard Newton iteration. + */ + newton = KIN_NONE, + /** + * Newton iteration with linesearch. + */ + linesearch = KIN_LINESEARCH, + /** + * Fixed point iteration. + */ + fixed_point = KIN_FP, + /** + * Picard iteration. + */ + picard = KIN_PICARD, + }; + + /** + * Initialization parameters for KINSOL. + * + * Global parameters: + * + * @param strategy Solution strategy + * @param maximum_non_linear_iterations Maximum number of nonlinear iterations + * @param function_tolerance Function norm stopping tolerance + * @param step_tolerance Scaled step stopping tolerance + * + * Newton parameters: + * + * @param no_init_setup No initial matrix setup + * @param maximum_setup_calls Maximum iterations without matrix setup + * @param maximum_newton_step Maximum allowable scaled length of the Newton step + * @param dq_relative_error Relative error for different quotient computation + * + * Linesearch parameters: + * + * @param maximum_beta_failures Maximum number of beta-condition failures + * + * Fixed point and Picard parameters: + * + * @param anderson_subspace_size Anderson acceleration subspace size + */ + AdditionalData( + // Global parameters + const SolutionStrategy &strategy = linesearch, + const unsigned int &maximum_non_linear_iterations = 200, + const double &function_tolerance = 0.0, + const double &step_tolerance = 0.0, + const bool &no_init_setup = false, + const unsigned int &maximum_setup_calls = 0, + const double &maximum_newton_step = 0.0, + const double &dq_relative_error = 0.0, + const unsigned int &maximum_beta_failures = 0, + const unsigned int &anderson_subspace_size = 0) : + strategy(strategy), + maximum_non_linear_iterations(maximum_non_linear_iterations), + function_tolerance(function_tolerance), + step_tolerance(step_tolerance), + no_init_setup(no_init_setup), + maximum_setup_calls(maximum_setup_calls), + maximum_newton_step(maximum_newton_step), + dq_relative_error(dq_relative_error), + maximum_beta_failures(maximum_beta_failures), + anderson_subspace_size(anderson_subspace_size) + {}; + + /** + * Add all AdditionalData() parameters to the given ParameterHandler + * object. When the parameters are parsed from a file, the internal + * parameters are automatically updated. + * + * The following parameters are declared: + * + * @code + * @endcode + * + * These are one-to-one with the options you can pass at construction time. + * + * The options you pass at construction time are set as default values in + * the ParameterHandler object `prm`. You can later modify them by parsing + * a parameter file using `prm`. The values of the parameter will be updated + * whenever the content of `prm` is updated. + * + * Make sure that this class lives longer than `prm`. Undefined behaviour + * will occurr if you destroy this class, and then parse a parameter file + * using `prm`. + */ + void add_parameters(ParameterHandler &prm) + { + static std::string strategy_str("newton"); + prm.add_parameter("Solution strategy", strategy_str, + "Choose among newton|linesearch|fixed_point|picard", + Patterns::Selection("newton|linesearch|fixed_point|picard")); + prm.add_action("Solution strategy", [&](const std::string &value) + { + if (value == "newton") + strategy = newton; + else if (value == "linesearch") + strategy = linesearch; + else if (value == "fixed_point") + strategy = fixed_point; + else if (value == "picard") + strategy = picard; + else + Assert(false, ExcInternalError()); + } + ); + prm.add_parameter("Maximum number of nonlinear iterations", + maximum_non_linear_iterations); + prm.add_parameter("Function norm stopping tolerance", + function_tolerance); + prm.add_parameter("Scaled step stopping tolerance", + step_tolerance); + + prm.enter_subsection("Newton parameters"); + prm.add_parameter("No initial matrix setup", + no_init_setup); + prm.add_parameter("Maximum iterations without matrix setup", + maximum_setup_calls); + prm.add_parameter("Maximum allowable scaled length of the Newton step", + maximum_newton_step); + prm.add_parameter("Relative error for different quotient computation", + dq_relative_error); + prm.leave_subsection(); + + prm.enter_subsection("Linesearch parameters"); + prm.add_parameter("Maximum number of beta-condition failures", + maximum_beta_failures); + prm.leave_subsection(); + + + prm.enter_subsection("Fixed point and Picard parameters"); + prm.add_parameter("Anderson acceleration subspace size", + anderson_subspace_size); + prm.leave_subsection(); + } + + /** + * The solution strategy to use. If you choose SolutionStrategy::newton + * or SolutionStrategy::linesearch, you have to provide also the function + * residual(). If you choose SolutionStrategy::picard or + * SolutionStrategy::fixed_point, you have to provide also the function + * iteration_function(). + */ + SolutionStrategy strategy; + + /** + * Maximum number of nonlinear iterations allowed. + */ + unsigned int maximum_non_linear_iterations; + + /** + * Specifies the scalar used as a stopping tolerance on the scaled + * maximum norm of the system function $F(u)$ or $G(u)$. + * + * Pass 0.0 to use KINSOL defaults. + */ + double function_tolerance; + + /** + * Specifies the scalar used as a stopping tolerance on the minimum + * scaled step length. + * + * Pass 0.0 to use KINSOL defaults. + */ + double step_tolerance; + + /** + * Specifies whether an initial call to the preconditioner or Jacobian + * setup function should be made or not. + * + * A call to this function is useful when solving a sequence of problems, + * in which the final preconditioner or Jacobian value from one problem + * is to be used initially for the next problem. + */ + bool no_init_setup; + + /** + * Specifies the maximum number of nonlinear iterations that can be + * performed between calls to the setup_jacobian() function. + * + * Pass 0.0 to use KINSOL defaults. + */ + unsigned int maximum_setup_calls; + + /** + * Specifies the maximum allowable scaled length of the Newton step. + * + * Pass 0.0 to use KINSOL defaults. + */ + double maximum_newton_step; + + /** + * Specifies the relative error in computing $F(u)$, which is used in the + * difference quotient approximation to the Jacobian matrix when the user + * does not supply a solve_jacobian_system_matrix() function. + * + * Pass 0.0 to use KINSOL defaults. + */ + double dq_relative_error; + + /** + * Specifies the maximum number of beta-condition failures in the + * linesearch algorithm. Only used if + * strategy==SolutionStrategy::linesearch. + */ + unsigned int maximum_beta_failures; + + /** + * Specifies the size of the subspace used with Anderson acceleration + * in conjunction with Picard or fixed-point iteration. + * + * If you set this to 0, no acceleration is used. + */ + unsigned int anderson_subspace_size; + }; + + /** + * Constructor. It is possible to fine tune the SUNDIALS KINSOL solver by + * passing an AdditionalData() object that sets all of the solver + * parameters. + * + * @param data KINSOL configuration data + * @param mpi_comm MPI communicator + */ + KINSOL(const AdditionalData &data=AdditionalData(), + const MPI_Comm mpi_comm = MPI_COMM_WORLD); + + /** + * Destructor. + */ + ~KINSOL(); + + /** + * Solve the non linear sytem. Return the number of nonlinear steps taken + * to converge. KINSOL uses the content of `solution` as initial guess, and + * stores the final solution in the same vector. + */ + unsigned int solve(VectorType &initial_guess_and_solution); + + /** + * A function object that users need to supply and that is intended to + * reinit the given vector. + */ + std::function reinit_vector; + + /** + * A function object that users should may 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: + * - 0: Success + * - >0: Recoverable error (KINSOLReinit will be called if this happens, and + * then last function will be attempted again + * - <0: Unrecoverable error the computation will be aborted and an assertion + * will be thrown. + */ + std::function residual; + + /** + * A function object that users may 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: + * - 0: Success + * - >0: Recoverable error (KINSOLReinit will be called if this happens, and + * then last function will be attempted again + * - <0: Unrecoverable error the computation will be aborted and an assertion + * will be thrown. + */ + std::function iteration_function; + + + /** + * 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. + * + * @param current_u Current value of u + * @param current_f Current value of F(u) or G(u) + * + * This function should return: + * - 0: Success + * - >0: Recoverable error (KINSOLReinit will be called if this happens, and + * then last function will be attempted again + * - <0: Unrecoverable error the computation will be aborted and an assertion + * will be thrown. + */ + std::function setup_jacobian; + + /** + * 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 coefficents 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 a fixed + * point iteration is used instead of a Newton method. Notice that this may + * not converge, or may converge very slowly. + * + * The jacobian $J$ should be (an approximation of) the system Jacobian + * \f[ + * J = M - \gamma \frac{\partial f_I}{\partial y} + * \f] + * evaluated at `t`, `ycur`. `fcur` is $f_I(t,ycur)$. + * + * 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 responsability to + * set up proper solvers and preconditioners inside this function. + * + * + * Arguments to the function are + * + * @param[in] t the current time + * @param[in] gamma the current factor to use in the jacobian computation + * @param[in] ycur is the current $y$ vector for the current KINSOL internal step + * @param[in] fcur is the current value of the implicit right-hand side at ycur, + * $f_I (t_n, ypred)$. + * + * + * This function should return: + * - 0: Success + * - >0: Recoverable error (KINSOLReinit will be called if this happens, and + * then last function will be attempted again + * - <0: Unrecoverable error the computation will be aborted and an assertion + * will be thrown. + */ + std::function solve_jacobian_system; + + /** + * 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. + */ + std::function 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 function evaluation away from the solution. The + * implementation of this function is optional, and it is used only if + * implemented. + */ + std::function get_function_scaling; + + /** + * Handle KINSOL exceptions. + */ + DeclException1(ExcKINSOLError, int, << "One of the SUNDIALS KINSOL internal functions " + << " returned a negative error code: " + << arg1 << ". Please consult SUNDIALS manual."); + + + private: + + /** + * Throw an exception when a function with the given name is not implemented. + */ + DeclException1(ExcFunctionNotProvided, std::string, + << "Please provide an implementation for the function \"" << arg1 << "\""); + + /** + * This function is executed at construction time to set the + * std::function above to trigger an assert if they are not + * implemented. + */ + void set_functions_to_trigger_an_assert(); + + /** + * KINSOL configuration data. + */ + AdditionalData data; + + /** + * KINSOL memory object. + */ + void *kinsol_mem; + + /** + * KINSOL solution vector. + */ + N_Vector solution; + + /** + * KINSOL solution scale. + */ + N_Vector u_scale; + + /** + * KINSOL f scale. + */ + N_Vector f_scale; + +#ifdef DEAL_II_WITH_MPI + /** + * MPI communicator. SUNDIALS solver runs happily in parallel. + */ + MPI_Comm communicator; +#endif + + /** + * Memory pool of vectors. + */ + GrowingVectorMemory mem; + }; + +} + + +DEAL_II_NAMESPACE_CLOSE +#endif + + +#endif diff --git a/source/sundials/CMakeLists.txt b/source/sundials/CMakeLists.txt index 238e7370d5..c21a9125a5 100644 --- a/source/sundials/CMakeLists.txt +++ b/source/sundials/CMakeLists.txt @@ -19,6 +19,7 @@ SET(_src arkode.cc ida.cc copy.cc + kinsol.cc ) SET(_inst diff --git a/source/sundials/kinsol.cc b/source/sundials/kinsol.cc new file mode 100644 index 0000000000..fbffd73efb --- /dev/null +++ b/source/sundials/kinsol.cc @@ -0,0 +1,336 @@ +//----------------------------------------------------------- +// +// Copyright (C) 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +//----------------------------------------------------------- + + +#include +#include + +#ifdef DEAL_II_WITH_SUNDIALS + +#include +#include +#ifdef DEAL_II_WITH_TRILINOS +#include +#include +#endif +#ifdef DEAL_II_WITH_PETSC +#include +#include +#endif +#include +#include + +#include + +#include +#include + +DEAL_II_NAMESPACE_OPEN + +namespace SUNDIALS +{ + using namespace internal; + + namespace + { + template + int t_kinsol_function(N_Vector yy, + N_Vector FF, + void *user_data) + { + KINSOL &solver = *static_cast *>(user_data); + GrowingVectorMemory mem; + + typename VectorMemory::Pointer src_yy(mem); + solver.reinit_vector(*src_yy); + + typename VectorMemory::Pointer dst_FF(mem); + solver.reinit_vector(*dst_FF); + + copy(*src_yy, yy); + + int err = 0; + if (solver.residual) + err = solver.residual(*src_yy, *dst_FF); + else if (solver.iteration_function) + err = solver.iteration_function(*src_yy, *dst_FF); + else + Assert(false, ExcInternalError()); + + copy(FF, *dst_FF); + + return err; + } + + + + template + int t_kinsol_setup_jacobian(KINMem kinsol_mem) + { + KINSOL &solver = *static_cast *>(kinsol_mem->kin_user_data); + GrowingVectorMemory mem; + + typename VectorMemory::Pointer src_ycur(mem); + solver.reinit_vector(*src_ycur); + + typename VectorMemory::Pointer src_fcur(mem); + solver.reinit_vector(*src_fcur); + + copy(*src_ycur, kinsol_mem->kin_uu); + copy(*src_fcur, kinsol_mem->kin_fval); + + int err = solver.setup_jacobian(*src_ycur, *src_fcur); + return err; + } + + + + template + int t_kinsol_solve_jacobian(KINMem kinsol_mem, + N_Vector x, + N_Vector b, + realtype *sJpnorm, + realtype *sFdotJp) + { + KINSOL &solver = *static_cast *>(kinsol_mem->kin_user_data); + GrowingVectorMemory mem; + + typename VectorMemory::Pointer src_ycur(mem); + solver.reinit_vector(*src_ycur); + + typename VectorMemory::Pointer src_fcur(mem); + solver.reinit_vector(*src_fcur); + + copy(*src_ycur, kinsol_mem->kin_uu); + copy(*src_fcur, kinsol_mem->kin_fval); + + typename VectorMemory::Pointer src(mem); + solver.reinit_vector(*src); + + typename VectorMemory::Pointer dst(mem); + solver.reinit_vector(*dst); + + copy(*src, b); + + int err = solver.solve_jacobian_system(*src_ycur, *src_fcur, + *src,*dst); + copy(x, *dst); + + *sJpnorm = N_VWL2Norm(b, kinsol_mem->kin_fscale); + N_VProd(b, kinsol_mem->kin_fscale, b); + N_VProd(b, kinsol_mem->kin_fscale, b); + *sFdotJp = N_VDotProd(kinsol_mem->kin_fval, b); + + return err; + } + } + + template + KINSOL::KINSOL(const AdditionalData &data, const MPI_Comm mpi_comm) : + data(data), + kinsol_mem(nullptr), + communicator(Utilities::MPI::duplicate_communicator(mpi_comm)) + { + set_functions_to_trigger_an_assert(); + } + + + + template + KINSOL::~KINSOL() + { + if (kinsol_mem) + KINFree(&kinsol_mem); + MPI_Comm_free(&communicator); + } + + + + template + unsigned int KINSOL::solve(VectorType &initial_guess_and_solution) + { + unsigned int system_size = initial_guess_and_solution.size(); + unsigned int local_system_size = system_size; + + // The solution is stored in + // solution. Here we take only a + // view of it. +#ifdef DEAL_II_WITH_MPI + if (is_serial_vector::value == false) + { + IndexSet is = initial_guess_and_solution.locally_owned_elements(); + local_system_size = is.n_elements(); + + solution = N_VNew_Parallel(communicator, + local_system_size, + system_size); + + u_scale = N_VNew_Parallel(communicator, + local_system_size, + system_size); + N_VConst_Parallel( 1.e0, u_scale ); + + f_scale = N_VNew_Parallel(communicator, + local_system_size, + system_size); + N_VConst_Parallel( 1.e0, f_scale ); + } + else +#endif + { + Assert(is_serial_vector::value, + ExcInternalError("Trying to use a serial code with a parallel vector.")); + solution = N_VNew_Serial(system_size); + u_scale = N_VNew_Serial(system_size); + N_VConst_Serial( 1.e0, u_scale ); + f_scale = N_VNew_Serial(system_size); + N_VConst_Serial( 1.e0, f_scale ); + } + + if (get_solution_scaling) + copy(u_scale, get_solution_scaling()); + + if (get_function_scaling) + copy(f_scale, get_function_scaling()); + + copy(solution, initial_guess_and_solution); + + if (kinsol_mem) + KINFree(&kinsol_mem); + + kinsol_mem = KINCreate(); + + int status = KINInit(kinsol_mem, t_kinsol_function , solution); + AssertKINSOL(status); + + status = KINSetUserData(kinsol_mem, (void *) this); + AssertKINSOL(status); + + status = KINSetNumMaxIters(kinsol_mem, data.maximum_non_linear_iterations); + AssertKINSOL(status); + + status = KINSetFuncNormTol(kinsol_mem, data.function_tolerance); + AssertKINSOL(status); + + status = KINSetScaledStepTol(kinsol_mem, data.step_tolerance); + AssertKINSOL(status); + + status = KINSetMaxSetupCalls(kinsol_mem, data.maximum_setup_calls); + AssertKINSOL(status); + + status = KINSetNoInitSetup(kinsol_mem, (int) data.no_init_setup); + AssertKINSOL(status); + + status = KINSetMaxNewtonStep(kinsol_mem, data.maximum_newton_step); + AssertKINSOL(status); + + status = KINSetMaxBetaFails(kinsol_mem, data.maximum_beta_failures); + AssertKINSOL(status); + + status = KINSetMAA(kinsol_mem, data.anderson_subspace_size); + AssertKINSOL(status); + + status = KINSetRelErrFunc(kinsol_mem, data.dq_relative_error); + AssertKINSOL(status); + + if (solve_jacobian_system) + { + KINMem KIN_mem = (KINMem) kinsol_mem; + KIN_mem->kin_lsolve = t_kinsol_solve_jacobian; + if (setup_jacobian) + { + KIN_mem->kin_lsetup = t_kinsol_setup_jacobian; + KIN_mem->kin_setupNonNull = true; + } + } + else + { + status = KINDense(kinsol_mem, system_size); + AssertKINSOL(status); + } + + if (data.strategy == AdditionalData::newton || + data.strategy == AdditionalData::linesearch) + Assert(residual, ExcFunctionNotProvided("residual")); + + if (data.strategy == AdditionalData::fixed_point || + data.strategy == AdditionalData::picard) + Assert(iteration_function, ExcFunctionNotProvided("iteration_function")); + + // call to KINSol + status = KINSol(kinsol_mem, solution, (int) data.strategy, u_scale, f_scale); + AssertKINSOL(status); + + copy(initial_guess_and_solution, solution ); + + // Free the vectors which are no longer used. +#ifdef DEAL_II_WITH_MPI + if (is_serial_vector::value == false) + { + N_VDestroy_Parallel(solution); + N_VDestroy_Parallel(u_scale); + N_VDestroy_Parallel(f_scale); + } + else +#endif + { + N_VDestroy_Serial(solution); + N_VDestroy_Serial(u_scale); + N_VDestroy_Serial(f_scale); + } + + long nniters; + status = KINGetNumNonlinSolvIters(kinsol_mem, &nniters); + AssertKINSOL(status); + + + KINFree(&kinsol_mem); + + return (unsigned int) nniters; + } + + template + void KINSOL::set_functions_to_trigger_an_assert() + { + + reinit_vector = [](VectorType &) + { + AssertThrow(false, ExcFunctionNotProvided("reinit_vector")); + }; + + } + + template class KINSOL >; + template class KINSOL >; + +#ifdef DEAL_II_WITH_MPI + +#ifdef DEAL_II_WITH_TRILINOS + template class KINSOL; + template class KINSOL; +#endif + +#ifdef DEAL_II_WITH_PETSC + template class KINSOL; + template class KINSOL; +#endif + +#endif + +} + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/tests/sundials/kinsol_01.cc b/tests/sundials/kinsol_01.cc new file mode 100644 index 0000000000..9ad613ab87 --- /dev/null +++ b/tests/sundials/kinsol_01.cc @@ -0,0 +1,82 @@ +//----------------------------------------------------------- +// +// Copyright (C) 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +//----------------------------------------------------------- + +#include "../tests.h" +#include +#include +#include +#include + +// provide only residual function, use internal solver. + +/** + * Solve the non linear problem + * + * F(u) = 0 , where f_i(u) = u_i^2 - i^2, 0 <= i < N + * + */ +int main (int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, numbers::invalid_unsigned_int); + + typedef Vector VectorType; + + SUNDIALS::KINSOL::AdditionalData data; + ParameterHandler prm; + data.add_parameters(prm); + + if (false) + { + std::ofstream ofile(SOURCE_DIR "/kinsol_01.prm"); + prm.print_parameters(ofile, ParameterHandler::ShortText); + ofile.close(); + } + + std::ifstream ifile(SOURCE_DIR "/kinsol_01.prm"); + prm.parse_input(ifile); + + // Size of the problem + unsigned int N=10; + + SUNDIALS::KINSOL kinsol(data); + + kinsol.reinit_vector = [N] (VectorType &v) + { + v.reinit(N); + }; + + kinsol.residual = [] (const VectorType &u, VectorType &F) -> int + { + for (unsigned int i=0; i int + { + for (unsigned int i=0; i