]> https://gitweb.dealii.org/ - dealii.git/commitdiff
PETScWrappers::TimeStepper fix backward compatibility to 3.7 15126/head
authorStefano Zampini <stefano.zampini@gmail.com>
Sun, 16 Apr 2023 13:35:46 +0000 (16:35 +0300)
committerStefano Zampini <stefano.zampini@gmail.com>
Fri, 21 Apr 2023 22:38:26 +0000 (01:38 +0300)
Align interface with SNES

include/deal.II/lac/petsc_ts.h
include/deal.II/lac/petsc_ts.templates.h
source/lac/petsc_ts.cc
tests/petsc/petsc_ts_00.cc
tests/petsc/petsc_ts_00.output
tests/petsc/petsc_ts_01.cc
tests/petsc/petsc_ts_02.cc

index a5f6153bae739212ca920a11ea4360be711885bb..3dc726407f56686bd52af5d0f03134e529bae55e 100644 (file)
@@ -1,6 +1,6 @@
 //-----------------------------------------------------------
 //
-//    Copyright (C) 2022 by the deal.II authors
+//    Copyright (C) 2023 by the deal.II authors
 //
 //    This file is part of the deal.II library.
 //
@@ -12,8 +12,6 @@
 //    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
@@ -53,8 +51,8 @@ namespace PETScWrappers
      *
      * 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 options_prefix The string indicating the options prefix for command line customization.
+     * @param ts_type 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.
@@ -76,13 +74,15 @@ namespace PETScWrappers
      * 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
+     * @note All parameters values specified here can be overriden by
      * command line choices.
+     *
+     * @ingroup PETScWrappers
      */
     TimeStepperData(
       // Running parameters
-      const std::string &opts_prefix       = "",
-      const std::string &tstype            = "",
+      const std::string &options_prefix    = "",
+      const std::string &ts_type           = "",
       const real_type    initial_time      = 0.0,
       const real_type    final_time        = 0.0,
       const real_type    initial_step_size = 0.0,
@@ -95,8 +95,8 @@ namespace PETScWrappers
       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)
+      : options_prefix(options_prefix)
+      , ts_type(ts_type)
       , initial_time(initial_time)
       , final_time(final_time)
       , initial_step_size(initial_step_size)
@@ -110,68 +110,21 @@ namespace PETScWrappers
       , ignore_algebraic_lte(ignore_algebraic_lte)
     {}
 
+    /**
+     * Import parameter values.
+     */
     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();
-    }
+    add_parameters(ParameterHandler &prm);
 
     /**
      * Options database prefix.
      */
-    std::string opts_prefix;
+    std::string options_prefix;
 
     /**
      * PETSc solver type.
      */
-    std::string tstype;
+    std::string ts_type;
 
     /**
      * Initial time for the DAE.
@@ -185,18 +138,20 @@ namespace PETScWrappers
 
     /**
      * Initial step size.
-     * Non-positive values are ignored.
+     *
+     * @note Non-positive values are ignored.
      */
     real_type initial_step_size;
 
     /**
      * Maximum number of steps to be taken.
-     * Negative values are ignored.
+     *
+     * @note Negative values are ignored.
      */
     int max_steps;
 
     /**
-     * Match final time requested?
+     * Flag to indicate to stop exactly at the requested final time.
      */
     bool match_step;
 
@@ -207,25 +162,29 @@ namespace PETScWrappers
 
     /**
      * Minimum allowed step size for adaptive time stepping.
-     * Non-positive values indicate to use PETSc's default.
+     *
+     * @note 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.
+     *
+     * @note 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.
+     *
+     * @note 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.
+     *
+     * @note Negative values indicate to use PETSc's default.
      */
     real_type relative_tolerance;
 
@@ -300,10 +259,10 @@ namespace PETScWrappers
    *
    * 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
+   *  - 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
+   *  - deal.II style using TimeStepper::setup_jacobian and
+   * TimeStepper::solve_with_jacobian
    *
    * In case both approaches are coded, the deal.II style
    * will be used.
