43 #ifndef PARTICLEITERATOR_HPP
44 #define PARTICLEITERATOR_HPP 1
51 #include <gsl/gsl_odeiv.h>
52 #include <gsl/gsl_poly.h>
99 return(
_x[0] < cd.
_x[0] );
111 const PP &x1,
const PP &x2 ) {
115 for(
size_t a = 0; a < PP::dim(); a++ ) {
117 int a1 = (int)floor( (x1[2*a+1]-mesh.
origo(a))/mesh.
h() );
118 int a2 = (int)floor( (x2[2*a+1]-mesh.
origo(a))/mesh.
h() );
125 for(
int b = a1+1; b <= a2; b++ ) {
128 double K = (b*mesh.
h() + mesh.
origo(a) - x1[2*a+1]) /
129 (x2[2*a+1] - x1[2*a+1]);
130 if( K < 0.0 ) K = 0.0;
131 else if( K > 1.0 ) K = 1.0;
134 if( x2[2*a+1] > x1[2*a+1] )
135 coldata.push_back(
ColData( x1 + (x2-x1)*K, a+1 ) );
137 coldata.push_back(
ColData( x1 + (x2-x1)*K, -a-1 ) );
142 sort( coldata.begin(), coldata.end() );
154 const PP &x1,
const PP &x2 ) {
156 #ifdef DEBUG_PARTICLE_ITERATOR
157 std::cout <<
"Building coldata using polynomial interpolation\n";
164 for(
size_t a = 0; a < PP::dim(); a++ ) {
166 x1[2*a+1], x1[2*a+2],
167 x2[2*a+1], x2[2*a+2] );
171 for(
size_t a = 0; a < PP::dim(); a++ ) {
174 int i = (int)floor( (x1[2*a+1]-mesh.
origo(a))/mesh.
h() );
177 for(
int dj = -1; dj <= 1; dj += 2 ) {
186 double val = mesh.
origo(a) + mesh.
h() * j;
187 if( val < mesh.
origo(a) )
189 else if( val > mesh.
max(a) )
192 #ifdef DEBUG_PARTICLE_ITERATOR
193 std::cout <<
" Searching intersections at coord(" << a <<
") = " << val <<
"\n";
195 Kcount = traj[a].
solve( K, val );
199 #ifdef DEBUG_PARTICLE_ITERATOR
200 std::cout <<
" Found " << Kcount <<
" valid roots: ";
201 for(
int p = 0; p < Kcount; p++ )
202 std::cout << K[p] <<
" ";
207 for(
int b = 0; b < Kcount; b++ ) {
210 xcol(0) = x1[0] + K[b]*(x2[0]-x1[0]);
211 for(
size_t c = 0; c < PP::dim(); c++ ) {
212 traj[c].
coord( x, v, K[b] );
220 xcol[5] = x1[5] + K[b]*(x2[5]-x1[5]);
221 if( xcol[2*a+2] >= 0.0 )
222 coldata.push_back(
ColData( xcol, a+1 ) );
224 coldata.push_back(
ColData( xcol, -a-1 ) );
233 sort( coldata.begin(), coldata.end() );
235 #ifdef DEBUG_PARTICLE_ITERATOR
236 std::cout <<
" Coldata built\n";
252 gsl_odeiv_system _system;
253 gsl_odeiv_step *_step;
254 gsl_odeiv_control *_control;
255 gsl_odeiv_evolve *_evolve;
273 pthread_mutex_t *_scharge_mutex;
277 std::vector<PP> _traj;
278 std::vector<ColData<PP> > _coldata;
288 void save_trajectory_point( PP x ) {
291 _traj.push_back( x );
292 }
catch( std::bad_alloc ) {
305 bool check_collision(
Particle<PP> &particle,
const PP &x1,
const PP &x2, PP &status_x ) {
308 Vec3D v2 = x2.location();
313 Vec3D v1 = x1.location();
317 for(
size_t a = 0; a < PP::size(); a++ )
318 status_x[a] = x2[a] + K*(x1[a]-x2[a]);
321 for(
size_t a = _traj.size()-1; a > 0; a-- ) {
322 if( _traj[a][0] > status_x[0] )
346 void handle_mirror(
size_t c,
int i[3],
size_t a,
int border, PP &x2 ) {
348 #ifdef DEBUG_PARTICLE_ITERATOR
349 std::cout <<
" handle_mirror( c = " << c
350 <<
", i = (" << i[0] <<
", " << i[1] <<
", " << i[2]
351 <<
"), a = " << a <<
", border = " << border
364 #ifdef DEBUG_PARTICLE_ITERATOR
365 std::cout <<
" xmirror = " << xmirror <<
"\n";
366 std::cout <<
" i = (" << i[0] <<
", " << i[1] <<
", " << i[2] <<
")\n";
367 std::cout <<
" xi = " << _xi <<
"\n";
371 bool caught_at_boundary =
false;
372 if( _coldata[c]._dir == border*((
int)a+1) &&
373 ( i[a] == 0 || i[a] == (
int)_pidata.
_geom->
size(a)-2 ) ) {
374 caught_at_boundary =
true;
375 #ifdef DEBUG_PARTICLE_ITERATOR
376 std::cout <<
" caught_at_boundary\n";
381 if( caught_at_boundary ) {
382 save_trajectory_point( _coldata[c]._x );
384 for(
int b = _traj.size()-1; b > 0; b-- ) {
385 if( _traj[b][0] >= _xi[0] ) {
387 #ifdef DEBUG_PARTICLE_ITERATOR
388 std::cout <<
" mirroring traj[" << b <<
"] = " << _traj[b] <<
"\n";
390 _traj[b][2*a+1] = 2.0*xmirror - _traj[b][2*a+1];
391 _traj[b][2*a+2] *= -1.0;
398 for(
size_t b = c; b < _coldata.size(); b++ ) {
399 if( (
size_t)abs(_coldata[b]._dir) == a+1 )
400 _coldata[b]._dir *= -1;
401 _coldata[b]._x[2*a+1] = 2.0*xmirror - _coldata[b]._x[2*a+1];
402 _coldata[b]._x[2*a+2] *= -1.0;
405 if( caught_at_boundary )
406 save_trajectory_point( _coldata[c]._x );
409 x2[2*a+1] = 2.0*xmirror - x2[2*a+1];
413 gsl_odeiv_step_reset( _step );
414 gsl_odeiv_evolve_reset( _evolve );
418 void handle_collision(
Particle<PP> &particle, uint32_t bound,
size_t c, PP &status_x ) {
420 #ifdef DEBUG_PARTICLE_ITERATOR
421 std::cout <<
" handle_collision()\n";
425 status_x = _coldata[c]._x;
439 bool handle_trajectory_advance(
Particle<PP> &particle,
size_t c,
int i[3], PP &x2 ) {
442 if( PP::dim() == 2 ) {
443 if( _coldata[c]._dir == -1 ) {
444 if( ( abs(_pidata.
_geom->
mesh(i[0], i[1] )) >= 7 ||
445 abs(_pidata.
_geom->
mesh(i[0], i[1]+1)) >= 7 ) &&
446 !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
449 }
else if( _coldata[c]._dir == +1 ) {
450 if( ( abs(_pidata.
_geom->
mesh(i[0]+1,i[1] )) >= 7 ||
451 abs(_pidata.
_geom->
mesh(i[0]+1,i[1]+1)) >= 7 ) &&
452 !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
455 }
else if( _coldata[c]._dir == -2 ) {
456 if( ( abs(_pidata.
_geom->
mesh(i[0], i[1] )) >= 7 ||
457 abs(_pidata.
_geom->
mesh(i[0]+1,i[1] )) >= 7 ) &&
458 !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
462 if( ( abs(_pidata.
_geom->
mesh(i[0], i[1]+1)) >= 7 ||
463 abs(_pidata.
_geom->
mesh(i[0]+1,i[1]+1)) >= 7 ) &&
464 !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
468 }
else if( PP::dim() == 3 ) {
469 if( _coldata[c]._dir == -1 ) {
470 if( ( abs(_pidata.
_geom->
mesh(i[0], i[1], i[2] )) >= 7 ||
471 abs(_pidata.
_geom->
mesh(i[0], i[1]+1,i[2] )) >= 7 ||
472 abs(_pidata.
_geom->
mesh(i[0], i[1], i[2]+1)) >= 7 ||
473 abs(_pidata.
_geom->
mesh(i[0], i[1]+1,i[2]+1)) >= 7 ) &&
474 !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
477 }
else if( _coldata[c]._dir == +1 ) {
478 if( ( abs(_pidata.
_geom->
mesh(i[0]+1,i[1], i[2] )) >= 7 ||
479 abs(_pidata.
_geom->
mesh(i[0]+1,i[1]+1,i[2] )) >= 7 ||
480 abs(_pidata.
_geom->
mesh(i[0]+1,i[1], i[2]+1)) >= 7 ||
481 abs(_pidata.
_geom->
mesh(i[0]+1,i[1]+1,i[2]+1)) >= 7 ) &&
482 !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
485 }
else if( _coldata[c]._dir == -2 ) {
486 if( ( abs(_pidata.
_geom->
mesh(i[0], i[1],i[2] )) >= 7 ||
487 abs(_pidata.
_geom->
mesh(i[0]+1,i[1],i[2] )) >= 7 ||
488 abs(_pidata.
_geom->
mesh(i[0], i[1],i[2]+1)) >= 7 ||
489 abs(_pidata.
_geom->
mesh(i[0]+1,i[1],i[2]+1)) >= 7 ) &&
490 !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
493 }
else if( _coldata[c]._dir == +2 ) {
494 if( ( abs(_pidata.
_geom->
mesh(i[0], i[1]+1,i[2] )) >= 7 ||
495 abs(_pidata.
_geom->
mesh(i[0]+1,i[1]+1,i[2] )) >= 7 ||
496 abs(_pidata.
_geom->
mesh(i[0], i[1]+1,i[2]+1)) >= 7 ||
497 abs(_pidata.
_geom->
mesh(i[0]+1,i[1]+1,i[2]+1)) >= 7 ) &&
498 !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
501 }
else if( _coldata[c]._dir == -3 ) {
502 if( ( abs(_pidata.
_geom->
mesh(i[0], i[1], i[2] )) >= 7 ||
503 abs(_pidata.
_geom->
mesh(i[0]+1,i[1], i[2] )) >= 7 ||
504 abs(_pidata.
_geom->
mesh(i[0], i[1]+1,i[2])) >= 7 ||
505 abs(_pidata.
_geom->
mesh(i[0]+1,i[1]+1,i[2])) >= 7 ) &&
506 !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
510 if( ( abs(_pidata.
_geom->
mesh(i[0], i[1], i[2]+1)) >= 7 ||
511 abs(_pidata.
_geom->
mesh(i[0]+1,i[1], i[2]+1)) >= 7 ||
512 abs(_pidata.
_geom->
mesh(i[0], i[1]+1,i[2]+1)) >= 7 ||
513 abs(_pidata.
_geom->
mesh(i[0]+1,i[1]+1,i[2]+1)) >= 7 ) &&
514 !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
524 for(
size_t a = 0; a < PP::dim(); a++ ) {
528 handle_mirror( c, i, a, -1, x2 );
530 handle_collision( particle, 1+2*a, c, x2 );
533 }
else if( i[a] >= (_pidata.
_geom->
size(a)-1) ) {
535 handle_mirror( c, i, a, +1, x2 );
537 handle_collision( particle, 2+2*a, c, x2 );
551 bool limit_trajectory_advance(
const PP &x1, PP &x2 ) {
553 bool touched =
false;
555 for(
size_t a = 0; a < PP::dim(); a++ ) {
562 if( x2[2*a+1] < lim1 ) {
564 double K = (lim1 - x1[2*a+1]) / (x2[2*a+1] - x1[2*a+1]);
567 #ifdef DEBUG_PARTICLE_ITERATOR
568 std::cout <<
"Limiting step to:\n";
569 std::cout <<
" x2: " << x2 <<
"\n";
571 }
else if(x2[2*a+1] > lim2 ) {
573 double K = (lim2 - x1[2*a+1]) / (x2[2*a+1] - x1[2*a+1]);
576 #ifdef DEBUG_PARTICLE_ITERATOR
577 std::cout <<
"Limiting step to:\n";
578 std::cout <<
" x2: " << x2 <<
"\n";
590 void build_coldata(
bool force_linear,
const PP &x1,
const PP &x2 ) {
593 if( _polyint && !force_linear )
597 }
catch( std::bad_alloc ) {
622 bool handle_trajectory(
Particle<PP> &particle,
const PP &x1, PP &x2,
623 bool force_linear,
bool first_step ) {
625 #ifdef DEBUG_PARTICLE_ITERATOR
626 std::cout <<
"Handle trajectory from x1 to x2:\n";
627 std::cout <<
" x1: " << x1 <<
"\n";
628 std::cout <<
" x2: " << x2 <<
"\n";
633 if( limit_trajectory_advance( x1, x2 ) )
636 build_coldata( force_linear, x1, x2 );
643 if( _coldata.size() == 0 ) {
644 #ifdef DEBUG_PARTICLE_ITERATOR
645 std::cout <<
"No coldata\n";
651 int i[3] = {0, 0, 0};
652 for(
size_t cdir = 0; cdir < PP::dim(); cdir++ )
653 i[cdir] = (
int)floor( (x1[2*cdir+1]-_pidata.
_geom->
origo(cdir))/_pidata.
_geom->
h() );
656 #ifdef DEBUG_PARTICLE_ITERATOR
657 std::cout <<
"Process coldata points:\n";
659 for(
size_t a = 0; a < _coldata.size(); a++ ) {
661 #ifdef DEBUG_PARTICLE_ITERATOR
662 std::cout <<
" Coldata " << std::setw(4) << a <<
": "
663 << _coldata[a]._x <<
", "
664 << std::setw(3) << i[0] <<
" "
665 << std::setw(3) << i[1] <<
" "
666 << std::setw(3) << i[2] <<
" "
667 << std::setw(3) << _coldata[a]._dir <<
"\n";
672 handle_trajectory_advance( particle, a, i, x2 );
677 _xi, _coldata[a]._x );
681 (*_thand_cb)( &particle, &_coldata[a]._x, &x2 );
683 #ifdef DEBUG_PARTICLE_ITERATOR
685 std::cout <<
" Particle out\n";
686 std::cout <<
" x = " << x2 <<
"\n";
688 std::cout <<
" Particle collided\n";
689 std::cout <<
" x = " << x2 <<
"\n";
694 save_trajectory_point( x2 );
700 _xi = _coldata[a]._x;
703 #ifdef DEBUG_PARTICLE_ITERATOR
704 std::cout <<
"Coldata done\n";
713 bool axis_mirror_required(
const PP &x2 ) {
716 x2[3] <= 0.01*_pidata.
_geom->
h() &&
717 x2[3]*fabs(x2[5]) <= 1.0e-9*fabs(x2[4]) );
727 bool handle_axis_mirror_step(
Particle<PP> &particle,
const PP &x1, PP &x2 ) {
731 PP::get_derivatives( x2[0], &x2[1], dxdt, (
void *)&_pidata );
735 double dt = -x2[3]/x2[4];
738 xc[1] = x2[1]+(x2[2]+0.5*dxdt[1]*dt)*dt;
740 xc[3] = x2[3]+x2[4]*dt;
750 #ifdef DEBUG_PARTICLE_ITERATOR
751 std::cout <<
"Particle mirror:\n";
752 std::cout <<
" x1: " << x1 <<
"\n";
753 std::cout <<
" x2: " << x2 <<
"\n";
754 std::cout <<
" xc: " << xc <<
"\n";
755 std::cout <<
" x3: " << x3 <<
"\n";
759 if( !handle_trajectory( particle, x2, x3,
true,
false ) )
763 save_trajectory_point( x2 );
764 save_trajectory_point( xc );
767 save_trajectory_point( xc );
771 gsl_odeiv_step_reset( _step );
772 gsl_odeiv_evolve_reset( _evolve );
792 bool check_particle_definition( PP &x ) {
794 #ifdef DEBUG_PARTICLE_ITERATOR
795 std::cout <<
"Particle defined at:\n";
796 std::cout <<
" x = " << x <<
"\n";
800 for(
size_t a = 0; a < PP::size(); a++ ) {
844 double calculate_dt(
const PP &x,
const double *dxdt ) {
846 double spd = 0.0, acc = 0.0;
848 for(
size_t a = 0; a < (PP::size()-1)/2; a++ ) {
850 spd += dxdt[2*a]*dxdt[2*a];
852 acc += dxdt[2*a+1]*dxdt[2*a+1];
857 spd += x[3]*x[3]*x[5]*x[5];
859 acc += x[3]*x[3]*dxdt[4]*dxdt[4];
863 spd = _pidata.
_geom->
h() / sqrt(spd);
864 acc = sqrt( 2.0*_pidata.
_geom->
h() / sqrt(acc) );
866 return( spd < acc ? spd : acc );
898 bool polyint, uint32_t maxsteps,
double maxt,
899 uint32_t trajdiv,
bool mirror[6],
ScalarField *scharge,
900 pthread_mutex_t *scharge_mutex,
903 : _type(type), _polyint(polyint), _epsabs(epsabs), _epsrel(epsrel), _maxsteps(maxsteps), _maxt(maxt),
904 _trajdiv(trajdiv), _pidata(scharge,efield,bfield,geom),
905 _thand_cb(0), _tend_cb(0), _bsup_cb(0), _pdb(0), _scharge_mutex(scharge_mutex),
906 _stat(geom->number_of_boundaries()) {
909 _mirror[0] = mirror[0];
910 _mirror[1] = mirror[1];
911 _mirror[2] = mirror[2];
912 _mirror[3] = mirror[3];
913 _mirror[4] = mirror[4];
914 _mirror[5] = mirror[5];
917 _system.jacobian = NULL;
918 _system.params = (
void *)&_pidata;
919 _system.function = PP::get_derivatives;
920 _system.dimension = PP::size()-1;
926 double scale_abs[PP::size()-1];
927 for( uint32_t a = 0; a < (uint32_t)PP::size()-2; a+=2 ) {
928 scale_abs[a+0] = 1.0;
929 scale_abs[a+1] = 1.0e6;
935 _step = gsl_odeiv_step_alloc( gsl_odeiv_step_rkck, _system.dimension );
937 _control = gsl_odeiv_control_scaled_new( _epsabs, _epsrel, 1.0, 1.0, scale_abs, PP::size()-1 );
938 _evolve = gsl_odeiv_evolve_alloc( _system.dimension );
945 gsl_odeiv_evolve_free( _evolve );
946 gsl_odeiv_control_free( _control );
947 gsl_odeiv_step_free( _step );
954 _thand_cb = thand_cb;
994 PP x = particle->
x();
997 if( !check_particle_definition( x ) ) {
1006 save_trajectory_point( x );
1007 #ifdef DEBUG_PARTICLE_ITERATOR
1008 std::cout << x[0] <<
" "
1014 _pidata.
_qm = particle->
qm();
1018 gsl_odeiv_step_reset( _step );
1019 gsl_odeiv_evolve_reset( _evolve );
1023 double dxdt[PP::size()-1];
1024 PP::get_derivatives( 0.0, &x[1], dxdt, (
void *)&_pidata );
1025 double dt = calculate_dt( x, dxdt );
1027 #ifdef DEBUG_PARTICLE_ITERATOR
1028 std::cout <<
"dxdt = ";
1029 for(
size_t a = 0; a < PP::size()-1; a++ )
1030 std::cout << dxdt[a] <<
" ";
1032 std::cout <<
"dt = " << dt <<
"\n";
1033 std::cout <<
"*** Starting iteration\n";
1040 while( nstp < _maxsteps && x[0] < _maxt ) {
1042 #ifdef DEBUG_PARTICLE_ITERATOR
1043 std::cout <<
"\n*** Step ***\n";
1044 std::cout <<
" x = " << x2 <<
"\n";
1045 std::cout <<
" dt = " << dt <<
" (proposed)\n";
1052 int retval = gsl_odeiv_evolve_apply( _evolve, _control, _step, &_system,
1053 &x2[0], _maxt, &dt, &x2[1] );
1055 #ifdef DEBUG_PARTICLE_ITERATOR
1056 std::cout <<
"Step rejected\n";
1057 std::cout <<
" x2 = " << x2 <<
"\n";
1058 std::cout <<
" dt = " << dt <<
"\n";
1067 }
else if( retval == GSL_SUCCESS ) {
1075 if( nstp >= _maxsteps )
1080 #ifdef DEBUG_PARTICLE_ITERATOR
1081 std::cout <<
"Step accepted from x1 to x2:\n";
1082 std::cout <<
" dt = " << dt <<
" (taken)\n";
1083 std::cout <<
" x1 = " << x <<
"\n";
1084 std::cout <<
" x2 = " << x2 <<
"\n";
1088 if( !handle_trajectory( *particle, x, x2,
false, nstp == 0 ) ) {
1095 if( axis_mirror_required( x2 ) ) {
1096 if( !handle_axis_mirror_step( *particle, x, x2 ) )
1104 save_trajectory_point( x2 );
1110 #ifdef DEBUG_PARTICLE_ITERATOR
1111 std::cout <<
"\n*** Done stepping ***\n";
1115 if( nstp == _maxsteps ) {
1118 }
else if( x[0] >= _maxt ) {
1127 if( _trajdiv != 0 && pi % _trajdiv == 0 )
1135 (*_tend_cb)( particle, _pdb );