arkode.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) 2017 - 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_h
#define dealii_sundials_arkode_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/arkode_stepper.h>
#  include <deal.II/sundials/n_vector.h>
#  include <deal.II/sundials/sundials_types.h>
#  include <deal.II/sundials/sunlinsol_wrapper.h>
 
#  include <boost/signals2.hpp>
 
#  include <sundials/sundials_linearsolver.h>
#  include <sundials/sundials_math.h>
 
#  include <exception>
#  include <memory>
 
#endif
 
DEAL_II_NAMESPACE_OPEN
 
#ifdef DEAL_II_WITH_SUNDIALS
namespace SUNDIALS
{
  template <typename VectorType = Vector<double>>
  class ARKode
  {
  private:
    template <class Fn>
    struct FunctionProxy
    {
      std::function<Fn> *target;
 
      FunctionProxy(std::function<Fn> *t = nullptr)
        : target(t)
      {}
 
      FunctionProxy &
      operator=(std::function<Fn> *t)
      {
        target = t;
        return *this;
      }
 
      FunctionProxy &
      operator=(std::function<Fn> f)
      {
        AssertThrow(
          target,
          ExcMessage(
            "Unable to assign the function since the target is not set in the "
            "proxy. Most probably, you explicitly provided a stepper to create "
            "the ARKode object and then tried to bind a callback to the latter. "
            "For such use cases, assign the user callbacks directly to the "
            "stepper. Assignment of callbacks to the ARKode object is permitted "
            "only if one of the constructors without a stepper parameter was "
            "used for instantiation of the ARKode object."));
 
        *target = std::move(f);
        return *this;
      }
    };
 
  public:
    class AdditionalData
    {
    public:
      DEAL_II_DEPRECATED_WITH_COMMENT(
        "Use another constructor and ARKStepper::AdditionalData to initialize "
        "with the ARKStepper relevant parameters instead.")
      AdditionalData(
        // Initial parameters
        const double initial_time,
        const double final_time,
        const double initial_step_size,
        const double output_period,
        // Running parameters
        const double       minimum_step_size,
        const unsigned int maximum_order                         = 5,
        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,
        // Error parameters
        const double absolute_tolerance = 1e-6,
        const double relative_tolerance = 1e-5);
 
      AdditionalData(
        // Initial parameters
        const double initial_time      = 0.0,
        const double final_time        = 1.0,
        const double initial_step_size = 1e-2,
        const double output_period     = 1e-1,
        // Running parameters
        const double minimum_step_size = 1e-6,
        // Error parameters
        const double absolute_tolerance = 1e-6,
        const double relative_tolerance = 1e-5);
 
      void
      add_parameters(ParameterHandler &prm);
 
      double initial_time;
 
      double final_time;
 
      double initial_step_size;
 
      double minimum_step_size;
 
      double absolute_tolerance;
 
      double relative_tolerance;
 
      DEAL_II_DEPRECATED_WITH_COMMENT(
        "Provide ARKStepper::AdditionalData::order instead.")
      unsigned int maximum_order;
 
      double output_period;
 
      DEAL_II_DEPRECATED_WITH_COMMENT(
        "Provide ARKStepper::AdditionalData::maximum_non_linear_iterations "
        "instead.")
      unsigned int maximum_non_linear_iterations;
 
      DEAL_II_DEPRECATED_WITH_COMMENT(
        "Provide ARKStepper::AdditionalData::implicit_function_is_linear "
        "instead.")
      bool implicit_function_is_linear;
 
      DEAL_II_DEPRECATED_WITH_COMMENT(
        "Provide ARKStepper::AdditionalData::implicit_function_is_time_independent "
        "instead.")
      bool implicit_function_is_time_independent;
 
      DEAL_II_DEPRECATED_WITH_COMMENT(
        "Provide ARKStepper::AdditionalData::mass_is_time_independent "
        "instead.")
      bool mass_is_time_independent;
 
      DEAL_II_DEPRECATED_WITH_COMMENT(
        "Provide ARKStepper::AdditionalData::anderson_acceleration_subspace "
        "instead.")
      int anderson_acceleration_subspace;
    };
 
    ARKode(const AdditionalData &data = AdditionalData());
 
    ARKode(const AdditionalData &data, const MPI_Comm mpi_comm);
 
    ARKode(ARKodeStepper<VectorType> &stepper,
           const AdditionalData      &data = AdditionalData());
 
    ARKode(ARKodeStepper<VectorType> &stepper,
           const AdditionalData      &data,
           const MPI_Comm             mpi_comm);
 
    ~ARKode();
 
    unsigned int
    solve_ode(VectorType &solution);
 
    unsigned int
    solve_ode_incrementally(VectorType  &solution,
                            const double intermediate_time,
                            const bool   reset_solver = false);
 
    void
    reset(const double t, const double h, const VectorType &y);
 
    void *
    get_arkode_memory() const;
 
    FunctionProxy<
      void(const double t, const VectorType &y, VectorType &explicit_f)>
      explicit_function;
 
    FunctionProxy<void(const double t, const VectorType &y, VectorType &res)>
      implicit_function;
 
    FunctionProxy<void(const double t, const VectorType &v, VectorType &Mv)>
      mass_times_vector;
 
    FunctionProxy<void(const double t)> mass_times_setup;
 
    FunctionProxy<void(const VectorType &v,
                       VectorType       &Jv,
                       const double      t,
                       const VectorType &y,
                       const VectorType &fy)>
      jacobian_times_vector;
 
