arkode.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2017 - 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_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/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>
 
 
DEAL_II_NAMESPACE_OPEN
 
 
// Shorthand notation for ARKODE error codes.
#  define AssertARKode(code) Assert(code >= 0, ExcARKodeError(code))
 
namespace SUNDIALS
{
  template <typename VectorType = Vector<double>>
  class ARKode
  {
  public:
    class AdditionalData
    {
    public:
      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,
        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);
 
      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;
 
      unsigned int maximum_order;
 
      double output_period;
 
      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;
    };
 
    ARKode(const AdditionalData &data = AdditionalData());
 
    ARKode(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;
 
    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_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_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(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
    setup_system_solver(const VectorType &solution);
 
    void
    setup_mass_solver(const VectorType &solution);
 
    void
    set_functions_to_trigger_an_assert();
 
    AdditionalData data;
 
    void *arkode_mem;
 
#  if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
    SUNContext arkode_ctx;
#  endif
 
    MPI_Comm mpi_communicator;
 
    double last_end_time;
 
    std::unique_ptr<internal::LinearSolverWrapper<VectorType>> linear_solver;
    std::unique_ptr<internal::LinearSolverWrapper<VectorType>> mass_solver;
 
    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
  };
 
 
  DeclException1(ExcARKodeError,
                 int,
                 << "One of the SUNDIALS ARKode internal functions "
                 << " returned a negative error code: " << arg1
                 << ". Please consult SUNDIALS manual.");
 
 
  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>
  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
 
DEAL_II_NAMESPACE_CLOSE
#endif
 
 
#endif