petsc_ts.h

Go to the documentation of this file.

//-----------------------------------------------------------
//
//    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 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>
 
#  include <exception>
 
#  if defined(DEAL_II_HAVE_CXX20)
#    include <concepts>
#  endif
 
DEAL_II_NAMESPACE_OPEN
 
namespace PETScWrappers
{
  class TimeStepperData
  {
  public:
    using real_type = PetscReal;
 
    TimeStepperData(
      // Running parameters
      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,
      const int          max_steps         = -1,
      const bool         match_step        = false,
      // Error parameters
      const std::string &ts_adapt_type        = "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)
      : options_prefix(options_prefix)
      , ts_type(ts_type)
      , initial_time(initial_time)
      , final_time(final_time)
      , initial_step_size(initial_step_size)
      , max_steps(max_steps)
      , match_step(match_step)
      , ts_adapt_type(ts_adapt_type)
      , 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);
 
    std::string options_prefix;
 
    std::string ts_type;
 
    real_type initial_time;
 
    real_type final_time;
 
    real_type initial_step_size;
 
    int max_steps;
 
    bool match_step;
 
    std::string ts_adapt_type;
 
    real_type minimum_step_size;
 
    real_type maximum_step_size;
 
    real_type absolute_tolerance;
 
    real_type relative_tolerance;
 
    bool ignore_algebraic_lte;
  };
 
  template <typename VectorType  = PETScWrappers::VectorBase,
            typename PMatrixType = PETScWrappers::MatrixBase,
            typename AMatrixType = PMatrixType>
#  if defined(DEAL_II_HAVE_CXX20) && \
    !defined(DEAL_II_DOXYGEN_DO_NOT_PARSE_REQUIRES_CLAUSES)
  requires(
    (concepts::is_dealii_petsc_vector_type<VectorType> ||
     std::constructible_from<
       VectorType,
       Vec>)&&(concepts::is_dealii_petsc_matrix_type<PMatrixType> ||
               std::constructible_from<
                 PMatrixType,
                 Mat>)&&(concepts::is_dealii_petsc_matrix_type<AMatrixType> ||
                         std::constructible_from<AMatrixType, Mat>))
#  endif
    class TimeStepper
  {
  public:
    using real_type = PetscReal;
 
    TimeStepper(const TimeStepperData &data     = TimeStepperData(),
                const MPI_Comm         mpi_comm = PETSC_COMM_WORLD);
 
    ~TimeStepper();
 
    operator TS() const;
 
    TS
    petsc_ts();
 
    MPI_Comm
    get_mpi_communicator() const;
 
    void
    reinit();
 
    void
    reinit(const TimeStepperData &data);
 
    void
    set_matrix(PMatrixType &P);
 
    void
    set_matrices(AMatrixType &A, PMatrixType &P);
 
    real_type
    get_time();
 
    real_type
    get_time_step();
 
    unsigned int
    solve(VectorType &y);
 
    unsigned int
    solve(VectorType &y, PMatrixType &P);
 
    unsigned int
    solve(VectorType &y, AMatrixType &A, PMatrixType &P);
 
    std::function<void(const real_type   t,
                       const VectorType &y,
                       const VectorType &y_dot,
                       VectorType &      res)>
      implicit_function;
 
    std::function<void(const real_type   t,
                       const VectorType &y,
                       const VectorType &y_dot,
                       const real_type   alpha,
                       AMatrixType &     A,
                       PMatrixType &     P)>
      implicit_jacobian;
 
    std::function<void(const real_type t, const VectorType &y, VectorType &res)>
      explicit_function;
 
    std::function<void(const real_type   t,
                       const VectorType &y,
                       AMatrixType &     A,
                       PMatrixType &     P)>
      explicit_jacobian;
 
    std::function<void(const real_type    t,
                       const VectorType & y,
                       const unsigned int step_number)>
      monitor;
 
    std::function<void(const real_type   t,
                       const VectorType &y,
                       const VectorType &ydot,
                       const real_type   alpha)>
      setup_jacobian;
 
    std::function<void(const VectorType &src, VectorType &dst)>
      solve_with_jacobian;
 
    std::function<IndexSet()> algebraic_components;
 
  private:
    TS ts;
 
    SmartPointer<AMatrixType, TimeStepper> A;
    SmartPointer<PMatrixType, TimeStepper> P;
 
    bool need_dae_tolerances;
 
    bool need_dummy_assemble;
 
    mutable std::exception_ptr pending_exception;
 
    bool error_in_function;
  };
 
} // namespace PETScWrappers
 
DEAL_II_NAMESPACE_CLOSE
 
#endif // DEAL_II_WITH_PETSC
 
#endif