1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
|
/*
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;
}
|
