On 07/22/2014 02:44 PM, Erik Schultheis wrote:
Thanks for the fast response.
Hi
On 07/22/2014 01:24 AM, Erik Schultheis wrote:
Dear all,
i am simulating particle trajectories, and depending on the data i want to generate i may or may not have to track additional data [boost.ublas c_matrix] along those trajectories. Right now i enable/disable that via a template parameter (so, using std::tuple(vec, vec) or tuple(vec, vec, c_matrix) and fusion_algebra), which has a few drawbacks. The algorithm is already templated to support varying dimensions via a template parameter, so the additional data practically doubles the number of template instantiations for system and all possible observers. This actually got so bad that gcc (mingw gcc 4.8) refuses to compile with debugging enabled. Therefore, i would like to make the switch between tracking the matrix and ignoring it a runtime choice instead of a template parameter.
I don't understand fully what you exactly want to do. Can you show us some more code? Can you put the tracking into an observer? I mean, do you need to solve an ODE to track the additional data?
Unfortunatly, the tracking is an ode itself, so the change depends on the position and the current matrix value, so it cannot be done inside of an observer.
Ok, very interesting problem. I had a project where we had a given fluid flow and a system of ODEs describing interacting particles moving in the fluid. This sounds similar to your project.
As for the current implementation: I have a class BasicTracer<S, I>, where S is the template parameter that decides the state type (currently one of std::array<vec<2>, 2>, std::array<vec<3>, 2>, std::tuple<vec<2>, <vec<2>, matrix<4>> and std::tuple<vec<3>, vec<3>, matrix<6>>) and I the integration startegy (either Integrator, that only operates on the vectors, or MonodromyIntegrator, that also operates on the matrix). The reason for this is that I have two different use cases for the program: 1) a human-viewable visualisation of the trajectories that does not need matrices, and 2) the actual data i want to generate. Since 1) needs far more trajectories to be calculated, it is really helpful that the computation time in this case is lower. With the mentioned limitation (gcc refuses to compile with -g, "tracer.o: too many sections (50892)" and "ccQF7l7o.s: Fatal error: can't close .obj\Debug\src\tracer.o: File too big") the code works as it is now.
For now i have only worked with small problem sizes, where the wait time for a result is in the order of a few minutes, however, this is going to change and therefore i want to parallelize the code. The differentiation between he use cases via templates make the code quite complex at times, and i hope to reduce that complexity before i add the parallelization. Therefore, I want to transform from a templated approach to a runtime approach. I am not sure if it is possible to make the code select the dimension at runtime without significant performance impact, but at least the second parameter should be removable without increasing runtime significantly.
I am not really sure which parts of the code to post to significantly inprove my description. The tracer itself is a type like
template class BasicTracer<std::tuple<vec<3>, vec<3>, monodromy_matrix<3>>, MonodromyIntegrator<std::tuple<vec<3>, vec<3>, monodromy_matrix<3>>>>;
or
template class BasicTracer<std::array<vec<2>, 2>, Integrator<std::array<vec<2>, 2>>>;
the integration is performed like
double starte = mEnergy(p);
try
{
mMasterObserver.startTrajectory(p);
perform_integration(mIntegrator, p, 2.0, 1.0 / mPotential.getSize(), mMasterObserver, mAbsErrorBound, mRelErrorBound);
} catch(int& i) {};
/// \todo this is the only reason for getLastState's existence
mEnergyError += std::pow( starte - mEnergy( mMasterObserver.getLastState() ), 2 );
where perform integration just sets parameters to the odeint function:
template<typename System, typename State, typename Observer> std::size_t perform_integration(const System& system, State& start, double t_end, double dt, Observer& obs, double abs_err, double rel_err) { using namespace boost::numeric::odeint;
// when range access is possible, try range algebra, otherwise, try fusion using algebra_type = typename std::conditional<boost::has_range_iterator<State>::value, range_algebra, fusion_algebra>::type;
typedef runge_kutta_cash_karp54<State, double, State, double, algebra_type> error_stepper_type; typedef controlled_runge_kutta<error_stepper_type> controlled_stepper_type;
typedef default_error_checker<typename controlled_stepper_type::value_type, typename controlled_stepper_type::algebra_type> error_checker_type;
controlled_stepper_type stepper(error_checker_type(abs_err, rel_err)); return integrate_const(stepper, system, start, 0.0, t_end, dt, std::ref(obs) ); }
and the operator() of the integrators look like
template<class S>
void Integrator<S>::operator()( const S& state, S& deriv, double/*t*/ ) const
{
using namespace system_state;
auto position = get_position(state);
auto p = position * mPotential.getSize();
double x = p[0];
double y = p[1];
// prevent out of bounds access
if(y < 0 || y > mPotential.getSize()-1 || x < 0 || x > mPotential.getSize()-1)
{
throw(1);
};
// calculate acceleration
auto& acceleration = get_velocity(deriv);
for(int i = 0; i < DIMENSION; ++i)
{
/// the compiler should be able to ompimize this
std::array<int, DIMENSION> didx;
didx.fill(0);
didx[i] = 1;
acceleration[i] = -linearInterpolate( mPotential.getDerivative(didx), p );
}
// change in position = current velocity
get_position( deriv ) = get_velocity( state );
}
template<class S> void MonodromyIntegrator<S>::operator()( const S& state, S& deriv, double/*t*/ ) const { using namespace system_state; auto position = get_position(state);
auto p = position * mPotential.getSize(); double x = p[0]; double y = p[1];
// prevent out of bounds access if(y < 0 || y > mPotential.getSize()-1 || x < 0 || x > mPotential.getSize()-1) { throw(1); };
// calculate acceleration auto& acceleration = get_velocity(deriv); for(int i = 0; i < DIMENSION; ++i) { std::array<int, DIMENSION> didx; didx.fill(0); didx[i] = 1; acceleration[i] = -linearInterpolate( mPotential.getDerivative(didx), p ); }
// change in position = current velocity get_position( deriv ) = get_velocity( state );
// change in monodromy boost::numeric::ublas::axpy_prod(getMonodromyCoeff( mPotential, p ), get_data(state), get_data(deriv), true); }
If there is some other part of the code that would be helpful, please let me know.
Is BasicTracer the "main" class and the function you showed here are methods of BasicTracer? In this case you could remove the dependency of the integrator into a boost::function or std::function object, like std::function< void( S const& , S const& , double > mIntegrator; This remove the integrator template parameter. If BasicTracer is not the "main" class, can you explain a bit more how BasicTracer works and what it should do?