]> https://gitweb.dealii.org/ - dealii.git/commitdiff
PETScWrappers: Add support for ODE solver 15066/head
authorStefano Zampini <stefano.zampini@gmail.com>
Tue, 17 May 2022 15:45:34 +0000 (18:45 +0300)
committerStefano Zampini <stefano.zampini@gmail.com>
Thu, 13 Apr 2023 08:11:20 +0000 (11:11 +0300)
12 files changed:
include/deal.II/base/exceptions.h
include/deal.II/lac/petsc_ts.h [new file with mode: 0644]
include/deal.II/lac/petsc_ts.templates.h [new file with mode: 0644]
source/lac/CMakeLists.txt
source/lac/petsc_ts.cc [new file with mode: 0644]
tests/petsc/petsc_ts_00.cc [new file with mode: 0644]
tests/petsc/petsc_ts_00.output [new file with mode: 0644]
tests/petsc/petsc_ts_00_in.prm [new file with mode: 0644]
tests/petsc/petsc_ts_01.cc [new file with mode: 0644]
tests/petsc/petsc_ts_01.output [new file with mode: 0644]
tests/petsc/petsc_ts_02.cc [new file with mode: 0644]
tests/petsc/petsc_ts_02.output [new file with mode: 0644]

index f59caf18523725842a7c0af6f8a9cdcf113a6a4d..28745c7fc93039f28af872713016a3db1cc21800 100644 (file)
@@ -860,6 +860,14 @@ namespace StandardExceptions
                  << "Please provide an implementation for the function \""
                  << arg1 << "\"");
 
+  /**
+   * This exception is used if some user function returns nonzero exit codes.
+   */
+  DeclException2(ExcFunctionNonzeroReturn,
+                 std::string,
+                 int,
+                 << "Function \"" << arg1 << "\" returned " << arg2);
+
   /**
    * This exception is used if some object is found uninitialized.
    */