@@ -315,7 +274,7 @@ namespace PETScWrappers
    * @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.
+   * still using TimeStepper::solve_with_jacobian 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
@@ -324,6 +283,8 @@ namespace PETScWrappers
    * @code
    * ./myApp -ksp_type cg -pc_type gamg
    * @endcode
+   *
+   * @ingroup PETScWrappers
    */
   template <typename VectorType  = PETScWrappers::VectorBase,
             typename PMatrixType = PETScWrappers::MatrixBase,
@@ -349,6 +310,13 @@ namespace PETScWrappers
      */
     ~TimeStepper();
 
+    /**
+     * Conversion operator to gain access to the underlying PETSc type. If you
+     * do this, you cut this class off some information it may need, so this
+     * conversion operator should only be used if you know what you do.
+     */
+    operator TS() const;
+
     /**
      * Return the PETSc TS object.
      */
@@ -369,7 +337,7 @@ namespace PETScWrappers
 
     /**
      * Reset solver.
-     * Change customization according to `data`.
+     * Change customization according to @p data.
      */
     void
     reinit(const TimeStepperData &data);
@@ -378,7 +346,7 @@ namespace PETScWrappers
      * Set the preconditioning matrix only.
      *
      * When used with TimeStepper::setup_jacobian and
-     * TimeStepper::solve_for_jacobian_system, PETSc will approximate the linear
+     * TimeStepper::solve_with_jacobian, PETSc will approximate the linear
      * system matrix-vector product using an internal matrix-free
      * representation.
      *
@@ -387,14 +355,14 @@ namespace PETScWrappers
      * preconditioning and matrix-vector products.
      */
     void
-    reinit_matrix(PMatrixType &P);
+    set_matrix(PMatrixType &P);
 
     /**
      * Set both the linear system matrix and the preconditioning matrix
      * that PETSc will use.
      */
     void
-    reinit_matrices(AMatrixType &A, PMatrixType &P);
+    set_matrices(AMatrixType &A, PMatrixType &P);
 
     /**
      * Return current time.
@@ -425,7 +393,7 @@ namespace PETScWrappers
      * 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.
+     * Here we also set the matrix to precondition the tangent system.
      */
     unsigned int
     solve(VectorType &y, PMatrixType &P);
@@ -437,7 +405,8 @@ namespace PETScWrappers
      * 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.
+     * Here we also set the matrices to describe and precondition
+     * the tangent system.
      */
     unsigned int
     solve(VectorType &y, AMatrixType &A, PMatrixType &P);
@@ -454,7 +423,7 @@ namespace PETScWrappers
     /**
      * Callback for the computation of the implicit Jacobian
      * $\dfrac{\partial F}{\partial y} + \alpha \dfrac{\partial F}{\partial \dot
-     * y}$
+     * y}$.
      *
      * All implicit solvers implementations are recast to use the above
      * linearization. The $\alpha$ parameter is time-step and solver-type
@@ -496,7 +465,7 @@ namespace PETScWrappers
       monitor;
 
     /**
-     * Set up Jacobian callback without matrices.
+     * Callback for the set up of the Jacobian system.
      *
      * This callback gives full control to users to set up the linearized
      * equations
@@ -507,7 +476,7 @@ namespace PETScWrappers
      * linearization. The $\alpha$ parameter is time-step and solver-type
      * specific.
      *
-     * Solvers must be provided via TimeStepper::solve_for_jacobian_system.
+     * Solvers must be provided via TimeStepper::solve_with_jacobian.
      */
     std::function<int(const real_type   t,
                       const VectorType &y,
@@ -516,15 +485,16 @@ namespace PETScWrappers
       setup_jacobian;
 
     /**
-     * Solution of the Jacobian system set up with TimeStepper::setup_jacobian.
+     * Callback for the solution of the tangent 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;
+      solve_with_jacobian;
 
     /**
-     * Return an index set containing the algebraic components.
+     * Callback to 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
@@ -533,55 +503,25 @@ namespace PETScWrappers
      */
     std::function<IndexSet()> algebraic_components;
 
-  protected:
+  private:
     /**
-     * The PETSc object
+     * 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));
-  }
+    /**
+     * This flag is set when changing the customization and used within solve.
+     */
+    bool need_dae_tolerances;
 
-  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
+    /**
+     * This flag is used to support versions of PETSc older than 3.13.
+     */
+    bool need_dummy_assemble;
+  };
 
 } // namespace PETScWrappers
 
