1 /*
2 Copyright 2011 Mario Mulansky
3 Copyright 2012 Karsten Ahnert
4
5 Distributed under the Boost Software License, Version 1.0.
6 (See accompanying file LICENSE_1_0.txt or
7 copy at http://www.boost.org/LICENSE_1_0.txt)
8 */
9
10
11 /*
12 * Example of a 2D simulation of nonlinearly coupled oscillators.
13 * Program just prints final energy which should be close to the initial energy (1.0).
14 * No parallelization is employed here.
15 * Run time on a 2.3GHz Intel Core-i5: about 10 seconds for 100 steps.
16 * Compile simply via bjam or directly:
17 * g++ -O3 -I${BOOST_ROOT} -I../../../../.. spreading.cpp
18 */
19
20
21 #include <iostream>
22 #include <fstream>
23 #include <vector>
24 #include <cstdlib>
25 #include <sys/time.h>
26
27 #include <boost/ref.hpp>
28 #include <boost/numeric/odeint/stepper/symplectic_rkn_sb3a_mclachlan.hpp>
29
30 // we use a vector< vector< double > > as state type,
31 // for that some functionality has to be added for odeint to work
32 #include "nested_range_algebra.hpp"
33 #include "vector_vector_resize.hpp"
34
35 // defines the rhs of our dynamical equation
36 #include "lattice2d.hpp"
37 /* dynamical equations (Hamiltonian structure):
38 dqdt_{i,j} = p_{i,j}
39 dpdt_{i,j} = - omega_{i,j}*q_{i,j} - \beta*[ (q_{i,j} - q_{i,j-1})^3
40 +(q_{i,j} - q_{i,j+1})^3
41 +(q_{i,j} - q_{i-1,j})^3
42 +(q_{i,j} - q_{i+1,j})^3 ]
43 */
44
45
46 using namespace std;
47
48 static const int MAX_N = 1024;//2048;
49
50 static const size_t KAPPA = 2;
51 static const size_t LAMBDA = 4;
52 static const double W = 1.0;
53 static const double gap = 0.0;
54 static const size_t steps = 100;
55 static const double dt = 0.1;
56
57 double initial_e = 1.0;
58 double beta = 1.0;
59 int realization_index = 0;
60
61 //the state type
62 typedef vector< vector< double > > state_type;
63
64 //the stepper, choose a symplectic one to account for hamiltonian structure
65 //use nested_range_algebra for calculations on vector< vector< ... > >
66 typedef boost::numeric::odeint::symplectic_rkn_sb3a_mclachlan<
67 state_type , state_type , double , state_type , state_type , double ,
68 nested_range_algebra< boost::numeric::odeint::range_algebra > ,
69 boost::numeric::odeint::default_operations > stepper_type;
70
time_diff_in_ms(timeval & t1,timeval & t2)71 double time_diff_in_ms( timeval &t1 , timeval &t2 )
72 { return (t2.tv_sec - t1.tv_sec)*1000.0 + (t2.tv_usec - t1.tv_usec)/1000.0 + 0.5; }
73
74
main(int argc,const char * argv[])75 int main( int argc, const char* argv[] ) {
76
77 srand( time(NULL) );
78
79 lattice2d< KAPPA , LAMBDA > lattice( beta );
80
81
82 lattice.generate_pot( W , gap , MAX_N );
83
84 state_type q( MAX_N , vector< double >( MAX_N , 0.0 ) );
85
86 state_type p( q );
87
88 state_type energy( q );
89
90 p[MAX_N/2][MAX_N/2] = sqrt( 0.5*initial_e );
91 p[MAX_N/2+1][MAX_N/2] = sqrt( 0.5*initial_e );
92 p[MAX_N/2][MAX_N/2+1] = sqrt( 0.5*initial_e );
93 p[MAX_N/2+1][MAX_N/2+1] = sqrt( 0.5*initial_e );
94
95 cout.precision(10);
96
97 lattice.local_energy( q , p , energy );
98 double e=0.0;
99 for( size_t i=0 ; i<energy.size() ; ++i )
100 for( size_t j=0 ; j<energy[i].size() ; ++j )
101 {
102 e += energy[i][j];
103 }
104
105 cout << "initial energy: " << lattice.energy( q , p ) << endl;
106
107 timeval elapsed_time_start , elapsed_time_end;
108 gettimeofday(&elapsed_time_start , NULL);
109
110 stepper_type stepper;
111
112 for( size_t step=0 ; step<=steps ; ++step )
113 {
114 stepper.do_step( boost::ref( lattice ) ,
115 make_pair( boost::ref( q ) , boost::ref( p ) ) ,
116 0.0 , 0.1 );
117 }
118
119 gettimeofday(&elapsed_time_end , NULL);
120 double elapsed_time = time_diff_in_ms( elapsed_time_start , elapsed_time_end );
121 cout << steps << " steps in " << elapsed_time/1000 << " s (energy: " << lattice.energy( q , p ) << ")" << endl;
122 }
123