arkode_stepper.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) 2026 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Detailed license information governing the source code and contributions
// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
//
// -----------------------------------------------------------------------------
 
 
#ifndef dealii_sundials_arkode_stepper_h
#define dealii_sundials_arkode_stepper_h
 
#include <deal.II/base/config.h>
 
 
#ifdef DEAL_II_WITH_SUNDIALS
 
#  include <deal.II/base/conditional_ostream.h>
#  include <deal.II/base/exceptions.h>
#  include <deal.II/base/mpi_stub.h>
#  include <deal.II/base/parameter_handler.h>
 
#  ifdef DEAL_II_WITH_PETSC
#    include <deal.II/lac/petsc_block_vector.h>
#    include <deal.II/lac/petsc_vector.h>
#  endif
#  include <deal.II/lac/vector.h>
#  include <deal.II/lac/vector_memory.h>
 
#  include <arkode/arkode.h>
#  include <nvector/nvector_serial.h>
#  ifdef DEAL_II_WITH_MPI
#    include <nvector/nvector_parallel.h>
#  endif
#  include <deal.II/base/discrete_time.h>
 
#  include <deal.II/sundials/invocation_context.h>
#  include <deal.II/sundials/n_vector.h>
#  include <deal.II/sundials/sundials_types.h>
#  include <deal.II/sundials/sunlinsol_wrapper.h>
#  include <deal.II/sundials/sunnonlinsol_wrapper.h>
 
#  include <boost/signals2.hpp>
 
#  include <sundials/sundials_linearsolver.h>
#  include <sundials/sundials_math.h>
#  include <sundials/sundials_nonlinearsolver.h>
 
#  include <complex>
#  include <exception>
#  include <memory>
#  include <optional>
#  include <string>
 
#endif // DEAL_II_WITH_SUNDIALS
 
DEAL_II_NAMESPACE_OPEN
 
#ifdef DEAL_II_WITH_SUNDIALS
 
namespace SUNDIALS
{
  template <typename VectorType>
  class ARKode;
 
  template <typename VectorType>
  class ARKodeStepper
  {
  public:
    friend class ARKode<VectorType>;
 
    virtual ~ARKodeStepper() = default;
 
    virtual void *
    get_arkode_memory() const = 0;
 
  private:
    void
    attach(std::exception_ptr &pending_exception
#  if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
           ,
           const internal::SharedSUNContext &context
#  endif
    )
    {
      Assert(pending_exception_ptr == nullptr,
             ExcMessage("This stepper is already attached to an ARKode "
                        "object."));
      pending_exception_ptr = &pending_exception;
#  if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
      sundials_context = context;
#  endif
    }
 
    virtual void
    reinit(double t0, const VectorType &y0) = 0;
 
    virtual void
    resize(double t0, const VectorType &y0, double hscale) = 0;
 
  protected:
    template <typename Stepper>
    struct CallbackContext
    {
      Stepper            *stepper           = nullptr;
      std::exception_ptr *pending_exception = nullptr;
    };
 
    using ARKodeMemoryPtr = std::unique_ptr<void, void (*)(void *)>;
 
    internal::InvocationContext
    invocation_context() const
    {
      Assert(pending_exception_ptr != nullptr,
             ExcMessage("This stepper is not attached to an ARKode object."));
      return internal::InvocationContext
      {
        *pending_exception_ptr
#  if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
          ,
          sundials_context
#  endif
      };
    }
 
  private:
    std::exception_ptr *pending_exception_ptr = nullptr;
 
#  if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
    internal::SharedSUNContext sundials_context;
#  endif
  };
 
  template <typename VectorType = Vector<double>>
  class ARKStepper : public ARKodeStepper<VectorType>
  {
  public:
    class AdditionalData
    {
    public:
      explicit AdditionalData(
        const unsigned int order                                 = 0,
        const unsigned int maximum_non_linear_iterations         = 10,
        const bool         implicit_function_is_linear           = false,
        const bool         implicit_function_is_time_independent = false,
        const bool         mass_is_time_independent              = false,
        const int          anderson_acceleration_subspace        = 3,
        const std::string &implicit_butcher_table                = "",
        const std::string &explicit_butcher_table                = "");
 
      void
      add_parameters(ParameterHandler &prm);
 
      unsigned int order;
 
      unsigned int maximum_non_linear_iterations;
 
      bool implicit_function_is_linear;
 
      bool implicit_function_is_time_independent;
 
      bool mass_is_time_independent;
 
      int anderson_acceleration_subspace;
 
