[Boost-commit] svn:boost r66708 - sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper

Subject: [Boost-commit] svn:boost r66708 - sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper
From: mario.mulansky_at_[hidden]
Date: 2010-11-24 04:47:23

Author: mariomulansky
Date: 2010-11-24 04:47:19 EST (Wed, 24 Nov 2010)
New Revision: 66708
URL: http://svn.boost.org/trac/boost/changeset/66708

Log:
generic runge kutta stepper with fusion vectors and boost arrays
Added:
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/rk_test.cpp   (contents, props changed)
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/runge_kutta_stepper.hpp   (contents, props changed)
Text files modified: 
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/Jamfile            |     4 ++++                                    
   sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/fusion_stepper.hpp |    26 +++++++++++++++++++++-----              
   2 files changed, 25 insertions(+), 5 deletions(-)

Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/Jamfile
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/Jamfile	(original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/Jamfile	2010-11-24 04:47:19 EST (Wed, 24 Nov 2010)
@@ -18,3 +18,7 @@
 exe test 
         : test.cpp
         ;
+
+exe rk_test
+    : rk_test.cpp
+    ;
\ No newline at end of file

Modified: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/fusion_stepper.hpp
==============================================================================
--- sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/fusion_stepper.hpp	(original)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/fusion_stepper.hpp	2010-11-24 04:47:19 EST (Wed, 24 Nov 2010)
@@ -195,8 +195,23 @@
 {
 
         typedef State state_type;
-	typedef vector< vector< double > > parameter_array_type2D;
-	typedef vector< double > parameter_array_type1D;
+
+
+	typedef typename fusion::result_of::as_vector
+    <
+        typename mpl::copy
+        <
+            mpl::range_c< size_t , 1 , N > ,
+            mpl::inserter
+            <
+                mpl::vector0<> ,
+                mpl::push_back< mpl::_1 , array_wrapper< double , mpl::_2 > >
+            >
+        >::type
+    >::type coef_a_type;
+	typedef boost::array< double , N > coef_b_type;
+	typedef boost::array< double , N > coef_c_type;
+
 
     template< class System >
     struct calculate_stage
@@ -252,9 +267,9 @@
     >::type stage_vector;
     typedef stepper_last_stage< State , stage_count > last_stage_type;
 
-    runge_kutta_stepper( const parameter_array_type2D &a ,
-							const parameter_array_type1D &b ,
-							const parameter_array_type1D &c )
+    runge_kutta_stepper( const coef_a_type &a ,
+							const coef_b_type &b ,
+							const coef_c_type &c )
     {
             fusion::for_each( m_stages , init_stage( a , c ) );
             m_last_stage.init( b , c );
@@ -267,6 +282,7 @@
         m_last_stage( system , m_x_tmp , x , m_k_vector , t , dt );
     }
 
+
 private:
 
         state_type m_k_vector[stage_count];

Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/rk_test.cpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/rk_test.cpp	2010-11-24 04:47:19 EST (Wed, 24 Nov 2010)
@@ -0,0 +1,57 @@
+/*
+ * rk_test.cpp
+ *
+ *  Created on: Nov 23, 2010
+ *      Author: mario
+ */
+
+#include <iostream>
+#include <vector>
+#include <tr1/array>
+#include <boost/array.hpp>
+#include <boost/fusion/container.hpp>
+
+#include "runge_kutta_stepper.hpp"
+
+using namespace std;
+
+typedef tr1::array< double , 3 > state_type;
+typedef runge_kutta_stepper< state_type , 1 > euler_stepper;
+typedef runge_kutta_stepper< state_type , 2 > midpoint_stepper;
+
+
+const double sigma = 10.0;
+const double R = 28.0;
+const double b = 8.0 / 3.0;
+
+void lorenz( const state_type &x , state_type &dxdt , double t )
+{
+    dxdt[0] = sigma * ( x[1] - x[0] );
+    dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
+    dxdt[2] = x[0]*x[1] - b * x[2];
+}
+
+const double dt = 0.001;
+
+int main( void )
+{
+
+    boost::array< double , 1 > b = {{ 1.0 }};
+    boost::array< double , 1 > c = {{ 0.0 }};
+
+    euler_stepper euler( fusion::vector0<>() , b , c  );
+    euler.print_vals();
+
+    boost::array< double , 1 > a2 = {{ 0.5 }};
+    boost::array< double , 2 > b2 = {{ 0.0 , 0.5 }};
+    boost::array< double , 2 > c2 = {{ 0.0 , 1.0 }};
+
+    cout << typeid(midpoint_stepper::coef_a_type).name() << endl;
+    cout << typeid(midpoint_stepper::coef_b_type).name() << endl;
+    cout << typeid(midpoint_stepper::coef_c_type).name() << endl;
+    cout << typeid(midpoint_stepper::stage_vector_base).name() << endl;
+
+//    midpoint_stepper midpoint( fusion::make_vector( a2 ) , b2 , c2  );
+//    midpoint.print_vals();
+
+}

Added: sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/runge_kutta_stepper.hpp
==============================================================================
--- (empty file)
+++ sandbox/odeint/branches/karsten/libs/numeric/odeint/ideas/generic_stepper/runge_kutta_stepper.hpp	2010-11-24 04:47:19 EST (Wed, 24 Nov 2010)
@@ -0,0 +1,177 @@
+
+#include <boost/mpl/vector.hpp>
+#include <boost/mpl/size.hpp>
+#include <boost/mpl/push_back.hpp>
+#include <boost/mpl/for_each.hpp>
+#include <boost/mpl/zip_view.hpp>
+#include <boost/mpl/range_c.hpp>
+#include <boost/mpl/int.hpp>
+#include <boost/mpl/at.hpp>
+#include <boost/mpl/copy.hpp>
+#include <boost/mpl/insert_range.hpp>
+
+#include <vector>
+#include <algorithm>
+#include <iostream>
+
+#include <boost/fusion/container.hpp>
+#include <boost/fusion/algorithm.hpp>
+#include <boost/fusion/include/mpl.hpp>
+
+#include <boost/array.hpp>
+
+#include <typeinfo>
+
+
+#include "algebra.hpp"
+
+namespace mpl = boost::mpl;
+namespace fusion = boost::fusion;
+
+using namespace std;
+
+
+template< class T , class Constant >
+struct array_wrapper
+{
+    typedef typename boost::array< T , Constant::value > type;
+};
+
+
+template< class T , class Constant >
+struct stage_wrapper
+{
+    typedef typename  fusion::result_of::as_vector< mpl::vector< size_t , T , boost::array< T , Constant::value > > >::type type;
+};
+
+struct print_stage
+{
+    template< typename Stage >
+    void operator() ( const Stage &stage ) const
+    {
+        std::cout<< fusion::at_c<0>( stage ) << " " << fusion::at_c<1>( stage ) << " " << std::endl;
+    }
+};
+
+
+/* Runge Kutta Stepper - consisting of several stages */
+
+template< typename State , size_t stage_count >
+class runge_kutta_stepper
+{
+public:
+    typedef State state_type;
+
+    typedef mpl::range_c< size_t , 1 , stage_count > stage_indices;
+
+    typedef typename fusion::result_of::as_vector
+    <
+        typename mpl::copy
+        <
+            stage_indices ,
+            mpl::inserter
+            <
+                mpl::vector0<> ,
+                mpl::push_back< mpl::_1 , array_wrapper< double , mpl::_2 > >
+            >
+        >::type
+    >::type coef_a_type;
+
+    typedef boost::array< double , stage_count > coef_b_type;
+
+    typedef boost::array< double , stage_count > coef_c_type;
+
+
+    typedef typename fusion::result_of::as_vector<
+        typename mpl::copy<
+            stage_indices,
+            mpl::inserter
+              <
+                mpl::vector0<> ,
+                mpl::push_back< mpl::_1 , stage_wrapper< double , mpl::_2 > >
+              >
+            >::type
+        >::type stage_vector_base;
+
+
+    struct stage_vector_type : public stage_vector_base
+    {
+        struct init_stage{
+
+            const coef_a_type &m_a;
+            const coef_b_type &m_c;
+
+            init_stage( const coef_a_type &a , const coef_b_type &b , const coef_c_type &c )
+                : m_a( a ) , m_c( c ) {}
+
+
+            template< typename Stage >
+            void operator() ( const Stage &stage ) const
+            {
+                const size_t n = fusion::at_c<2>( stage ).size();
+                fusion::at_c<0>( stage ) = n;
+                fusion::at_c<1>( stage ) = m_c[n];
+                fusion::at_c<2>( stage ) = m_a[n];
+            }
+
+        };
+
+
+        stage_vector_type( const coef_a_type &a , const coef_b_type &b , const coef_c_type &c )
+        {
+            fusion::for_each( static_cast< stage_vector_base& >( *this ) , init_stage( a , b , c ) );
+        }
+    };
+
+
+//    template< class System >
+//    struct calculate_stage
+//    {
+//        System &system;
+//        state_type &x , &x_tmp;
+//        state_type *k_vector;
+//        const double t;
+//        const double dt;
+//
+//        calculate_stage( System &_system , state_type &_x , state_type &_x_tmp , state_type *_k_vector ,
+//                            const double _t , const double _dt )
+//        : system( _system ) , x( _x ) , x_tmp( _x_tmp ) , k_vector( _k_vector ) , t( _t ) , dt( _dt )
+//        {}
+//
+//        template< class Stage >
+//        void operator()( Stage const &stage ) const
+//        {
+//            stage( system , x_tmp , x , k_vector , t , dt );
+//        }
+//
+//    };
+
+
+public:
+
+    runge_kutta_stepper( const coef_a_type &a ,
+                            const coef_b_type &b ,
+                            const coef_c_type &c )
+        : m_stages( a , b , c )
+    { }
+
+    template< class System >
+    void do_step( System &system , state_type &x , double t , const double dt )
+    {
+//        fusion::for_each( m_stages , calculate_stage< System >( system , x , m_x_tmp , m_k_vector , t , dt ) );
+//        m_last_stage( system , m_x_tmp , x , m_k_vector , t , dt );
+    }
+
+    void print_vals()
+    {
+        fusion::for_each( m_stages , print_stage() );
+    }
+
+
+private:
+
+    state_type m_k_vector[stage_count];
+    state_type m_x_tmp;
+    const stage_vector_type m_stages;
+//    last_stage_type m_last_stage;
+};