diff --git a/include/deal.II/lac/petsc_ts.h b/include/deal.II/lac/petsc_ts.h
new file mode 100644 (file)
index 0000000..a5f6153
--- /dev/null
@@ -0,0 +1,592 @@
+//-----------------------------------------------------------
+//
+//    Copyright (C) 2022 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.md at
+//    the top level directory of deal.II.
+//
+//---------------------------------------------------------------
+//
+// Author: Stefano Zampini, King Abdullah University of Science and Technology.
+
+#ifndef dealii_petsc_ts_h
+#define dealii_petsc_ts_h
+
+#include <deal.II/base/config.h>
+
+#ifdef DEAL_II_WITH_PETSC
+#  include <deal.II/base/mpi.h>
+#  include <deal.II/base/parameter_handler.h>
+#  include <deal.II/base/smartpointer.h>
+
+#  include <deal.II/lac/petsc_matrix_base.h>
+#  include <deal.II/lac/petsc_precondition.h>
+#  include <deal.II/lac/petsc_vector_base.h>
+
+#  include <petscts.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace PETScWrappers
+{
+  /**
+   * Additional parameters that can be passed to the TimeStepper class.
+   */
+  class TimeStepperData
+  {
+  public:
+    /**
+     * Type that holds real-valued numbers.
+     *
+     * Used to represent time and norms tolerances.
+     */
+    using real_type = PetscReal;
+
+    /**
+     * Initialization parameters for TimeStepper.
+     *
+     * Running parameters:
+     *
+     * @param opts_prefix The string indicating the options prefix for command line customization.
+     * @param tstype The string indicating the PETSc solver type.
+     * @param initial_time Initial simulation time.
+     * @param final_time Final simulation time.
+     * @param initial_step_size Initial step size.
+     * @param max_steps Maximum number of steps allowed.
+     * @param match_step Whether or not to exactly stop at final time or step over it.
+     *
+     * Error parameters:
+     *
+     * @param tsadapttype The string indicating the PETSc time step adaptor type.
+     * @param minimum_step_size Minimum step size allowed.
+     * @param maximum_step_size Maximum step size allowed.
+     * @param absolute_tolerance Absolute error tolerance.
+     * @param relative_tolerance Relative error tolerance.
+     * @param ignore_algebraic_lte Ignore algebraic terms for
+     * error computations
+     *
+     * Note that one between `final_time` or `max_steps` must
+     * be specified by the user, otherwise PETSc will complain.
+     * Adaptive time stepping is disabled by default.
+     * Negative values indicate using PETSc's default.
+     *
+     * Note that all parameters values specified here can be overriden by
+     * command line choices.
+     */
+    TimeStepperData(
+      // Running parameters
+      const std::string &opts_prefix       = "",
+      const std::string &tstype            = "",
+      const real_type    initial_time      = 0.0,
+      const real_type    final_time        = 0.0,
+      const real_type    initial_step_size = 0.0,
+      const int          max_steps         = -1,
+      const bool         match_step        = false,
+      // Error parameters
+      const std::string &tsadapttype          = "none",
+      const real_type    minimum_step_size    = -1.0,
+      const real_type    maximum_step_size    = -1.0,
+      const real_type    absolute_tolerance   = -1.0,
+      const real_type    relative_tolerance   = -1.0,
+      const bool         ignore_algebraic_lte = true)
+      : opts_prefix(opts_prefix)
+      , tstype(tstype)
+      , initial_time(initial_time)
+      , final_time(final_time)
+      , initial_step_size(initial_step_size)
+      , max_steps(max_steps)
+      , match_step(match_step)
+      , tsadapttype(tsadapttype)
+      , minimum_step_size(minimum_step_size)
+      , maximum_step_size(maximum_step_size)
+      , absolute_tolerance(absolute_tolerance)
+      , relative_tolerance(relative_tolerance)
+      , ignore_algebraic_lte(ignore_algebraic_lte)
+    {}
+
+    void
+    add_parameters(ParameterHandler &prm)
+    {
+      prm.enter_subsection("Running parameters");
+      prm.add_parameter(
+        "options prefix",
+        opts_prefix,
+        "The string indicating the options prefix for command line customization.");
+      prm.add_parameter("solver type",
+                        tstype,
+                        "The string indicating the PETSc TS type.");
+      prm.add_parameter("initial time",
+                        initial_time,
+                        "The value for the initial time.");
+      prm.add_parameter("final time",
+                        final_time,
+                        "The value for the final time.");
+      prm.add_parameter("initial step size",
+                        initial_step_size,
+                        "The value for the initial time step.");
+      prm.add_parameter("maximum number of steps",
+                        max_steps,
+                        "Maximum number of time steps allowed.");
+      prm.add_parameter(
+        "match final time",
+        match_step,
+        "Whether or not to exactly stop at final time or step over it.");
+      prm.leave_subsection();
+
+      prm.enter_subsection("Error control");
+      prm.add_parameter("adaptor type",
+                        tsadapttype,
+                        "The string for the TSAdapt type.");
+      prm.add_parameter("minimum step size",
+                        minimum_step_size,
+                        "Minimum time step size allowed.");
+      prm.add_parameter("maximum step size",
+                        maximum_step_size,
+                        "Maximum time step size allowed.");
+      prm.add_parameter("absolute error tolerance",
+                        absolute_tolerance,
+                        "Absolute error tolerance.");
+      prm.add_parameter("relative error tolerance",
+                        relative_tolerance,
+                        "Absolute error tolerance.");
+      prm.add_parameter(
+        "ignore algebraic lte",
+        ignore_algebraic_lte,
+        "Indicate whether or not to suppress algebraic variables "
+        "in the local truncation error test.");
+      prm.leave_subsection();
+    }
+
+    /**
+     * Options database prefix.
+     */
+    std::string opts_prefix;
+
+    /**
+     * PETSc solver type.
+     */
+    std::string tstype;
+
+    /**
+     * Initial time for the DAE.
+     */
+    real_type initial_time;
+
+    /**
+     * Final time.
+     */
+    real_type final_time;
+
+    /**
+     * Initial step size.
+     * Non-positive values are ignored.
+     */
+    real_type initial_step_size;
+
+    /**
+     * Maximum number of steps to be taken.
+     * Negative values are ignored.
+     */
+    int max_steps;
+
+    /**
+     * Match final time requested?
+     */
+    bool match_step;
+
+    /**
+     * PETSc time step adaptor type.
+     */
+    std::string tsadapttype;
+
+    /**
+     * Minimum allowed step size for adaptive time stepping.
+     * Non-positive values indicate to use PETSc's default.
+     */
+    real_type minimum_step_size;
+
+    /**
+     * Maximum allowed step size for adaptive time stepping.
+     * Non-positive values indicate to use PETSc's default.
+     */
+    real_type maximum_step_size;
+
+    /**
+     * Absolute error tolerance for adaptive time stepping.
+     * Negative values indicate to use PETSc's default.
+     */
+    real_type absolute_tolerance;
+
+    /**
+     * Relative error tolerance for adaptive time stepping.
+     * Negative values indicate to use PETSc's default.
+     */
+    real_type relative_tolerance;
+
+    /**
+     * Ignore algebraic terms for local truncation error.
+     */
+    bool ignore_algebraic_lte;
+  };
+
+  /**
+   * Interface to the PETSc TS solver for Ordinary Differential Equations
+   * and Differential-Algebraic Equations. The TS solver is described in the
+   * [PETSc manual](https://petsc.org/release/manual/ts/).
+   *
+   * This class supports two kinds of formulations.
+   * The explicit formulation:
+   * \f[
+   *   \begin{cases}
+   *       \dot y = G(t,y)\, , \\
+   *       y(t_0) = y_0\, , \\
+   *   \end{cases}
+   * \f]
+   * and the implicit formulation:
+   * \f[
+   *   \begin{cases}
+   *       F(t,y,\dot y) = 0\, , \\
+   *       y(t_0) = y_0\, . \\
+   *   \end{cases}
+   * \f]
+   *
+   * The interface to PETSc is realized by means of std::function callbacks
+   * like in the SUNDIALS::IDA and SUNDIALS::ARKode classes.
+   *
+   * TimeStepper supports any vector and matrix type having constructors and
+   * methods:
+   *
+   * @code
+   * class VectorType : public Subscriptor
+   *    ...
+   *    explicit VectorType(Vec);
+   *    ...
+   *    Vec & petsc_vector();
+   *    ...
+   * @endcode
+   *
+   * @code
+   * class MatrixType : public Subscriptor
+   *    ...
+   *    explicit MatrixType(Mat);
+   *    ...
+   *    Mat & petsc_matrix();
+   *    ...
+   * @endcode
+   *
+   * In particular, the supported types are the ones that can *wrap*
+   * PETSc's native vector and matrix classes, that are able to modify
+   * them in place, and that can return PETSc native types when requested.
+   *
+   * To use explicit solvers (like for example explicit Runge-Kutta methods),
+   * the user only needs to provide the implementation of $G$ via the
+   * TimeStepper::explicit_function. For implicit solvers, users have also the
+   * alternative of providing the $F$ function via
+   * TimeStepper::implicit_function. IMEX methods are also supported by
+   * providing both callbacks.
+   *
+   * The default linearization procedure of an implicit solver instantiated with
+   * this class consists in using Jacobian-Free-Newton-Krylov; the action of
+   * tangent matrices inside a linear solver process are approximated via
+   * matrix-free finite-differencing of the nonlinear residual equations that
+   * are ODE-solver specific. For details, consult the [PETSc
+   * manual](https://petsc.org/release/manual/snes/#sec-nlmatrixfree).
+   *
+   * In alternative, users can also provide the implementations of the
+   * *Jacobians*. This can be accomplished in two ways:
+   *  - a-la-PETSc style using TimeStepper::implicit_jacobian
+   *    and TimeStepper::explicit_jacobian.
+   *  - a-la-Deal.II style using TimeStepper::setup_jacobian and
+   * TimeStepper::solve_for_jacobian_system
+   *
+   * In case both approaches are coded, the deal.II style
+   * will be used.
+   *
+   * The second approach is more in style with the deal.II philosophy
+   * and it also allows command line customization, like for example,
+   * @code
+   * ./myApp -snes_type newtontr -ksp_type cg
+   * @endcode
+   * in case the user wants to change the default nonlinear solver to
+   * a trust region solver and iterate on the tangent system with CG,
+   * still using TimeStepper::solve_for_jacobian_system as a preconditioner.
+   *
+   * The first approach has instead the advantage that only the matrix assembly
+   * procedure has to be coded, thus allowing quicker implementations and
+   * faster turnaround for experimenting with linear solver preconditioning
+   * configurations via command line customizations, like for example,
+   * @code
+   * ./myApp -ksp_type cg -pc_type gamg
+   * @endcode
+   */
+  template <typename VectorType  = PETScWrappers::VectorBase,
+            typename PMatrixType = PETScWrappers::MatrixBase,
+            typename AMatrixType = PMatrixType>
+  class TimeStepper
+  {
+  public:
+    /**
+     * Type that holds real-valued numbers.
+     *
+     * Used to represent time and norms tolerances.
+     */
+    using real_type = PetscReal;
+
+    /**
+     * Constructor.
+     */
+    TimeStepper(const TimeStepperData &data     = TimeStepperData(),
+                const MPI_Comm &       mpi_comm = PETSC_COMM_WORLD);
+
+    /**
+     * Destructor.
+     */
+    ~TimeStepper();
+
+    /**
+     * Return the PETSc TS object.
+     */
+    TS
+    petsc_ts();
+
+    /**
+     * Return the underlying MPI communicator.
+     */
+    MPI_Comm
+    get_mpi_communicator() const;
+
+    /**
+     * Reset the solver, it does not change the customization.
+     */
+    void
+    reinit();
+
+    /**
+     * Reset solver.
+     * Change customization according to `data`.
+     */
+    void
+    reinit(const TimeStepperData &data);
+
+    /**
+     * Set the preconditioning matrix only.
+     *
+     * When used with TimeStepper::setup_jacobian and
+     * TimeStepper::solve_for_jacobian_system, PETSc will approximate the linear
+     * system matrix-vector product using an internal matrix-free
+     * representation.
+     *
+     * When used with TimeStepper::implicit_jacobian or
+     * TimeStepper::explicit_jacobian, PETSc will use the same matrix for both
+     * preconditioning and matrix-vector products.
+     */
+    void
+    reinit_matrix(PMatrixType &P);
+
+    /**
+     * Set both the linear system matrix and the preconditioning matrix
+     * that PETSc will use.
+     */
+    void
+    reinit_matrices(AMatrixType &A, PMatrixType &P);
+
+    /**
+     * Return current time.
+     */
+    real_type
+    get_time();
+
+    /**
+     * Return current time step.
+     */
+    real_type
+    get_time_step();
+
+    /**
+     * Integrate the differential-algebraic equations starting from @p y.
+     *
+     * This function returns the final number of computed steps.
+     * Upon returning, the @p y vector contains the solution of the DAE at
+     * the end time.
+     */
+    unsigned int
+    solve(VectorType &y);
+
+    /**
+     * Integrate the differential-algebraic equations starting from @p y.
+     *
+     * This function returns the final number of computed steps.
+     * Upon returning, the @p y vector contains the solution of the DAE at
+     * the end time.
+     *
+     * Here we also pass the matrix to handle the Jacobian.
+     */
+    unsigned int
+    solve(VectorType &y, PMatrixType &P);
+
+    /**
+     * Integrate the differential-algebraic equations starting from @p y.
+     *
+     * This function returns the final number of computed steps.
+     * Upon returning, the @p y vector contains the solution of the DAE at
+     * the end time.
+     *
+     * Here we also pass the matrices to handle Jacobians.
+     */
+    unsigned int
+    solve(VectorType &y, AMatrixType &A, PMatrixType &P);
+
+    /**
+     * Callback for the computation of the implicit residual $F(t, y, \dot y)$.
+     */
+    std::function<int(const real_type   t,
+                      const VectorType &y,
+                      const VectorType &y_dot,
+                      VectorType &      res)>
+      implicit_function;
+
+    /**
+     * Callback for the computation of the implicit Jacobian
+     * $\dfrac{\partial F}{\partial y} + \alpha \dfrac{\partial F}{\partial \dot
+     * y}$
+     *
+     * All implicit solvers implementations are recast to use the above
+     * linearization. The $\alpha$ parameter is time-step and solver-type
+     * specific.
+     */
+    std::function<int(const real_type   t,
+                      const VectorType &y,
+                      const VectorType &y_dot,
+                      const real_type   alpha,
+                      AMatrixType &     A,
+                      PMatrixType &     P)>
+      implicit_jacobian;
+
+    /**
+     * Callback for the computation of the explicit residual $G(t, y)$.
+     */
+    std::function<int(const real_type t, const VectorType &y, VectorType &res)>
+      explicit_function;
+
+    /**
+     * Callback for the computation of the explicit Jacobian $\dfrac{\partial
+     * G}{\partial y}$.
+     */
+    std::function<int(const real_type   t,
+                      const VectorType &y,
+                      AMatrixType &     A,
+                      PMatrixType &     P)>
+      explicit_jacobian;
+
+    /**
+     * Callback for monitoring the solution process.
+     *
+     * This function is called by TimeStepper at the beginning
+     * of each time step.
+     */
+    std::function<int(const real_type    t,
+                      const VectorType & y,
+                      const unsigned int step_number)>
+      monitor;
+
+    /**
+     * Set up Jacobian callback without matrices.
+     *
+     * This callback gives full control to users to set up the linearized
+     * equations
+     * $\dfrac{\partial F}{\partial y} + \alpha \dfrac{\partial F}{\partial \dot
+     * y}$.
+     *
+     * All implicit solvers implementations are recast to use the above
+     * linearization. The $\alpha$ parameter is time-step and solver-type
+     * specific.
+     *
+     * Solvers must be provided via TimeStepper::solve_for_jacobian_system.
+     */
+    std::function<int(const real_type   t,
+                      const VectorType &y,
+                      const VectorType &ydot,
+                      const real_type   alpha)>
+      setup_jacobian;
+
+    /**
+     * Solution of the Jacobian system set up with TimeStepper::setup_jacobian.
+     *
+     * This is used as a preconditioner inside the Krylov process.
+     */
+    std::function<int(const VectorType &src, VectorType &dst)>
+      solve_for_jacobian_system;
+
+    /**
+     * Return an index set containing the algebraic components.
+     *
+     * Implementation of this function is optional. If your equation is also
+     * algebraic (i.e., it contains algebraic constraints, or Lagrange
+     * multipliers), you should implement this function in order to return only
+     * these components of your system.
+     */
+    std::function<IndexSet()> algebraic_components;
+
+  protected:
+    /**
+     * The PETSc object
+     */
+    TS ts;
+
+    SmartPointer<AMatrixType, TimeStepper> A;
+    SmartPointer<PMatrixType, TimeStepper> P;
+    bool                                   need_dae_tolerances;
+  };
+
+#  ifndef DOXYGEN
+  /* Inline functions */
+
+  template <typename VectorType, typename PMatrixType, typename AMatrixType>
+  inline TS
+  TimeStepper<VectorType, PMatrixType, AMatrixType>::petsc_ts()
+  {
+    return ts;
+  }
+
+  template <typename VectorType, typename PMatrixType, typename AMatrixType>
+  inline MPI_Comm
+  TimeStepper<VectorType, PMatrixType, AMatrixType>::get_mpi_communicator()
+    const
+  {
+    return PetscObjectComm(reinterpret_cast<PetscObject>(ts));
+  }
+
+  template <typename VectorType, typename PMatrixType, typename AMatrixType>
+  inline typename TimeStepper<VectorType, PMatrixType, AMatrixType>::real_type
+  TimeStepper<VectorType, PMatrixType, AMatrixType>::get_time()
+  {
+    PetscReal      t;
+    PetscErrorCode ierr = TSGetTime(ts, &t);
+    AssertThrow(ierr == 0, ExcPETScError(ierr));
+    return t;
+  }
+
+  template <typename VectorType, typename PMatrixType, typename AMatrixType>
+  inline typename TimeStepper<VectorType, PMatrixType, AMatrixType>::real_type
+  TimeStepper<VectorType, PMatrixType, AMatrixType>::get_time_step()
+  {
+    PetscReal      dt;
+    PetscErrorCode ierr = TSGetTimeStep(ts, &dt);
+    AssertThrow(ierr == 0, ExcPETScError(ierr));
+    return dt;
+  }
+#  endif
+
+} // namespace PETScWrappers
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_WITH_PETSC
+
+#endif
diff --git a/include/deal.II/lac/petsc_ts.templates.h b/include/deal.II/lac/petsc_ts.templates.h
new file mode 100644 (file)
index 0000000..93c0fb0
--- /dev/null
@@ -0,0 +1,483 @@
+//-----------------------------------------------------------
+//
+//    Copyright (C) 2022 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.md at
+//    the top level directory of deal.II.
+//
+//---------------------------------------------------------------
+//
+// Author: Stefano Zampini, King Abdullah University of Science and Technology.
+
+#ifndef dealii_petsc_ts_templates_h
+#define dealii_petsc_ts_templates_h
+
+#include <deal.II/base/config.h>
+
+#ifdef DEAL_II_WITH_PETSC
+#  include <deal.II/base/exceptions.h>
+
+#  include <deal.II/lac/petsc_precondition.h>
+#  include <deal.II/lac/petsc_ts.h>
+
+#  include <petscdm.h>
+#  include <petscts.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+// Shorthand notation for PETSc error codes.
+#  define AssertTS(code)                             \
+    do                                               \
+      {                                              \
+        PetscErrorCode ierr = (code);                \
+        AssertThrow(ierr == 0, ExcPETScError(ierr)); \
+      }                                              \
+    while (0)
+
+// Shorthand notation for User error codes.
+#  define AssertUser(code, name)                                               \
+    do                                                                         \
+      {                                                                        \
+        int ierr = (code);                                                     \
+        AssertThrow(ierr == 0,                                                 \
+                    StandardExceptions::ExcFunctionNonzeroReturn(name, ierr)); \
+      }                                                                        \
+    while (0)
+
+namespace PETScWrappers
+{
+  template <typename VectorType, typename PMatrixType, typename AMatrixType>
+  TimeStepper<VectorType, PMatrixType, AMatrixType>::TimeStepper(
+    const TimeStepperData &data,
+    const MPI_Comm &       mpi_comm)
+  {
+    AssertTS(TSCreate(mpi_comm, &ts));
+    AssertTS(TSSetApplicationContext(ts, this));
+    reinit(data);
+  }
+
+
+
+  template <typename VectorType, typename PMatrixType, typename AMatrixType>
+  TimeStepper<VectorType, PMatrixType, AMatrixType>::~TimeStepper()
+  {
+    AssertTS(TSDestroy(&ts));
+  }
+
+
+
+  template <typename VectorType, typename PMatrixType, typename AMatrixType>
+  void
+  TimeStepper<VectorType, PMatrixType, AMatrixType>::reinit()
+  {
+    AssertTS(TSReset(ts));
+  }
+
+
+
+  template <typename VectorType, typename PMatrixType, typename AMatrixType>
+  void
+  TimeStepper<VectorType, PMatrixType, AMatrixType>::reinit(
+    const TimeStepperData &data)
+  {
+    reinit();
+
+    // Solver type
+    if (data.tstype.size())
+      AssertTS(TSSetType(ts, data.tstype.c_str()));
+
+    // Options prefix
+    if (data.opts_prefix.size())
+      AssertTS(TSSetOptionsPrefix(ts, data.opts_prefix.c_str()));
+
+    // Time and steps limits
+    AssertTS(TSSetTime(ts, data.initial_time));
+    if (data.final_time > data.initial_time)
+      AssertTS(TSSetMaxTime(ts, data.final_time));
+    if (data.initial_step_size > 0.0)
+      AssertTS(TSSetTimeStep(ts, data.initial_step_size));
+    if (data.max_steps >= 0)
+      AssertTS(TSSetMaxSteps(ts, data.max_steps));
+
+    // Decide how to end the integration. Either stepover the final time or
+    // match it.
+    AssertTS(TSSetExactFinalTime(ts,
+                                 data.match_step ? TS_EXACTFINALTIME_MATCHSTEP :
+                                                   TS_EXACTFINALTIME_STEPOVER));
+
+    // Adaptive tolerances
+    const PetscReal atol = data.absolute_tolerance > 0.0 ?
+                             data.absolute_tolerance :
+                             static_cast<PetscReal>(PETSC_DEFAULT);
+    const PetscReal rtol = data.relative_tolerance > 0.0 ?
+                             data.relative_tolerance :
+                             static_cast<PetscReal>(PETSC_DEFAULT);
+    AssertTS(TSSetTolerances(ts, atol, nullptr, rtol, nullptr));
+
+    // At this point we do not know the problem size so we cannot
+    // set variable tolerances for differential and algebratic equations
+    // Store this value and use it during solve.
+    this->need_dae_tolerances = data.ignore_algebraic_lte;
+
+    // Adaptive time stepping
+    TSAdapt tsadapt;
+    AssertTS(TSGetAdapt(ts, &tsadapt));
+    AssertTS(TSAdaptSetType(tsadapt, data.tsadapttype.c_str()));
+
+    // As of 3.19, PETSc does not propagate options prefixes to the
+    // adaptors.
+    if (data.opts_prefix.size())
+      AssertTS(TSAdaptSetOptionsPrefix(tsadapt, data.opts_prefix.c_str()));
+
+    // Time step limits
+    const PetscReal hmin = data.minimum_step_size > 0.0 ?
+                             data.minimum_step_size :
+                             static_cast<PetscReal>(PETSC_DEFAULT);
+    const PetscReal hmax = data.maximum_step_size > 0.0 ?
+                             data.maximum_step_size :
+                             static_cast<PetscReal>(PETSC_DEFAULT);
+    AssertTS(TSAdaptSetStepLimits(tsadapt, hmin, hmax));
+  }
+
+
+
+  template <typename VectorType, typename PMatrixType, typename AMatrixType>
+  void
+  TimeStepper<VectorType, PMatrixType, AMatrixType>::reinit_matrix(
+    PMatrixType &P)
+  {
+    this->A = nullptr;
+    this->P = &P;
+  }
+
+
+
+  template <typename VectorType, typename PMatrixType, typename AMatrixType>
+  void
+  TimeStepper<VectorType, PMatrixType, AMatrixType>::reinit_matrices(
+    AMatrixType &A,
+    PMatrixType &P)
+  {
+    this->A = &A;
+    this->P = &P;
+  }
+
+
+
+  template <typename VectorType, typename PMatrixType, typename AMatrixType>
+  unsigned int
+  TimeStepper<VectorType, PMatrixType, AMatrixType>::solve(VectorType &y)
+  {
+    auto ts_ifunction_ =
+      [](TS ts, PetscReal t, Vec x, Vec xdot, Vec f, void *ctx)
+      -> PetscErrorCode {
+      PetscFunctionBeginUser;
+      (void)ts;
+      auto user = static_cast<TimeStepper *>(ctx);
+
+      VectorType xdealii(x);
+      VectorType xdotdealii(xdot);
+      VectorType fdealii(f);
+      AssertUser(user->implicit_function(t, xdealii, xdotdealii, fdealii),
+                 "implicit_function");
+      PetscFunctionReturn(0);
+    };
+
+    auto ts_ijacobian_ = [](TS        ts,
+                            PetscReal t,
+                            Vec       x,
+                            Vec       xdot,
+                            PetscReal s,
+                            Mat       A,
+                            Mat       P,
+                            void *    ctx) -> PetscErrorCode {
+      PetscFunctionBeginUser;
+      (void)ts;
+      auto user = static_cast<TimeStepper *>(ctx);
+
+      VectorType  xdealii(x);
+      VectorType  xdotdealii(xdot);
+      AMatrixType Adealii(A);
+      PMatrixType Pdealii(P);
+
+      AssertUser(
+        user->implicit_jacobian(t, xdealii, xdotdealii, s, Adealii, Pdealii),
+        "implicit_jacobian");
+
+      // Handle the Jacobian-free case
+      // This call allow to resample the linearization point
+      // of the MFFD tangent operator
+      PetscBool flg;
+      AssertTS(PetscObjectTypeCompare((PetscObject)A, MATMFFD, &flg));
+      if (flg)
+        {
+          AssertTS(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
+          AssertTS(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
+        }
+      PetscFunctionReturn(0);
+    };
+
+    auto ts_ijacobian_with_setup_ = [](TS        ts,
+                                       PetscReal t,
+                                       Vec       x,
+                                       Vec       xdot,
+                                       PetscReal s,
+                                       Mat       A,
+                                       Mat       P,
+                                       void *    ctx) -> PetscErrorCode {
+      PetscFunctionBeginUser;
+      (void)ts;
+      auto user = static_cast<TimeStepper *>(ctx);
+
+      VectorType  xdealii(x);
+      VectorType  xdotdealii(xdot);
+      AMatrixType Adealii(A);
+      PMatrixType Pdealii(P);
+
+      user->A = &Adealii;
+      user->P = &Pdealii;
+      AssertUser(user->setup_jacobian(t, xdealii, xdotdealii, s),
+                 "setup_jacobian");
+
+      // Handle the Jacobian-free case
+      // This call allow to resample the linearization point
+      // of the MFFD tangent operator
+      PetscBool flg;
+      AssertTS(PetscObjectTypeCompare((PetscObject)A, MATMFFD, &flg));
+      if (flg)
+        {
+          AssertTS(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
+          AssertTS(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
+        }
+      PetscFunctionReturn(0);
+    };
+
+    auto ts_rhsfunction_ =
+      [](TS ts, PetscReal t, Vec x, Vec f, void *ctx) -> PetscErrorCode {
+      PetscFunctionBeginUser;
+      (void)ts;
+      auto user = static_cast<TimeStepper *>(ctx);
+
+      VectorType xdealii(x);
+      VectorType fdealii(f);
+
+      AssertUser(user->explicit_function(t, xdealii, fdealii),
+                 "explicit_function");
+      PetscFunctionReturn(0);
+    };
+
+    auto ts_rhsjacobian_ =
+      [](TS ts, PetscReal t, Vec x, Mat A, Mat P, void *ctx) -> PetscErrorCode {
+      PetscFunctionBeginUser;
+      (void)ts;
+      auto user = static_cast<TimeStepper *>(ctx);
+
+      VectorType  xdealii(x);
+      AMatrixType Adealii(A);
+      PMatrixType Pdealii(P);
+
+      AssertUser(user->explicit_jacobian(t, xdealii, Adealii, Pdealii),
+                 "explicit_jacobian");
+
+      // Handle the Jacobian-free case
+      // This call allow to resample the linearization point
+      // of the MFFD tangent operator
+      PetscBool flg;
+      AssertTS(PetscObjectTypeCompare((PetscObject)A, MATMFFD, &flg));
+      if (flg)
+        {
+          AssertTS(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
+          AssertTS(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
+        }
+      PetscFunctionReturn(0);
+    };
+
+    auto ts_monitor_ =
+      [](TS ts, PetscInt it, PetscReal t, Vec x, void *ctx) -> PetscErrorCode {
+      PetscFunctionBeginUser;
+      (void)ts;
+      auto user = static_cast<TimeStepper *>(ctx);
+
+      VectorType xdealii(x);
+      AssertUser(user->monitor(t, xdealii, it), "monitor");
+      PetscFunctionReturn(0);
+    };
+
+    AssertThrow(explicit_function || implicit_function,
+                StandardExceptions::ExcFunctionNotProvided(
+                  "explicit_function || implicit_function"));
+
+    AssertTS(TSSetSolution(ts, y.petsc_vector()));
+
+    if (explicit_function)
+      AssertTS(TSSetRHSFunction(ts, nullptr, ts_rhsfunction_, this));
+
+    if (implicit_function)
+      AssertTS(TSSetIFunction(ts, nullptr, ts_ifunction_, this));
+
+    if (setup_jacobian)
+      {
+        AssertTS(TSSetIJacobian(ts,
+                                A ? A->petsc_matrix() : nullptr,
+                                P ? P->petsc_matrix() : nullptr,
+                                ts_ijacobian_with_setup_,
+                                this));
+
+        // Tell PETSc to setup a MFFD operator for the linear system matrix
+        SNES snes;
+        AssertTS(TSGetSNES(ts, &snes));
+        if (!A)
+          AssertTS(SNESSetUseMatrixFree(snes, PETSC_TRUE, PETSC_FALSE));
+
+        // Do not waste memory by creating a dummy AIJ matrix inside PETSc.
+        if (!P)
+          {
+            DM dm;
+            AssertTS(SNESGetDM(snes, &dm));
+            AssertTS(DMSetMatType(dm, MATSHELL));
+          }
+      }
+    else
+      {
+        if (explicit_jacobian)
+          {
+            AssertTS(TSSetRHSJacobian(ts,
+                                      A ? A->petsc_matrix() :
+                                          (P ? P->petsc_matrix() : nullptr),
+                                      P ? P->petsc_matrix() : nullptr,
+                                      ts_rhsjacobian_,
+                                      this));
+          }
+
+        if (implicit_jacobian)
+          {
+            AssertTS(TSSetIJacobian(ts,
+                                    A ? A->petsc_matrix() :
+                                        (P ? P->petsc_matrix() : nullptr),
+                                    P ? P->petsc_matrix() : nullptr,
+                                    ts_ijacobian_,
+                                    this));
+          }
+
+        // The user did not set any Jacobian callback. PETSc default in this
+        // case is to use FD and thus assemble a dense operator by finite
+        // differencing the residual callbacks. Here instead we decide to
+        // use a full matrix-free approach by default. This choice can always
+        // be overriden from command line.
+        if (!explicit_jacobian && !implicit_jacobian)
+          {
+            SNES snes;
+            AssertTS(TSGetSNES(ts, &snes));
+            AssertTS(SNESSetUseMatrixFree(snes, PETSC_FALSE, PETSC_TRUE));
+          }
+      }
+
+    // In case solve_for_jacobian_system is provided, create a shell
+    // preconditioner wrapping the user call. The default internal Krylov
+    // solver only applies the preconditioner. This choice
+    // can be overriden by command line and users can use any other
+    // Krylov method if their solve is not accurate enough.
+    // Using solve_for_jacobian_system as a preconditioner allow users
+    // to provide approximate solvers and possibly iterate on a matrix-free
+    // approximation of the tangent operator.
+    PreconditionShell precond(
+      PetscObjectComm(reinterpret_cast<PetscObject>(ts)));
+    if (solve_for_jacobian_system)
+      {
+        precond.vmult = [&](VectorBase &indst, const VectorBase &insrc) -> int {
+          VectorType       dst(static_cast<const Vec &>(indst));
+          const VectorType src(static_cast<const Vec &>(insrc));
+          return solve_for_jacobian_system(src, dst);
+        };
+
+        // Default Krylov solver (preconditioner only)
+        SNES snes;
+        KSP  ksp;
+        AssertTS(TSGetSNES(ts, &snes));
+        AssertTS(SNESGetKSP(snes, &ksp));
+        AssertTS(KSPSetType(ksp, KSPPREONLY));
+        AssertTS(KSPSetPC(ksp, precond.get_pc()));
+      }
+
+    // Attach user monitoring routine.
+    if (monitor)
+      AssertTS(TSMonitorSet(ts, ts_monitor_, this, nullptr));
+
+    // Allow command line customization.
+    AssertTS(TSSetFromOptions(ts));
+
+    // Handle algebraic components.
+    if (algebraic_components && need_dae_tolerances)
+      {
+        PetscReal atol, rtol;
+        AssertTS(TSGetTolerances(ts, &atol, nullptr, &rtol, nullptr));
+
+        Vec av, rv;
+        AssertTS(VecDuplicate(y.petsc_vector(), &av));
+        AssertTS(VecDuplicate(y.petsc_vector(), &rv));
+
+        VectorBase avdealii(av);
+        VectorBase rvdealii(rv);
+        avdealii = atol;
+        rvdealii = rtol;
+        for (auto i : algebraic_components())
+          {
+            avdealii[i] = -1.0;
+            rvdealii[i] = -1.0;
+          }
+        avdealii.compress(VectorOperation::insert);
+        rvdealii.compress(VectorOperation::insert);
+        AssertTS(TSSetTolerances(ts, atol, av, rtol, rv));
+        AssertTS(VecDestroy(&av));
+        AssertTS(VecDestroy(&rv));
+      }
+
+    // Having set everything up, now do the actual work
+    // and let PETSc do the time stepping.
+    AssertTS(TSSolve(ts, nullptr));
+
+    // Return the number of steps taken.
+    PetscInt nt;
+    AssertTS(TSGetStepNumber(ts, &nt));
+    return nt;
+  }
+
+
+
+  template <typename VectorType, typename PMatrixType, typename AMatrixType>
+  unsigned int
+  TimeStepper<VectorType, PMatrixType, AMatrixType>::solve(VectorType & y,
+                                                           PMatrixType &P)
+  {
+    reinit_matrix(P);
+    return solve(y);
+  }
+
+
+
+  template <typename VectorType, typename PMatrixType, typename AMatrixType>
+  unsigned int
+  TimeStepper<VectorType, PMatrixType, AMatrixType>::solve(VectorType & y,
+                                                           AMatrixType &A,
+                                                           PMatrixType &P)
+  {
+    reinit_matrices(A, P);
+    return solve(y);
+  }
+
+} // namespace PETScWrappers
+
+#  undef AssertTS
+#  undef AssertUser
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_WITH_PETSC
+
+#endif
index 1ada9bea314190ccf979fc31d7274e28d4beb329..8b6cdd0ded92a17a336df4d011accf7aeec8ec82 100644 (file)
@@ -98,6 +98,7 @@ set(_inst
 if(DEAL_II_WITH_PETSC)
   set(_unity_include_src
     ${_unity_include_src}
+    petsc_communication_pattern.cc
     petsc_full_matrix.cc
     petsc_matrix_base.cc
     petsc_matrix_free.cc
@@ -107,7 +108,7 @@ if(DEAL_II_WITH_PETSC)
     petsc_parallel_vector.cc
     petsc_precondition.cc
     petsc_solver.cc
-    petsc_communication_pattern.cc
+    petsc_ts.cc
     petsc_sparse_matrix.cc
     petsc_vector_base.cc
   )
diff --git a/source/lac/petsc_ts.cc b/source/lac/petsc_ts.cc
new file mode 100644 (file)
index 0000000..d701808
--- /dev/null
@@ -0,0 +1,41 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2023 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#ifndef DOXYGEN
+
+#  include <deal.II/base/config.h>
+
+#  ifdef DEAL_II_WITH_PETSC
+
+#    include <deal.II/lac/petsc_block_sparse_matrix.h>
+#    include <deal.II/lac/petsc_ts.templates.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+template class PETScWrappers::TimeStepper<>;
+template class PETScWrappers::TimeStepper<PETScWrappers::MPI::Vector>;
+template class PETScWrappers::TimeStepper<PETScWrappers::MPI::BlockVector>;
+template class PETScWrappers::TimeStepper<PETScWrappers::MPI::Vector,
+                                          PETScWrappers::MPI::SparseMatrix>;
+template class PETScWrappers::TimeStepper<
+  PETScWrappers::MPI::BlockVector,
+  PETScWrappers::MPI::BlockSparseMatrix>;
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#  endif // DEAL_II_WITH_PETSC
+#endif
diff --git a/tests/petsc/petsc_ts_00.cc b/tests/petsc/petsc_ts_00.cc
new file mode 100644 (file)
index 0000000..d95ae8a
--- /dev/null
@@ -0,0 +1,264 @@
+//-----------------------------------------------------------
+//
+//    Copyright (C) 2022 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.md at
+//    the top level directory of deal.II.
+//
+//-----------------------------------------------------------
+//
+// Author: Stefano Zampini, King Abdullah University of Science and Technology.
+
+#include <deal.II/base/parameter_handler.h>
+
+#include <deal.II/lac/petsc_matrix_base.h>
+#include <deal.II/lac/petsc_ts.h>
+#include <deal.II/lac/petsc_vector.h>
+
+#include <cmath>
+
+#include "../tests.h"
+
+/**
+ * Solves the harmonic oscillator problem
+ *
+ * u'' = -k^2 u
+ * u (t0) = sin(k * t0)
+ * u'(t0) = k cos (k* t0)
+ *
+ * using the PETScWrappers::TimeStepper class
+ * that interfaces PETSc TS ODE solver object.
+ * The ODE can be an EXPLICIT first order ode:
+ *
+ * y[0]' =       y[1]   \
+ *                       -> y' = G(y,t)
+ * y[1]' = - k^2 y[0]   /
+ *
+ * or specified in IMPLICIT form:
+ *
+ * y[0]' -     y[1]  = 0 \
+ *                        -> F(y',y,t) = 0
+ * y[1]' + k^2 y[0]  = 0 /
+ *
+ * The exact solution is
+ *
+ * y[0](t) = sin(k t)
+ * y[1](t) = k cos(k t)
+ *
+ * We use the same class to test both formulations.
+ * The class also supports a third approach where
+ * control for linear systems generation and
+ * solution is completely on users. In this case
+ * users are in charge of solving for the
+ * implicit Jacobian.
+ * The model used to wrap PETSc's TS is the same
+ * used for the nonlinear solver SNES. Check
+ * petsc_snes_00.cc for additional information.
+ *
+ */
+using VectorType  = PETScWrappers::MPI::Vector;
+using MatrixType  = PETScWrappers::MatrixBase;
+using TimeStepper = PETScWrappers::TimeStepper<VectorType, MatrixType>;
+using real_type   = TimeStepper::real_type;
+
+class HarmonicOscillator
+{
+public:
+  HarmonicOscillator(real_type                                      _kappa,
+                     const typename PETScWrappers::TimeStepperData &data,
+                     bool                                           setjac,
+                     bool                                           implicit,
+                     bool                                           user,
+                     std::ostream &                                 _out)
+    : time_stepper(data)
+    , out(_out)
+    , kappa(_kappa)
+  {
+    // In this case we use the implicit form
+    if (implicit)
+      {
+        time_stepper.implicit_function = [&](const real_type   t,
+                                             const VectorType &y,
+                                             const VectorType &y_dot,
+                                             VectorType &      res) -> int {
+          res(0) = y_dot(0) - y(1);
+          res(1) = y_dot(1) + kappa * kappa * y(0);
+          res.compress(VectorOperation::insert);
+          return 0;
+        };
+
+        // We either have the possibility of using PETSc standard
+        // functionalities for Jacobian setup and solve, or take full control of
+        // the Jacobian solutions. The latter aligns more with deal.II tutorials
+        // For our testing class, 'user' discriminate between those two
+        // approaches.
+        if (!user && setjac)
+          {
+            // Callback for Jacobian evaluation
+            // This callback populates the P matrix with
+            //     shift * dF/dydot + dF/dy
+            // A and P are the same matrix if not otherwise specified
+            // PETSc can compute a Jacobian by finite-differencing the residual
+            // evaluation, and it can be selected at command line.
+            time_stepper.implicit_jacobian = [&](const real_type   t,
+                                                 const VectorType &y,
+                                                 const VectorType &y_dot,
+                                                 const real_type   shift,
+                                                 MatrixType &      A,
+                                                 MatrixType &      P) -> int {
+              P.set(0, 0, shift);
+              P.set(0, 1, -1);
+              P.set(1, 0, kappa * kappa);
+              P.set(1, 1, shift);
+              P.compress(VectorOperation::insert);
+              return 0;
+            };
+          }
+        else if (user)
+          {
+            // In this code block, users are completely in charge of
+            // setting up the Jacobian system and solve for it.
+            // In this example we only store the solver shift
+            // during setup.
+            time_stepper.setup_jacobian = [&](const real_type   t,
+                                              const VectorType &y,
+                                              const VectorType &y_dot,
+                                              const real_type   shift) -> int {
+              myshift = shift;
+              return 0;
+            };
+
+            // In the solve phase we se the stored shift to solve
+            // for the implicit Jacobian system
+            time_stepper.solve_for_jacobian_system =
+              [&](const VectorType &src, VectorType &dst) -> int {
+              auto sf = 1. / (kappa * kappa + myshift * myshift);
+              dst(0)  = sf * (myshift * src(0) + src(1));
+              dst(1)  = sf * (-kappa * kappa * src(0) + myshift * src(1));
+              dst.compress(VectorOperation::insert);
+              return 0;
+            };
+          }
+      }
+    else
+      { // Here we instead use the explicit form
+        // This is the only function one would populate in case an explicit
+        // solver is used.
+        time_stepper.explicit_function =
+          [&](const real_type t, const VectorType &y, VectorType &res) -> int {
+          res(0) = y(1);
+          res(1) = -kappa * kappa * y(0);
+          res.compress(VectorOperation::insert);
+          return 0;
+        };
+
+        // The explicit Jacobian callback is not needed in case
+        // an explicit solver is used. We need it when moving to
+        // an implicit solver like in the parameter file used for
+        // this test
+        if (setjac)
+          {
+            time_stepper.explicit_jacobian = [&](const real_type   t,
+                                                 const VectorType &y,
+                                                 MatrixType &      A,
+                                                 MatrixType &      P) -> int {
+              P.set(0, 0, 0);
+              P.set(0, 1, 1);
+              P.set(1, 0, -kappa * kappa);
+              P.set(1, 1, 0);
+              P.compress(VectorOperation::insert);
+              return 0;
+            };
+          }
+      }
+
+    // Monitoring routine. Here we print diagnostic for the exact
+    // solution to the log file.
+    time_stepper.monitor = [&](const real_type    t,
+                               const VectorType & y,
+                               const unsigned int step_number) -> int {
+      std::vector<real_type> exact(2);
+      exact[0] = std::sin(kappa * t);
+      exact[1] = kappa * std::cos(kappa * t);
+      out << t << " " << y(0) << " (" << exact[0] << ")"
+          << " " << y(1) << " (" << exact[1] << ")" << std::endl;
+      return 0;
+    };
+  }
+
+  void
+  run()
+  {
+    // Set initial conditions
+    // This test uses a nonzero initial time
+    // We can access the time value with the
+    // get_time method
+    auto t = time_stepper.get_time();
+
+    VectorType y(MPI_COMM_SELF, 2, 2);
+    y[0] = std::sin(kappa * t);
+    y[1] = kappa * std::cos(kappa * t);
+    y.compress(VectorOperation::insert);
+
+    // Integrate the ODE.
+    time_stepper.solve(y);
+  }
+
+private:
+  TimeStepper   time_stepper;
+  std::ostream &out;     // Used by the monitoring routine
+  real_type     kappa;   // Defines the oscillator
+  real_type     myshift; // Used by the user solve
+};
+
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+  std::ofstream out("output");
+
+  PETScWrappers::TimeStepperData data;
+  ParameterHandler               prm;
+
+  data.add_parameters(prm);
+  out << "# Default Parameters" << std::endl;
+  prm.print_parameters(out, ParameterHandler::ShortText);
+
+  std::ifstream ifile(SOURCE_DIR "/petsc_ts_00_in.prm");
+  prm.parse_input(ifile);
+
+  out << "# Testing Parameters" << std::endl;
+  prm.print_parameters(out, ParameterHandler::ShortText);
+
+  for (int setjaci = 0; setjaci < 2; setjaci++)
+    {
+      bool setjac = setjaci ? true : false;
+
+      {
+        out << "# Test explicit interface (J " << setjac << ")" << std::endl;
+        HarmonicOscillator ode_expl(1.0, data, setjac, false, false, out);
+        ode_expl.run();
+      }
+
+      {
+        out << "# Test implicit interface (J " << setjac << ")" << std::endl;
+        HarmonicOscillator ode_impl(1.0, data, setjac, true, false, out);
+        ode_impl.run();
+      }
+    }
+
+  {
+    out << "# Test user interface" << std::endl;
+    HarmonicOscillator ode_user(1.0, data, true, true, true, out);
+    ode_user.run();
+  }
+  return 0;
+}
diff --git a/tests/petsc/petsc_ts_00.output b/tests/petsc/petsc_ts_00.output
new file mode 100644 (file)
index 0000000..3eba927
--- /dev/null
@@ -0,0 +1,446 @@
+# Default Parameters
+subsection Error control
+  set absolute error tolerance = -1
+  set adaptor type             = none
+  set ignore algebraic lte     = true
+  set maximum step size        = -1
+  set minimum step size        = -1
+  set relative error tolerance = -1
+end
+subsection Running parameters
+  set final time              = 0
+  set initial step size       = 0
+  set initial time            = 0
+  set match final time        = false
+  set maximum number of steps = -1
+  set options prefix          = 
+  set solver type             = 
+end
+# Testing Parameters
+subsection Error control
+  set absolute error tolerance = -1
+  set adaptor type             = none
+  set ignore algebraic lte     = true
+  set maximum step size        = 1e+7
+  set minimum step size        = 1e-7
+  set relative error tolerance = -1
+end
+subsection Running parameters
+  set final time              = 6.0
+  set initial step size       = 0.0625
+  set initial time            = 1.0
+  set match final time        = false
+  set maximum number of steps = -1
+  set options prefix          = 
+  set solver type             = bdf
+end
+# Test explicit interface (J 0)
+1 0.841471 (0.841471) 0.540302 (0.540302)
+1.0625 0.873 (0.873575) 0.486385 (0.48669)
+1.125 0.901462 (0.902268) 0.430852 (0.431177)
+1.1875 0.926521 (0.927437) 0.373735 (0.37398)
+1.25 0.948008 (0.948985) 0.315196 (0.315322)
+1.3125 0.965813 (0.966827) 0.255441 (0.255434)
+1.375 0.979861 (0.980893) 0.194695 (0.194548)
+1.4375 0.990093 (0.991129) 0.133193 (0.132902)
+1.5 0.996469 (0.997495) 0.0711722 (0.0707372)
+1.5625 0.998964 (0.999966) 0.00887518 (0.00829623)
+1.625 0.997569 (0.998531) -0.0534558 (-0.0541771)
+1.6875 0.992288 (0.993198) -0.115578 (-0.116439)
+1.75 0.983142 (0.983986) -0.17725 (-0.178246)
+1.8125 0.970167 (0.970932) -0.238231 (-0.239357)
+1.875 0.953414 (0.954086) -0.298283 (-0.299534)
+1.9375 0.932947 (0.933514) -0.357174 (-0.35854)
+2 0.908847 (0.909297) -0.414672 (-0.416147)
+2.0625 0.881207 (0.88153) -0.470556 (-0.472128)
+2.125 0.850135 (0.85032) -0.524606 (-0.526266)
+2.1875 0.815752 (0.815789) -0.576612 (-0.578349)
+2.25 0.778193 (0.778073) -0.626372 (-0.628174)
+2.3125 0.737602 (0.737319) -0.673692 (-0.675545)
+2.375 0.694139 (0.693685) -0.718388 (-0.720278)
+2.4375 0.647973 (0.647343) -0.760286 (-0.762199)
+2.5 0.599283 (0.598472) -0.799221 (-0.801144)
+2.5625 0.54826 (0.547265) -0.835044 (-0.83696)
+2.625 0.495102 (0.49392) -0.867614 (-0.869507)
+2.6875 0.440015 (0.438647) -0.896805 (-0.898659)
+2.75 0.383215 (0.381661) -0.922503 (-0.924302)
+2.8125 0.324924 (0.323185) -0.944607 (-0.946336)
+2.875 0.265367 (0.263446) -0.963032 (-0.964674)
+2.9375 0.204777 (0.202679) -0.977706 (-0.979245)
+3 0.143389 (0.14112) -0.988572 (-0.989992)
+3.0625 0.0814441 (0.0790102) -0.995588 (-0.996874)
+3.125 0.0191822 (0.0165919) -0.998725 (-0.999862)
+3.1875 -0.043154 (-0.0458912) -0.997973 (-0.998946)
+3.25 -0.105322 (-0.108195) -0.993334 (-0.99413)
+3.3125 -0.167079 (-0.170077) -0.984826 (-0.985431)
+3.375 -0.228184 (-0.231294) -0.972482 (-0.972884)
+3.4375 -0.288401 (-0.291608) -0.956351 (-0.956538)
+3.5 -0.347494 (-0.350783) -0.936495 (-0.936457)
+3.5625 -0.405233 (-0.408589) -0.912992 (-0.912719)
+3.625 -0.461393 (-0.464799) -0.885933 (-0.885416)
+3.6875 -0.515755 (-0.519194) -0.855423 (-0.854657)
+3.75 -0.568109 (-0.571561) -0.821582 (-0.820559)
+3.8125 -0.618249 (-0.621697) -0.784541 (-0.783258)
+3.875 -0.665981 (-0.669405) -0.744445 (-0.742898)
+3.9375 -0.711119 (-0.714499) -0.70145 (-0.699637)
+4 -0.753487 (-0.756802) -0.655723 (-0.653644)
+4.0625 -0.792919 (-0.796151) -0.607442 (-0.605098)
+4.125 -0.829263 (-0.832391) -0.556796 (-0.55419)
+4.1875 -0.862377 (-0.86538) -0.503981 (-0.501117)
+4.25 -0.892132 (-0.894989) -0.449204 (-0.446087)
+4.3125 -0.918412 (-0.921104) -0.392678 (-0.389316)
+4.375 -0.941115 (-0.943622) -0.334623 (-0.331024)
+4.4375 -0.960152 (-0.962455) -0.275265 (-0.27144)
+4.5 -0.975449 (-0.97753) -0.214835 (-0.210796)
+4.5625 -0.986946 (-0.988788) -0.153569 (-0.149328)
+4.625 -0.9946 (-0.996184) -0.0917056 (-0.0872778)
+4.6875 -0.99838 (-0.99969) -0.0294852 (-0.0248864)
+4.75 -0.998271 (-0.999293) 0.0328495 (0.0376022)
+4.8125 -0.994274 (-0.994993) 0.0950559 (0.0999439)
+4.875 -0.986405 (-0.986808) 0.156891 (0.161895)
+4.9375 -0.974693 (-0.974769) 0.218116 (0.223215)
+5 -0.959186 (-0.958924) 0.27849 (0.283662)
+5.0625 -0.939943 (-0.939335) 0.337779 (0.343002)
+5.125 -0.917038 (-0.916077) 0.395752 (0.401003)
+5.1875 -0.890563 (-0.889242) 0.452183 (0.457437)
+5.25 -0.860619 (-0.858934) 0.506852 (0.512085)
+5.3125 -0.827323 (-0.825273) 0.559547 (0.564734)
+5.375 -0.790805 (-0.788389) 0.610062 (0.615177)
+5.4375 -0.751207 (-0.748426) 0.658201 (0.663219)
+5.5 -0.708684 (-0.70554) 0.703776 (0.70867)
+5.5625 -0.663401 (-0.6599) 0.746609 (0.751354)
+5.625 -0.615534 (-0.611682) 0.786534 (0.791104)
+5.6875 -0.565271 (-0.561076) 0.823395 (0.827764)
+5.75 -0.512806 (-0.508279) 0.857049 (0.861192)
+5.8125 -0.458344 (-0.453497) 0.887365 (0.891258)
+5.875 -0.402097 (-0.396944) 0.914224 (0.917843)
+5.9375 -0.344285 (-0.338842) 0.937523 (0.940843)
+6 -0.285132 (-0.279415) 0.957169 (0.96017)
+# Test implicit interface (J 0)
+1 0.841471 (0.841471) 0.540302 (0.540302)
+1.0625 0.873 (0.873575) 0.486385 (0.48669)
+1.125 0.901462 (0.902268) 0.430852 (0.431177)
+1.1875 0.926521 (0.927437) 0.373735 (0.37398)
+1.25 0.948008 (0.948985) 0.315196 (0.315322)
+1.3125 0.965813 (0.966827) 0.255441 (0.255434)
+1.375 0.979861 (0.980893) 0.194695 (0.194548)
+1.4375 0.990093 (0.991129) 0.133193 (0.132902)
+1.5 0.996469 (0.997495) 0.0711722 (0.0707372)
+1.5625 0.998964 (0.999966) 0.00887518 (0.00829623)
+1.625 0.997569 (0.998531) -0.0534558 (-0.0541771)
+1.6875 0.992288 (0.993198) -0.115578 (-0.116439)
+1.75 0.983142 (0.983986) -0.17725 (-0.178246)
+1.8125 0.970167 (0.970932) -0.238231 (-0.239357)
+1.875 0.953414 (0.954086) -0.298283 (-0.299534)
+1.9375 0.932947 (0.933514) -0.357174 (-0.35854)
+2 0.908847 (0.909297) -0.414672 (-0.416147)
+2.0625 0.881207 (0.88153) -0.470556 (-0.472128)
+2.125 0.850135 (0.85032) -0.524606 (-0.526266)
+2.1875 0.815752 (0.815789) -0.576612 (-0.578349)
+2.25 0.778193 (0.778073) -0.626372 (-0.628174)
+2.3125 0.737602 (0.737319) -0.673692 (-0.675545)
+2.375 0.694139 (0.693685) -0.718388 (-0.720278)
+2.4375 0.647973 (0.647343) -0.760286 (-0.762199)
+2.5 0.599283 (0.598472) -0.799221 (-0.801144)
+2.5625 0.54826 (0.547265) -0.835044 (-0.83696)
+2.625 0.495102 (0.49392) -0.867614 (-0.869507)
+2.6875 0.440015 (0.438647) -0.896805 (-0.898659)
+2.75 0.383215 (0.381661) -0.922503 (-0.924302)
+2.8125 0.324924 (0.323185) -0.944607 (-0.946336)
+2.875 0.265367 (0.263446) -0.963032 (-0.964674)
+2.9375 0.204777 (0.202679) -0.977706 (-0.979245)
+3 0.143389 (0.14112) -0.988572 (-0.989992)
+3.0625 0.0814441 (0.0790102) -0.995588 (-0.996874)
+3.125 0.0191822 (0.0165919) -0.998725 (-0.999862)
+3.1875 -0.043154 (-0.0458912) -0.997973 (-0.998946)
+3.25 -0.105322 (-0.108195) -0.993334 (-0.99413)
+3.3125 -0.167079 (-0.170077) -0.984826 (-0.985431)
+3.375 -0.228184 (-0.231294) -0.972482 (-0.972884)
+3.4375 -0.288401 (-0.291608) -0.956351 (-0.956538)
+3.5 -0.347494 (-0.350783) -0.936495 (-0.936457)
+3.5625 -0.405233 (-0.408589) -0.912992 (-0.912719)
+3.625 -0.461393 (-0.464799) -0.885933 (-0.885416)
+3.6875 -0.515755 (-0.519194) -0.855423 (-0.854657)
+3.75 -0.568109 (-0.571561) -0.821582 (-0.820559)
+3.8125 -0.618249 (-0.621697) -0.784541 (-0.783258)
+3.875 -0.665981 (-0.669405) -0.744445 (-0.742898)
+3.9375 -0.711119 (-0.714499) -0.70145 (-0.699637)
+4 -0.753487 (-0.756802) -0.655723 (-0.653644)
+4.0625 -0.792919 (-0.796151) -0.607442 (-0.605098)
+4.125 -0.829263 (-0.832391) -0.556796 (-0.55419)
+4.1875 -0.862377 (-0.86538) -0.503981 (-0.501117)
+4.25 -0.892132 (-0.894989) -0.449204 (-0.446087)
+4.3125 -0.918412 (-0.921104) -0.392678 (-0.389316)
+4.375 -0.941115 (-0.943622) -0.334623 (-0.331024)
+4.4375 -0.960152 (-0.962455) -0.275265 (-0.27144)
+4.5 -0.975449 (-0.97753) -0.214835 (-0.210796)
+4.5625 -0.986946 (-0.988788) -0.153569 (-0.149328)
+4.625 -0.9946 (-0.996184) -0.0917056 (-0.0872778)
+4.6875 -0.99838 (-0.99969) -0.0294852 (-0.0248864)
+4.75 -0.998271 (-0.999293) 0.0328495 (0.0376022)
+4.8125 -0.994274 (-0.994993) 0.0950559 (0.0999439)
+4.875 -0.986405 (-0.986808) 0.156891 (0.161895)
+4.9375 -0.974693 (-0.974769) 0.218116 (0.223215)
+5 -0.959186 (-0.958924) 0.27849 (0.283662)
+5.0625 -0.939943 (-0.939335) 0.337779 (0.343002)
+5.125 -0.917038 (-0.916077) 0.395752 (0.401003)
+5.1875 -0.890563 (-0.889242) 0.452183 (0.457437)
+5.25 -0.860619 (-0.858934) 0.506852 (0.512085)
+5.3125 -0.827323 (-0.825273) 0.559547 (0.564734)
+5.375 -0.790805 (-0.788389) 0.610062 (0.615177)
+5.4375 -0.751207 (-0.748426) 0.658201 (0.663219)
+5.5 -0.708684 (-0.70554) 0.703776 (0.70867)
+5.5625 -0.663401 (-0.6599) 0.746609 (0.751354)
+5.625 -0.615534 (-0.611682) 0.786534 (0.791104)
+5.6875 -0.565271 (-0.561076) 0.823395 (0.827764)
+5.75 -0.512806 (-0.508279) 0.857049 (0.861192)
+5.8125 -0.458344 (-0.453497) 0.887365 (0.891258)
+5.875 -0.402097 (-0.396944) 0.914224 (0.917843)
+5.9375 -0.344285 (-0.338842) 0.937523 (0.940843)
+6 -0.285132 (-0.279415) 0.957169 (0.96017)
+# Test explicit interface (J 1)
+1 0.841471 (0.841471) 0.540302 (0.540302)
+1.0625 0.873 (0.873575) 0.486385 (0.48669)
+1.125 0.901462 (0.902268) 0.430852 (0.431177)
+1.1875 0.926521 (0.927437) 0.373735 (0.37398)
+1.25 0.948008 (0.948985) 0.315196 (0.315322)
+1.3125 0.965813 (0.966827) 0.255441 (0.255434)
+1.375 0.979861 (0.980893) 0.194695 (0.194548)
+1.4375 0.990093 (0.991129) 0.133193 (0.132902)
+1.5 0.996469 (0.997495) 0.0711722 (0.0707372)
+1.5625 0.998964 (0.999966) 0.00887518 (0.00829623)
+1.625 0.997569 (0.998531) -0.0534558 (-0.0541771)
+1.6875 0.992288 (0.993198) -0.115578 (-0.116439)
+1.75 0.983142 (0.983986) -0.17725 (-0.178246)
+1.8125 0.970167 (0.970932) -0.238231 (-0.239357)
+1.875 0.953414 (0.954086) -0.298283 (-0.299534)
+1.9375 0.932947 (0.933514) -0.357174 (-0.35854)
+2 0.908847 (0.909297) -0.414672 (-0.416147)
+2.0625 0.881207 (0.88153) -0.470556 (-0.472128)
+2.125 0.850135 (0.85032) -0.524606 (-0.526266)
+2.1875 0.815752 (0.815789) -0.576612 (-0.578349)
+2.25 0.778193 (0.778073) -0.626372 (-0.628174)
+2.3125 0.737602 (0.737319) -0.673692 (-0.675545)
+2.375 0.694139 (0.693685) -0.718388 (-0.720278)
+2.4375 0.647973 (0.647343) -0.760286 (-0.762199)
+2.5 0.599283 (0.598472) -0.799221 (-0.801144)
+2.5625 0.54826 (0.547265) -0.835044 (-0.83696)
+2.625 0.495102 (0.49392) -0.867614 (-0.869507)
+2.6875 0.440015 (0.438647) -0.896805 (-0.898659)
+2.75 0.383215 (0.381661) -0.922503 (-0.924302)
+2.8125 0.324924 (0.323185) -0.944607 (-0.946336)
+2.875 0.265367 (0.263446) -0.963032 (-0.964674)
+2.9375 0.204777 (0.202679) -0.977706 (-0.979245)
+3 0.143389 (0.14112) -0.988572 (-0.989992)
+3.0625 0.0814441 (0.0790102) -0.995588 (-0.996874)
+3.125 0.0191822 (0.0165919) -0.998725 (-0.999862)
+3.1875 -0.043154 (-0.0458912) -0.997973 (-0.998946)
+3.25 -0.105322 (-0.108195) -0.993334 (-0.99413)
+3.3125 -0.167079 (-0.170077) -0.984826 (-0.985431)
+3.375 -0.228184 (-0.231294) -0.972482 (-0.972884)
+3.4375 -0.288401 (-0.291608) -0.956351 (-0.956538)
+3.5 -0.347494 (-0.350783) -0.936495 (-0.936457)
+3.5625 -0.405233 (-0.408589) -0.912992 (-0.912719)
+3.625 -0.461393 (-0.464799) -0.885933 (-0.885416)
+3.6875 -0.515755 (-0.519194) -0.855423 (-0.854657)
+3.75 -0.568109 (-0.571561) -0.821582 (-0.820559)
+3.8125 -0.618249 (-0.621697) -0.784541 (-0.783258)
+3.875 -0.665981 (-0.669405) -0.744445 (-0.742898)
+3.9375 -0.711119 (-0.714499) -0.70145 (-0.699637)
+4 -0.753487 (-0.756802) -0.655723 (-0.653644)
+4.0625 -0.792919 (-0.796151) -0.607442 (-0.605098)
+4.125 -0.829263 (-0.832391) -0.556796 (-0.55419)
+4.1875 -0.862377 (-0.86538) -0.503981 (-0.501117)
+4.25 -0.892132 (-0.894989) -0.449204 (-0.446087)
+4.3125 -0.918412 (-0.921104) -0.392678 (-0.389316)
+4.375 -0.941115 (-0.943622) -0.334623 (-0.331024)
+4.4375 -0.960152 (-0.962455) -0.275265 (-0.27144)
+4.5 -0.975449 (-0.97753) -0.214835 (-0.210796)
+4.5625 -0.986946 (-0.988788) -0.153569 (-0.149328)
+4.625 -0.9946 (-0.996184) -0.0917056 (-0.0872778)
+4.6875 -0.99838 (-0.99969) -0.0294852 (-0.0248864)
+4.75 -0.998271 (-0.999293) 0.0328495 (0.0376022)
+4.8125 -0.994274 (-0.994993) 0.0950559 (0.0999439)
+4.875 -0.986405 (-0.986808) 0.156891 (0.161895)
+4.9375 -0.974693 (-0.974769) 0.218116 (0.223215)
+5 -0.959186 (-0.958924) 0.27849 (0.283662)
+5.0625 -0.939943 (-0.939335) 0.337779 (0.343002)
+5.125 -0.917038 (-0.916077) 0.395752 (0.401003)
+5.1875 -0.890563 (-0.889242) 0.452183 (0.457437)
+5.25 -0.860619 (-0.858934) 0.506852 (0.512085)
+5.3125 -0.827323 (-0.825273) 0.559547 (0.564734)
+5.375 -0.790805 (-0.788389) 0.610062 (0.615177)
+5.4375 -0.751207 (-0.748426) 0.658201 (0.663219)
+5.5 -0.708684 (-0.70554) 0.703776 (0.70867)
+5.5625 -0.663401 (-0.6599) 0.746609 (0.751354)
+5.625 -0.615534 (-0.611682) 0.786534 (0.791104)
+5.6875 -0.565271 (-0.561076) 0.823395 (0.827764)
+5.75 -0.512806 (-0.508279) 0.857049 (0.861192)
+5.8125 -0.458344 (-0.453497) 0.887365 (0.891258)
+5.875 -0.402097 (-0.396944) 0.914224 (0.917843)
+5.9375 -0.344285 (-0.338842) 0.937523 (0.940843)
+6 -0.285132 (-0.279415) 0.957169 (0.96017)
+# Test implicit interface (J 1)
+1 0.841471 (0.841471) 0.540302 (0.540302)
+1.0625 0.873 (0.873575) 0.486385 (0.48669)
+1.125 0.901462 (0.902268) 0.430852 (0.431177)
+1.1875 0.926521 (0.927437) 0.373735 (0.37398)
+1.25 0.948008 (0.948985) 0.315196 (0.315322)
+1.3125 0.965813 (0.966827) 0.255441 (0.255434)
+1.375 0.979861 (0.980893) 0.194695 (0.194548)
+1.4375 0.990093 (0.991129) 0.133193 (0.132902)
+1.5 0.996469 (0.997495) 0.0711722 (0.0707372)
+1.5625 0.998964 (0.999966) 0.00887518 (0.00829623)
+1.625 0.997569 (0.998531) -0.0534558 (-0.0541771)
+1.6875 0.992288 (0.993198) -0.115578 (-0.116439)
+1.75 0.983142 (0.983986) -0.17725 (-0.178246)
+1.8125 0.970167 (0.970932) -0.238231 (-0.239357)
+1.875 0.953414 (0.954086) -0.298283 (-0.299534)
+1.9375 0.932947 (0.933514) -0.357174 (-0.35854)
+2 0.908847 (0.909297) -0.414672 (-0.416147)
+2.0625 0.881207 (0.88153) -0.470556 (-0.472128)
+2.125 0.850135 (0.85032) -0.524606 (-0.526266)
+2.1875 0.815752 (0.815789) -0.576612 (-0.578349)
+2.25 0.778193 (0.778073) -0.626372 (-0.628174)
+2.3125 0.737602 (0.737319) -0.673692 (-0.675545)
+2.375 0.694139 (0.693685) -0.718388 (-0.720278)
+2.4375 0.647973 (0.647343) -0.760286 (-0.762199)
+2.5 0.599283 (0.598472) -0.799221 (-0.801144)
+2.5625 0.54826 (0.547265) -0.835044 (-0.83696)
+2.625 0.495102 (0.49392) -0.867614 (-0.869507)
+2.6875 0.440015 (0.438647) -0.896805 (-0.898659)
+2.75 0.383215 (0.381661) -0.922503 (-0.924302)
+2.8125 0.324924 (0.323185) -0.944607 (-0.946336)
+2.875 0.265367 (0.263446) -0.963032 (-0.964674)
+2.9375 0.204777 (0.202679) -0.977706 (-0.979245)
+3 0.143389 (0.14112) -0.988572 (-0.989992)
+3.0625 0.0814441 (0.0790102) -0.995588 (-0.996874)
+3.125 0.0191822 (0.0165919) -0.998725 (-0.999862)
+3.1875 -0.043154 (-0.0458912) -0.997973 (-0.998946)
+3.25 -0.105322 (-0.108195) -0.993334 (-0.99413)
+3.3125 -0.167079 (-0.170077) -0.984826 (-0.985431)
+3.375 -0.228184 (-0.231294) -0.972482 (-0.972884)
+3.4375 -0.288401 (-0.291608) -0.956351 (-0.956538)
+3.5 -0.347494 (-0.350783) -0.936495 (-0.936457)
+3.5625 -0.405233 (-0.408589) -0.912992 (-0.912719)
+3.625 -0.461393 (-0.464799) -0.885933 (-0.885416)
+3.6875 -0.515755 (-0.519194) -0.855423 (-0.854657)
+3.75 -0.568109 (-0.571561) -0.821582 (-0.820559)
+3.8125 -0.618249 (-0.621697) -0.784541 (-0.783258)
+3.875 -0.665981 (-0.669405) -0.744445 (-0.742898)
+3.9375 -0.711119 (-0.714499) -0.70145 (-0.699637)
+4 -0.753487 (-0.756802) -0.655723 (-0.653644)
+4.0625 -0.792919 (-0.796151) -0.607442 (-0.605098)
+4.125 -0.829263 (-0.832391) -0.556796 (-0.55419)
+4.1875 -0.862377 (-0.86538) -0.503981 (-0.501117)
+4.25 -0.892132 (-0.894989) -0.449204 (-0.446087)
+4.3125 -0.918412 (-0.921104) -0.392678 (-0.389316)
+4.375 -0.941115 (-0.943622) -0.334623 (-0.331024)
+4.4375 -0.960152 (-0.962455) -0.275265 (-0.27144)
+4.5 -0.975449 (-0.97753) -0.214835 (-0.210796)
+4.5625 -0.986946 (-0.988788) -0.153569 (-0.149328)
+4.625 -0.9946 (-0.996184) -0.0917056 (-0.0872778)
+4.6875 -0.99838 (-0.99969) -0.0294852 (-0.0248864)
+4.75 -0.998271 (-0.999293) 0.0328495 (0.0376022)
+4.8125 -0.994274 (-0.994993) 0.0950559 (0.0999439)
+4.875 -0.986405 (-0.986808) 0.156891 (0.161895)
+4.9375 -0.974693 (-0.974769) 0.218116 (0.223215)
+5 -0.959186 (-0.958924) 0.27849 (0.283662)
+5.0625 -0.939943 (-0.939335) 0.337779 (0.343002)
+5.125 -0.917038 (-0.916077) 0.395752 (0.401003)
+5.1875 -0.890563 (-0.889242) 0.452183 (0.457437)
+5.25 -0.860619 (-0.858934) 0.506852 (0.512085)
+5.3125 -0.827323 (-0.825273) 0.559547 (0.564734)
+5.375 -0.790805 (-0.788389) 0.610062 (0.615177)
+5.4375 -0.751207 (-0.748426) 0.658201 (0.663219)
+5.5 -0.708684 (-0.70554) 0.703776 (0.70867)
+5.5625 -0.663401 (-0.6599) 0.746609 (0.751354)
+5.625 -0.615534 (-0.611682) 0.786534 (0.791104)
+5.6875 -0.565271 (-0.561076) 0.823395 (0.827764)
+5.75 -0.512806 (-0.508279) 0.857049 (0.861192)
+5.8125 -0.458344 (-0.453497) 0.887365 (0.891258)
+5.875 -0.402097 (-0.396944) 0.914224 (0.917843)
+5.9375 -0.344285 (-0.338842) 0.937523 (0.940843)
+6 -0.285132 (-0.279415) 0.957169 (0.96017)
+# Test user interface
+1 0.841471 (0.841471) 0.540302 (0.540302)
+1.0625 0.873 (0.873575) 0.486385 (0.48669)
+1.125 0.901462 (0.902268) 0.430852 (0.431177)
+1.1875 0.926521 (0.927437) 0.373735 (0.37398)
+1.25 0.948008 (0.948985) 0.315196 (0.315322)
+1.3125 0.965813 (0.966827) 0.255441 (0.255434)
+1.375 0.979861 (0.980893) 0.194695 (0.194548)
+1.4375 0.990093 (0.991129) 0.133193 (0.132902)
+1.5 0.996469 (0.997495) 0.0711722 (0.0707372)
+1.5625 0.998964 (0.999966) 0.00887518 (0.00829623)
+1.625 0.997569 (0.998531) -0.0534558 (-0.0541771)
+1.6875 0.992288 (0.993198) -0.115578 (-0.116439)
+1.75 0.983142 (0.983986) -0.17725 (-0.178246)
+1.8125 0.970167 (0.970932) -0.238231 (-0.239357)
+1.875 0.953414 (0.954086) -0.298283 (-0.299534)
+1.9375 0.932947 (0.933514) -0.357174 (-0.35854)
+2 0.908847 (0.909297) -0.414672 (-0.416147)
+2.0625 0.881207 (0.88153) -0.470556 (-0.472128)
+2.125 0.850135 (0.85032) -0.524606 (-0.526266)
+2.1875 0.815752 (0.815789) -0.576612 (-0.578349)
+2.25 0.778193 (0.778073) -0.626372 (-0.628174)
+2.3125 0.737602 (0.737319) -0.673692 (-0.675545)
+2.375 0.694139 (0.693685) -0.718388 (-0.720278)
+2.4375 0.647973 (0.647343) -0.760286 (-0.762199)
+2.5 0.599283 (0.598472) -0.799221 (-0.801144)
+2.5625 0.54826 (0.547265) -0.835044 (-0.83696)
+2.625 0.495102 (0.49392) -0.867614 (-0.869507)
+2.6875 0.440015 (0.438647) -0.896805 (-0.898659)
+2.75 0.383215 (0.381661) -0.922503 (-0.924302)
+2.8125 0.324924 (0.323185) -0.944607 (-0.946336)
+2.875 0.265367 (0.263446) -0.963032 (-0.964674)
+2.9375 0.204777 (0.202679) -0.977706 (-0.979245)
+3 0.143389 (0.14112) -0.988572 (-0.989992)
+3.0625 0.0814441 (0.0790102) -0.995588 (-0.996874)
+3.125 0.0191822 (0.0165919) -0.998725 (-0.999862)
+3.1875 -0.043154 (-0.0458912) -0.997973 (-0.998946)
+3.25 -0.105322 (-0.108195) -0.993334 (-0.99413)
+3.3125 -0.167079 (-0.170077) -0.984826 (-0.985431)
+3.375 -0.228184 (-0.231294) -0.972482 (-0.972884)
+3.4375 -0.288401 (-0.291608) -0.956351 (-0.956538)
+3.5 -0.347494 (-0.350783) -0.936495 (-0.936457)
+3.5625 -0.405233 (-0.408589) -0.912992 (-0.912719)
+3.625 -0.461393 (-0.464799) -0.885933 (-0.885416)
+3.6875 -0.515755 (-0.519194) -0.855423 (-0.854657)
+3.75 -0.568109 (-0.571561) -0.821582 (-0.820559)
+3.8125 -0.618249 (-0.621697) -0.784541 (-0.783258)
+3.875 -0.665981 (-0.669405) -0.744445 (-0.742898)
+3.9375 -0.711119 (-0.714499) -0.70145 (-0.699637)
+4 -0.753487 (-0.756802) -0.655723 (-0.653644)
+4.0625 -0.792919 (-0.796151) -0.607442 (-0.605098)
+4.125 -0.829263 (-0.832391) -0.556796 (-0.55419)
+4.1875 -0.862377 (-0.86538) -0.503981 (-0.501117)
+4.25 -0.892132 (-0.894989) -0.449204 (-0.446087)
+4.3125 -0.918412 (-0.921104) -0.392678 (-0.389316)
+4.375 -0.941115 (-0.943622) -0.334623 (-0.331024)
+4.4375 -0.960152 (-0.962455) -0.275265 (-0.27144)
+4.5 -0.975449 (-0.97753) -0.214835 (-0.210796)
+4.5625 -0.986946 (-0.988788) -0.153569 (-0.149328)
+4.625 -0.9946 (-0.996184) -0.0917056 (-0.0872778)
+4.6875 -0.99838 (-0.99969) -0.0294852 (-0.0248864)
+4.75 -0.998271 (-0.999293) 0.0328495 (0.0376022)
+4.8125 -0.994274 (-0.994993) 0.0950559 (0.0999439)
+4.875 -0.986405 (-0.986808) 0.156891 (0.161895)
+4.9375 -0.974693 (-0.974769) 0.218116 (0.223215)
+5 -0.959186 (-0.958924) 0.27849 (0.283662)
+5.0625 -0.939943 (-0.939335) 0.337779 (0.343002)
+5.125 -0.917038 (-0.916077) 0.395752 (0.401003)
+5.1875 -0.890563 (-0.889242) 0.452183 (0.457437)
+5.25 -0.860619 (-0.858934) 0.506852 (0.512085)
+5.3125 -0.827323 (-0.825273) 0.559547 (0.564734)
+5.375 -0.790805 (-0.788389) 0.610062 (0.615177)
+5.4375 -0.751207 (-0.748426) 0.658201 (0.663219)
+5.5 -0.708684 (-0.70554) 0.703776 (0.70867)
+5.5625 -0.663401 (-0.6599) 0.746609 (0.751354)
+5.625 -0.615534 (-0.611682) 0.786534 (0.791104)
+5.6875 -0.565271 (-0.561076) 0.823395 (0.827764)
+5.75 -0.512806 (-0.508279) 0.857049 (0.861192)
+5.8125 -0.458344 (-0.453497) 0.887365 (0.891258)
+5.875 -0.402097 (-0.396944) 0.914224 (0.917843)
+5.9375 -0.344285 (-0.338842) 0.937523 (0.940843)
+6 -0.285132 (-0.279415) 0.957169 (0.96017)
diff --git a/tests/petsc/petsc_ts_00_in.prm b/tests/petsc/petsc_ts_00_in.prm
new file mode 100644 (file)
index 0000000..5127d37
--- /dev/null
@@ -0,0 +1,10 @@
+subsection Running parameters
+  set initial time      = 1.0
+  set final time        = 6.0
+  set solver type       = bdf 
+  set initial step size = 0.0625
+end
+subsection Error control
+  set minimum step size = 1e-7
+  set maximum step size = 1e+7
+end
diff --git a/tests/petsc/petsc_ts_01.cc b/tests/petsc/petsc_ts_01.cc
new file mode 100644 (file)
index 0000000..b2aee68
--- /dev/null
@@ -0,0 +1,88 @@
+//-----------------------------------------------------------
+//
+//    Copyright (C) 2023 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.md at
+//    the top level directory of deal.II.
+//
+//-----------------------------------------------------------
+//
+// Author: Stefano Zampini, King Abdullah University of Science and Technology.
+
+#include <deal.II/lac/exceptions.h>
+#include <deal.II/lac/petsc_ts.templates.h>
+
+#include "../tests.h"
+
+/**
+ * Tests user defined Vector and Matrix types
+ * and exceptions handling for PETSCWrappers::TimeStepper.
+ */
+
+class VectorType : public Subscriptor
+{
+public:
+  explicit VectorType(Vec v)
+    : v(v)
+  {}
+
+  Vec &
+  petsc_vector()
+  {
+    return v;
+  }
+
+private:
+  Vec v;
+};
+
+class MatrixType : public Subscriptor
+{
+public:
+  explicit MatrixType(Mat A)
+    : A(A)
+  {}
+
+  Mat &
+  petsc_matrix()
+  {
+    return A;
+  }
+
+private:
+  Mat A;
+};
+
+using TimeStepper = PETScWrappers::TimeStepper<VectorType, MatrixType>;
+
+int
+main(int argc, char **argv)
+{
+  initlog();
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+  VectorType v(nullptr);
+  MatrixType A(nullptr);
+
+  TimeStepper myode;
+
+  try
+    {
+      auto ts = myode.petsc_ts();
+      auto t0 = myode.get_time();
+      auto dt = myode.get_time_step();
+      myode.solve(v, A);
+    }
+  catch (std::exception &exc)
+    {
+      deallog << exc.what() << std::endl;
+    }
+  return 0;
+}
diff --git a/tests/petsc/petsc_ts_01.output b/tests/petsc/petsc_ts_01.output
new file mode 100644 (file)
index 0000000..761c766
--- /dev/null
@@ -0,0 +1,12 @@
+
+DEAL::
+--------------------------------------------------------
+An error occurred in file <petsc_ts.templates.h> in function
+    unsigned int dealii::PETScWrappers::TimeStepper<VectorType, PMatrixType, AMatrixType>::solve(VectorType&) [with VectorType = VectorType; PMatrixType = MatrixType; AMatrixType = MatrixType]
+The violated condition was: 
+    explicit_function || implicit_function
+Additional information: 
+    Please provide an implementation for the function "explicit_function
+    || implicit_function"
+--------------------------------------------------------
+
diff --git a/tests/petsc/petsc_ts_02.cc b/tests/petsc/petsc_ts_02.cc
new file mode 100644 (file)
index 0000000..30ea46b
--- /dev/null
@@ -0,0 +1,92 @@
+//-----------------------------------------------------------
+//
+//    Copyright (C) 2023 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.md at
+//    the top level directory of deal.II.
+//
+//-----------------------------------------------------------
+//
+// Author: Stefano Zampini, King Abdullah University of Science and Technology.
+
+//#include <deal.II/base/parameter_handler.h>
+
+#include <deal.II/lac/petsc_ts.h>
+#include <deal.II/lac/petsc_vector.h>
+
+#include "../tests.h"
+
+/**
+ * Solves a dummy index-1 DAE to test the solver.
+ */
+using VectorType  = PETScWrappers::MPI::Vector;
+using TimeStepper = PETScWrappers::TimeStepper<VectorType>;
+using real_type   = TimeStepper::real_type;
+
+
+class Dummy : public TimeStepper
+{
+public:
+  Dummy()
+  {
+    // Customize solver: use BDF and basic adaptive time stepping
+    PETScWrappers::TimeStepperData data;
+    data.tstype      = "bdf";
+    data.final_time  = 1.0;
+    data.tsadapttype = "basic";
+    reinit(data);
+
+    // Here we solve the two variables system:
+    //   x_dot = x
+    //   z = x
+    // with y = (x,z)
+    implicit_function = [&](const real_type   t,
+                            const VectorType &y,
+                            const VectorType &y_dot,
+                            VectorType &      res) -> int {
+      res(0) = y_dot(0) - y(0);
+      res(1) = y(1) - y(0);
+      res.compress(VectorOperation::insert);
+      return 0;
+    };
+
+    // Return the index set representing the
+    // algebraic components
+    algebraic_components = [&]() -> IndexSet {
+      IndexSet algebraic_is;
+
+      algebraic_is.set_size(2);
+      algebraic_is.add_index(1);
+      return algebraic_is;
+    };
+  }
+
+  unsigned int
+  solve()
+  {
+    VectorType y(MPI_COMM_SELF, 2, 2);
+    y(0) = 1.0;
+    y(1) = y(0);
+    y.compress(VectorOperation::insert);
+
+    return TimeStepper::solve(y);
+  }
+};
+
+int
+main(int argc, char **argv)
+{
+  initlog();
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+  Dummy dae;
+  dae.solve();
+  return 0;
+}
diff --git a/tests/petsc/petsc_ts_02.output b/tests/petsc/petsc_ts_02.output
new file mode 100644 (file)
index 0000000..8b13789
--- /dev/null
@@ -0,0 +1 @@
+

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.