      std::string implicit_butcher_table;
 
      std::string explicit_butcher_table;
    };
 
    ARKStepper(const AdditionalData &data = AdditionalData());
 
    void *
    get_arkode_memory() const override;
 
    std::function<
      void(const double t, const VectorType &y, VectorType &explicit_f)>
      explicit_function;
 
    std::function<void(const double t, const VectorType &y, VectorType &res)>
      implicit_function;
 
    std::function<void(const double t, const VectorType &v, VectorType &Mv)>
      mass_times_vector;
 
    std::function<void(const double t)> mass_times_vector_setup;
 
    std::function<void(const VectorType &v,
                       VectorType       &Jv,
                       const double      t,
                       const VectorType &y,
                       const VectorType &fy)>
      jacobian_times_vector;
 
    std::function<
      void(const double t, const VectorType &y, const VectorType &fy)>
      jacobian_times_vector_setup;
 
    LinearSolveFunction<VectorType> solve_linearized_system;
 
    LinearSolveFunction<VectorType> solve_mass;
 
    std::function<void(const double      t,
                       const VectorType &y,
                       const VectorType &fy,
                       const VectorType &r,
                       VectorType       &z,
                       const double      gamma,
                       const double      tol,
                       const int         lr)>
      jacobian_preconditioner_solve;
 
    std::function<void(const double      t,
                       const VectorType &y,
                       const VectorType &fy,
                       const int         jok,
                       int              &jcur,
                       const double      gamma)>
      jacobian_preconditioner_setup;
 
    std::function<void(const double      t,
                       const VectorType &r,
                       VectorType       &z,
                       const double      tol,
                       const int         lr)>
      mass_preconditioner_solve;
 
    std::function<void(const double t)> mass_preconditioner_setup;
 
    std::function<void(void *arkode_mem)> custom_setup;
 
  private:
    using CallbackContext = typename ARKodeStepper<
      VectorType>::template CallbackContext<ARKStepper<VectorType>>;
 
    using ARKodeMemoryPtr = typename ARKodeStepper<VectorType>::ARKodeMemoryPtr;
 
    void
    reinit(double t0, const VectorType &y0) override;
 
    void
    resize(double t0, const VectorType &y0, double hscale) override;
 
    void
    setup_system_solver(const VectorType &solution);
 
    void
    setup_mass_solver(const VectorType &solution);
 
    AdditionalData data;
 
    std::optional<internal::LinearSolverWrapper<VectorType>> linear_solver;
 
    std::optional<internal::LinearSolverWrapper<VectorType>> mass_solver;
 
    std::optional<internal::SundialsLinearSolver> default_linear_solver;
 
    std::optional<internal::SundialsLinearSolver> default_mass_solver;
 
    std::optional<internal::SundialsNonlinearSolver> fixed_point_solver;
 
    ARKodeMemoryPtr arkode_mem;
 
    CallbackContext callback_ctx;
  };
 
 
  template <typename VectorType>
  ARKStepper<VectorType>::AdditionalData::AdditionalData(
    const unsigned int order,
    const unsigned int maximum_non_linear_iterations,
    const bool         implicit_function_is_linear,
    const bool         implicit_function_is_time_independent,
    const bool         mass_is_time_independent,
    const int          anderson_acceleration_subspace,
    const std::string &implicit_butcher_table,
    const std::string &explicit_butcher_table)
    : order(order)
    , maximum_non_linear_iterations(maximum_non_linear_iterations)
    , implicit_function_is_linear(implicit_function_is_linear)
    , implicit_function_is_time_independent(
        implicit_function_is_time_independent)
    , mass_is_time_independent(mass_is_time_independent)
    , anderson_acceleration_subspace(anderson_acceleration_subspace)
    , implicit_butcher_table(implicit_butcher_table)
    , explicit_butcher_table(explicit_butcher_table)
  {}
 
 
 
  template <typename VectorType>
  void
  ARKStepper<VectorType>::AdditionalData::add_parameters(ParameterHandler &prm)
  {
    prm.add_parameter("Integration accuracy order", order);
    prm.add_parameter("Maximum number of nonlinear iterations",
                      maximum_non_linear_iterations);
    prm.add_parameter("Implicit function is linear",
                      implicit_function_is_linear);
    prm.add_parameter("Implicit function is time independent",
                      implicit_function_is_time_independent);
    prm.add_parameter("Mass is time independent", mass_is_time_independent);
    prm.add_parameter("Anderson-acceleration subspace",
                      anderson_acceleration_subspace);
    prm.add_parameter("Implicit Butcher table", implicit_butcher_table);
    prm.add_parameter("Explicit Butcher table", explicit_butcher_table);
  }
 
 
  template <typename VectorType = Vector<double>>
  class ERKStepper : public ARKodeStepper<VectorType>
  {
  public:
    class AdditionalData
    {
    public:
      AdditionalData(const unsigned int order                  = 0,
                     const std::string &explicit_butcher_table = "");
 
