...one of the most highly
regarded and expertly designed C++ library projects in the
world.
— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards
Home | Libraries | People | FAQ | More |
odeint supports changes of the state size during integration if a state_type
is used which can be resized, like std::vector
.
The adjustment of the state's size has to be done from outside and the stepper
has to be instantiated with always_resizer
as the template argument for the resizer_type
.
In this configuration, the stepper checks for changes in the state size and
adjust it's internal storage accordingly.
We show this for a Hamiltonian system of nonlinear, disordered oscillators with nonlinear nearest neighbor coupling.
The system function is implemented in terms of a class that also provides functions for calculating the energy. Note, that this class stores the random potential internally which is not resized, but rather a start index is kept which should be changed whenever the states' size change.
typedef vector< double > coord_type; typedef pair< coord_type , coord_type > state_type; struct compacton_lattice { const int m_max_N; const double m_beta; int m_pot_start_index; vector< double > m_pot; compacton_lattice( int max_N , double beta , int pot_start_index ) : m_max_N( max_N ) , m_beta( beta ) , m_pot_start_index( pot_start_index ) , m_pot( max_N ) { srand( time( NULL ) ); // fill random potential with iid values from [0,1] boost::mt19937 rng; boost::uniform_real<> unif( 0.0 , 1.0 ); boost::variate_generator< boost::mt19937&, boost::uniform_real<> > gen( rng , unif ); generate( m_pot.begin() , m_pot.end() , gen ); } void operator()( const coord_type &q , coord_type &dpdt ) { // calculate dpdt = -dH/dq of this hamiltonian system // dp_i/dt = - V_i * q_i^3 - beta*(q_i - q_{i-1})^3 + beta*(q_{i+1} - q_i)^3 const int N = q.size(); double diff = q[0] - q[N-1]; for( int i=0 ; i<N ; ++i ) { dpdt[i] = - m_pot[m_pot_start_index+i] * q[i]*q[i]*q[i] - m_beta * diff*diff*diff; diff = q[(i+1) % N] - q[i]; dpdt[i] += m_beta * diff*diff*diff; } } void energy_distribution( const coord_type &q , const coord_type &p , coord_type &energies ) { // computes the energy per lattice site normalized by total energy const size_t N = q.size(); double en = 0.0; for( size_t i=0 ; i<N ; i++ ) { const double diff = q[(i+1) % N] - q[i]; energies[i] = p[i]*p[i]/2.0 + m_pot[m_pot_start_index+i]*q[i]*q[i]*q[i]*q[i]/4.0 + m_beta/4.0 * diff*diff*diff*diff; en += energies[i]; } en = 1.0/en; for( size_t i=0 ; i<N ; i++ ) { energies[i] *= en; } } double energy( const coord_type &q , const coord_type &p ) { // calculates the total energy of the excitation const size_t N = q.size(); double en = 0.0; for( size_t i=0 ; i<N ; i++ ) { const double diff = q[(i+1) % N] - q[i]; en += p[i]*p[i]/2.0 + m_pot[m_pot_start_index+i]*q[i]*q[i]*q[i]*q[i] / 4.0 + m_beta/4.0 * diff*diff*diff*diff; } return en; } void change_pot_start( const int delta ) { m_pot_start_index += delta; } };
The total size we allow is 1024 and we start with an initial state size of 60.
//start with 60 sites const int N_start = 60; coord_type q( N_start , 0.0 ); q.reserve( max_N ); coord_type p( N_start , 0.0 ); p.reserve( max_N ); // start with uniform momentum distribution over 20 sites fill( p.begin()+20 , p.end()-20 , 1.0/sqrt(20.0) ); coord_type distr( N_start , 0.0 ); distr.reserve( max_N ); // create the system compacton_lattice lattice( max_N , beta , (max_N-N_start)/2 ); //create the stepper, note that we use an always_resizer because state size might change during steps typedef symplectic_rkn_sb3a_mclachlan< coord_type , coord_type , double , coord_type , coord_type , double , range_algebra , default_operations , always_resizer > hamiltonian_stepper; hamiltonian_stepper stepper; hamiltonian_stepper::state_type state = make_pair( q , p );
The lattice gets resized whenever the energy distribution comes close to
the borders distr[10] >
1E-150
, distr[distr.size()-10] >
1E-150
. If we increase to the left,
q
and p
have to be rotated because their resize function always appends at the end.
Additionally, the start index of the potential changes in this case.
double t = 0.0; const double dt = 0.1; const int steps = 10000; for( int step = 0 ; step < steps ; ++step ) { stepper.do_step( boost::ref(lattice) , state , t , dt ); lattice.energy_distribution( state.first , state.second , distr ); if( distr[10] > 1E-150 ) { do_resize( state.first , state.second , distr , state.first.size()+20 ); rotate( state.first.begin() , state.first.end()-20 , state.first.end() ); rotate( state.second.begin() , state.second.end()-20 , state.second.end() ); lattice.change_pot_start( -20 ); cout << t << ": resized left to " << distr.size() << ", energy = " << lattice.energy( state.first , state.second ) << endl; } if( distr[distr.size()-10] > 1E-150 ) { do_resize( state.first , state.second , distr , state.first.size()+20 ); cout << t << ": resized right to " << distr.size() << ", energy = " << lattice.energy( state.first , state.second ) << endl; } t += dt; }
The do_resize
function simply
calls vector.resize
of q
, p
and distr
.
void do_resize( coord_type &q , coord_type &p , coord_type &distr , const int N ) { q.resize( N ); p.resize( N ); distr.resize( N ); }
The full example can be found in resizing_lattice.cpp