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