ida.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_ida_h
#define dealii_sundials_ida_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>
 
#  ifdef DEAL_II_SUNDIALS_WITH_IDAS
#    include <idas/idas.h>
#  else
#    include <ida/ida.h>
#  endif
 
#  include <deal.II/sundials/sundials_types.h>
#  include <deal.II/sundials/sunlinsol_wrapper.h>
 
#  include <boost/signals2.hpp>
 
#  include <nvector/nvector_serial.h>
#  include <sundials/sundials_config.h>
#  include <sundials/sundials_math.h>
 
#  include <cmath>
#  include <limits>
#  include <memory>
 
#endif // DEAL_II_WITH_SUNDIALS
 
DEAL_II_NAMESPACE_OPEN
 
#ifdef DEAL_II_WITH_SUNDIALS
namespace SUNDIALS
{
  template <typename VectorType = Vector<double>>
  class IDA
  {
  public:
    class AdditionalData
    {
    public:
      enum InitialConditionCorrection
      {
        none = 0,
 
        use_y_diff = 1,
 
        use_y_dot = 2
      };
 
      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-15,
        const unsigned int maximum_order                 = 5,
        const unsigned int maximum_non_linear_iterations = 10,
        const double       ls_norm_factor                = 0,
        // Error parameters
        const double absolute_tolerance                = 1e-6,
        const double relative_tolerance                = 1e-5,
        const bool   ignore_algebraic_terms_for_errors = true,
        // Initial conditions parameters
        const InitialConditionCorrection &ic_type    = use_y_diff,
        const InitialConditionCorrection &reset_type = use_y_diff,
        const unsigned int                maximum_non_linear_iterations_ic = 5,
        // Additional running and nonlinear solver parameters
        const double       maximum_step_size                 = 0.0,
        const long int     maximum_number_of_steps           = 500,
        const unsigned int maximum_error_test_failures       = 10,
        const unsigned int maximum_convergence_failures      = 10,
        const double       nonlinear_convergence_coefficient = 0.33,
        // Additional initial-condition correction parameters
        const unsigned int maximum_number_of_steps_ic      = 5,
        const unsigned int maximum_number_of_jacobians_ic  = 4,
        const unsigned int maximum_number_of_backtracks_ic = 100,
        const bool         line_search_enabled_ic          = true,
        const double       step_tolerance_ic =
          std::pow(std::numeric_limits<double>::epsilon(), 2.0 / 3.0),
        const double nonlinear_convergence_coefficient_ic = 0.0033,
        // Additional linear solver parameters
        const double linear_tolerance_factor = 0.05,
        const double increment_factor        = 1.0,
        const bool   linear_solution_scaling = true,
        const double delta_cj_lsetup         = 0.25,
        // Step adaptivity parameters
        const double       eta_min                     = 0.5,
        const double       eta_max                     = 2.0,
        const double       eta_low                     = 0.9,
        const double       eta_min_err_fail            = 0.25,
        const double       eta_conv_fail               = 0.25,
        const double       eta_min_fixed               = 1.0,
        const double       eta_max_fixed               = 2.0,
        const unsigned int maximum_constraint_failures = 10)
        : initial_time(initial_time)
        , final_time(final_time)
        , initial_step_size(initial_step_size)
        , minimum_step_size(minimum_step_size)
        , maximum_step_size(maximum_step_size)
        , maximum_number_of_steps(maximum_number_of_steps)
        , maximum_error_test_failures(maximum_error_test_failures)
        , maximum_convergence_failures(maximum_convergence_failures)
        , nonlinear_convergence_coefficient(nonlinear_convergence_coefficient)
        , maximum_number_of_steps_ic(maximum_number_of_steps_ic)
        , maximum_number_of_jacobians_ic(maximum_number_of_jacobians_ic)
        , maximum_number_of_backtracks_ic(maximum_number_of_backtracks_ic)
        , line_search_enabled_ic(line_search_enabled_ic)
        , step_tolerance_ic(step_tolerance_ic)
        , nonlinear_convergence_coefficient_ic(
            nonlinear_convergence_coefficient_ic)
        , linear_tolerance_factor(linear_tolerance_factor)
        , increment_factor(increment_factor)
        , linear_solution_scaling(linear_solution_scaling)
        , delta_cj_lsetup(delta_cj_lsetup)
        , eta_min(eta_min)
        , eta_max(eta_max)
        , eta_low(eta_low)
        , eta_min_err_fail(eta_min_err_fail)
        , eta_conv_fail(eta_conv_fail)
        , eta_min_fixed(eta_min_fixed)
        , eta_max_fixed(eta_max_fixed)
        , maximum_constraint_failures(maximum_constraint_failures)
        , absolute_tolerance(absolute_tolerance)
        , relative_tolerance(relative_tolerance)
        , maximum_order(maximum_order)
        , output_period(output_period)
        , ignore_algebraic_terms_for_errors(ignore_algebraic_terms_for_errors)
        , ic_type(ic_type)
        , reset_type(reset_type)
        , maximum_non_linear_iterations_ic(maximum_non_linear_iterations_ic)
        , maximum_non_linear_iterations(maximum_non_linear_iterations)
        , ls_norm_factor(ls_norm_factor)
      {}
 
