]> https://gitweb.dealii.org/ - dealii.git/commitdiff
ARKode interface.
authorLuca Heltai <luca.heltai@sissa.it>
Tue, 26 Sep 2017 14:45:49 +0000 (16:45 +0200)
committerLuca Heltai <luca.heltai@sissa.it>
Wed, 27 Sep 2017 12:34:41 +0000 (14:34 +0200)
include/deal.II/sundials/arkode.h [new file with mode: 0644]
source/sundials/CMakeLists.txt
source/sundials/arkode.cc [new file with mode: 0644]

diff --git a/include/deal.II/sundials/arkode.h b/include/deal.II/sundials/arkode.h
new file mode 100644 (file)
index 0000000..0ed5c85
--- /dev/null
@@ -0,0 +1,869 @@
+//-----------------------------------------------------------
+//
+//    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_arkode_h
+#define dealii_sundials_arkode_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 <arkode/arkode.h>
+#include <arkode/arkode_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 ARKODE error codes.
+#define AssertARKode(code) Assert(code >= 0, ExcARKodeError(code))
+
+namespace SUNDIALS
+{
+  /**
+   * Interface to SUNDIALS additive Runge-Kutta methods (ARKode).
+   *
+   * The class ARKode is a wrapper to SUNDIALS variable-step, embedded,
+   * additive Runge-Kutta solver which is a general purpose solver for systems
+   * of ordinary differential equations characterized by the presence of both
+   * fast and slow dynamics.
+   *
+   * Fast dynamics are treated implicitly, and slow dynamics are treated
+   * explicitly, using nested families of implicit and explicit Runge-Kutta
+   * solvers.
+   *
+   * Citing directly from ARKode documentation:
+   *
+   * ARKode solves ODE initial value problems (IVPs) in $R^N$. These problems
+   * should be posed in explicit form as
+   *
+   * \f[
+   *   M\dot y = f_E(t, y) + f_I (t, y), \qquad y(t_0) = y_0.
+   * \f]
+   *
+   * Here, $t$ is the independent variable (e.g. time), and the dependent
+   * variables are given by $y \in R^N$, and we use notation $\dot y$ to denote
+   * $dy/dt$. $M$ is a user-supplied nonsingular operator from $R^N \to R^N$.
+   * This operator may depend on $t$ but not on $y$.
+   *
+   * For standard systems of ordinary differential equations and for problems
+   * arising from the spatial semi-discretization of partial differential
+   * equations using finite difference or finite volume methods, $M$ is
+   * typically the identity matrix, $I$. For PDEs using a finite-element
+   * spatial semi-discretization $M$ is typically a well-conditioned mass
+   * matrix.
+   *
+   * The two right-hand side functions may be described as:
+   *
+   * - $f_E(t, y)$: contains the "slow" time scale components of the system.
+   *                This will be integrated using explicit methods.
+   * - $f_I(t, y)$: contains the "fast" time scale components of the system.
+   *                This will be integrated using implicit methods.
+   *
+   * ARKode may be used to solve stiff, nonstiff and multi-rate problems.
+   * Roughly speaking, stiffness is characterized by the presence of at least
+   * one rapidly damped mode, whose time constant is small compared to the time
+   * scale of the solution itself. In the implicit/explicit (ImEx) splitting
+   * above, these stiff components should be included in the right-hand side
+   * function $f_I (t, y)$.
+   *
+   * For multi-rate problems, a user should provide both of the functions $f_E$
+   * and $f_I$ that define the IVP system.
+   *
+   * For nonstiff problems, only $f_E$ should be provided, and $f_I$ is assumed
+   * to be zero, i.e. the system reduces to the non-split IVP:
+   *
+   * \f[
+   *   M\dot y = f_E(t, y), \qquad y(t_0) = y_0.
+   * \f]
+   *
+   * In this scenario, the ARK methods reduce to classical explicit Runge-Kutta
+   * methods (ERK). For these classes of methods, ARKode allows orders of
+   * accuracy $q = \{2, 3, 4, 5, 6, 8\}$, with embeddings of orders $p = \{1,
+   * 2, 3, 4, 5, 7\}$. These default to the Heun-Euler-2-1-2,
+   * Bogacki-Shampine-4-2-3, Zonneveld-5-3-4, Cash-Karp-6-4-5, Verner-8-5-6 and
+   * Fehlberg-13-7-8 methods, respectively.
+   *
+   * Finally, for stiff (linear or nonlinear) problems the user may provide only
+   * $f_I$, implying that $f_E = 0$, so that the system reduces to the non-split
+   * IVP
+   *
+   * \f[
+   *   M\dot y = f_I(t, y), \qquad y(t_0) = y_0.
+   * \f]
+   *
+   * Similarly to ERK methods, in this scenario the ARK methods reduce to
+   * classical diagonally-implicit Runge-Kutta methods (DIRK). For these
+   * classes of methods, ARKode allows orders of accuracy $q = \{2, 3, 4, 5\}$,
+   * with embeddings of orders $p = \{1, 2, 3, 4\}$. These default to the
+   * SDIRK-2-1-2, ARK-4-2-3 (implicit), SDIRK-5-3-4 and ARK-8-4-5 (implicit)
+   * methods, respectively.
+   *
+   * For both DIRK and ARK methods, an implicit system of the form
+   * \f[
+   *  G(z_i) := M z_i − h_n A^I_{i,i} f_I (t^I_{n,i}, z_i) − a_i = 0
+   * \f]
+   * must be solved for each stage $z_i , i = 1, \ldot, s$, where
+   * we have the data
+   * \[
+   *  a_i :=  M y_{n−1} + h_n \sum_{j=1}^{i−1} [ A^E_{i,j} f_E(t^E_{n,j}, z_j)
+   *  + A^I_{i,j} f_I (t^I_{n,j}, z_j)]
+   * \]
+   * for the ARK methods, or
+   * \[
+   *  a_i :=  M y_{n−1} + h_n \sum_{j=1}^{i−1} A^I_{i,j} f_I (t^I_{n,j}, z_j)
+   * \]
+   * for the DIRK methods. Here $A^I_{i,j}$ and $A^E_{i,j}$ are the Butcher's
+   * tables for the chosen solver.
+   *
+   * If $f_I(t,y)$ depends nonlinearly on $y$ then the systems above correspond
+   * to a nonlinear system of equations; if $f_I (t, y)$ depends linearly on
+   * $y$ then this is a linear system of equations. By specifying the flag
+   * `implicit_function_is_linear`, ARKode takes some shortcuts that allow a
+   * faster solution process.
+   *
+   * For systems of either type, ARKode allows a choice of solution strategy.
+   * The default solver choice is a variant of Newton’s method,
+   * \[
+   *  z_i^{m+1} = z_i^m +\delta^{m+1},
+   * \]
+   * where $m$ is the Newton index, and the Newton update $\delta^{m+1}$
+   * requires the solution of the linear Newton system
+   * \[
+   *  N(z_i^m) \delta^{m+1} = -G(z_i^m),
+   * \]
+   * where
+   * \[
+   * N := M - \gamma J, \quad J := \frac{\partial f_I}{\partial y},
+   * \qquad \gamma:= h_n A^I_{i,i}.
+   * \]
+   *
+   * As an alternate to Newton’s method, ARKode may solve for each stage $z_i ,i
+   * = 1, \ldots , s$ using an Anderson-accelerated fixed point iteration
+   * \[
+   * z_i^{m+1} = g(z_i^{m}), m=0,1,\ldots.
+   * \]
+   *
+   * Unlike with Newton’s method, this option does not require the solution of
+   * a linear system at each iteration, instead opting for solution of a
+   * low-dimensional least-squares solution to construct the nonlinear update.
+   *
+   * Finally, if the user specifies `implicit_function_is_linear`, i.e.,
+   * $f_I(t, y)$ depends linearly on $y$, and if the Newton-based nonlinear
+   * solver is chosen, then the system will be solved using only a single
+   * Newton iteration. Notice that in order for the Newton solver to be used,
+   * at least the solve_jacobian_system() function should be supplied. If this
+   * function is not supplied, then only the fixed-point iteration will be
+   * supported, and the `implicit_function_is_linear` setting is ignored.
+   *
+   * The optimal solver (Newton vs fixed-point) is highly problem-dependent.
+   * Since fixed-point solvers do not require the solution of any linear
+   * systems, each iteration may be significantly less costly than their Newton
+   * counterparts. However, this can come at the cost of slower convergence (or
+   * even divergence) in comparison with Newton-like methods. These fixed-point
+   * solvers do allow for user specification of the Anderson-accelerated
+   * subspace size, $m_k$. While the required amount of solver memory grows
+   * proportionately to $m_k N$, larger values of $m_k$ may result in faster
+   * convergence.
+   *
+   * This improvement may be significant even for "small" values, e.g. $1 \leq
+   * m_k \leq 5$, and convergence may not improve (or even deteriorate) for
+   * larger values of $m_k$. While ARKode uses a Newton-based iteration as its
+   * default solver due to its increased robustness on very stiff problems, it
+   * is highly recommended that users also consider the fixed-point solver for
+   * their cases when attempting a new problem.
+   *
+   * For either the Newton or fixed-point solvers, it is well-known that both
+   * the efficiency and robustness of the algorithm intimately depends on the
+   * choice of a good initial guess. In ARKode, the initial guess for either
+   * nonlinear solution method is a predicted value $z_i(0)$ that is computed
+   * explicitly from the previously-computed data (e.g. $y_{n−2}, y_{n−1}$, and
+   * $z_j$ where $j < i$). Additional information on the specific predictor
+   * algorithms implemented in ARKode is provided in ARKode documentation.
+   *
+   * The user has to provide the implementation of the following std::function:
+   *  - reinit_vector;
+   * and either one or both of
+   *  - implicit_function;
+   *  - explicit_function;
+   *
+   * If the mass matrix is different from the identity, the user should supply
+   *  - solve_mass_system;
+   * and, optionally,
+   *  - setup_mass;
+   *
+   * If the use of a Newton method is desired, then the user should also supply
+   *  - solve_jacobian_system;
+   * and optionally
+   *  - setup_jacobian;
+   *
+   * Also the following functions could be rewritten. By default
+   * they do nothing, or are not required.
+   *  - solver_should_restart;
+   *  - get_local_tolerances;
+   *
+   * To produce output at fixed steps, overload the function
+   *  - output_step;
+   *
+   *
+   * To provide a simple example, consider the harmonic oscillator problem:
+   * \f[
+   * \begin{split}
+   *   u'' & = -k^2 u \\
+   *   u (0) & = 0 \\
+   *   u'(0) & = k
+   * \end{split}
+   * \f]
+   *
+   * We write it in terms of a first order ode:
+   *\f[
+   * \begin{matrix}
+   *   y_0' & =  y_1 \\
+   *   y_1' & = - k^2 y_0
+   * \end{matrix}
+   * \f]
+   *
+   * That is $y' = A y$
+   * where
+   * \f[
+   * A:=
+   * \begin{matrix}
+   * 0 & 1 \\
+   * -k^2 &0
+   * \end{matrix}
+   * \f]
+   * and $y(0)=(0, k)$.
+   *
+   * The exact solution is $y_0(t) = \sin(k t)$, $y_1(t) = y_0'(t) = k \cos(k t)$,
+   * $y_1'(t) = -k^2 \sin(k t)$.
+   *
+   * A minimal implementation, using only explicit RK methods, is given by the
+   * following code snippet:
+   *
+   * @code
+   * typedef Vector<double> VectorType;
+   *
+   * SUNDIALS::ARKode<VectorType> ode;
+   *
+   * ode.reinit_vector = [] (VectorType&v)
+   * {
+   *   v.reinit(2);
+   * };
+   *
+   * double kappa = 1.0;
+   *
+   * ode.explicit_function = [kappa] (double,
+   *                                  const VectorType &y,
+   *                                  VectorType &ydot) -> int
+   * {
+   *   ydot[0] = y[1];
+   *   ydot[1] = -kappa*kappa*y[0];
+   *   return 0;
+   * };
+   *
+   * Vector<double> y(2);
+   * y[1] = kappa;
+   *
+   * ode.solve_ode(y);
+   * @endcode
+   *
+   * @author Luca Heltai, 2017.
+   */
+  template<typename VectorType=Vector<double> >
+  class ARKode
+  {
+  public:
+
+    /**
+     * Additional parameters that can be passed to the ARKode class.
+     */
+    class AdditionalData
+    {
+    public:
+      /**
+       * Initialization parameters for ARKode.
+       *
+       * Global parameters:
+       *
+       * @param initial_time Initial time
+       * @param final_time Final time
+       * @param initial_step_size Initial step size
+       * @param output_period Time interval between each output
+       *
+       * Running parameters:
+       *
+       * @param minimum_step_size Minimum step size
+       * @param maximum_order Maximum ARK order
+       * @param maximum_non_linear_iterations Maximum number of nonlinear iterations
+       * @param implicit_function_is_linear Specifies that the implicit portion of the problem is linear
+       * @param implicit_function_is_time_independent Specifies that the implicit portion of the problem
+       *        is linear and time independent
+       *
+       * Error parameters:
+       *
+       * @param absolute_tolerance Absolute error tolerance
+       * @param relative_tolerance Relative error tolerance
+       */
+      AdditionalData(
+        // Initial parameters
+        const double &initial_time = 0.0,
+        const double &final_time = 1.0,
+        const double &initial_step_size = 1e-2,
+        const double &output_period = 1e-1,
+        // Running parameters
+        const double &minimum_step_size = 1e-6,
+        const unsigned int &maximum_order = 5,
+        const unsigned int &maximum_non_linear_iterations = 10,
+        const bool implicit_function_is_linear = false,
+        const bool implicit_function_is_time_independent = false,
+        // Error parameters
+        const double &absolute_tolerance = 1e-6,
+        const double &relative_tolerance = 1e-5) :
+        initial_time(initial_time),
+        final_time(final_time),
+        initial_step_size(initial_step_size),
+        minimum_step_size(minimum_step_size),
+        absolute_tolerance(absolute_tolerance),
+        relative_tolerance(relative_tolerance),
+        maximum_order(maximum_order),
+        output_period(output_period),
+        maximum_non_linear_iterations(maximum_non_linear_iterations),
+        implicit_function_is_linear(implicit_function_is_linear),
+        implicit_function_is_time_independent(implicit_function_is_time_independent)
+      {};
+
+      /**
+       * 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
+       * set Final time                        = 1.000000
+       * set Initial time                      = 0.000000
+       * set Time interval between each output = 0.2
+       * subsection Error control
+       *   set Absolute error tolerance                      = 0.000001
+       *   set Ignore algebraic terms for error computations = true
+       *   set Relative error tolerance                      = 0.00001
+       *   set Use local tolerances                          = false
+       * end
+       * subsection Initial condition correction parameters
+       *   set Correction type at initial time        = none
+       *   set Correction type after restart          = none
+       *   set Maximum number of nonlinear iterations = 5
+       * end
+       * subsection Running parameters
+       *   set Initial step size                      = 0.1
+       *   set Maximum number of nonlinear iterations = 10
+       *   set Maximum order of ARK                   = 5
+       *   set Minimum step size                      = 0.000001
+       * end
+       * @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)
+      {
+        prm.add_parameter("Initial time", initial_time);
+        prm.add_parameter("Final time", final_time);
+        prm.add_parameter("Time interval between each output", output_period);
+
+        prm.enter_subsection("Running parameters");
+        prm.add_parameter("Initial step size",initial_step_size);
+        prm.add_parameter("Minimum step size", minimum_step_size);
+        prm.add_parameter("Maximum order of ARK", maximum_order);
+        prm.add_parameter("Maximum number of nonlinear iterations", maximum_non_linear_iterations);
+        prm.add_parameter("Implicit function is linear", implicit_function_is_linear);
+        prm.add_parameter("Implicit function is time independent", implicit_function_is_time_independent);
+        prm.leave_subsection();
+
+        prm.enter_subsection("Error control");
+        prm.add_parameter("Absolute error tolerance", absolute_tolerance);
+        prm.add_parameter("Relative error tolerance", relative_tolerance);
+        prm.leave_subsection();
+      }
+
+      /**
+       * Initial time for the DAE.
+       */
+      double initial_time;
+
+      /**
+       * Final time.
+       */
+      double final_time;
+
+      /**
+       * Initial step size.
+       */
+      double initial_step_size;
+
+      /**
+       * Minimum step size.
+       */
+      double minimum_step_size;
+
+      /**
+       * Absolute error tolerance for adaptive time stepping.
+       */
+      double absolute_tolerance;
+
+      /**
+       * Relative error tolerance for adaptive time stepping.
+       */
+      double relative_tolerance;
+
+      /**
+       * Maximum order of ARK.
+       */
+      unsigned int maximum_order;
+
+      /**
+       * Time period between each output.
+       */
+      double output_period;
+
+      /**
+       * Maximum number of iterations for Newton or fixed point method during
+       * time advancement.
+       */
+      unsigned int maximum_non_linear_iterations;
+
+      /**
+       * Specifies that the implicit portion of the problem is linear.
+       */
+      bool implicit_function_is_linear;
+
+      /**
+       * Specifies that the implicit portion of the problem is linear and time
+       * independent.
+       */
+      bool implicit_function_is_time_independent;
+    };
+
+    /**
+     * Constructor. It is possible to fine tune the SUNDIALS ARKode solver by
+     * passing an AdditionalData() object that sets all of the solver
+     * parameters.
+     *
+     * @param mpi_comm MPI communicator
+     * @param data ARKode configuration data
+     */
+    ARKode(const MPI_Comm mpi_comm = MPI_COMM_WORLD,
+           const AdditionalData &data=AdditionalData());
+
+    /**
+     * Destructor.
+     */
+    ~ARKode();
+
+    /**
+     * Integrate the initial value problem. This function returns the final
+     * number of computed steps.
+     */
+    unsigned int solve_ode(VectorType &solution);
+
+    /**
+     * Clear internal memory and start with clean objects. This function is
+     * called when the simulation starts and when the user returns true to a
+     * call to solver_should_restart().
+     *
+     * By default solver_should_restart() returns false. If the user needs to
+     * implement, for example, local adaptivity in space, he or she may assign
+     * a different function to solver_should_restart() that performs all mesh
+     * changes, transfers the solution to the new mesh, and returns true.
+     *
+     * @param[in] t  The new starting time
+     * @param[in] h  The new starting time step
+     * @param[in,out] y   The new initial solution
+     */
+    void reset(const double &t,
+               const double &h,
+               const VectorType &y);
+
+    /**
+     * 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 may supply and that is intended to compute
+     * the explicit part of the IVP right hand side. Sets $explicit_f = f_E(t,
+     * y)$.
+     *
+     * At least one of explicit_function() or implicit_function() must be
+     * provided. According to which one is provided, explicit, implicit, or
+     * mixed RK methods are used.
+     *
+     * This function should return:
+     * - 0: Success
+     * - >0: Recoverable error (ARKodeReinit 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 double t,
+                      const VectorType &y,
+                      VectorType &explicit_f)> explicit_function;
+
+    /**
+     * A function object that users may supply and that is intended to compute
+     * the implicit part of the IVP right hand side. Sets $implicit_f = f_I(t,
+     * y)$.
+     *
+     * At least one of explicit_function() or implicit_function() must be
+     * provided. According to which one is provided, explicit, implicit, or
+     * mixed RK methods are used.
+     *
+     * This function should return:
+     * - 0: Success
+     * - >0: Recoverable error (ARKodeReinit 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 double t,
+                      const VectorType &y,
+                      VectorType &res)> implicit_function;
+
+    /**
+     * A function object that users may supply and that is intended to
+     * prepare the linear solver for subsequent calls to
+     * solve_jacobian_system().
+     *
+     * Make sure that after a call to this function, we know how to compute
+     * solutions of systems $A x = b$, where $A$ is some approximation to the
+     * Newton matrix, $M − \gamma \partial f_I/\partial y$. This function is
+     * optional. If the user does not provide it, then solve_jacobian_system()
+     * is assumed to also perform the setup internally.
+     *
+     * The setup_jacobian() function may call a user-supplied function to
+     * compute needed data related to the Jacobian matrix. Alterntively, it may
+     * choose to retrieve and use stored values of this data. In either case,
+     * setup_jacobian() 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
+     * stage solve (or even every time step), 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 time steps.
+     *
+     * If the user uses a matrix based computation of the Jacobian, then this
+     * is the right place where an assembly routine shoulde be called to
+     * assemble both a matrix and a preconditioner for the Jacobian system.
+     * Subsequent calls (possibly more than one) to solve_jacobian_system() can
+     * assume that this function has been called at least once.
+     *
+     * Notice that no assumption is made by this interface on what the user
+     * should do in this function. ARKode only assumes that after a call to
+     * setup_jacobian() it is possible to call solve_jacobian_system(), to
+     * obtain a solution $x$ to the system $J x = b$. If this function is not
+     * provided, then it is never called.
+     *
+     * 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] ypred  is the predicted $y$ vector for the current ARKode internal step
+     * @param[in] fpred  is the value of the implicit right-hand side at ypred,
+     *        $f_I (t_n, ypred)$.
+     *
+     * @param[in] convfail – an input flag used to indicate any problem that occurred
+     *   during the solution of the nonlinear equation on the current time step
+     *   for which the linear solver is being used. This flag can be used to help
+     *   decide whether the Jacobian data kept by a linear solver needs to be
+     *   updated or not. Its possible values are:
+     *
+     *   - ARK_NO_FAILURES: this value is passed if either this is the first call
+     *     for this step, or the local error test failed on the previous attempt at
+     *     this step (but the Newton iteration converged).
+     *
+     *   - ARK_FAIL_BAD_J: this value is passed if (a) the previous Newton
+     *     corrector iteration did not converge and the linear solver's setup
+     *     function indicated that its Jacobian-related data is not current, or (b)
+     *     during the previous Newton corrector iteration, the linear solver's
+     *     solve function failed in a recoverable manner and the linear solver's
+     *     setup function indicated that its Jacobian-related data is not current.
+     *
+     *   - ARK_FAIL_OTHER: this value is passed if during the current internal
+     *     step try, the previous Newton iteration failed to converge even though
+     *     the linear solver was using current Jacobian-related data.
+     *
+     * @param[out] j_is_current: a boolean to be filled in by setup_jacobian(). The value
+     * should be set to `true` if the Jacobian data is current after the call,
+     * and should be set set to `false` if its Jacobian data is not current. If
+     * setup_jacobian() calls for re-evaluation of Jacobian data (based on
+     * convfail and ARKode state data), then it should set `j_is_current` to
+     * `true` unconditionally, otherwise an infinite loop can result.
+     *
+     * This function should return:
+     * - 0: Success
+     * - >0: Recoverable error (ARKodeReinit 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 int convfail,
+                      const double t,
+                      const double gamma,
+                      const VectorType &ypred,
+                      const VectorType &fpred,
+                      bool &j_is_current)> 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 ARKode
+     * (possibly several times) after setup_jacobian() has been called at least
+     * once. ARKode tries to do its best to call setup_jacobian() the minimum
+     * amount of times. If convergence can be achieved without updating the
+     * Jacobian, then ARKode does not call setup_jacobian() again. If, on the
+     * contrary, internal ARKode convergence tests fail, then ARKode 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 ARKode 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 (ARKodeReinit 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 double t,
+                      const double gamma,
+                      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 setup
+     * the mass matrix. This function is called by ARKode any time a mass
+     * matrix update is required. The user should compute the mass matrix (or
+     * update all the variables that allow the application of the mass matrix).
+     * This function is called by ARKode once, before any call to
+     * solve_mass_system().
+     *
+     * ARKode supports the case where the mass matrix may depend on time, but
+     * not the case where the mass matrix depends on the solution itself.
+     *
+     * If the user does not provide a solve_mass_matrix() function, then the
+     * identity is used. If the setup_mass() function is not provided, then
+     * solve_mass_system() should do all the work by itself.
+     *
+     * If the user uses a matrix based computation of the mass matrix, then
+     * this is the right place where an assembly routine shoulde be called to
+     * assemble both a matrix and a preconditioner. Subsequent calls (possibly
+     * more than one) to solve_mass_system() can assume that this function
+     * has been called at least once.
+     *
+     * Notice that no assumption is made by this interface on what the user
+     * should do in this function. ARKode only assumes that after a call to
+     * setup_mass() it is possible to call solve_mass_system(), to
+     * obtain a solution $x$ to the system $M x = b$.
+     *
+     * This function should return:
+     * - 0: Success
+     * - >0: Recoverable error (ARKodeReinit 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 double t)> setup_mass;
+
+    /**
+     * A function object that users may supply and that is intended to solve
+     * the mass matrix linear system. This function will be called by ARKode
+     * (possibly several times) after setup_mass() has been called at least
+     * once. ARKode tries to do its best to call setup_mass() the minimum
+     * amount of times.
+     *
+     * A call to this function should store in `dst` the result of $M^{-1}$
+     * applied to `src`, i.e., `M*dst = src`. It is the users responsability to
+     * set up proper solvers and preconditioners inside this function.
+     *
+     * This function should return:
+     * - 0: Success
+     * - >0: Recoverable error (ARKodeReinit 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 &rhs, VectorType &dst)> solve_mass_system;
+
+    /**
+     * A function object that users may supply and that is intended to
+     * postprocess the solution. This function is called by ARKode at fixed
+     * time increments (every `output_period` seconds), and it is passed a
+     * polynomial interpolation of the solution, computed using the current ARK
+     * order and the (internally stored) previously computed solution steps.
+     *
+     * Notice that it is well possible that internally ARKode computes a time step
+     * which is much larger than the `output_period` step, and therefore calls
+     * this function consecutively several times by simply performing all
+     * intermediate interpolations. There is no relationship between how many
+     * times this function is called and how many time steps have actually been
+     * computed.
+     */
+    std::function<void(const double t,
+                       const VectorType &sol,
+                       const unsigned int step_number)> output_step;
+
+    /**
+     * A function object that users may supply and that is intended to evaluate
+     * wether the solver should be restarted (for example because the number of
+     * degrees of freedom has changed).
+     *
+     * This function is supposed to perform all operations that are necessary
+     * in `sol` to make sure that the resulting vectors are consistent, and of
+     * the correct final size.
+     *
+     * For example, one may decide that a local refinement is necessary at time
+     * t. This function should then return true, and change the dimension of
+     * `sol` to reflect the new dimension. Since ARKode does not know about the
+     * new dimension, an internal reset is necessary.
+     *
+     * The default implementation simply returns `false`, i.e., no restart is
+     * performed during the evolution.
+     */
+    std::function<bool (const double t,
+                        VectorType &sol)> solver_should_restart;
+
+    /**
+     * A function object that users may supply and that is intended to return a
+     * vector whose components are the weights used by ARKode to compute the
+     * vector norm. The implementation of this function is optional, and it is
+     * used only if implemented.
+     */
+    std::function<VectorType&()> get_local_tolerances;
+
+    /**
+     * Handle ARKode exceptions.
+     */
+    DeclException1(ExcARKodeError, int, << "One of the SUNDIALS ARKode 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();
+
+    /**
+     * ARKode configuration data.
+     */
+    AdditionalData data;
+
+    /**
+     * ARKode memory object.
+     */
+    void *arkode_mem;
+
+    /**
+     * ARKode solution vector.
+     */
+    N_Vector yy;
+
+    /**
+     * ARKode absolute tolerances vector.
+     */
+    N_Vector abs_tolls;
+
+#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 437bf07466dc06b9f4909939c0c6d0d0d6302f6f..238e7370d562690f8e56284f47cb702fbc80bfd6 100644 (file)
@@ -16,6 +16,7 @@
 INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR})
 
 SET(_src
+  arkode.cc
   ida.cc
   copy.cc
   )
diff --git a/source/sundials/arkode.cc b/source/sundials/arkode.cc
new file mode 100644 (file)
index 0000000..759bfe8
--- /dev/null
@@ -0,0 +1,488 @@
+//-----------------------------------------------------------
+//
+//    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/arkode.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 <iostream>
+#include <iomanip>
+#include <arkode/arkode_impl.h>
+
+// Make sure we know how to call sundials own ARKode() function
+const auto &SundialsARKode = ARKode;
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace SUNDIALS
+{
+  using namespace internal;
+
+  namespace
+  {
+
+    template<typename VectorType>
+    int t_arkode_explicit_function(realtype tt,
+                                   N_Vector yy,
+                                   N_Vector yp,
+                                   void *user_data)
+    {
+      ARKode<VectorType> &solver = *static_cast<ARKode<VectorType> *>(user_data);
+      GrowingVectorMemory<VectorType> mem;
+
+      typename VectorMemory<VectorType>::Pointer src_yy(mem);
+      solver.reinit_vector(*src_yy);
+
+      typename VectorMemory<VectorType>::Pointer dst_yp(mem);
+      solver.reinit_vector(*dst_yp);
+
+      copy(*src_yy, yy);
+
+      int err = solver.explicit_function(tt, *src_yy, *dst_yp);
+
+      copy(yp, *dst_yp);
+
+      return err;
+    }
+
+
+
+    template<typename VectorType>
+    int t_arkode_implicit_function(realtype tt,
+                                   N_Vector yy,
+                                   N_Vector yp,
+                                   void *user_data)
+    {
+      ARKode<VectorType> &solver = *static_cast<ARKode<VectorType> *>(user_data);
+      GrowingVectorMemory<VectorType> mem;
+
+      typename VectorMemory<VectorType>::Pointer src_yy(mem);
+      solver.reinit_vector(*src_yy);
+
+      typename VectorMemory<VectorType>::Pointer dst_yp(mem);
+      solver.reinit_vector(*dst_yp);
+
+      copy(*src_yy, yy);
+
+      int err = solver.implicit_function(tt, *src_yy, *dst_yp);
+
+      copy(yp, *dst_yp);
+
+      return err;
+    }
+
+
+
+    template<typename VectorType>
+    int t_arkode_setup_jacobian(ARKodeMem arkode_mem,
+                                int convfail,
+                                N_Vector ypred,
+                                N_Vector fpred,
+                                booleantype *jcurPtr,
+                                N_Vector,
+                                N_Vector,
+                                N_Vector)
+    {
+      ARKode<VectorType> &solver = *static_cast<ARKode<VectorType> *>(arkode_mem->ark_user_data);
+      GrowingVectorMemory<VectorType> mem;
+
+      typename VectorMemory<VectorType>::Pointer src_ypred(mem);
+      solver.reinit_vector(*src_ypred);
+
+      typename VectorMemory<VectorType>::Pointer src_fpred(mem);
+      solver.reinit_vector(*src_fpred);
+
+      copy(*src_ypred, ypred);
+      copy(*src_fpred, fpred);
+
+      int err = solver.setup_jacobian(convfail,
+                                      arkode_mem->ark_tn,
+                                      arkode_mem->ark_gamma,
+                                      *src_ypred,
+                                      *src_fpred,
+                                      (bool &)*jcurPtr);
+
+      return err;
+    }
+
+
+
+    template<typename VectorType>
+    int t_arkode_solve_jacobian(ARKodeMem arkode_mem,
+                                N_Vector b,
+                                N_Vector,
+                                N_Vector ycur,
+                                N_Vector fcur)
+    {
+      ARKode<VectorType> &solver = *static_cast<ARKode<VectorType> *>(arkode_mem->ark_user_data);
+      GrowingVectorMemory<VectorType> mem;
+
+      typename VectorMemory<VectorType>::Pointer src(mem);
+      solver.reinit_vector(*src);
+
+      typename VectorMemory<VectorType>::Pointer src_ycur(mem);
+      solver.reinit_vector(*src_ycur);
+
+      typename VectorMemory<VectorType>::Pointer src_fcur(mem);
+      solver.reinit_vector(*src_fcur);
+
+      typename VectorMemory<VectorType>::Pointer dst(mem);
+      solver.reinit_vector(*dst);
+
+      copy(*src, b);
+      copy(*src_ycur, ycur);
+      copy(*src_fcur, fcur);
+
+      int err = solver.solve_jacobian_system(arkode_mem->ark_tn,
+                                             arkode_mem->ark_gamma,
+                                             *src_ycur, *src_fcur,
+                                             *src,*dst);
+      copy(b, *dst);
+
+      return err;
+    }
+
+
+
+    template<typename VectorType>
+    int t_arkode_setup_mass(ARKodeMem arkode_mem,
+                            N_Vector,
+                            N_Vector,
+                            N_Vector)
+    {
+      ARKode<VectorType> &solver = *static_cast<ARKode<VectorType> *>(arkode_mem->ark_user_data);
+      int err = solver.setup_mass(arkode_mem->ark_tn);
+      return err;
+    }
+
+
+
+    template<typename VectorType>
+    int t_arkode_solve_mass(ARKodeMem arkode_mem,
+                            N_Vector b,
+                            N_Vector)
+    {
+      ARKode<VectorType> &solver = *static_cast<ARKode<VectorType> *>(arkode_mem->ark_user_data);
+      GrowingVectorMemory<VectorType> mem;
+
+      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_mass_system(*src,*dst);
+      copy(b, *dst);
+
+      return err;
+    }
+  }
+
+  template <typename VectorType>
+  ARKode<VectorType>::ARKode(const MPI_Comm mpi_comm, const AdditionalData &data) :
+    data(data),
+    arkode_mem(nullptr),
+    communicator(Utilities::MPI::duplicate_communicator(mpi_comm))
+  {
+    set_functions_to_trigger_an_assert();
+  }
+
+  template <typename VectorType>
+  ARKode<VectorType>::~ARKode()
+  {
+    if (arkode_mem)
+      ARKodeFree(&arkode_mem);
+    MPI_Comm_free(&communicator);
+  }
+
+
+
+  template <typename VectorType>
+  unsigned int ARKode<VectorType>::solve_ode(VectorType &solution)
+  {
+
+    unsigned int system_size = solution.size();
+    unsigned int local_system_size = system_size;
+
+    double t = data.initial_time;
+    double h = data.initial_step_size;
+    unsigned int step_number = 0;
+
+    int status;
+    (void)status;
+
+    // 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 = solution.locally_owned_elements();
+        local_system_size = is.n_elements();
+
+        yy        = N_VNew_Parallel(communicator,
+                                    local_system_size,
+                                    system_size);
+
+        abs_tolls = N_VNew_Parallel(communicator,
+                                    local_system_size,
+                                    system_size);
+      }
+    else
+#endif
+      {
+        Assert(is_serial_vector<VectorType>::value,
+               ExcInternalError("Trying to use a serial code with a parallel vector."));
+        yy        = N_VNew_Serial(system_size);
+        abs_tolls = N_VNew_Serial(system_size);
+      }
+    reset(data.initial_time,
+          data.initial_step_size,
+          solution);
+
+    double next_time = data.initial_time;
+
+    output_step( 0, solution, 0);
+
+    while (t<data.final_time)
+      {
+
+        next_time += data.output_period;
+
+        status = SundialsARKode(arkode_mem, next_time, yy, &t, ARK_NORMAL);
+
+        AssertARKode(status);
+
+        status = ARKodeGetLastStep(arkode_mem, &h);
+        AssertARKode(status);
+
+        copy(solution, yy);
+
+        while (solver_should_restart(t, solution))
+          reset(t, h, solution);
+
+        step_number++;
+
+        if (output_step)
+          output_step(t, solution, step_number);
+      }
+
+    // Free the vectors which are no longer used.
+#ifdef DEAL_II_WITH_MPI
+    if (is_serial_vector<VectorType>::value == false)
+      {
+        N_VDestroy_Parallel(yy);
+        N_VDestroy_Parallel(abs_tolls);
+      }
+    else
+#endif
+      {
+        N_VDestroy_Serial(yy);
+        N_VDestroy_Serial(abs_tolls);
+      }
+
+    return step_number;
+  }
+
+  template <typename VectorType>
+  void ARKode<VectorType>::reset(const double &current_time,
+                                 const double &current_time_step,
+                                 const VectorType &solution)
+  {
+
+    unsigned int system_size;
+    unsigned int local_system_size;
+
+    if (arkode_mem)
+      ARKodeFree(&arkode_mem);
+
+    arkode_mem = ARKodeCreate();
+
+    // Free the vectors which are no longer used.
+    if (yy)
+      {
+#ifdef DEAL_II_WITH_MPI
+        if (is_serial_vector<VectorType>::value == false)
+          {
+            N_VDestroy_Parallel(yy);
+            N_VDestroy_Parallel(abs_tolls);
+          }
+        else
+#endif
+          {
+            N_VDestroy_Serial(yy);
+            N_VDestroy_Serial(abs_tolls);
+          }
+      }
+
+    int status;
+    (void)status;
+    system_size = solution.size();
+#ifdef DEAL_II_WITH_MPI
+    if (is_serial_vector<VectorType>::value == false)
+      {
+        IndexSet is = solution.locally_owned_elements();
+        local_system_size = is.n_elements();
+
+        yy        = N_VNew_Parallel(communicator,
+                                    local_system_size,
+                                    system_size);
+
+        abs_tolls = N_VNew_Parallel(communicator,
+                                    local_system_size,
+                                    system_size);
+
+      }
+    else
+#endif
+      {
+        yy        = N_VNew_Serial(system_size);
+        abs_tolls = N_VNew_Serial(system_size);
+      }
+
+    copy(yy, solution);
+
+    Assert(explicit_function || implicit_function,
+           ExcFunctionNotProvided("explicit_function || implicit_function"));
+
+    status = ARKodeInit(arkode_mem,
+                        explicit_function ? t_arkode_explicit_function<VectorType> : nullptr,
+                        implicit_function ? t_arkode_implicit_function<VectorType> : nullptr,
+                        current_time, yy);
+    AssertARKode(status);
+
+    if (get_local_tolerances)
+      {
+        copy(abs_tolls, get_local_tolerances());
+        status = ARKodeSVtolerances(arkode_mem, data.relative_tolerance, abs_tolls);
+        AssertARKode(status);
+      }
+    else
+      {
+        status = ARKodeSStolerances(arkode_mem, data.relative_tolerance, data.absolute_tolerance);
+        AssertARKode(status);
+      }
+
+    status = ARKodeSetInitStep(arkode_mem, current_time_step);
+    AssertARKode(status);
+
+    status = ARKodeSetUserData(arkode_mem, (void *) this);
+    AssertARKode(status);
+
+    status = ARKodeSetStopTime(arkode_mem, data.final_time);
+    AssertARKode(status);
+
+    status = ARKodeSetMaxNonlinIters(arkode_mem, data.maximum_non_linear_iterations);
+    AssertARKode(status);
+
+    // Initialize solver
+    ARKodeMem ARKode_mem = (ARKodeMem) arkode_mem;
+
+    if (solve_jacobian_system)
+      {
+        status = ARKodeSetNewton(arkode_mem);
+        AssertARKode(status);
+        if (data.implicit_function_is_linear)
+          {
+            status = ARKodeSetLinear(arkode_mem,
+                                     data.implicit_function_is_time_independent ? 0 : 1);
+            AssertARKode(status);
+          }
+
+
+        ARKode_mem->ark_lsolve = t_arkode_solve_jacobian<VectorType>;
+        if (setup_jacobian)
+          {
+            ARKode_mem->ark_lsetup = t_arkode_setup_jacobian<VectorType>;
+            ARKode_mem->ark_setupNonNull = true;
+          }
+      }
+    else
+      {
+        status = ARKodeSetFixedPoint(arkode_mem, data.maximum_non_linear_iterations);
+        AssertARKode(status);
+      }
+
+
+    if (solve_mass_system)
+      {
+        ARKode_mem->ark_msolve = t_arkode_solve_mass<VectorType>;
+
+        if (setup_mass)
+          {
+            ARKode_mem->ark_msetup = t_arkode_setup_mass<VectorType>;
+            ARKode_mem->ark_MassSetupNonNull = true;
+          }
+      }
+
+    status = ARKodeSetOrder(arkode_mem, data.maximum_order);
+    AssertARKode(status);
+  }
+
+  template<typename VectorType>
+  void ARKode<VectorType>::set_functions_to_trigger_an_assert()
+  {
+
+    reinit_vector = [](VectorType &)
+    {
+      AssertThrow(false, ExcFunctionNotProvided("reinit_vector"));
+    };
+
+    solver_should_restart = [](const double,
+                               VectorType &) ->bool
+    {
+      return false;
+    };
+  }
+
+  template class ARKode<Vector<double> >;
+  template class ARKode<BlockVector<double> >;
+
+#ifdef DEAL_II_WITH_MPI
+
+#ifdef DEAL_II_WITH_TRILINOS
+  template class ARKode<TrilinosWrappers::MPI::Vector>;
+  template class ARKode<TrilinosWrappers::MPI::BlockVector>;
+#endif
+
+#ifdef DEAL_II_WITH_PETSC
+  template class ARKode<PETScWrappers::MPI::Vector>;
+  template class ARKode<PETScWrappers::MPI::BlockVector>;
+#endif
+
+#endif
+
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif

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.