perf: change return type of some function to keep expressions

ponio/examples/belousov_zhabotinsky_pirock.cpp

@@ -207,9 +207,7 @@ main( int argc, char** argv )
 
     auto pb = ponio::make_imex_operator_problem( fd, fr, fr_t );
 
-    auto pirock = ponio::runge_kutta::pirock::pirock<1>( ponio::runge_kutta::pirock::beta_0<double>(),
-        eigmax_computer,
-        ponio::shampine_trick::shampine_trick<decltype( u_ini )>() );
+    auto pirock = ponio::runge_kutta::pirock::pirock<1>( ponio::runge_kutta::pirock::beta_0<double>(), eigmax_computer );
 
     // time loop  -------------------------------------------------------------
     auto sol_range = ponio::make_solver_range( pb, pirock, u_ini, t_span, dt );

ponio/examples/heat_rock.cpp

@@ -115,7 +115,7 @@ main()
         } );
 
     auto pb_heat         = heat_model( dx, xmin, xmax );
-    auto eigmax_computer = [=]( heat_model&, double, std::valarray<double>&, double )
+    auto eigmax_computer = [=]( auto&, double, std::valarray<double>&, double )
     {
         return 4. / ( dx * dx );
     };

ponio/examples/heat_samurai.py

@@ -16,12 +16,12 @@
 name = "heat_samurai"
 data_dir = name+"_data"
 
-# make = subprocess.Popen(["make", name])
-# make.wait()
+make = subprocess.Popen(["make", name])
+make.wait()
 
-# args = [os.path.join(".", name)]
-# process = subprocess.Popen(args)
-# process.wait()
+args = [os.path.join(".", name)]
+process = subprocess.Popen(args)
+process.wait()
 
 
 def read_file(filename: str):

ponio/examples/lorenz.cpp

@@ -11,6 +11,13 @@
 #include <ponio/solver.hpp>
 #include <ponio/time_span.hpp>
 