      void
      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,
                          "Smallest allowed integration step.");
        prm.add_parameter(
          "Maximum step size",
          maximum_step_size,
          "Largest allowed integration step; zero is unlimited.");
        prm.add_parameter("Maximum number of steps",
                          maximum_number_of_steps,
                          "Maximum internal steps between output times.");
        prm.add_parameter(
          "Maximum error test failures",
          maximum_error_test_failures,
          "Stops a solve after repeated local error test rejections.");
        prm.add_parameter(
          "Maximum nonlinear convergence failures",
          maximum_convergence_failures,
          "Stops a solve after repeated nonlinear solver failures.");
        prm.add_parameter("Maximum order of BDF", maximum_order);
        prm.add_parameter("Maximum number of nonlinear iterations",
                          maximum_non_linear_iterations);
        prm.leave_subsection();
 
        prm.enter_subsection("Nonlinear solver");
        prm.add_parameter(
          "Convergence coefficient",
          nonlinear_convergence_coefficient,
          "Controls when Newton corrections count as converged.");
        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.add_parameter(
          "Ignore algebraic terms for error computations",
          ignore_algebraic_terms_for_errors,
          "Indicate whether or not to suppress algebraic variables "
          "in the local error test.");
        prm.leave_subsection();
 
        prm.enter_subsection("Initial condition correction parameters");
        prm.add_parameter(
          "Correction type at initial time",
          ic_type,
          "This is one of the following three options for the "
          "initial condition calculation. \n"
          " none: do not try to make initial conditions consistent. \n"
          " use_y_diff: compute the algebraic components of y and differential\n"
          "    components of y_dot, given the differential components of y. \n"
          "    This option requires that the user specifies differential and \n"
          "    algebraic components in the function differential_components().\n"
          " use_y_dot: compute all components of y, given y_dot.");
        prm.add_parameter(
          "Correction type after restart",
          reset_type,
          "This is one of the following three options for the "
          "initial condition calculation. \n"
          " none: do not try to make initial conditions consistent. \n"
          " use_y_diff: compute the algebraic components of y and differential\n"
          "    components of y_dot, given the differential components of y. \n"
          "    This option requires that the user specifies differential and \n"
          "    algebraic components in the function differential_components().\n"
          " use_y_dot: compute all components of y, given y_dot.");
        prm.add_parameter("Maximum number of nonlinear iterations",
                          maximum_non_linear_iterations_ic);
        prm.add_parameter(
          "Factor to use when converting from the integrator tolerance to the linear solver tolerance",
          ls_norm_factor);
        prm.add_parameter(
          "Maximum number of steps",
          maximum_number_of_steps_ic,
          "Bounds work while making the initial values consistent.");
        prm.add_parameter(
          "Maximum number of Jacobians",
          maximum_number_of_jacobians_ic,
          "Bounds Jacobian refresh attempts during IC correction.");
        prm.add_parameter("Maximum number of backtracks",
                          maximum_number_of_backtracks_ic,
                          "Limits trial reductions in the IC line search.");
        prm.add_parameter("Line search enabled",
                          line_search_enabled_ic,
                          "Enable backtracking when adjusting initial values.");
        prm.add_parameter(
          "Step tolerance",
          step_tolerance_ic,
          "Treats a sufficiently small Newton correction as convergence.");
        prm.add_parameter(
          "Nonlinear convergence coefficient",
          nonlinear_convergence_coefficient_ic,
          "Sets the Newton convergence test used by IDACalcIC.");
        prm.leave_subsection();
 
        prm.enter_subsection("Linear solver");
        prm.add_parameter(
          "Linear tolerance factor",
          linear_tolerance_factor,
          "Relates linear solve accuracy to the nonlinear tolerance.");
        prm.add_parameter(
          "Difference quotient increment factor",
          increment_factor,
          "Scales finite-difference perturbations in Jacobian products.");
        prm.add_parameter(
          "Linear solution scaling",
          linear_solution_scaling,
          "Compensates for changes in Cj between linear solves.");
        prm.add_parameter(
          "Cj setup difference factor",
          delta_cj_lsetup,
          "Rebuilds linear setup when Cj changes by this fraction.");
        prm.leave_subsection();
 
