1 /*
2 * resizing_lattice.cpp
3 *
4 * Demonstrates the usage of resizing of the state type during integration.
5 * Examplary system is a strongly nonlinear, disordered Hamiltonian lattice
6 * where the spreading of energy is investigated
7 *
8 * Copyright 2011-2012 Mario Mulansky
9 * Copyright 2012-2013 Karsten Ahnert
10 * Distributed under the Boost Software License, Version 1.0. (See
11 * accompanying file LICENSE_1_0.txt or copy at
12 * http://www.boost.org/LICENSE_1_0.txt)
13 *
14 */
15
16 #include <iostream>
17 #include <utility>
18
19 #include <boost/numeric/odeint.hpp>
20
21 #include <boost/ref.hpp>
22 #include <boost/random.hpp>
23
24 using namespace std;
25 using namespace boost::numeric::odeint;
26
27 //[ resizing_lattice_system_class
28 typedef vector< double > coord_type;
29 typedef pair< coord_type , coord_type > state_type;
30
31 struct compacton_lattice
32 {
33 const int m_max_N;
34 const double m_beta;
35 int m_pot_start_index;
36 vector< double > m_pot;
37
compacton_latticecompacton_lattice38 compacton_lattice( int max_N , double beta , int pot_start_index )
39 : m_max_N( max_N ) , m_beta( beta ) , m_pot_start_index( pot_start_index ) , m_pot( max_N )
40 {
41 srand( time( NULL ) );
42 // fill random potential with iid values from [0,1]
43 boost::mt19937 rng;
44 boost::uniform_real<> unif( 0.0 , 1.0 );
45 boost::variate_generator< boost::mt19937&, boost::uniform_real<> > gen( rng , unif );
46 generate( m_pot.begin() , m_pot.end() , gen );
47 }
48
operator ()compacton_lattice49 void operator()( const coord_type &q , coord_type &dpdt )
50 {
51 // calculate dpdt = -dH/dq of this hamiltonian system
52 // dp_i/dt = - V_i * q_i^3 - beta*(q_i - q_{i-1})^3 + beta*(q_{i+1} - q_i)^3
53 const int N = q.size();
54 double diff = q[0] - q[N-1];
55 for( int i=0 ; i<N ; ++i )
56 {
57 dpdt[i] = - m_pot[m_pot_start_index+i] * q[i]*q[i]*q[i] -
58 m_beta * diff*diff*diff;
59 diff = q[(i+1) % N] - q[i];
60 dpdt[i] += m_beta * diff*diff*diff;
61 }
62 }
63
energy_distributioncompacton_lattice64 void energy_distribution( const coord_type &q , const coord_type &p , coord_type &energies )
65 {
66 // computes the energy per lattice site normalized by total energy
67 const size_t N = q.size();
68 double en = 0.0;
69 for( size_t i=0 ; i<N ; i++ )
70 {
71 const double diff = q[(i+1) % N] - q[i];
72 energies[i] = p[i]*p[i]/2.0
73 + m_pot[m_pot_start_index+i]*q[i]*q[i]*q[i]*q[i]/4.0
74 + m_beta/4.0 * diff*diff*diff*diff;
75 en += energies[i];
76 }
77 en = 1.0/en;
78 for( size_t i=0 ; i<N ; i++ )
79 {
80 energies[i] *= en;
81 }
82 }
83
energycompacton_lattice84 double energy( const coord_type &q , const coord_type &p )
85 {
86 // calculates the total energy of the excitation
87 const size_t N = q.size();
88 double en = 0.0;
89 for( size_t i=0 ; i<N ; i++ )
90 {
91 const double diff = q[(i+1) % N] - q[i];
92 en += p[i]*p[i]/2.0
93 + m_pot[m_pot_start_index+i]*q[i]*q[i]*q[i]*q[i] / 4.0
94 + m_beta/4.0 * diff*diff*diff*diff;
95 }
96 return en;
97 }
98
change_pot_startcompacton_lattice99 void change_pot_start( const int delta )
100 {
101 m_pot_start_index += delta;
102 }
103 };
104 //]
105
106 //[ resizing_lattice_resize_function
do_resize(coord_type & q,coord_type & p,coord_type & distr,const int N)107 void do_resize( coord_type &q , coord_type &p , coord_type &distr , const int N )
108 {
109 q.resize( N );
110 p.resize( N );
111 distr.resize( N );
112 }
113 //]
114
115 const int max_N = 1024;
116 const double beta = 1.0;
117
main()118 int main()
119 {
120 //[ resizing_lattice_initialize
121 //start with 60 sites
122 const int N_start = 60;
123 coord_type q( N_start , 0.0 );
124 q.reserve( max_N );
125 coord_type p( N_start , 0.0 );
126 p.reserve( max_N );
127 // start with uniform momentum distribution over 20 sites
128 fill( p.begin()+20 , p.end()-20 , 1.0/sqrt(20.0) );
129
130 coord_type distr( N_start , 0.0 );
131 distr.reserve( max_N );
132
133 // create the system
134 compacton_lattice lattice( max_N , beta , (max_N-N_start)/2 );
135
136 //create the stepper, note that we use an always_resizer because state size might change during steps
137 typedef symplectic_rkn_sb3a_mclachlan< coord_type , coord_type , double , coord_type , coord_type , double ,
138 range_algebra , default_operations , always_resizer > hamiltonian_stepper;
139 hamiltonian_stepper stepper;
140 hamiltonian_stepper::state_type state = make_pair( q , p );
141 //]
142
143 //[ resizing_lattice_steps_loop
144 double t = 0.0;
145 const double dt = 0.1;
146 const int steps = 10000;
147 for( int step = 0 ; step < steps ; ++step )
148 {
149 stepper.do_step( boost::ref(lattice) , state , t , dt );
150 lattice.energy_distribution( state.first , state.second , distr );
151 if( distr[10] > 1E-150 )
152 {
153 do_resize( state.first , state.second , distr , state.first.size()+20 );
154 rotate( state.first.begin() , state.first.end()-20 , state.first.end() );
155 rotate( state.second.begin() , state.second.end()-20 , state.second.end() );
156 lattice.change_pot_start( -20 );
157 cout << t << ": resized left to " << distr.size() << ", energy = " << lattice.energy( state.first , state.second ) << endl;
158 }
159 if( distr[distr.size()-10] > 1E-150 )
160 {
161 do_resize( state.first , state.second , distr , state.first.size()+20 );
162 cout << t << ": resized right to " << distr.size() << ", energy = " << lattice.energy( state.first , state.second ) << endl;
163 }
164 t += dt;
165 }
166 //]
167
168 cout << "final lattice size: " << distr.size() << ", final energy: " << lattice.energy( state.first , state.second ) << endl;
169 }
170