/* 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;
}