        prm.enter_subsection("Step adaptivity");
        prm.add_parameter(
          "Minimum step growth factor",
          eta_min,
          "Prevents overly aggressive reductions in step size.");
        prm.add_parameter("Maximum step growth factor",
                          eta_max,
                          "Caps step growth after a successful step.");
        prm.add_parameter("Low step growth factor",
                          eta_low,
                          "Applies when IDA takes a lower-order step.");
        prm.add_parameter("Minimum growth factor after error failure",
                          eta_min_err_fail,
                          "Limits step reduction after a failed error test.");
        prm.add_parameter(
          "Growth factor after convergence failure",
          eta_conv_fail,
          "Sets the next-step reduction after nonlinear failure.");
        prm.add_parameter(
          "Minimum fixed step growth factor",
          eta_min_fixed,
          "Lower step-size ratio bound for retaining the current step size.");
        prm.add_parameter(
          "Maximum fixed step growth factor",
          eta_max_fixed,
          "Upper step-size ratio bound for retaining the current step size.");
        prm.leave_subsection();
        prm.enter_subsection("Constraints");
        prm.add_parameter(
          "Maximum constraint failures",
          maximum_constraint_failures,
          "Rejects a step after this many constraint violations.");
        prm.leave_subsection();
      }
 
      double initial_time;
 
      double final_time;
 
      double initial_step_size;
 
      double minimum_step_size;
 
      double maximum_step_size = 0.0;
 
      long int maximum_number_of_steps = 500;
 
      unsigned int maximum_error_test_failures = 10;
 
      unsigned int maximum_convergence_failures = 10;
 
      double nonlinear_convergence_coefficient = 0.33;
 
      unsigned int maximum_number_of_steps_ic = 5;
 
      unsigned int maximum_number_of_jacobians_ic = 4;
 
      unsigned int maximum_number_of_backtracks_ic = 100;
 
      bool line_search_enabled_ic = true;
 
      double step_tolerance_ic =
        std::pow(std::numeric_limits<double>::epsilon(), 2.0 / 3.0);
 
      double nonlinear_convergence_coefficient_ic = 0.0033;
 
      double linear_tolerance_factor = 0.05;
 
      double increment_factor = 1.0;
 
      bool linear_solution_scaling = true;
 
      double delta_cj_lsetup = 0.25;
 
      double eta_min = 0.5;
 
      double eta_max = 2.0;
 
      double eta_low = 0.9;
 
      double eta_min_err_fail = 0.25;
 
      double eta_conv_fail = 0.25;
 
      double eta_min_fixed = 1.0;
 
      double eta_max_fixed = 2.0;
 
      unsigned int maximum_constraint_failures = 10;
 
      double absolute_tolerance;
 
      double relative_tolerance;
 
      unsigned int maximum_order;
 
      double output_period;
 
      bool ignore_algebraic_terms_for_errors;
 
      InitialConditionCorrection ic_type;
 
      InitialConditionCorrection reset_type;
 
      unsigned maximum_non_linear_iterations_ic;
 
      unsigned int maximum_non_linear_iterations;
 
      double ls_norm_factor;
    };
 
    IDA(const AdditionalData &data = AdditionalData());
 
    IDA(const AdditionalData &data, const MPI_Comm mpi_comm);
 
    ~IDA();
 
    unsigned int
    solve_dae(VectorType &solution, VectorType &solution_dot);
 
    void
    reset(const double t, const double h, VectorType &y, VectorType &yp);
 
    std::function<void(VectorType &)> reinit_vector;
 
    std::function<void(const double      t,
                       const VectorType &y,
                       const VectorType &y_dot,
                       VectorType       &res)>
      residual;
 
    std::function<void(const double      t,
                       const VectorType &y,
                       const VectorType &y_dot,
                       const double      alpha)>
      setup_jacobian;
 
    std::function<
      void(const VectorType &rhs, VectorType &dst, const double tolerance)>
      solve_with_jacobian;
 
    std::function<void(const double       t,
                       const VectorType  &sol,
                       const VectorType  &sol_dot,
                       const unsigned int step_number)>
      output_step;
 
    std::function<bool(const double t, VectorType &sol, VectorType &sol_dot)>
      solver_should_restart;
 
    std::function<IndexSet()> differential_components;
 
    std::function<VectorType &()> get_local_tolerances;
 
    std::function<VectorType &()> get_constraint_vector;
 
    void *
    get_ida_memory() const;
 
    std::function<void(void *ida_mem)> custom_setup;
 
  private:
    DeclException1(ExcFunctionNotProvided,
                   std::string,
                   << "Please provide an implementation for the function \""
                   << arg1 << "\"");
 
    void
    set_functions_to_trigger_an_assert();
 
    const AdditionalData data;
 
    void *ida_mem;
 
#  if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
    SUNContext ida_ctx;
#  endif
 
    MPI_Comm mpi_communicator;
 
    GrowingVectorMemory<VectorType> mem;
 
    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(ExcIDAError,
                 int,
                 << "One of SUNDIALS IDA's internal functions "
                 << "returned an error code: " << arg1
                 << ". Please consult SUNDIALS manual.");
 
} // namespace SUNDIALS
 
#endif // DEAL_II_WITH_SUNDIALS
 
DEAL_II_NAMESPACE_CLOSE
 
#endif