    FunctionProxy<
      void(const double t, const VectorType &y, const VectorType &fy)>
      jacobian_times_setup;
 
    FunctionProxy<void(SundialsOperator<VectorType>       &op,
                       SundialsPreconditioner<VectorType> &prec,
                       VectorType                         &x,
                       const VectorType                   &b,
                       double                              tol)>
      solve_linearized_system;
 
    FunctionProxy<void(SundialsOperator<VectorType>       &op,
                       SundialsPreconditioner<VectorType> &prec,
                       VectorType                         &x,
                       const VectorType                   &b,
                       double                              tol)>
      solve_mass;
 
    FunctionProxy<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;
 
    FunctionProxy<void(const double      t,
                       const VectorType &y,
                       const VectorType &fy,
                       const int         jok,
                       int              &jcur,
                       const double      gamma)>
      jacobian_preconditioner_setup;
 
    FunctionProxy<void(const double      t,
                       const VectorType &r,
                       VectorType       &z,
                       const double      tol,
                       const int         lr)>
      mass_preconditioner_solve;
 
    FunctionProxy<void(const double t)> mass_preconditioner_setup;
 
    std::function<void(const double       t,
                       const VectorType  &sol,
                       const unsigned int step_number)>
      output_step;
 
    std::function<bool(const double t, VectorType &sol)> solver_should_restart;
 
    std::function<VectorType &()> get_local_tolerances;
 
    std::function<void(void *arkode_mem)> custom_setup;
 
  private:
    DeclException1(ExcFunctionNotProvided,
                   std::string,
                   << "Please provide an implementation for the function \""
                   << arg1 << "\"");
 
    unsigned int
    do_evolve_time(VectorType           &solution,
                   ::DiscreteTime &time,
                   const bool            do_reset);
 
    void
    set_functions_to_trigger_an_assert();
 
    void
    initialize_context();
 
    AdditionalData data;
 
    std::unique_ptr<ARKStepper<VectorType>> ark_stepper_storage{nullptr};
 
    ARKodeStepper<VectorType> &stepper;
 
#  if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
    SUNContext arkode_ctx;
#  endif
 
    MPI_Comm mpi_communicator;
 
    double last_end_time;
 
    mutable std::exception_ptr pending_exception;
 
#  ifdef DEAL_II_WITH_PETSC
#    ifdef PETSC_USE_COMPLEX
    static_assert(!std::is_same_v<VectorType, PETScWrappers::MPI::Vector>,
                  "Sundials does not support complex scalar types, "
                  "but PETSc is configured to use a complex scalar type!");
 
    static_assert(!std::is_same_v<VectorType, PETScWrappers::MPI::BlockVector>,
                  "Sundials does not support complex scalar types, "
                  "but PETSc is configured to use a complex scalar type!");
#    endif // PETSC_USE_COMPLEX
#  endif   // DEAL_II_WITH_PETSC
  };
 
 
 
  template <typename VectorType>
  ARKode<VectorType>::AdditionalData::AdditionalData(
    // Initial parameters
    const double initial_time,
    const double final_time,
    const double initial_step_size,
    const double output_period,
    // Running parameters
    const double       minimum_step_size,
    const unsigned int maximum_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,
    // Error parameters
    const double absolute_tolerance,
    const double relative_tolerance)
    : initial_time(initial_time)
    , final_time(final_time)
    , initial_step_size(initial_step_size)
    , minimum_step_size(minimum_step_size)
    , absolute_tolerance(absolute_tolerance)
    , relative_tolerance(relative_tolerance)
    , maximum_order(maximum_order)
    , output_period(output_period)
    , 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)
  {}
 
 
 
  template <typename VectorType>
  ARKode<VectorType>::AdditionalData::AdditionalData(
    // Initial parameters
    const double initial_time,
    const double final_time,
    const double initial_step_size,
    const double output_period,
    // Running parameters
    const double minimum_step_size,
    // Error parameters
    const double absolute_tolerance,
    const double relative_tolerance)
    : initial_time(initial_time)
    , final_time(final_time)
    , initial_step_size(initial_step_size)
    , minimum_step_size(minimum_step_size)
    , absolute_tolerance(absolute_tolerance)
    , relative_tolerance(relative_tolerance)
    , maximum_order(5)
    , output_period(output_period)
    , maximum_non_linear_iterations(10)
    , implicit_function_is_linear(false)
    , implicit_function_is_time_independent(false)
    , mass_is_time_independent(false)
    , anderson_acceleration_subspace(3)
  {}
 
 
 
  template <typename VectorType>
  void
  ARKode<VectorType>::AdditionalData::add_parameters(ParameterHandler &prm)
  {
    prm.add_parameter("Initial time", initial_time);
    prm.add_parameter("Final time", final_time);
    prm.add_parameter("Time interval between each output", output_period);
    prm.enter_subsection("Running parameters");
    prm.add_parameter("Initial step size", initial_step_size);
    prm.add_parameter("Minimum step size", minimum_step_size);
    prm.add_parameter("Maximum order of ARK", maximum_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.leave_subsection();
    prm.enter_subsection("Error control");
    prm.add_parameter("Absolute error tolerance", absolute_tolerance);
    prm.add_parameter("Relative error tolerance", relative_tolerance);
    prm.leave_subsection();
  }
 
} // namespace SUNDIALS
 
#endif
 
DEAL_II_NAMESPACE_CLOSE
 
 
#endif