+template <std::size_t I, typename array_coeff_t>
+auto
+f( std::valarray<double> const& x, std::valarray<double> const& y, double const& a, array_coeff_t const& coeff )
+{
+    return x + a * coeff * y;
+}
+
 int
 main( int, char** )
 {
@@ -24,7 +31,7 @@ main( int, char** )
     double const rho   = 28.;
     double const beta  = 8. / 3.;
 
-    auto lorenz = [=]( double, state_t const& u ) -> state_t
+    auto lorenz = [=]( double, auto const& u ) -> state_t
     {
         auto du1 = sigma * ( u[1] - u[0] );
         auto du2 = rho * u[0] - u[1] - u[0] * u[2];

ponio/include/ponio/detail.hpp

@@ -59,14 +59,14 @@ namespace detail
 
     /* tpl_inner_product */
     template <typename state_t, typename value_t, typename ArrayA_t, typename ArrayB_t, std::size_t... Is>
-    constexpr state_t
+    inline constexpr auto
     tpl_inner_product_impl( ArrayA_t const& a,
         ArrayB_t const& b,
         state_t const& init,
-        [[maybe_unused]] value_t mul_coeff,
+        [[maybe_unused]] value_t const& mul_coeff,
         std::index_sequence<Is...> )
     {
-        return ( init + ... + ( mul_coeff * a[Is] * b[Is] ) );
+        return ( init + ... + ( mul_coeff * ( a[Is] * b[Is] ) ) );
     }
 
     /**
@@ -86,8 +86,8 @@ namespace detail
      * @details This function compute \f$\texttt{init} + \sum_{i=0}^N \texttt{mul_coeff}a_ib_i\f$ without loop thanks to template.
      */
     template <std::size_t N, typename state_t, typename value_t, typename ArrayA_t, typename ArrayB_t>
-    constexpr state_t
-    tpl_inner_product( ArrayA_t const& a, ArrayB_t const& b, state_t const& init, value_t mul_coeff = value_t{ 1.0 } )
+    inline constexpr auto
+    tpl_inner_product( ArrayA_t const& a, ArrayB_t const& b, state_t const& init, value_t const& mul_coeff )
     {
         return tpl_inner_product_impl( a, b, init, mul_coeff, std::make_index_sequence<N>() );
     }

ponio/include/ponio/method.hpp

@@ -58,7 +58,7 @@ namespace ponio
          */
         method( Algorithm_t const& alg_, state_t const& shadow_of_u0 )
             : alg( alg_ )
-            , kis( ::detail::init_fill_array<std::tuple_size<step_storage_t>::value>( shadow_of_u0 ) )
+            , kis( ::detail::init_fill_array<std::tuple_size_v<step_storage_t>>( shadow_of_u0 ) )
         {
         }
 
@@ -207,7 +207,7 @@ namespace ponio
 
         method( Algorithm_t const& alg_, state_t const& shadow_of_u0 )
             : alg( alg_ )
-            , kis( ::detail::init_fill_array<std::tuple_size<step_storage_t>::value>( shadow_of_u0 ) )
+            , kis( ::detail::init_fill_array<std::tuple_size_v<step_storage_t>>( shadow_of_u0 ) )
         {
         }
 

ponio/include/ponio/problem.hpp

@@ -26,11 +26,11 @@ namespace ponio
         simple_problem( Callable_t&& f_ );
 
         template <typename state_t, typename value_t>
-        state_t
+        auto
         operator()( value_t t, state_t&& u );
 
         template <typename state_t, typename value_t>
-        state_t
+        auto
         operator()( value_t t, state_t& u );
     };
 
@@ -52,15 +52,15 @@ namespace ponio
      */
     template <typename Callable_t>
     template <typename state_t, typename value_t>
-    inline state_t
+    inline auto
     simple_problem<Callable_t>::operator()( value_t t, state_t&& u )
     {
         return f( t, std::forward<state_t>( u ) );
     }
 
     template <typename Callable_t>
     template <typename state_t, typename value_t>
-    inline state_t
+    inline auto
     simple_problem<Callable_t>::operator()( value_t t, state_t& u )
     {
         return f( t, u );
@@ -175,14 +175,14 @@ namespace ponio
         imex_problem( Callable_explicit_t&& f_explicit, Implicit_problem_t&& pb_implicit );
 
         template <typename state_t, typename value_t>
-        state_t
+        auto
         operator()( value_t t, state_t&& u )
         {
             return explicit_part( t, std::forward<state_t>( u ) ) + implicit_part( t, std::forward<state_t>( u ) );
         }
 
         template <typename state_t, typename value_t>
-        state_t
+        auto
         operator()( value_t t, state_t& u )
         {
             return explicit_part( t, u ) + implicit_part( t, u );
@@ -254,7 +254,7 @@ namespace ponio
         lawson_problem( Linear_t&& l_, Nonlinear_t&& n_ );
 
         template <typename state_t, typename value_t>
-        state_t
+        auto
         operator()( value_t t, state_t&& u );
     };
 
@@ -278,7 +278,7 @@ namespace ponio
      */
     template <typename Linear_t, typename Nonlinear_t>
     template <typename state_t, typename value_t>
-    state_t
+    auto
     lawson_problem<Linear_t, Nonlinear_t>::operator()( value_t t, state_t&& u )
     {
         return l * u + n( t, std::forward<state_t>( u ) );
@@ -311,35 +311,35 @@ namespace ponio
         problem( Callables_t... args );
 
         template <typename value_t, typename state_t, std::size_t... Is>
-        state_t
+        auto
         _sum_components_impl( value_t t, state_t&& u, std::index_sequence<Is...> );
 
         template <typename value_t, typename state_t, std::size_t... Is>
-        state_t
+        auto
         _sum_components_impl( value_t t, state_t& u, std::index_sequence<Is...> );
 
         template <typename value_t, typename state_t>
-        state_t
+        auto
         operator()( value_t t, state_t&& u );
 
         template <typename value_t, typename state_t>
-        state_t
+        auto
         operator()( value_t t, state_t& u );
 
         template <std::size_t I, typename value_t, typename state_t>
-        state_t
+        auto
         operator()( std::integral_constant<std::size_t, I>, value_t t, state_t&& u );
 
         template <std::size_t I, typename value_t, typename state_t>
-        state_t
+        auto
         operator()( std::integral_constant<std::size_t, I>, value_t t, state_t& u );
 
         template <std::size_t Index, typename value_t, typename state_t>
-        state_t
+        auto
         call( value_t t, state_t&& u );
 
         template <std::size_t Index, typename value_t, typename state_t>
-        state_t
+        auto
         call( value_t t, state_t& u );
     };
 
@@ -362,15 +362,15 @@ namespace ponio
      */
     template <typename... Callables_t>
     template <typename value_t, typename state_t, std::size_t... Is>
-    inline state_t
+    inline auto
     problem<Callables_t...>::_sum_components_impl( value_t t, state_t&& u, std::index_sequence<Is...> )
     {
         return ( call<Is>( t, std::forward<state_t>( u ) ) + ... );
     }
 
     template <typename... Callables_t>
     template <typename value_t, typename state_t, std::size_t... Is>
-    inline state_t
+    inline auto
     problem<Callables_t...>::_sum_components_impl( value_t t, state_t& u, std::index_sequence<Is...> )
     {
         return ( call<Is>( t, u ) + ... );
@@ -384,15 +384,15 @@ namespace ponio
      */
     template <typename... Callables_t>
     template <typename value_t, typename state_t>
-    inline state_t
+    inline auto
     problem<Callables_t...>::operator()( value_t t, state_t&& u )
     {
         return _sum_components_impl( t, std::forward<state_t>( u ), std::make_index_sequence<size>{} );
     }
 
     template <typename... Callables_t>
     template <typename value_t, typename state_t>
-    inline state_t
+    inline auto
     problem<Callables_t...>::operator()( value_t t, state_t& u )
     {
         return _sum_components_impl( t, u, std::make_index_sequence<size>{} );
@@ -408,15 +408,15 @@ namespace ponio
      */
     template <typename... Callables_t>
     template <std::size_t I, typename value_t, typename state_t>
-    inline state_t
+    inline auto
     problem<Callables_t...>::operator()( std::integral_constant<std::size_t, I>, value_t t, state_t&& u )
     {
         return call<I>( t, std::forward<state_t>( u ) );
     }
 
     template <typename... Callables_t>
     template <std::size_t I, typename value_t, typename state_t>
-    inline state_t
+    inline auto
     problem<Callables_t...>::operator()( std::integral_constant<std::size_t, I>, value_t t, state_t& u )
     {
         return call<I>( t, u );
@@ -431,15 +431,15 @@ namespace ponio
      */
     template <typename... Callables_t>
     template <std::size_t Index, typename value_t, typename state_t>
-    inline state_t
+    inline auto
     problem<Callables_t...>::call( value_t t, state_t&& u )
     {
         return std::get<Index>( system )( t, std::forward<state_t>( u ) );
     }
 
     template <typename... Callables_t>
     template <std::size_t Index, typename value_t, typename state_t>
-    inline state_t
+    inline auto
     problem<Callables_t...>::call( value_t t, state_t& u )
     {
         return std::get<Index>( system )( t, u );

ponio/include/ponio/solver.hpp

@@ -121,7 +121,7 @@
             , pb( rhs.pb )
             , t_span( rhs.t_span )
             , it_next_time( std::begin( t_span ) + std::ranges::distance( std::begin( rhs.t_span ), rhs.it_next_time ) )
             , dt_reference( rhs.dt_reference )
         {
         }
 
@@ -205,7 +205,7 @@
         difference_type
         next_time() const
         {
             return sol.time + sol.time_step;
         }
 
         /**
@@ -527,7 +527,7 @@
 
         value_t last_time = t_span.back();
 
-        auto meth = make_method<value_t>( std::forward<Algorithm_t>( algo ), un );
+        auto meth = make_method<value_t>( std::forward<Algorithm_t>( algo ), u0 );
 
         obs( current_time, un, dt );
 

ponio/test/compute_order.hpp

@@ -141,20 +141,10 @@ namespace explicit_method
         T gamma = 1.;
         T delta = 1.;
 
-        // make a problem that can be use also for splitting (in 3 parts) methods
-        auto pb = ponio::make_problem(
-            [=]( T, state_t&& u ) -> state_t
-            {
-                return { alpha * u[0] - beta * u[0] * u[1], 0. };
-            },
-            [=]( T, state_t&& u ) -> state_t
-            {
-                return { 0., delta * u[0] * u[1] };
-            },
-            [=]( T, state_t&& u ) -> state_t
-            {
-                return { 0., -gamma * u[1] };
-            } );
+        auto pb = [=]( T, auto&& u ) -> state_t
+        {
+            return { alpha * u[0] - beta * u[0] * u[1], delta * u[0] * u[1] - gamma * u[1] };
+        };
 
         // invariant calculator
         auto V = [=]( state_t const& u ) -> T
@@ -338,3 +328,67 @@ namespace RDA_method
     }
 
 } // namespace RDA_method
+
+namespace splitting_method
+{
+
+    template <typename Algorithm_t, typename T = double>
+    auto
+    long_time_check_order( Algorithm_t& algo )
+    {
+        using state_t = std::valarray<T>;
+
+        T alpha = 2. / 3.;
+        T beta  = 4. / 3.;
+        T gamma = 1.;
+        T delta = 1.;
+
+        auto pb = ponio::make_problem(
+            [=]( T, auto&& u ) -> state_t
+            {
+                return { alpha * u[0] - beta * u[0] * u[1], 0. };
+            },
+            [=]( T, auto&& u ) -> state_t
+            {
+                return { 0., delta * u[0] * u[1] };
+            },
+            [=]( T, auto&& u ) -> state_t
+            {
+                return { 0., -gamma * u[1] };
+            } );
+
+        // invariant calculator
+        auto V = [=]( state_t const& u ) -> T
+        {
+            return delta * u[0] - gamma * std::log( u[0] ) + beta * u[1] - alpha * std::log( u[1] );
+        };
+
+        T x0          = 1.9;
+        state_t u_ini = { x0, x0 };
+        T V_ini       = V( u_ini );
+
+        ponio::time_span<T> t_span = { 0., 1000. };
+        std::vector<T> dts         = { 0.25, 0.125, 0.1, 0.075, 0.05 }; // find a way to adapt this range to the method
+        std::vector<T> relative_errors;
+        std::vector<T> log_dts;
+
+        for ( auto dt : dts )
+        {
+            state_t u_end = ::ponio::solve( pb, algo, u_ini, t_span, dt, []( T, state_t const&, T ) {} );
+            relative_errors.push_back( std::log10( relative_error( V_ini, V( u_end ) ) ) );
+            log_dts.push_back( std::log10( dt ) );
+        }
+
+        auto [a, b] = mayor_method( log_dts, relative_errors );
+
+        return a;
+    }
+
+    template <typename Algorithm_t, typename T = double>
+    T
+    check_order( Algorithm_t algo = Algorithm_t() )
+    {
+        return long_time_check_order( algo );
+    }
+
+} // splitting_method