      void
      add_parameters(ParameterHandler &prm);
 
      unsigned int order;
 
      std::string explicit_butcher_table;
    };
 
    ERKStepper(const AdditionalData &data = AdditionalData());
 
    void *
    get_arkode_memory() const override;
 
    std::function<
      void(const double t, const VectorType &y, VectorType &explicit_f)>
      explicit_function;
 
    std::function<void(void *arkode_mem)> custom_setup;
 
  private:
    using CallbackContext = typename ARKodeStepper<
      VectorType>::template CallbackContext<ERKStepper<VectorType>>;
 
    using ARKodeMemoryPtr = typename ARKodeStepper<VectorType>::ARKodeMemoryPtr;
 
    void
    reinit(double t0, const VectorType &y0) override;
 
    void
    resize(double t0, const VectorType &y0, double hscale) override;
 
    AdditionalData data;
 
    ARKodeMemoryPtr arkode_mem;
 
    CallbackContext callback_ctx;
  };
 
 
  template <typename VectorType>
  ERKStepper<VectorType>::AdditionalData::AdditionalData(
    const unsigned int order,
    const std::string &explicit_butcher_table)
    : order(order)
    , explicit_butcher_table(explicit_butcher_table)
  {}
 
 
 
  template <typename VectorType>
  void
  ERKStepper<VectorType>::AdditionalData::add_parameters(ParameterHandler &prm)
  {
    prm.add_parameter("Integration accuracy order", order);
    prm.add_parameter("Explicit Butcher table", explicit_butcher_table);
  }
 
 
#  if DEAL_II_SUNDIALS_VERSION_GTE(7, 2, 0)
 
  template <typename VectorType = Vector<double>>
  class LSRKStepperSTS : public ARKodeStepper<VectorType>
  {
  public:
    class AdditionalData
    {
    public:
      AdditionalData(
        const std::string &method_name       = "",
        const unsigned int dom_eig_frequency = numbers::invalid_unsigned_int,
        const unsigned int max_num_stages    = 0,
        const unsigned int dom_eig_estimator_max_iters = 0,
        const double       dom_eig_estimator_rel_tol   = 0,
        const unsigned int dom_eig_estimator_num_warmups =
          numbers::invalid_unsigned_int)
        : method_name(method_name)
        , dom_eig_frequency(dom_eig_frequency)
        , max_num_stages(max_num_stages)
        , dom_eig_estimator_max_iters(dom_eig_estimator_max_iters)
        , dom_eig_estimator_rel_tol(dom_eig_estimator_rel_tol)
        , dom_eig_estimator_num_warmups(dom_eig_estimator_num_warmups)
      {}
 
      void
      add_parameters(ParameterHandler &prm);
 
      std::string method_name;
 
      unsigned int dom_eig_frequency;
 
      unsigned int max_num_stages;
 
      unsigned int dom_eig_estimator_max_iters;
 
      double dom_eig_estimator_rel_tol;
 
      unsigned int dom_eig_estimator_num_warmups;
    };
 
    LSRKStepperSTS(const AdditionalData &data = AdditionalData());
 
    void *
    get_arkode_memory() const override;
 
    std::function<
      void(const double t, const VectorType &y, VectorType &explicit_f)>
      explicit_function;
 
    std::function<std::complex<double>(const double      t,
                                       const VectorType &y,
                                       const VectorType &f)>
      dominant_eigenvalue_function;
 
    std::function<void(void *arkode_mem)> custom_setup;
 
  private:
    using CallbackContext = typename ARKodeStepper<
      VectorType>::template CallbackContext<LSRKStepperSTS<VectorType>>;
 
    using ARKodeMemoryPtr = typename ARKodeStepper<VectorType>::ARKodeMemoryPtr;
 
    void
    reinit(double t0, const VectorType &y0) override;
 
    void
    resize(double t0, const VectorType &y0, double hscale) override;
 
#    if DEAL_II_SUNDIALS_VERSION_GTE(7, 5, 0)
    void
    attach_default_dom_eig_estimator(N_Vector y, SUNContext context);
#    endif
 