index 93c0fb09c275cafa855a0c7494b347868c1c9824..91d3268c64aef26b14604ec6f89f61d4b0f87919 100644 (file)
@@ -1,6 +1,6 @@
 //-----------------------------------------------------------
 //
-//    Copyright (C) 2022 by the deal.II authors
+//    Copyright (C) 2023 by the deal.II authors
 //
 //    This file is part of the deal.II library.
 //
@@ -12,8 +12,6 @@
 //    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
@@ -32,7 +30,7 @@
 DEAL_II_NAMESPACE_OPEN
 
 // Shorthand notation for PETSc error codes.
-#  define AssertTS(code)                             \
+#  define AssertPETSc(code)                          \
     do                                               \
       {                                              \
         PetscErrorCode ierr = (code);                \
@@ -57,8 +55,8 @@ namespace PETScWrappers
     const TimeStepperData &data,
     const MPI_Comm &       mpi_comm)
   {
-    AssertTS(TSCreate(mpi_comm, &ts));
-    AssertTS(TSSetApplicationContext(ts, this));
+    AssertPETSc(TSCreate(mpi_comm, &ts));
+    AssertPETSc(TSSetApplicationContext(ts, this));
     reinit(data);
   }
 
@@ -67,7 +65,58 @@ namespace PETScWrappers
   template <typename VectorType, typename PMatrixType, typename AMatrixType>
   TimeStepper<VectorType, PMatrixType, AMatrixType>::~TimeStepper()
   {
-    AssertTS(TSDestroy(&ts));
+    AssertPETSc(TSDestroy(&ts));
+  }
+
+
+
+  template <typename VectorType, typename PMatrixType, typename AMatrixType>
+  TimeStepper<VectorType, PMatrixType, AMatrixType>::operator TS() const
+  {
+    return ts;
+  }
+
+
+
+  template <typename VectorType, typename PMatrixType, typename AMatrixType>
+  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>
+  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>
+  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;
   }
 
 
@@ -76,7 +125,7 @@ namespace PETScWrappers
   void
   TimeStepper<VectorType, PMatrixType, AMatrixType>::reinit()
   {
-    AssertTS(TSReset(ts));
+    AssertPETSc(TSReset(ts));
   }
 
 
@@ -89,27 +138,28 @@ namespace PETScWrappers
     reinit();
 
     // Solver type
-    if (data.tstype.size())
-      AssertTS(TSSetType(ts, data.tstype.c_str()));
+    if (data.ts_type.size())
+      AssertPETSc(TSSetType(ts, data.ts_type.c_str()));
 
     // Options prefix
-    if (data.opts_prefix.size())
-      AssertTS(TSSetOptionsPrefix(ts, data.opts_prefix.c_str()));
+    if (data.options_prefix.size())
+      AssertPETSc(TSSetOptionsPrefix(ts, data.options_prefix.c_str()));
 
     // Time and steps limits
-    AssertTS(TSSetTime(ts, data.initial_time));
+    AssertPETSc(TSSetTime(ts, data.initial_time));
     if (data.final_time > data.initial_time)
-      AssertTS(TSSetMaxTime(ts, data.final_time));
+      ts_set_max_time(ts, data.final_time);
     if (data.initial_step_size > 0.0)
