]> https://gitweb.dealii.org/ - dealii.git/commitdiff
KINSOL support.
authorLuca Heltai <luca.heltai@sissa.it>
Thu, 28 Sep 2017 14:48:04 +0000 (16:48 +0200)
committerLuca Heltai <luca.heltai@sissa.it>
Thu, 28 Sep 2017 14:48:15 +0000 (16:48 +0200)
doc/news/changes/major/20170928LucaHeltai [new file with mode: 0644]
include/deal.II/sundials/kinsol.h [new file with mode: 0644]
source/sundials/CMakeLists.txt
source/sundials/kinsol.cc [new file with mode: 0644]
tests/sundials/kinsol_01.cc [new file with mode: 0644]
tests/sundials/kinsol_01.output [new file with mode: 0644]
tests/sundials/kinsol_01.prm [new file with mode: 0644]

diff --git a/doc/news/changes/major/20170928LucaHeltai b/doc/news/changes/major/20170928LucaHeltai
new file mode 100644 (file)
index 0000000..f37bcb2
--- /dev/null
@@ -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. 
+<br>
+(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 (file)
index 0000000..1e2b7a3
--- /dev/null
@@ -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 <deal.II/base/config.h>
+#ifdef DEAL_II_WITH_SUNDIALS
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/parameter_handler.h>
+#include <deal.II/base/conditional_ostream.h>
+#include <deal.II/base/mpi.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/vector_view.h>
+#include <deal.II/lac/vector_memory.h>
+
+#include <kinsol/kinsol.h>
+#include <kinsol/kinsol_impl.h>
+#include <nvector/nvector_serial.h>
+#include <sundials/sundials_math.h>
+#include <sundials/sundials_types.h>
+
+#include <boost/signals2.hpp>
+#include <memory>
+
+
+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<typename VectorType=Vector<double> >
+  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<void(VectorType &)> 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<int(const VectorType &src,
+                      VectorType &dst)> 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<int(const VectorType &src,
+                      VectorType &dst)> 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<int(const VectorType &current_u,
+                      const VectorType &current_f)> 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<int(const VectorType &ycur,
+                      const VectorType &fcur,
+                      const VectorType &rhs,
+                      VectorType &dst)> 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<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 function evaluation away from the solution. The
+     * implementation of this function is optional, and it is used only if
+     * implemented.
+     */
+    std::function<VectorType&()> 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<VectorType> mem;
+  };
+
+}
+
+
+DEAL_II_NAMESPACE_CLOSE
+#endif
+
+
+#endif
index 238e7370d562690f8e56284f47cb702fbc80bfd6..c21a9125a5505862520b3b9a69e34da401d8f2a4 100644 (file)
@@ -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 (file)
index 0000000..fbffd73
--- /dev/null
@@ -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 <deal.II/sundials/kinsol.h>
+#include <deal.II/base/config.h>
+
+#ifdef DEAL_II_WITH_SUNDIALS
+
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/block_vector.h>
+#ifdef DEAL_II_WITH_TRILINOS
+#include <deal.II/lac/trilinos_parallel_block_vector.h>
+#include <deal.II/lac/trilinos_vector.h>
+#endif
+#ifdef DEAL_II_WITH_PETSC
+#include <deal.II/lac/petsc_parallel_block_vector.h>
+#include <deal.II/lac/petsc_parallel_vector.h>
+#endif
+#include <deal.II/base/utilities.h>
+#include <deal.II/sundials/copy.h>
+
+#include <kinsol/kinsol_dense.h>
+
+#include <iostream>
+#include <iomanip>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace SUNDIALS
+{
+  using namespace internal;
+
+  namespace
+  {
+    template<typename VectorType>
+    int t_kinsol_function(N_Vector yy,
+                          N_Vector FF,
+                          void *user_data)
+    {
+      KINSOL<VectorType> &solver = *static_cast<KINSOL<VectorType> *>(user_data);
+      GrowingVectorMemory<VectorType> mem;
+
+      typename VectorMemory<VectorType>::Pointer src_yy(mem);
+      solver.reinit_vector(*src_yy);
+
+      typename VectorMemory<VectorType>::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<typename VectorType>
+    int t_kinsol_setup_jacobian(KINMem kinsol_mem)
+    {
+      KINSOL<VectorType> &solver = *static_cast<KINSOL<VectorType> *>(kinsol_mem->kin_user_data);
+      GrowingVectorMemory<VectorType> mem;
+
+      typename VectorMemory<VectorType>::Pointer src_ycur(mem);
+      solver.reinit_vector(*src_ycur);
+
+      typename VectorMemory<VectorType>::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<typename VectorType>
+    int t_kinsol_solve_jacobian(KINMem kinsol_mem,
+                                N_Vector x,
+                                N_Vector b,
+                                realtype *sJpnorm,
+                                realtype *sFdotJp)
+    {
+      KINSOL<VectorType> &solver = *static_cast<KINSOL<VectorType> *>(kinsol_mem->kin_user_data);
+      GrowingVectorMemory<VectorType> mem;
+
+      typename VectorMemory<VectorType>::Pointer src_ycur(mem);
+      solver.reinit_vector(*src_ycur);
+
+      typename VectorMemory<VectorType>::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<VectorType>::Pointer src(mem);
+      solver.reinit_vector(*src);
+
+      typename VectorMemory<VectorType>::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 <typename VectorType>
+  KINSOL<VectorType>::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 <typename VectorType>
+  KINSOL<VectorType>::~KINSOL()
+  {
+    if (kinsol_mem)
+      KINFree(&kinsol_mem);
+    MPI_Comm_free(&communicator);
+  }
+
+
+
+  template <typename VectorType>
+  unsigned int KINSOL<VectorType>::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<VectorType>::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<VectorType>::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<VectorType> , 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<VectorType>;
+        if (setup_jacobian)
+          {
+            KIN_mem->kin_lsetup = t_kinsol_setup_jacobian<VectorType>;
+            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<VectorType>::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<typename VectorType>
+  void KINSOL<VectorType>::set_functions_to_trigger_an_assert()
+  {
+
+    reinit_vector = [](VectorType &)
+    {
+      AssertThrow(false, ExcFunctionNotProvided("reinit_vector"));
+    };
+
+  }
+
+  template class KINSOL<Vector<double> >;
+  template class KINSOL<BlockVector<double> >;
+
+#ifdef DEAL_II_WITH_MPI
+
+#ifdef DEAL_II_WITH_TRILINOS
+  template class KINSOL<TrilinosWrappers::MPI::Vector>;
+  template class KINSOL<TrilinosWrappers::MPI::BlockVector>;
+#endif
+
+#ifdef DEAL_II_WITH_PETSC
+  template class KINSOL<PETScWrappers::MPI::Vector>;
+  template class KINSOL<PETScWrappers::MPI::BlockVector>;
+#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 (file)
index 0000000..9ad613a
--- /dev/null
@@ -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 <deal.II/sundials/kinsol.h>
+#include <deal.II/base/parameter_handler.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+// 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<double> VectorType;
+
+  SUNDIALS::KINSOL<VectorType>::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<VectorType> 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<u.size(); ++i)
+      F[i] = u[i]*u[i]-(i+1)*(i+1);
+    return 0;
+  };
+
+
+  kinsol.iteration_function = [] (const VectorType &u, VectorType &F) -> int
+  {
+    for (unsigned int i=0; i<u.size(); ++i)
+      F[i] = u[i]*u[i]-i*i - u[i];
+    return 0;
+  };
+
+  VectorType v(N);
+  v = 1.0;
+  auto niter  = kinsol.solve(v);
+  deallog << v << std::endl;
+  deallog << "Converged in " << niter << " iterations." << std::endl;
+}
diff --git a/tests/sundials/kinsol_01.output b/tests/sundials/kinsol_01.output
new file mode 100644 (file)
index 0000000..5fd595f
--- /dev/null
@@ -0,0 +1,4 @@
+
+DEAL::1.00000 2.00000 3.00000 4.00000 5.00000 6.00000 7.00000 8.00000 9.00000 10.0000 
+DEAL::
+DEAL::Converged in 23 iterations.
diff --git a/tests/sundials/kinsol_01.prm b/tests/sundials/kinsol_01.prm
new file mode 100644 (file)
index 0000000..a37daf6
--- /dev/null
@@ -0,0 +1,16 @@
+set Function norm stopping tolerance       = 0.000000
+set Maximum number of nonlinear iterations = 200
+set Scaled step stopping tolerance         = 0.000000
+set Solution strategy                      = linesearch
+subsection Fixed point and Picard parameters
+  set Anderson acceleration subspace size = 5
+end
+subsection Linesearch parameters
+  set Maximum number of beta-condition failures = 0
+end
+subsection Newton parameters
+  set Maximum allowable scaled length of the Newton step = 0.000000
+  set Maximum iterations without matrix setup            = 0
+  set No initial matrix setup                            = false
+  set Relative error for different quotient computation  = 0.000000
+end

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.