    AdditionalData data;
 
#    if DEAL_II_SUNDIALS_VERSION_GTE(7, 5, 0)
    std::unique_ptr<void, void (*)(void *)> dom_eig_estimator;
#    endif
 
    ARKodeMemoryPtr arkode_mem;
 
    CallbackContext callback_ctx;
  };
 
 
  template <typename VectorType>
  void
  LSRKStepperSTS<VectorType>::AdditionalData::add_parameters(
    ParameterHandler &prm)
  {
    prm.add_parameter(
      "Method name",
      method_name,
      "STS method name passed to LSRKStepSetSTSMethodByName. "
      "If empty, SUNDIALS uses the default (ARKODE_LSRK_RKC_2).");
    prm.add_parameter(
      "Dominant eigenvalue recomputation frequency",
      dom_eig_frequency,
      "Number of successful steps after which the dominant eigenvalue "
      "estimate is recomputed. 0 means constant dominant eigenvalue "
      "(no change throughout the simulation). Any negative value leaves "
      "this option at the SUNDIALS default of 25.");
    prm.add_parameter("Maximum number of stages",
                      max_num_stages,
                      "Maximum number of polynomial stages per time step. "
                      "0 means no limit.");
    prm.add_parameter(
      "Dominant eigenvalue estimator maximum iterations",
      dom_eig_estimator_max_iters,
      "Maximum number of power iterations performed by the built-in "
      "dominant-eigenvalue estimator. Only used when no dominant eigenvalue "
      "function is provided.");
    prm.add_parameter(
      "Dominant eigenvalue estimator relative tolerance",
      dom_eig_estimator_rel_tol,
      "Relative convergence tolerance of the power iteration used by the "
      "built-in dominant-eigenvalue estimator. Only used when no dominant "
      "eigenvalue function is provided.");
    prm.add_parameter(
      "Dominant eigenvalue estimator warmup iterations",
      dom_eig_estimator_num_warmups,
      "Number of preprocessing warmup iterations performed once when the "
      "built-in dominant-eigenvalue estimator is first initialized. Only used "
      "when no dominant eigenvalue function is provided.");
  }
 
 
  template <typename VectorType = Vector<double>>
  class LSRKStepperSSP : public ARKodeStepper<VectorType>
  {
  public:
    class AdditionalData
    {
    public:
      AdditionalData(const std::string &method_name = "",
                     const unsigned int num_stages  = 0)
        : method_name(method_name)
        , num_stages(num_stages)
      {}
 
      void
      add_parameters(ParameterHandler &prm);
 
      std::string method_name;
 
      unsigned int num_stages;
    };
 
    LSRKStepperSSP(const AdditionalData &data = AdditionalData());
 
    void *
    get_arkode_memory() const override;
 
    std::function<
      void(const double t, const VectorType &y, VectorType &explicit_f)>
      explicit_function;
 
    std::function<void(void *arkode_mem)> custom_setup;
 
  private:
    using CallbackContext = typename ARKodeStepper<
      VectorType>::template CallbackContext<LSRKStepperSSP<VectorType>>;
 
    using ARKodeMemoryPtr = typename ARKodeStepper<VectorType>::ARKodeMemoryPtr;
 
    void
    reinit(double t0, const VectorType &y0) override;
 
    void
    resize(double t0, const VectorType &y0, double hscale) override;
 
    AdditionalData data;
 
    ARKodeMemoryPtr arkode_mem;
 
    CallbackContext callback_ctx;
  };
 
 
  template <typename VectorType>
  void
  LSRKStepperSSP<VectorType>::AdditionalData::add_parameters(
    ParameterHandler &prm)
  {
    prm.add_parameter(
      "Method name",
      method_name,
      "SSP method name passed to LSRKStepSetSSPMethodByName. "
      "If empty, SUNDIALS uses the default (ARKODE_LSRK_SSP_10_4).");
    prm.add_parameter("Number of SSP stages",
                      num_stages,
                      "Number of stages for variable-stage SSP methods "
                      "(ARKODE_LSRK_SSP_S_2, ARKODE_LSRK_SSP_S_3). "
                      "0 means SUNDIALS default.");
  }
 
#  endif // DEAL_II_SUNDIALS_VERSION_GTE(7, 2, 0)
 
} // namespace SUNDIALS
 
#endif // DEAL_II_WITH_SUNDIALS
 
DEAL_II_NAMESPACE_CLOSE
 
#endif