Update Sundials and fix boost mpi include issue

CMakeLists.txt

@@ -140,7 +140,7 @@ else()
     sundials
     DOWNLOAD_EXTRACT_TIMESTAMP ON
     GIT_REPOSITORY https://github.com/LLNL/sundials
-    GIT_TAG        v6.1.1
+    GIT_TAG        v7.1.1
       # adjust this to the version you need
   )
   FetchContent_GetProperties(sundials)
@@ -150,7 +150,7 @@ else()
   endif()
 endif()
 
-
+set(BOOST_NUMERIC_ODEINT_NO_ADAPTORS ON)
 if (STAN_USE_SYSTEM_BOOST)
   find_package(Boost REQUIRED)
 else()
@@ -162,13 +162,13 @@ else()
   set(BOOST_ENABLE_PYTHON OFF)
   set(BOOST_BUILD_TESTS OFF)
   set(CMAKE_BUILD_TYPE Release)
-  set(BOOST_ENABLE_MPI ON)
+  set(BOOST_ENABLE_MPI ${STAN_MPI})
 
   set(BOOST_INCLUDE_LIBRARIES math numeric/odeint lexical_cast optional random mpi serialization)
   FetchContent_Declare(
     Boost
     DOWNLOAD_EXTRACT_TIMESTAMP ON
-    URL https://github.com/boostorg/boost/releases/download/boost-1.85.0/boost-1.85.0-cmake.tar.xz
+    URL https://github.com/boostorg/boost/releases/download/boost-1.86.0/boost-1.86.0-cmake.tar.xz
   )
   FetchContent_MakeAvailable(Boost)
 endif()

stan/math/rev/functor/algebra_solver_fp.hpp

@@ -11,7 +11,7 @@
 #include <stan/math/prim/fun/to_array_1d.hpp>
 #include <stan/math/prim/fun/to_vector.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
-#include <sundials/sundials_context.h>
+#include <sundials/sundials_context.hpp>
 #include <kinsol/kinsol.h>
 #include <sunmatrix/sunmatrix_dense.h>
 #include <sunlinsol/sunlinsol_dense.h>

stan/math/rev/functor/cvodes_integrator.hpp

@@ -8,7 +8,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/functor/apply.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
-#include <sundials/sundials_context.h>
+#include <sundials/sundials_context.hpp>
 #include <cvodes/cvodes.h>
 #include <nvector/nvector_serial.h>
 #include <sunlinsol/sunlinsol_dense.h>
