petsc_ts.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2023 - 2024 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.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,
      const bool         restart_if_remesh = 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)
      , restart_if_remesh(restart_if_remesh)
      , 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;
 
    bool restart_if_remesh;
 
    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
    get_step_number();
 
    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;
 
    DEAL_II_DEPRECATED_EARLY
    std::function<void(const real_type t, VectorType &y)> distribute;
 
    std::function<void(const real_type t, VectorType &y)>
      update_constrained_components;
 
    DEAL_II_DEPRECATED_EARLY
    std::function<void(const real_type    t,
                       const unsigned int step,
                       const VectorType  &y,
                       bool              &resize)>
      decide_for_coarsening_and_refinement;
 
    std::function<
      bool(const real_type t, const unsigned int step, const VectorType &y)>
      decide_and_prepare_for_remeshing;
 
    DEAL_II_DEPRECATED_EARLY
    std::function<void(const std::vector<VectorType> &all_in,
                       std::vector<VectorType>       &all_out)>
      interpolate;
 
    std::function<void(const real_type                t,
                       const std::vector<VectorType> &all_in,
                       std::vector<VectorType>       &all_out)>
      transfer_solution_vectors_to_new_mesh;
 
  private:
    TS ts;
 
    SmartPointer<AMatrixType, TimeStepper> A;
    SmartPointer<PMatrixType, TimeStepper> P;
 
    PreconditionShell solve_with_jacobian_pc;
 
    bool restart_if_remesh;
 
    bool need_dae_tolerances;
 
    bool need_dummy_assemble;
 
    mutable std::exception_ptr pending_exception;
 
    bool error_in_function;
 
    void
    setup_callbacks();
 
    void
    setup_algebraic_constraints(const VectorType &y);
  };
 
} // namespace PETScWrappers
 
DEAL_II_NAMESPACE_CLOSE
 
#endif // DEAL_II_WITH_PETSC
 
#endif