123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147 |
- /* Boost lorenz.cpp test file
- Copyright 2012 Karsten Ahnert
- Copyright 2012 Mario Mulansky
- This file tests the odeint library with the vexcl types
- 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)
- */
- #define BOOST_TEST_MODULE odeint_vexcl
- #include <vector>
- #include <iostream>
- #include <vexcl/vexcl.hpp>
- #include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
- #include <boost/numeric/odeint/stepper/runge_kutta_cash_karp54.hpp>
- #include <boost/numeric/odeint/stepper/runge_kutta_dopri5.hpp>
- #include <boost/numeric/odeint/stepper/controlled_runge_kutta.hpp>
- #include <boost/numeric/odeint/stepper/dense_output_runge_kutta.hpp>
- #include <boost/numeric/odeint/algebra/vector_space_algebra.hpp>
- #include <boost/numeric/odeint/integrate/integrate_const.hpp>
- #include <boost/numeric/odeint/external/vexcl/vexcl.hpp>
- #include <boost/test/unit_test.hpp>
- using namespace boost::unit_test;
- namespace odeint = boost::numeric::odeint;
- typedef vex::vector< double > vector_type;
- typedef vex::multivector< double, 3 > state_type;
- const double sigma = 10.0;
- const double b = 8.0 / 3.0;
- struct sys_func
- {
- const vector_type &R;
- sys_func( const vector_type &_R ) : R( _R ) { }
- void operator()( const state_type &x , state_type &dxdt , double t ) const
- {
- dxdt(0) = -sigma * ( x(0) - x(1) );
- dxdt(1) = R * x(0) - x(1) - x(0) * x(2);
- dxdt(2) = - b * x(2) + x(0) * x(1);
- }
- };
- const size_t n = 1024 ;
- const double dt = 0.01;
- const double t_max = 1.0;
- BOOST_AUTO_TEST_CASE( integrate_const_rk4 )
- {
- vex::Context ctx( vex::Filter::Type(CL_DEVICE_TYPE_GPU) );
- double Rmin = 0.1 , Rmax = 50.0 , dR = ( Rmax - Rmin ) / double( n - 1 );
- std::vector<double> x( n * 3 ) , r( n );
- for( size_t i=0 ; i<n ; ++i ) r[i] = Rmin + dR * double( i );
- state_type X( ctx.queue() , n );
- X(0) = 10.0;
- X(1) = 10.0;
- X(2) = 10.0;
- vector_type R( ctx.queue() , r );
- odeint::runge_kutta4<
- state_type , double , state_type , double ,
- odeint::vector_space_algebra , odeint::default_operations
- > stepper;
- odeint::integrate_const( stepper , sys_func( R ) , X , 0.0 , t_max , dt );
- std::vector< double > res( 3 * n );
- vex::copy( X(0) , res );
- std::cout << res[0] << std::endl;
- }
- BOOST_AUTO_TEST_CASE( integrate_const_controlled_rk54 )
- {
- vex::Context ctx( vex::Filter::Type(CL_DEVICE_TYPE_GPU) );
- double Rmin = 0.1 , Rmax = 50.0 , dR = ( Rmax - Rmin ) / double( n - 1 );
- std::vector<double> x( n * 3 ) , r( n );
- for( size_t i=0 ; i<n ; ++i ) r[i] = Rmin + dR * double( i );
- state_type X( ctx.queue() , n );
- X(0) = 10.0;
- X(1) = 10.0;
- X(2) = 10.0;
- vector_type R( ctx.queue() , r );
- typedef odeint::runge_kutta_cash_karp54<
- state_type , double , state_type , double ,
- odeint::vector_space_algebra , odeint::default_operations
- > stepper_type;
- typedef odeint::controlled_runge_kutta< stepper_type > controlled_stepper_type;
- odeint::integrate_const( controlled_stepper_type() , sys_func( R ) , X , 0.0 , t_max , dt );
- std::vector< double > res( 3 * n );
- vex::copy( X(0) , res );
- std::cout << res[0] << std::endl;
- }
- BOOST_AUTO_TEST_CASE( integrate_const_dense_output_dopri5 )
- {
- vex::Context ctx( vex::Filter::Type(CL_DEVICE_TYPE_GPU) );
- double Rmin = 0.1 , Rmax = 50.0 , dR = ( Rmax - Rmin ) / double( n - 1 );
- std::vector<double> x( n * 3 ) , r( n );
- for( size_t i=0 ; i<n ; ++i ) r[i] = Rmin + dR * double( i );
- state_type X( ctx.queue() , n );
- X(0) = 10.0;
- X(1) = 10.0;
- X(2) = 10.0;
- vector_type R( ctx.queue() , r );
- typedef odeint::runge_kutta_dopri5<
- state_type , double , state_type , double ,
- odeint::vector_space_algebra , odeint::default_operations
- > stepper_type;
- typedef odeint::controlled_runge_kutta< stepper_type > controlled_stepper_type;
- typedef odeint::dense_output_runge_kutta< controlled_stepper_type > dense_output_stepper_type;
- odeint::integrate_const( dense_output_stepper_type() , sys_func( R ) , X , 0.0 , t_max , dt );
- std::vector< double > res( 3 * n );
- vex::copy( X(0) , res );
- std::cout << res[0] << std::endl;
- }
|