/* Copyright 2011 Mario Mulansky Copyright 2012 Karsten Ahnert 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) */ /* * 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 <iostream> #include <fstream> #include <vector> #include <cstdlib> #include <sys/time.h> #include <boost/ref.hpp> #include <boost/numeric/odeint/stepper/symplectic_rkn_sb3a_mclachlan.hpp> // 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<energy.size() ; ++i ) for( size_t j=0 ; j<energy[i].size() ; ++j ) { e += energy[i][j]; } cout << "initial energy: " << lattice.energy( q , p ) << endl; timeval elapsed_time_start , elapsed_time_end; gettimeofday(&elapsed_time_start , NULL); stepper_type stepper; for( size_t step=0 ; step<=steps ; ++step ) { stepper.do_step( boost::ref( lattice ) , make_pair( boost::ref( q ) , boost::ref( p ) ) , 0.0 , 0.1 ); } gettimeofday(&elapsed_time_end , NULL); double elapsed_time = time_diff_in_ms( elapsed_time_start , elapsed_time_end ); cout << steps << " steps in " << elapsed_time/1000 << " s (energy: " << lattice.energy( q , p ) << ")" << endl; }