From bd519120263530dc75a407bb6483855e52fe04a4 Mon Sep 17 00:00:00 2001 From: Mario Mulansky Date: Tue, 20 Dec 2011 17:07:01 +0100 Subject: [PATCH] 2d example --- .../numeric/odeint/stepper/bulirsch_stoer.hpp | 10 +- .../odeint/examples/2d_lattice/Jamfile | 14 ++ .../odeint/examples/2d_lattice/lattice2d.hpp | 155 ++++++++++++++++++ .../2d_lattice/nested_range_algebra.hpp | 36 ++++ .../odeint/examples/2d_lattice/spreading.cpp | 112 +++++++++++++ .../2d_lattice/vector_vector_resize.hpp | 95 +++++++++++ libs/numeric/odeint/examples/Jamfile | 1 + libs/numeric/odeint/examples/bs_bug.cpp | 39 +++++ 8 files changed, 458 insertions(+), 4 deletions(-) create mode 100644 libs/numeric/odeint/examples/2d_lattice/Jamfile create mode 100644 libs/numeric/odeint/examples/2d_lattice/lattice2d.hpp create mode 100644 libs/numeric/odeint/examples/2d_lattice/nested_range_algebra.hpp create mode 100644 libs/numeric/odeint/examples/2d_lattice/spreading.cpp create mode 100644 libs/numeric/odeint/examples/2d_lattice/vector_vector_resize.hpp create mode 100644 libs/numeric/odeint/examples/bs_bug.cpp diff --git a/boost/numeric/odeint/stepper/bulirsch_stoer.hpp b/boost/numeric/odeint/stepper/bulirsch_stoer.hpp index 8367bdd9..45f0d730 100644 --- a/boost/numeric/odeint/stepper/bulirsch_stoer.hpp +++ b/boost/numeric/odeint/stepper/bulirsch_stoer.hpp @@ -197,13 +197,13 @@ public: value_vector h_opt( m_k_max+1 ); value_vector work( m_k_max+1 ); - //std::cout << "t=" << t <<", dt=" << dt << "(" << m_dt_last << ")" << ", k_opt=" << m_current_k_opt << std::endl; + std::cout << "t=" << t <<", dt=" << dt << "(" << m_dt_last << ")" << ", k_opt=" << m_current_k_opt << std::endl; time_type new_h = dt; for( size_t k = 0 ; k <= m_current_k_opt+1 ; k++ ) { - //std::cout << "k=" << k <<": " << ", first: " << m_first << std::endl; + std::cout << " k=" << k; //<<": " << ", first: " << m_first << std::endl; m_midpoint.set_steps( m_interval_sequence[k] ); if( k == 0 ) { @@ -220,7 +220,7 @@ public: h_opt[k] = calc_h_opt( dt , error , k ); work[k] = m_cost[k]/h_opt[k]; //std::cout << '\t' << "h_opt=" << h_opt[k] << ", work=" << work[k] << std::endl; - //std::cout << '\t' << "error: " << error << std::endl; + std::cout << '\t' << "error: " << error << std::endl; if( (k == m_current_k_opt-1) || m_first ) { // convergence before k_opt ? @@ -285,7 +285,6 @@ public: new_h = h_opt[m_current_k_opt]; } else { - //std::cout << "REJECT!" << std::endl; reject = true; new_h = h_opt[m_current_k_opt]; } @@ -295,7 +294,10 @@ public: } if( !reject ) + { t += dt; + } else + std::cout << "REJECT!" << std::endl; if( !m_last_step_rejected || (new_h < dt) ) { diff --git a/libs/numeric/odeint/examples/2d_lattice/Jamfile b/libs/numeric/odeint/examples/2d_lattice/Jamfile new file mode 100644 index 00000000..81e56e67 --- /dev/null +++ b/libs/numeric/odeint/examples/2d_lattice/Jamfile @@ -0,0 +1,14 @@ +# Copyright 2011 Karsten Ahnert and Mario Mulansky. +# Distributed under the Boost Software License, Version 1.0. (See +# accompanying file LICENSE_1_0.txt or copy at +# http://www.boost.org/LICENSE_1_0.txt) + +project + : requirements + ../../../../.. + BOOST_ALL_NO_LIB=1 + : build-dir . + : default-build release + ; + +exe spreading : spreading.cpp ; \ No newline at end of file diff --git a/libs/numeric/odeint/examples/2d_lattice/lattice2d.hpp b/libs/numeric/odeint/examples/2d_lattice/lattice2d.hpp new file mode 100644 index 00000000..2c3f7e07 --- /dev/null +++ b/libs/numeric/odeint/examples/2d_lattice/lattice2d.hpp @@ -0,0 +1,155 @@ +/* stronlgy nonlinear hamiltonian lattice in 2d */ + +#ifndef LATTICE2D_HPP +#define LATTICE2D_HPP + +#include + +#include + +using boost::math::pow; + +template< int Kappa , int Lambda > +struct lattice2d { + + const double m_beta; + std::vector< std::vector< double > > m_omega; + + lattice2d( const double beta ) + : m_beta( beta ) + { } + + template< class StateIn , class StateOut > + void operator()( const StateIn &q , StateOut &dpdt ) + { + // q and dpdt are 2d + const int N = q.size(); + + int i; + for( i = 0 ; i < N ; ++i ) + { + const int i_l = (i-1+N) % N; + const int i_r = (i+1) % N; + for( int j = 0 ; j < N ; ++j ) + { + const int j_l = (j-1+N) % N; + const int j_r = (j+1) % N; + dpdt[i][j] = - m_omega[i][j] * pow( q[i][j] ) + - m_beta * pow( q[i][j] - q[i][j_l] ) + - m_beta * pow( q[i][j] - q[i][j_r] ) + - m_beta * pow( q[i][j] - q[i_l][j] ) + - m_beta * pow( q[i][j] - q[i_r][j] ); + } + } + } + + template< class StateIn > + double energy( const StateIn &q , const StateIn &p ) + { + // q and dpdt are 2d + const int N = q.size(); + double energy = 0.0; + int i; + for( i = 0 ; i < N ; ++i ) + { + const int i_l = (i-1+N) % N; + const int i_r = (i+1) % N; + for( int j = 0 ; j < N ; ++j ) + { + const int j_l = (j-1+N) % N; + const int j_r = (j+1) % N; + energy += p[i][j]*p[i][j] / 2.0 + + m_omega[i][j] * pow( q[i][j] ) / Kappa + + m_beta * pow( q[i][j] - q[i][j_l] ) / Lambda / 2 + + m_beta * pow( q[i][j] - q[i][j_r] ) / Lambda / 2 + + m_beta * pow( q[i][j] - q[i_l][j] ) / Lambda / 2 + + m_beta * pow( q[i][j] - q[i_r][j] ) / Lambda / 2; + } + } + return energy; + } + + + template< class StateIn , class StateOut > + double local_energy( const StateIn &q , const StateIn &p , StateOut &energy ) + { + // q and dpdt are 2d + const int N = q.size(); + double e = 0.0; + int i; + for( i = 0 ; i < N ; ++i ) + { + const int i_l = (i-1+N) % N; + const int i_r = (i+1) % N; + for( int j = 0 ; j < N ; ++j ) + { + const int j_l = (j-1+N) % N; + const int j_r = (j+1) % N; + energy[i][j] = p[i][j]*p[i][j] / 2.0 + + m_omega[i][j] * pow( q[i][j] ) / Kappa + + m_beta * pow( q[i][j] - q[i][j_l] ) / Lambda / 2 + + m_beta * pow( q[i][j] - q[i][j_r] ) / Lambda / 2 + + m_beta * pow( q[i][j] - q[i_l][j] ) / Lambda / 2 + + m_beta * pow( q[i][j] - q[i_r][j] ) / Lambda / 2; + e += energy[i][j]; + } + } + //rescale + e = 1.0/e; + for( i = 0 ; i < N ; ++i ) + for( int j = 0 ; j < N ; ++j ) + energy[i][j] *= e; + return 1.0/e; + } + + void load_pot( const char* filename , const double W , const double gap , + const size_t dim ) + { + std::ifstream in( filename , std::ios::in | std::ios::binary ); + if( !in.is_open() ) { + std::cerr << "pot file not found: " << filename << std::endl; + exit(0); + } else { + std::cout << "using pot file: " << filename << std::endl; + } + + m_omega.resize( dim ); + for( int i = 0 ; i < dim ; ++i ) + { + m_omega[i].resize( dim ); + for( size_t j = 0 ; j < dim ; ++j ) + { + if( !in.good() ) + { + std::cerr << "I/O Error: " << filename << std::endl; + exit(0); + } + double d; + in.read( (char*) &d , sizeof(d) ); + if( (d < 0) || (d > 1.0) ) + { + std::cerr << "ERROR: " << d << std::endl; + exit(0); + } + m_omega[i][j] = W*d + gap; + } + } + + } + + void generate_pot( const double W , const double gap , const size_t dim ) + { + m_omega.resize( dim ); + for( size_t i = 0 ; i < dim ; ++i ) + { + m_omega[i].resize( dim ); + for( size_t j = 0 ; j < dim ; ++j ) + { + m_omega[i][j] = W*static_cast(rand())/RAND_MAX + gap; + } + } + } + +}; + +#endif diff --git a/libs/numeric/odeint/examples/2d_lattice/nested_range_algebra.hpp b/libs/numeric/odeint/examples/2d_lattice/nested_range_algebra.hpp new file mode 100644 index 00000000..c094c796 --- /dev/null +++ b/libs/numeric/odeint/examples/2d_lattice/nested_range_algebra.hpp @@ -0,0 +1,36 @@ +/* nested range algebra */ + +#ifndef NESTED_RANGE_ALGEBRA +#define NESTED_RANGE_ALGEBRA + +namespace detail { + + template< class Iterator1 , class Iterator2 , class Iterator3 , class Operation , class Algebra > + void for_each3( Iterator1 first1 , Iterator1 last1 , Iterator2 first2 , Iterator3 first3, Operation op , Algebra &algebra ) +{ + for( ; first1 != last1 ; ) + algebra.for_each3( *first1++ , *first2++ , *first3++ , op ); +} +} + + +template< class InnerAlgebra > +struct nested_range_algebra +{ + + nested_range_algebra() + : m_inner_algebra() + { } + + template< class S1 , class S2 , class S3 , class Op > + void for_each3( S1 &s1 , S2 &s2 , S3 &s3 , Op op ) + { + detail::for_each3( boost::begin( s1 ) , boost::end( s1 ) , boost::begin( s2 ) , boost::begin( s3 ) , op , m_inner_algebra ); + } + + +private: + InnerAlgebra m_inner_algebra; +}; + +#endif diff --git a/libs/numeric/odeint/examples/2d_lattice/spreading.cpp b/libs/numeric/odeint/examples/2d_lattice/spreading.cpp new file mode 100644 index 00000000..465952c9 --- /dev/null +++ b/libs/numeric/odeint/examples/2d_lattice/spreading.cpp @@ -0,0 +1,112 @@ +/* + * Example of a 2D simulation of nonlinearly coupled oscillators. + * Program just prints final energy which should be close to the initial energy (1.0). + * No parallelization is employed here. + * Run time on a 2.3GHz Intel Core-i5: about 10 seconds for 100 steps. + * Compile simply via bjam or directly: + * g++ -O3 -I${BOOST_ROOT} -I../../../../.. spreading.cpp + */ + + +#include +#include +#include +#include +#include + +#include +#include + +// we use a vector< vector< double > > as state type, +// for that some functionality has to be added for odeint to work +#include "nested_range_algebra.hpp" +#include "vector_vector_resize.hpp" + +// defines the rhs of our dynamical equation +#include "lattice2d.hpp" +/* dynamical equations (Hamiltonian structure): +dqdt_{i,j} = p_{i,j} +dpdt_{i,j} = - omega_{i,j}*q_{i,j} - \beta*[ (q_{i,j} - q_{i,j-1})^3 + +(q_{i,j} - q_{i,j+1})^3 + +(q_{i,j} - q_{i-1,j})^3 + +(q_{i,j} - q_{i+1,j})^3 ] +*/ + + +using namespace std; + +static const int MAX_N = 1024;//2048; + +static const size_t KAPPA = 2; +static const size_t LAMBDA = 4; +static const double W = 1.0; +static const double gap = 0.0; +static const size_t steps = 100; +static const double dt = 0.1; + +double initial_e = 1.0; +double beta = 1.0; +int realization_index = 0; + +//the state type +typedef vector< vector< double > > state_type; + +//the stepper, choose a symplectic one to account for hamiltonian structure +//use nested_range_algebra for calculations on vector< vector< ... > > +typedef boost::numeric::odeint::symplectic_rkn_sb3a_mclachlan< + state_type , state_type , double , state_type , state_type , double , + nested_range_algebra< boost::numeric::odeint::range_algebra > , + boost::numeric::odeint::default_operations > stepper_type; + +double time_diff_in_ms( timeval &t1 , timeval &t2 ) +{ return (t2.tv_sec - t1.tv_sec)*1000.0 + (t2.tv_usec - t1.tv_usec)/1000.0 + 0.5; } + + +int main( int argc, const char* argv[] ) { + + srand( time(NULL) ); + + lattice2d< KAPPA , LAMBDA > lattice( beta ); + + + lattice.generate_pot( W , gap , MAX_N ); + + state_type q( MAX_N , vector< double >( MAX_N , 0.0 ) ); + + state_type p( q ); + + state_type energy( q ); + + p[MAX_N/2][MAX_N/2] = sqrt( 0.5*initial_e ); + p[MAX_N/2+1][MAX_N/2] = sqrt( 0.5*initial_e ); + p[MAX_N/2][MAX_N/2+1] = sqrt( 0.5*initial_e ); + p[MAX_N/2+1][MAX_N/2+1] = sqrt( 0.5*initial_e ); + + cout.precision(10); + + lattice.local_energy( q , p , energy ); + double e=0.0; + for( size_t i=0 ; i + +#include + +namespace boost { namespace numeric { namespace odeint { + +template<> +struct is_resizeable< std::vector< std::vector< double > > > +{ + typedef boost::true_type type; + const static bool value = type::value; +}; + +template<> +struct same_size_impl< std::vector< std::vector< double > > , std::vector< std::vector< double > > > +{ + typedef std::vector< std::vector< double > > state_type; + + static bool same_size( const state_type &x1 , + const state_type &x2 ) + { + bool same = ( boost::size( x1 ) == boost::size( x2 ) ); + if( !same ) + return false; + typename state_type::const_iterator begin1 = boost::begin( x1 ); + typename state_type::const_iterator begin2 = boost::begin( x2 ); + while( begin1 != boost::end( x1 ) ) + same &= ( boost::size( *begin1++ ) == boost::size( *begin2++ ) ); + return same; + } +}; + +template<> +struct resize_impl< std::vector< std::vector< double > > , std::vector< std::vector< double > > > +{ + typedef std::vector< std::vector< double > > state_type; + + static void resize( state_type &x1 , const state_type &x2 ) + { + x1.resize( boost::size( x2 ) ); + typename state_type::iterator begin1 = boost::begin( x1 ); + typename state_type::const_iterator begin2 = boost::begin( x2 ); + while( begin1 != boost::end( x1 ) ) + (*begin1++).resize( boost::size( *begin2++ ) ); + } +}; + +template<> +struct state_wrapper< std::vector< std::vector< double > > , true > +{ + typedef std::vector< std::vector< double > > state_type; + typedef state_wrapper< state_type > state_wrapper_type; + typedef boost::true_type is_resizeable; + + state_type m_v; + + template< class State > + bool same_size( const State &x ) + { + bool same = ( boost::size( m_v ) == boost::size( x ) ); + if( !same ) + return false; + typename state_type::iterator begin1 = boost::begin( m_v ); + typename State::const_iterator begin2 = boost::begin( x ); + while( begin1 != boost::end( m_v ) ) + same &= ( boost::size( *begin1++ ) == boost::size( *begin2++ ) ); + return same; + } + + template< class State > + bool resize( const State &x ) + { + if( !same_size( x ) ) + { + m_v.resize( boost::size( x ) ); + typename state_type::iterator begin1 = boost::begin( m_v ); + typename State::const_iterator begin2 = boost::begin( x ); + while( begin1 != boost::end( m_v ) ) + (*begin1++).resize( boost::size( *begin2++ ) ); + + return true; + } else + return false; + } + +}; + +} } } + +#endif diff --git a/libs/numeric/odeint/examples/Jamfile b/libs/numeric/odeint/examples/Jamfile index 3c2c2180..d97e528e 100644 --- a/libs/numeric/odeint/examples/Jamfile +++ b/libs/numeric/odeint/examples/Jamfile @@ -28,6 +28,7 @@ exe list_lattice : list_lattice.cpp ; exe stepper_details : stepper_details.cpp ; exe my_vector : my_vector.cpp ; exe lorenz_point : lorenz_point.cpp ; +exe bs_bug : bs_bug.cpp ; # build-project mtl ; # build-project ublas ; diff --git a/libs/numeric/odeint/examples/bs_bug.cpp b/libs/numeric/odeint/examples/bs_bug.cpp new file mode 100644 index 00000000..939ea3eb --- /dev/null +++ b/libs/numeric/odeint/examples/bs_bug.cpp @@ -0,0 +1,39 @@ +#include +#include +#include + +using namespace std; +using namespace boost::numeric::odeint; + +typedef boost::array< double , 3 > state_type; + +typedef boost::numeric::odeint::bulirsch_stoer< state_type > bulirsch_stoer_type; + +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]; +} + +int main() { + + const double eps_abs = 10.0;//1E-4; + const double eps_rel = 1E-4; + const double t_end = 10.0; + + bulirsch_stoer_type bulirsch_stoer( eps_abs , eps_rel ); + state_type x = {{ 10.0 , 10.0 , 10.0 }}; + + cout.precision(17); + + size_t steps = integrate_adaptive( bulirsch_stoer , lorenz , x , 0.0 , t_end , 0.1 ); + + cout << steps << " steps" << endl; + + return 0; +}