-      AssertTS(TSSetTimeStep(ts, data.initial_step_size));
+      AssertPETSc(TSSetTimeStep(ts, data.initial_step_size));
     if (data.max_steps >= 0)
-      AssertTS(TSSetMaxSteps(ts, data.max_steps));
+      ts_set_max_steps(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));
+    AssertPETSc(TSSetExactFinalTime(ts,
+                                    data.match_step ?
+                                      TS_EXACTFINALTIME_MATCHSTEP :
+                                      TS_EXACTFINALTIME_STEPOVER));
 
     // Adaptive tolerances
     const PetscReal atol = data.absolute_tolerance > 0.0 ?
@@ -118,7 +168,7 @@ namespace PETScWrappers
     const PetscReal rtol = data.relative_tolerance > 0.0 ?
                              data.relative_tolerance :
                              static_cast<PetscReal>(PETSC_DEFAULT);
-    AssertTS(TSSetTolerances(ts, atol, nullptr, rtol, nullptr));
+    AssertPETSc(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
@@ -127,13 +177,14 @@ namespace PETScWrappers
 
     // Adaptive time stepping
     TSAdapt tsadapt;
-    AssertTS(TSGetAdapt(ts, &tsadapt));
-    AssertTS(TSAdaptSetType(tsadapt, data.tsadapttype.c_str()));
+    AssertPETSc(TSGetAdapt(ts, &tsadapt));
+    AssertPETSc(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()));
+    if (data.options_prefix.size())
+      AssertPETSc(
+        TSAdaptSetOptionsPrefix(tsadapt, data.options_prefix.c_str()));
 
     // Time step limits
     const PetscReal hmin = data.minimum_step_size > 0.0 ?
@@ -142,15 +193,14 @@ namespace PETScWrappers
     const PetscReal hmax = data.maximum_step_size > 0.0 ?
                              data.maximum_step_size :
                              static_cast<PetscReal>(PETSC_DEFAULT);
-    AssertTS(TSAdaptSetStepLimits(tsadapt, hmin, hmax));
+    AssertPETSc(TSAdaptSetStepLimits(tsadapt, hmin, hmax));
   }
 
 
 
   template <typename VectorType, typename PMatrixType, typename AMatrixType>
   void
