some other fixes

Added a concept to make sure that the forcing term has the right
signature

Examples/src/RKFSolver/ButcherRKF.hpp

@@ -39,7 +39,7 @@ template <unsigned int NSTAGES> struct ButcherArray
     std::transform(A.begin(), A.end(), c.begin(), [](auto const &row) {
       return std::accumulate(row.begin(), row.end(), 0.0);
     });
-    set_implicit<NSTAGES>();
+    implicit_ = set_implicit<NSTAGES>();
   }
   /* I store the full array even if only the part below the main diagonal is
    * different from zero. For simplicity
@@ -82,13 +82,18 @@ template <unsigned int NSTAGES> struct ButcherArray
 protected:
   bool implicit_ = false;
   template <unsigned int N>
-  constexpr void
+  constexpr bool
   set_implicit()
   {
+    /*   if constexpr(N <= 1u)
+         implicit_ = implicit_ or (A[0][0] != 0);
+       else
+         implicit_ = implicit_ or (A[N - 1][N - 1] != 0);
+     */
     if constexpr(N <= 1u)
-      implicit_ = implicit_ or (A[0][0] != 0);
+      return (A[0][0] != 0);
     else
-      implicit_ = implicit_ or (A[N - 1][N - 1] != 0);
+      return (A[N - 1][N - 1] != 0) or set_implicit<N - 1>();
   }
 };
 
@@ -156,13 +161,13 @@ namespace RKFScheme
     {}
   };
 
-  //! ESDIRK12 scheme
+  //! ESDIRK12 scheme (Lobatto III A)
   struct ESDIRK12_t : public ButcherArray<2>
   {
     constexpr ESDIRK12_t()
       : ButcherArray<2>{{{
                           {{0., 0.}},
-                          {{0., 1.}},
+                          {{0., 1.0}},
                         }},
                         {{0., 1.}},   // 1st order
                         {{0.5, 0.5}}, // 2nd order

Examples/src/RKFSolver/RKF.hpp

@@ -12,6 +12,7 @@
 #include "Newton.hpp"
 #include "RKFTraits.hpp"
 #include <cmath>
+#include <concepts>
 #include <functional>
 #include <iostream>
 #include <limits>
@@ -42,7 +43,8 @@ template <RKFKind KIND> struct RKFResult
 };
 
 /*!
- * A calss for explicit Runge-Kutta Fehlberg type solution of ODEs
+ * A class for explicit ir diagonally implicit Runge-Kutta Fehlberg type
+ * solution of ODEs
  * @tparam B The Butcher table of the scheme. Must be defined following the
  * scheme shown in ButcherRKF.hpp
  * @tparam KIND The type of traits to be used: SCALAR, VECTOR, MATRIX
@@ -54,11 +56,15 @@ class RKF : public RKFTraits<KIND>
   using VariableType = typename RKFTraits<KIND>::VariableType;
   using Function = typename RKFTraits<KIND>::ForcingTermType;
   //! Constructor just taking the function
-  template <class F = Function> RKF(F &&f) : M_f{std::forward<F>(f)} {};
+  template <class F = Function>
+    requires std::convertible_to<F, Function>
+  RKF(F &&f) : M_f{std::forward<F>(f)} {};
   // Constructor passing butcher table and forcing function
-  // Note: eliminated since butcher table is now a constant expression
-  // this should allow to sped up the code
-  template <class F = Function> RKF(B const &bt, F &&f) : RKF{f} {};
+  // Butcher table is now a constant expression
+  // this constructor is there only to activate template parameter deduction
+  template <class F = Function>
+    requires std::convertible_to<F, Function>
+  RKF(B const &bt, F &&f) : RKF{f} {};
 
   //! Default constructor
   RKF() = default;
@@ -99,7 +105,8 @@ class RKF : public RKFTraits<KIND>
 private:
   Function           M_f;
   static constexpr B ButcherTable = B{};
-  // The default options for the quasi newton solver
+  //! The default options for the quasi newton solver
+  //! @todo make it a member of the class. You should be able to change them!
   static constexpr apsc::NewtonOptions newtonOptions{
     1.e-10, // tolerance on step \f$||x_{new}-x_{old}||<tolerance\$
     1.e-10, // tolerance on residual
@@ -181,7 +188,7 @@ RKF<B, KIND>::operator()(const double &T0, const double &T,
   //
   // Now I need a factor to specify when I can enlarge the time step
   // to avoid reducing and expanding the time step repeatedly
-  // I need to take into account the order of the scheme is >2
+  // I need to take into account the order of the scheme may be >2
   double factor_contraction = 1. / (ButcherTable.order);
   double factor_expansion = 1. / (ButcherTable.order + 1.0);
 

Examples/src/RKFSolver/main_rk45.cpp

@@ -18,6 +18,7 @@ main()
   double errorDesired = 1.e-4;
   {
     apsc::RKF<RKFScheme::RK45_t, RKFKind::SCALAR> solver{fun};
+
     // RKF<RKFScheme::RK23_t> solver{fun};
     // RKF<RKFScheme::RK12_t> solver{fun};
     auto   solution = solver(t0, T, y0, h_init, errorDesired);
@@ -49,17 +50,17 @@ main()
       return out;
     };
     // apsc::RKF<RKFScheme::RK45_t, RKFKind::VECTOR> solver{fun};
-    //  apsc::RKF<RKFScheme::ESDIRK34_t, RKFKind::VECTOR> solver{fun};
-    //  std::cout << "ESDIRK34 is implicit? " << std::boolalpha
-    //            << apsc::RKFScheme::ESDIRK34_t{}.implicit() << std::endl;
+    // apsc::RKF<RKFScheme::ESDIRK34_t, RKFKind::VECTOR> solver{fun};
+    // std::cout << "ESDIRK12 is implicit? " << std::boolalpha
+    //          << apsc::RKFScheme::ESDIRK12_t{}.implicit() << std::endl;
     apsc::RKF<RKFScheme::ESDIRK12_t, RKFKind::VECTOR> solver{fun};
     t0 = 0;
     T = 100.;
     Eigen::VectorXd y0(2);
     y0[0] = 1.;
     y0[1] = 1.;
     int maxSteps = 2000;
-    errorDesired = 1.e-1;
+    errorDesired = 1.e-2;
     auto solution = solver(t0, T, y0, h_init, errorDesired, maxSteps);
     std::cout << " Desired max error " << errorDesired
               << " Failed:" << solution.failed << " Estimated Error "