@@ -66,7 +66,7 @@ class cvodes_integrator {
    * Implements the function of type CVRhsFn which is the user-defined
    * ODE RHS passed to CVODES.
    */
-  static int cv_rhs(realtype t, N_Vector y, N_Vector ydot, void* user_data) {
+  static int cv_rhs(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) {
     cvodes_integrator* integrator = static_cast<cvodes_integrator*>(user_data);
     integrator->rhs(t, NV_DATA_S(y), NV_DATA_S(ydot));
     return 0;
@@ -76,7 +76,7 @@ class cvodes_integrator {
    * Implements the function of type CVSensRhsFn which is the
    * RHS of the sensitivity ODE system.
    */
-  static int cv_rhs_sens(int Ns, realtype t, N_Vector y, N_Vector ydot,
+  static int cv_rhs_sens(int Ns, sunrealtype t, N_Vector y, N_Vector ydot,
                          N_Vector* yS, N_Vector* ySdot, void* user_data,
                          N_Vector tmp1, N_Vector tmp2) {
     cvodes_integrator* integrator = static_cast<cvodes_integrator*>(user_data);
@@ -90,7 +90,7 @@ class cvodes_integrator {
    * ode_rhs wrt to the states y. The jacobian is stored in column
    * major format.
    */
-  static int cv_jacobian_states(realtype t, N_Vector y, N_Vector fy,
+  static int cv_jacobian_states(sunrealtype t, N_Vector y, N_Vector fy,
                                 SUNMatrix J, void* user_data, N_Vector tmp1,
                                 N_Vector tmp2, N_Vector tmp3) {
     cvodes_integrator* integrator = static_cast<cvodes_integrator*>(user_data);

stan/math/rev/functor/cvodes_integrator_adjoint.hpp

@@ -9,7 +9,7 @@
 #include <stan/math/prim/fun/value_of.hpp>
 #include <stan/math/prim/functor/apply.hpp>
 #include <stan/math/prim/functor/for_each.hpp>
-#include <sundials/sundials_context.h>
+#include <sundials/sundials_context.hpp>
 #include <cvodes/cvodes.h>
 #include <nvector/nvector_serial.h>
 #include <sunmatrix/sunmatrix_dense.h>
@@ -576,7 +576,7 @@ class cvodes_integrator_adjoint_vari : public vari_base {
    * Implements the function of type CVRhsFn which is the user-defined
    * ODE RHS passed to CVODES.
    */
-  constexpr static int cv_rhs(realtype t, N_Vector y, N_Vector ydot,
+  constexpr static int cv_rhs(sunrealtype t, N_Vector y, N_Vector ydot,
                               void* user_data) {
     return cast_to_self(user_data)->rhs(t, NV_DATA_S(y), NV_DATA_S(ydot));
   }
@@ -611,7 +611,7 @@ class cvodes_integrator_adjoint_vari : public vari_base {
    * Implements the function of type CVRhsFnB which is the
    * RHS of the backward ODE system.
    */
-  constexpr static int cv_rhs_adj(realtype t, N_Vector y, N_Vector yB,
+  constexpr static int cv_rhs_adj(sunrealtype t, N_Vector y, N_Vector yB,
                                   N_Vector yBdot, void* user_data) {
     return cast_to_self(user_data)->rhs_adj(t, y, yB, yBdot);
   }
@@ -654,7 +654,7 @@ class cvodes_integrator_adjoint_vari : public vari_base {
    * Implements the function of type CVQuadRhsFnB which is the
    * RHS of the backward ODE system's quadrature.
    */
-  constexpr static int cv_quad_rhs_adj(realtype t, N_Vector y, N_Vector yB,
+  constexpr static int cv_quad_rhs_adj(sunrealtype t, N_Vector y, N_Vector yB,
                                        N_Vector qBdot, void* user_data) {
     return cast_to_self(user_data)->quad_rhs_adj(t, y, yB, qBdot);
   }
@@ -693,7 +693,7 @@ class cvodes_integrator_adjoint_vari : public vari_base {
    * ode_rhs wrt to the states y. The jacobian is stored in column
    * major format.
    */
-  constexpr static int cv_jacobian_rhs_states(realtype t, N_Vector y,
+  constexpr static int cv_jacobian_rhs_states(sunrealtype t, N_Vector y,
                                               N_Vector fy, SUNMatrix J,
                                               void* user_data, N_Vector tmp1,
                                               N_Vector tmp2, N_Vector tmp3) {
@@ -722,7 +722,7 @@ class cvodes_integrator_adjoint_vari : public vari_base {
    * Implements the CVLsJacFnB function for evaluating the jacobian of
    * the adjoint problem wrt to the backward states.
    */
-  constexpr static int cv_jacobian_rhs_adj_states(realtype t, N_Vector y,
+  constexpr static int cv_jacobian_rhs_adj_states(sunrealtype t, N_Vector y,
                                                   N_Vector yB, N_Vector fyB,
                                                   SUNMatrix J, void* user_data,
                                                   N_Vector tmp1, N_Vector tmp2,

stan/math/rev/functor/idas_service.hpp

@@ -9,7 +9,7 @@
 #include <nvector/nvector_serial.h>
 #include <sunmatrix/sunmatrix_dense.h>
 #include <sunlinsol/sunlinsol_dense.h>
-#include <sundials/sundials_context.h>
+#include <sundials/sundials_context.hpp>
 #include <ostream>
 #include <vector>
 #include <algorithm>
@@ -106,8 +106,8 @@ struct idas_service {
           "Failed to allocate yps N_Vectors for sensitivities.");
     }
     for (size_t is = 0; is < ns; ++is) {
-      N_VConst(RCONST(0.0), yys[is]);
-      N_VConst(RCONST(0.0), yps[is]);
+      N_VConst(SUN_RCONST(0.0), yys[is]);
+      N_VConst(SUN_RCONST(0.0), yps[is]);
     }
     set_init_sens(yys, yps, n);
     // std::cout << "Calling IDASensInit with ns=" << ns << " inner ns: " << ns_

stan/math/rev/functor/kinsol_data.hpp

@@ -6,7 +6,7 @@
 #include <stan/math/prim/fun/to_array_1d.hpp>
 #include <stan/math/prim/fun/to_vector.hpp>
 #include <stan/math/prim/functor/apply.hpp>
-#include <sundials/sundials_context.h>
+#include <sundials/sundials_context.hpp>
 #include <kinsol/kinsol.h>
 #include <sunmatrix/sunmatrix_dense.h>
 #include <sunlinsol/sunlinsol_dense.h>

stan/math/rev/functor/kinsol_solve.hpp

@@ -6,7 +6,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/to_array_1d.hpp>
 #include <stan/math/prim/fun/to_vector.hpp>
-#include <sundials/sundials_context.h>
+#include <sundials/sundials_context.hpp>
 #include <kinsol/kinsol.h>
 #include <sunmatrix/sunmatrix_dense.h>
 #include <sunlinsol/sunlinsol_dense.h>