-  TimeStepper<VectorType, PMatrixType, AMatrixType>::reinit_matrix(
-    PMatrixType &P)
+  TimeStepper<VectorType, PMatrixType, AMatrixType>::set_matrix(PMatrixType &P)
   {
     this->A = nullptr;
     this->P = &P;
@@ -160,7 +210,7 @@ namespace PETScWrappers
 
   template <typename VectorType, typename PMatrixType, typename AMatrixType>
   void
-  TimeStepper<VectorType, PMatrixType, AMatrixType>::reinit_matrices(
+  TimeStepper<VectorType, PMatrixType, AMatrixType>::set_matrices(
     AMatrixType &A,
     PMatrixType &P)
   {
@@ -174,11 +224,9 @@ namespace PETScWrappers
   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 {
+    const auto ts_ifunction =
+      [](TS, PetscReal t, Vec x, Vec xdot, Vec f, void *ctx) -> PetscErrorCode {
       PetscFunctionBeginUser;
-      (void)ts;
       auto user = static_cast<TimeStepper *>(ctx);
 
       VectorType xdealii(x);
@@ -186,19 +234,14 @@ namespace PETScWrappers
       VectorType fdealii(f);
       AssertUser(user->implicit_function(t, xdealii, xdotdealii, fdealii),
                  "implicit_function");
+      petsc_increment_state_counter(f);
       PetscFunctionReturn(0);
     };
 
-    auto ts_ijacobian_ = [](TS        ts,
-                            PetscReal t,
-                            Vec       x,
-                            Vec       xdot,
-                            PetscReal s,
-                            Mat       A,
-                            Mat       P,
-                            void *    ctx) -> PetscErrorCode {
+    const auto ts_ijacobian =
+      [](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);
@@ -209,30 +252,27 @@ namespace PETScWrappers
       AssertUser(
         user->implicit_jacobian(t, xdealii, xdotdealii, s, Adealii, Pdealii),
         "implicit_jacobian");
+      petsc_increment_state_counter(P);
 
       // 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));
+      AssertPETSc(PetscObjectTypeCompare((PetscObject)A, MATMFFD, &flg));
       if (flg)
         {
-          AssertTS(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
-          AssertTS(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
+          AssertPETSc(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
+          AssertPETSc(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
         }
+      else
+        petsc_increment_state_counter(A);
       PetscFunctionReturn(0);
     };
 
-    auto ts_ijacobian_with_setup_ = [](TS        ts,
-                                       PetscReal t,
-                                       Vec       x,
-                                       Vec       xdot,
-                                       PetscReal s,
-                                       Mat       A,
-                                       Mat       P,
-                                       void *    ctx) -> PetscErrorCode {
+    const auto ts_ijacobian_with_setup =
+      [](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);
@@ -244,24 +284,38 @@ namespace PETScWrappers
       user->P = &Pdealii;
       AssertUser(user->setup_jacobian(t, xdealii, xdotdealii, s),
                  "setup_jacobian");
+      petsc_increment_state_counter(P);
+
+      // Handle older versions of PETSc for which we cannot pass a MATSHELL
+      // matrix to DMSetMatType. This has been fixed from 3.13 on.
+      // What we need to do for older version of PETSc is instead to have
+      // a zero matrix with all diagonal entries present.
+      if (user->need_dummy_assemble)
+        {
+          AssertPETSc(MatZeroEntries(P));
+          AssertPETSc(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
+          AssertPETSc(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
+          AssertPETSc(MatShift(P, 0.0));
+        }
 
       // 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));
+      AssertPETSc(PetscObjectTypeCompare((PetscObject)A, MATMFFD, &flg));
       if (flg)
         {
-          AssertTS(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
-          AssertTS(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
+          AssertPETSc(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
+          AssertPETSc(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
         }
+      else
+        petsc_increment_state_counter(A);
       PetscFunctionReturn(0);
     };
 
-    auto ts_rhsfunction_ =
-      [](TS ts, PetscReal t, Vec x, Vec f, void *ctx) -> PetscErrorCode {
+    const auto ts_rhsfunction =
+      [](TS, PetscReal t, Vec x, Vec f, void *ctx) -> PetscErrorCode {
       PetscFunctionBeginUser;
-      (void)ts;
       auto user = static_cast<TimeStepper *>(ctx);
 
       VectorType xdealii(x);
@@ -269,13 +323,13 @@ namespace PETScWrappers
 
       AssertUser(user->explicit_function(t, xdealii, fdealii),
                  "explicit_function");
+      petsc_increment_state_counter(f);
       PetscFunctionReturn(0);
     };
 
-    auto ts_rhsjacobian_ =
-      [](TS ts, PetscReal t, Vec x, Mat A, Mat P, void *ctx) -> PetscErrorCode {
+    const auto ts_rhsjacobian =
+      [](TS, PetscReal t, Vec x, Mat A, Mat P, void *ctx) -> PetscErrorCode {
       PetscFunctionBeginUser;
-      (void)ts;
       auto user = static_cast<TimeStepper *>(ctx);
 
       VectorType  xdealii(x);
@@ -284,24 +338,36 @@ namespace PETScWrappers
 
       AssertUser(user->explicit_jacobian(t, xdealii, Adealii, Pdealii),
                  "explicit_jacobian");
+      petsc_increment_state_counter(P);
+
+      // Handle older versions of PETSc for which we cannot pass a MATSHELL
+      // matrix to DMSetMatType
+      if (user->need_dummy_assemble)
+        {
+          AssertPETSc(MatZeroEntries(P));
+          AssertPETSc(MatAssemblyBegin(P, MAT_FINAL_ASSEMBLY));
+          AssertPETSc(MatAssemblyEnd(P, MAT_FINAL_ASSEMBLY));
+          AssertPETSc(MatShift(P, 0.0));
+        }
 
       // 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));
+      AssertPETSc(PetscObjectTypeCompare((PetscObject)A, MATMFFD, &flg));
       if (flg)
         {
-          AssertTS(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
-          AssertTS(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
+          AssertPETSc(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
+          AssertPETSc(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
         }
+      else
+        petsc_increment_state_counter(A);
       PetscFunctionReturn(0);
     };
 
-    auto ts_monitor_ =
-      [](TS ts, PetscInt it, PetscReal t, Vec x, void *ctx) -> PetscErrorCode {
+    const auto ts_monitor =
+      [](TS, PetscInt it, PetscReal t, Vec x, void *ctx) -> PetscErrorCode {
       PetscFunctionBeginUser;
-      (void)ts;
       auto user = static_cast<TimeStepper *>(ctx);
 
       VectorType xdealii(x);
@@ -313,56 +379,61 @@ namespace PETScWrappers
                 StandardExceptions::ExcFunctionNotProvided(
                   "explicit_function || implicit_function"));
 
-    AssertTS(TSSetSolution(ts, y.petsc_vector()));
+    AssertPETSc(TSSetSolution(ts, y.petsc_vector()));
 
     if (explicit_function)
-      AssertTS(TSSetRHSFunction(ts, nullptr, ts_rhsfunction_, this));
+      AssertPETSc(TSSetRHSFunction(ts, nullptr, ts_rhsfunction, this));
 
     if (implicit_function)
-      AssertTS(TSSetIFunction(ts, nullptr, ts_ifunction_, this));
+      AssertPETSc(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));
+        AssertPETSc(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));
+          set_use_matrix_free(ts, true, false);
 
         // Do not waste memory by creating a dummy AIJ matrix inside PETSc.
+        this->need_dummy_assemble = false;
         if (!P)
           {
-            DM dm;
-            AssertTS(SNESGetDM(snes, &dm));
-            AssertTS(DMSetMatType(dm, MATSHELL));
+#  if DEAL_II_PETSC_VERSION_GTE(3, 13, 0)
+            DM   dm;
+            SNES snes;
+            AssertPETSc(TSGetSNES(ts, &snes));
+            AssertPETSc(SNESGetDM(snes, &dm));
+            AssertPETSc(DMSetMatType(dm, MATSHELL));
+#  else
+            this->need_dummy_assemble = true;
+#  endif
           }
       }
     else
       {
         if (explicit_jacobian)
           {
-            AssertTS(TSSetRHSJacobian(ts,
-                                      A ? A->petsc_matrix() :
-                                          (P ? P->petsc_matrix() : nullptr),
-                                      P ? P->petsc_matrix() : nullptr,
-                                      ts_rhsjacobian_,
-                                      this));
+            AssertPETSc(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));
+            AssertPETSc(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
@@ -372,55 +443,53 @@ namespace PETScWrappers
         // be overriden from command line.
         if (!explicit_jacobian && !implicit_jacobian)
           {
-            SNES snes;
-            AssertTS(TSGetSNES(ts, &snes));
-            AssertTS(SNESSetUseMatrixFree(snes, PETSC_FALSE, PETSC_TRUE));
+            set_use_matrix_free(ts, false, true);
           }
       }
 
-    // In case solve_for_jacobian_system is provided, create a shell
+    // In case solve_with_jacobian 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
+    // Using solve_with_jacobian 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)
+    if (solve_with_jacobian)
       {
         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);
+          return solve_with_jacobian(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()));
+        AssertPETSc(TSGetSNES(ts, &snes));
+        AssertPETSc(SNESGetKSP(snes, &ksp));
+        AssertPETSc(KSPSetType(ksp, KSPPREONLY));
+        AssertPETSc(KSPSetPC(ksp, precond.get_pc()));
       }
 
     // Attach user monitoring routine.
     if (monitor)
-      AssertTS(TSMonitorSet(ts, ts_monitor_, this, nullptr));
+      AssertPETSc(TSMonitorSet(ts, ts_monitor, this, nullptr));
 
     // Allow command line customization.
-    AssertTS(TSSetFromOptions(ts));
+    AssertPETSc(TSSetFromOptions(ts));
 
     // Handle algebraic components.
     if (algebraic_components && need_dae_tolerances)
       {
         PetscReal atol, rtol;
-        AssertTS(TSGetTolerances(ts, &atol, nullptr, &rtol, nullptr));
+        AssertPETSc(TSGetTolerances(ts, &atol, nullptr, &rtol, nullptr));
 
         Vec av, rv;
-        AssertTS(VecDuplicate(y.petsc_vector(), &av));
-        AssertTS(VecDuplicate(y.petsc_vector(), &rv));
+        AssertPETSc(VecDuplicate(y.petsc_vector(), &av));
+        AssertPETSc(VecDuplicate(y.petsc_vector(), &rv));
 
         VectorBase avdealii(av);
         VectorBase rvdealii(rv);
@@ -433,18 +502,27 @@ namespace PETScWrappers
           }
         avdealii.compress(VectorOperation::insert);
         rvdealii.compress(VectorOperation::insert);
-        AssertTS(TSSetTolerances(ts, atol, av, rtol, rv));
-        AssertTS(VecDestroy(&av));
-        AssertTS(VecDestroy(&rv));
+        AssertPETSc(TSSetTolerances(ts, atol, av, rtol, rv));
+        AssertPETSc(VecDestroy(&av));
+        AssertPETSc(VecDestroy(&rv));
       }
 
     // Having set everything up, now do the actual work
     // and let PETSc do the time stepping.
-    AssertTS(TSSolve(ts, nullptr));
+    AssertPETSc(TSSolve(ts, nullptr));
+
+    // Get the number of steps taken.
+    auto nt = ts_get_step_number(ts);
+
+    // Raise an exception if the solver has not converged
+    TSConvergedReason reason;
+    AssertPETSc(TSGetConvergedReason(ts, &reason));
+    AssertThrow(reason > 0,
+                ExcMessage("TS solver did not converge after " +
+                           std::to_string(nt) + " iterations with reason " +
+                           TSConvergedReasons[reason]));
 
-    // Return the number of steps taken.
-    PetscInt nt;
-    AssertTS(TSGetStepNumber(ts, &nt));
+    // Finally return
     return nt;
   }
 
@@ -455,7 +533,7 @@ namespace PETScWrappers
   TimeStepper<VectorType, PMatrixType, AMatrixType>::solve(VectorType & y,
                                                            PMatrixType &P)
   {
-    reinit_matrix(P);
+    set_matrix(P);
     return solve(y);
   }
 
@@ -467,13 +545,13 @@ namespace PETScWrappers
                                                            AMatrixType &A,
                                                            PMatrixType &P)
   {
-    reinit_matrices(A, P);
+    set_matrices(A, P);
     return solve(y);
   }
 
 } // namespace PETScWrappers
 
-#  undef AssertTS
+#  undef AssertPETSc
 #  undef AssertUser
 
 DEAL_II_NAMESPACE_CLOSE
index d7018087d2e8036e13a24cb097df5bbe706c1af9..f54cdc0b2075e8c5aecf711ed59a7de0cbc84efd 100644 (file)
 
 DEAL_II_NAMESPACE_OPEN
 
+namespace PETScWrappers
+{
+  void
+  TimeStepperData::add_parameters(ParameterHandler &prm)
+  {
+    prm.enter_subsection("Running parameters");
+    prm.add_parameter(
+      "options prefix",
+      options_prefix,
+      "The string indicating the options prefix for command line customization.");
+    prm.add_parameter("solver type",
+                      ts_type,
+                      "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();
+  }
+
+} // namespace PETScWrappers
 
 template class PETScWrappers::TimeStepper<>;
 template class PETScWrappers::TimeStepper<PETScWrappers::MPI::Vector>;
index d95ae8aad624fa3306a11368bbd74c02daad7adb..98f197bc0d6bc00fbbda8a631ef3c01f49a94348 100644 (file)
@@ -1,6 +1,6 @@
 //-----------------------------------------------------------
 //
-//    Copyright (C) 2022 by the deal.II authors
+//    Copyright (C) 2023 by the deal.II authors
 //
 //    This file is part of the deal.II library.
 //
@@ -12,8 +12,6 @@
 //    the top level directory of deal.II.
 //
 //-----------------------------------------------------------
-//
-// Author: Stefano Zampini, King Abdullah University of Science and Technology.
 
 #include <deal.II/base/parameter_handler.h>
 
@@ -136,8 +134,8 @@ public:
 
             // 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 {
+            time_stepper.solve_with_jacobian = [&](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));
@@ -207,7 +205,8 @@ public:
     y.compress(VectorOperation::insert);
 
     // Integrate the ODE.
-    time_stepper.solve(y);
+    auto nt = time_stepper.solve(y);
+    out << "# Number of steps taken: " << nt << std::endl;
   }
 
 private:
index 3eba9270de3fecbfcef3f13607c4203911a42a80..c6233771c3ac5f329028f6282e05b8c7707c1973 100644 (file)
@@ -116,6 +116,7 @@ end
 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)
+# Number of steps taken: 80
 # Test implicit interface (J 0)
 1 0.841471 (0.841471) 0.540302 (0.540302)
 1.0625 0.873 (0.873575) 0.486385 (0.48669)
@@ -198,6 +199,7 @@ end
 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)
+# Number of steps taken: 80
 # Test explicit interface (J 1)
 1 0.841471 (0.841471) 0.540302 (0.540302)
 1.0625 0.873 (0.873575) 0.486385 (0.48669)
@@ -280,6 +282,7 @@ end
 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)
+# Number of steps taken: 80
 # Test implicit interface (J 1)
 1 0.841471 (0.841471) 0.540302 (0.540302)
 1.0625 0.873 (0.873575) 0.486385 (0.48669)
@@ -362,6 +365,7 @@ end
 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)
+# Number of steps taken: 80
 # Test user interface
 1 0.841471 (0.841471) 0.540302 (0.540302)
 1.0625 0.873 (0.873575) 0.486385 (0.48669)
@@ -444,3 +448,4 @@ end
 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)
+# Number of steps taken: 80
index b2aee68d48179d30b3c34cb61daf7fe4921c2200..a20c22d47ccd2e237a7fdfda3b20af20499b63f3 100644 (file)
@@ -12,8 +12,6 @@
 //    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>
@@ -78,6 +76,7 @@ main(int argc, char **argv)
       auto ts = myode.petsc_ts();
       auto t0 = myode.get_time();
       auto dt = myode.get_time_step();
+      AssertThrow(ts == static_cast<TS>(myode), ExcInternalError());
       myode.solve(v, A);
     }
   catch (std::exception &exc)
index 30ea46b1db963975020b54f530078277f614b0e7..43315dc92f63d1d6c1fc4151c55edcbde0315d4b 100644 (file)
@@ -12,8 +12,6 @@
 //    the top level directory of deal.II.
 //
 //-----------------------------------------------------------
-//
-// Author: Stefano Zampini, King Abdullah University of Science and Technology.
 
 //#include <deal.II/base/parameter_handler.h>
 
@@ -37,7 +35,7 @@ public:
   {
     // Customize solver: use BDF and basic adaptive time stepping
     PETScWrappers::TimeStepperData data;
-    data.tstype      = "bdf";
+    data.ts_type     = "bdf";
     data.final_time  = 1.0;
     data.tsadapttype = "basic";
     reinit(data);

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.