00001
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043 #ifndef PARTICLEITERATOR_HPP
00044 #define PARTICLEITERATOR_HPP 1
00045
00046
00047 #include <vector>
00048 #include <iostream>
00049 #include <algorithm>
00050 #include <iomanip>
00051 #include <gsl/gsl_odeiv.h>
00052 #include <gsl/gsl_poly.h>
00053 #include "geometry.hpp"
00054 #include "compmath.hpp"
00055 #include "trajectory.hpp"
00056 #include "particles.hpp"
00057 #include "vectorfield.hpp"
00058 #include "scalarfield.hpp"
00059 #include "scharge.hpp"
00060 #include "scheduler.hpp"
00061 #include "polysolver.hpp"
00062 #include "particledatabase.hpp"
00063
00064
00065
00066
00067
00070 enum particle_iterator_type_e {
00071 PARTICLE_ITERATOR_ADAPTIVE = 0,
00072 PARTICLE_ITERATOR_FIXED_STEP_LEN
00073 };
00074
00075
00084 template <class PP> class ColData {
00085 public:
00086 PP _x;
00087 int _dir;
00092 ColData( PP x, int dir ) : _x(x), _dir(dir) {}
00093
00098 bool operator<( const ColData &cd ) const {
00099 return( _x[0] < cd._x[0] );
00100 }
00101
00110 static void build_coldata_linear( std::vector<ColData> &coldata, const Mesh &mesh,
00111 const PP &x1, const PP &x2 ) {
00112
00113 coldata.clear();
00114
00115 for( size_t a = 0; a < PP::dim(); a++ ) {
00116
00117 int a1 = (int)floor( (x1[2*a+1]-mesh.origo(a))/mesh.h() );
00118 int a2 = (int)floor( (x2[2*a+1]-mesh.origo(a))/mesh.h() );
00119 if( a1 > a2 ) {
00120 int a = a2;
00121 a2 = a1;
00122 a1 = a;
00123 }
00124
00125 for( int b = a1+1; b <= a2; b++ ) {
00126
00127
00128 double K = (b*mesh.h() + mesh.origo(a) - x1[2*a+1]) /
00129 (x2[2*a+1] - x1[2*a+1]);
00130 if( K < 0.0 ) K = 0.0;
00131 else if( K > 1.0 ) K = 1.0;
00132
00133
00134 if( x2[2*a+1] > x1[2*a+1] )
00135 coldata.push_back( ColData( x1 + (x2-x1)*K, a+1 ) );
00136 else
00137 coldata.push_back( ColData( x1 + (x2-x1)*K, -a-1 ) );
00138 }
00139 }
00140
00141
00142 sort( coldata.begin(), coldata.end() );
00143 }
00144
00153 static void build_coldata_poly( std::vector<ColData> &coldata, const Mesh &mesh,
00154 const PP &x1, const PP &x2 ) {
00155
00156 #ifdef DEBUG_PARTICLE_ITERATOR
00157 std::cout << "Building coldata using polynomial interpolation\n";
00158 #endif
00159
00160 coldata.clear();
00161
00162
00163 TrajectoryRep1D traj[PP::dim()];
00164 for( size_t a = 0; a < PP::dim(); a++ ) {
00165 traj[a].construct( x2[0]-x1[0],
00166 x1[2*a+1], x1[2*a+2],
00167 x2[2*a+1], x2[2*a+2] );
00168 }
00169
00170
00171 for( size_t a = 0; a < PP::dim(); a++ ) {
00172
00173
00174 int i = (int)floor( (x1[2*a+1]-mesh.origo(a))/mesh.h() );
00175
00176
00177 for( int dj = -1; dj <= 1; dj += 2 ) {
00178 int j = i;
00179 if( dj == +1 )
00180 j = i+1;
00181 int Kcount;
00182 double K[3];
00183 while( 1 ) {
00184
00185
00186 double val = mesh.origo(a) + mesh.h() * j;
00187 if( val < mesh.origo(a) )
00188 break;
00189 else if( val > mesh.max(a) )
00190 break;
00191
00192 #ifdef DEBUG_PARTICLE_ITERATOR
00193 std::cout << " Searching intersections at coord(" << a << ") = " << val << "\n";
00194 #endif
00195 Kcount = traj[a].solve( K, val );
00196 if( Kcount == 0 )
00197 break;
00198
00199 #ifdef DEBUG_PARTICLE_ITERATOR
00200 std::cout << " Found " << Kcount << " valid roots: ";
00201 for( int p = 0; p < Kcount; p++ )
00202 std::cout << K[p] << " ";
00203 std::cout << "\n";
00204 #endif
00205
00206
00207 for( int b = 0; b < Kcount; b++ ) {
00208 PP xcol;
00209 double x, v;
00210 xcol(0) = x1[0] + K[b]*(x2[0]-x1[0]);
00211 for( size_t c = 0; c < PP::dim(); c++ ) {
00212 traj[c].coord( x, v, K[b] );
00213 if( a == c )
00214 xcol[2*c+1] = val;
00215 else
00216 xcol[2*c+1] = x;
00217 xcol[2*c+2] = v;
00218 }
00219 if( mesh.geom_mode() == MODE_CYL )
00220 xcol[5] = x1[5] + K[b]*(x2[5]-x1[5]);
00221 if( xcol[2*a+2] >= 0.0 )
00222 coldata.push_back( ColData( xcol, a+1 ) );
00223 else
00224 coldata.push_back( ColData( xcol, -a-1 ) );
00225 }
00226
00227 j += dj;
00228 }
00229 }
00230 }
00231
00232
00233 sort( coldata.begin(), coldata.end() );
00234
00235 #ifdef DEBUG_PARTICLE_ITERATOR
00236 std::cout << " Coldata built\n";
00237 #endif
00238 }
00239
00240 };
00241
00242
00250 template <class PP> class ParticleIterator {
00251
00252 gsl_odeiv_system _system;
00253 gsl_odeiv_step *_step;
00254 gsl_odeiv_control *_control;
00255 gsl_odeiv_evolve *_evolve;
00257 particle_iterator_type_e _type;
00259 bool _polyint;
00260 double _epsabs;
00261 double _epsrel;
00262 uint32_t _maxsteps;
00263 double _maxt;
00264 uint32_t _trajdiv;
00266 bool _mirror[6];
00268 ParticleIteratorData _pidata;
00269 const TrajectoryHandlerCallback *_thand_cb;
00270 const TrajectoryEndCallback *_tend_cb;
00271 const TrajectoryEndCallback *_bsup_cb;
00272 ParticleDataBase *_pdb;
00273 pthread_mutex_t *_scharge_mutex;
00275 PP _xi;
00277 std::vector<PP> _traj;
00278 std::vector<ColData<PP> > _coldata;
00280 ParticleStatistics _stat;
00288 void save_trajectory_point( PP x ) {
00289
00290 try {
00291 _traj.push_back( x );
00292 } catch( std::bad_alloc ) {
00293 throw( ErrorNoMem( ERROR_LOCATION, "Out of memory saving trajectory" ) );
00294 }
00295 }
00296
00305 bool check_collision( Particle<PP> &particle, const PP &x1, const PP &x2, PP &status_x ) {
00306
00307
00308 Vec3D v2 = x2.location();
00309 int32_t bound = _pidata._geom->inside( v2 );
00310 if( bound < 7 )
00311 return( true );
00312 Vec3D vc;
00313 Vec3D v1 = x1.location();
00314 double K = _pidata._geom->bracket_surface( bound, v2, v1, vc );
00315
00316
00317 for( size_t a = 0; a < PP::size(); a++ )
00318 status_x[a] = x2[a] + K*(x1[a]-x2[a]);
00319
00320
00321 for( size_t a = _traj.size()-1; a > 0; a-- ) {
00322 if( _traj[a][0] > status_x[0] )
00323 _traj.pop_back();
00324 else
00325 break;
00326 }
00327
00328
00329
00330 particle.set_status( PARTICLE_COLL );
00331
00332
00333 _stat.add_bound_collision( bound, particle.IQ() );
00334
00335 return( false );
00336 }
00337
00338
00346 void handle_mirror( size_t c, int i[3], size_t a, int border, PP &x2 ) {
00347
00348 #ifdef DEBUG_PARTICLE_ITERATOR
00349 std::cout << " handle_mirror( c = " << c
00350 << ", i = (" << i[0] << ", " << i[1] << ", " << i[2]
00351 << "), a = " << a << ", border = " << border
00352 << ")\n";
00353 #endif
00354
00355 double xmirror;
00356 if( border < 0 ) {
00357 xmirror = _pidata._geom->origo(a);
00358 i[a] = -i[a]-1;
00359 } else {
00360 xmirror = _pidata._geom->max(a);
00361 i[a] = 2*_pidata._geom->size(a)-i[a]-3;
00362 }
00363
00364 #ifdef DEBUG_PARTICLE_ITERATOR
00365 std::cout << " xmirror = " << xmirror << "\n";
00366 std::cout << " i = (" << i[0] << ", " << i[1] << ", " << i[2] << ")\n";
00367 std::cout << " xi = " << _xi << "\n";
00368 #endif
00369
00370
00371 bool caught_at_boundary = false;
00372 if( _coldata[c]._dir == border*((int)a+1) &&
00373 ( i[a] == 0 || i[a] == (int)_pidata._geom->size(a)-2 ) ) {
00374 caught_at_boundary = true;
00375 #ifdef DEBUG_PARTICLE_ITERATOR
00376 std::cout << " caught_at_boundary\n";
00377 #endif
00378 }
00379
00380
00381 if( caught_at_boundary ) {
00382 save_trajectory_point( _coldata[c]._x );
00383 } else {
00384 for( int b = _traj.size()-1; b > 0; b-- ) {
00385 if( _traj[b][0] >= _xi[0] ) {
00386
00387 #ifdef DEBUG_PARTICLE_ITERATOR
00388 std::cout << " mirroring traj[" << b << "] = " << _traj[b] << "\n";
00389 #endif
00390 _traj[b][2*a+1] = 2.0*xmirror - _traj[b][2*a+1];
00391 _traj[b][2*a+2] *= -1.0;
00392 } else
00393 break;
00394 }
00395 }
00396
00397
00398 for( size_t b = c; b < _coldata.size(); b++ ) {
00399 if( (size_t)abs(_coldata[b]._dir) == a+1 )
00400 _coldata[b]._dir *= -1;
00401 _coldata[b]._x[2*a+1] = 2.0*xmirror - _coldata[b]._x[2*a+1];
00402 _coldata[b]._x[2*a+2] *= -1.0;
00403 }
00404
00405 if( caught_at_boundary )
00406 save_trajectory_point( _coldata[c]._x );
00407
00408
00409 x2[2*a+1] = 2.0*xmirror - x2[2*a+1];
00410 x2[2*a+2] *= -1.0;
00411
00412
00413 gsl_odeiv_step_reset( _step );
00414 gsl_odeiv_evolve_reset( _evolve );
00415 }
00416
00417
00418 void handle_collision( Particle<PP> &particle, uint32_t bound, size_t c, PP &status_x ) {
00419
00420 #ifdef DEBUG_PARTICLE_ITERATOR
00421 std::cout << " handle_collision()\n";
00422 #endif
00423
00424
00425 status_x = _coldata[c]._x;
00426 particle.set_status( PARTICLE_OUT );
00427 _stat.add_bound_collision( bound, particle.IQ() );
00428 }
00429
00430
00439 bool handle_trajectory_advance( Particle<PP> &particle, size_t c, int i[3], PP &x2 ) {
00440
00441
00442 if( PP::dim() == 2 ) {
00443 if( _coldata[c]._dir == -1 ) {
00444 if( ( abs(_pidata._geom->mesh(i[0], i[1] )) >= 7 ||
00445 abs(_pidata._geom->mesh(i[0], i[1]+1)) >= 7 ) &&
00446 !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
00447 return( false );
00448 i[0]--;
00449 } else if( _coldata[c]._dir == +1 ) {
00450 if( ( abs(_pidata._geom->mesh(i[0]+1,i[1] )) >= 7 ||
00451 abs(_pidata._geom->mesh(i[0]+1,i[1]+1)) >= 7 ) &&
00452 !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
00453 return( false );
00454 i[0]++;
00455 } else if( _coldata[c]._dir == -2 ) {
00456 if( ( abs(_pidata._geom->mesh(i[0], i[1] )) >= 7 ||
00457 abs(_pidata._geom->mesh(i[0]+1,i[1] )) >= 7 ) &&
00458 !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
00459 return( false );
00460 i[1]--;
00461 } else {
00462 if( ( abs(_pidata._geom->mesh(i[0], i[1]+1)) >= 7 ||
00463 abs(_pidata._geom->mesh(i[0]+1,i[1]+1)) >= 7 ) &&
00464 !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
00465 return( false );
00466 i[1]++;
00467 }
00468 } else if( PP::dim() == 3 ) {
00469 if( _coldata[c]._dir == -1 ) {
00470 if( ( abs(_pidata._geom->mesh(i[0], i[1], i[2] )) >= 7 ||
00471 abs(_pidata._geom->mesh(i[0], i[1]+1,i[2] )) >= 7 ||
00472 abs(_pidata._geom->mesh(i[0], i[1], i[2]+1)) >= 7 ||
00473 abs(_pidata._geom->mesh(i[0], i[1]+1,i[2]+1)) >= 7 ) &&
00474 !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
00475 return( false );
00476 i[0]--;
00477 } else if( _coldata[c]._dir == +1 ) {
00478 if( ( abs(_pidata._geom->mesh(i[0]+1,i[1], i[2] )) >= 7 ||
00479 abs(_pidata._geom->mesh(i[0]+1,i[1]+1,i[2] )) >= 7 ||
00480 abs(_pidata._geom->mesh(i[0]+1,i[1], i[2]+1)) >= 7 ||
00481 abs(_pidata._geom->mesh(i[0]+1,i[1]+1,i[2]+1)) >= 7 ) &&
00482 !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
00483 return( false );
00484 i[0]++;
00485 } else if( _coldata[c]._dir == -2 ) {
00486 if( ( abs(_pidata._geom->mesh(i[0], i[1],i[2] )) >= 7 ||
00487 abs(_pidata._geom->mesh(i[0]+1,i[1],i[2] )) >= 7 ||
00488 abs(_pidata._geom->mesh(i[0], i[1],i[2]+1)) >= 7 ||
00489 abs(_pidata._geom->mesh(i[0]+1,i[1],i[2]+1)) >= 7 ) &&
00490 !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
00491 return( false );
00492 i[1]--;
00493 } else if( _coldata[c]._dir == +2 ) {
00494 if( ( abs(_pidata._geom->mesh(i[0], i[1]+1,i[2] )) >= 7 ||
00495 abs(_pidata._geom->mesh(i[0]+1,i[1]+1,i[2] )) >= 7 ||
00496 abs(_pidata._geom->mesh(i[0], i[1]+1,i[2]+1)) >= 7 ||
00497 abs(_pidata._geom->mesh(i[0]+1,i[1]+1,i[2]+1)) >= 7 ) &&
00498 !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
00499 return( false );
00500 i[1]++;
00501 } else if( _coldata[c]._dir == -3 ) {
00502 if( ( abs(_pidata._geom->mesh(i[0], i[1], i[2] )) >= 7 ||
00503 abs(_pidata._geom->mesh(i[0]+1,i[1], i[2] )) >= 7 ||
00504 abs(_pidata._geom->mesh(i[0], i[1]+1,i[2])) >= 7 ||
00505 abs(_pidata._geom->mesh(i[0]+1,i[1]+1,i[2])) >= 7 ) &&
00506 !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
00507 return( false );
00508 i[2]--;
00509 } else {
00510 if( ( abs(_pidata._geom->mesh(i[0], i[1], i[2]+1)) >= 7 ||
00511 abs(_pidata._geom->mesh(i[0]+1,i[1], i[2]+1)) >= 7 ||
00512 abs(_pidata._geom->mesh(i[0], i[1]+1,i[2]+1)) >= 7 ||
00513 abs(_pidata._geom->mesh(i[0]+1,i[1]+1,i[2]+1)) >= 7 ) &&
00514 !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
00515 return( false );
00516 i[2]++;
00517 }
00518 } else {
00519 throw( Error( ERROR_LOCATION, "unsupported dimension number" ) );
00520 }
00521
00522
00523
00524 for( size_t a = 0; a < PP::dim(); a++ ) {
00525
00526 if( i[a] < 0 ) {
00527 if( _mirror[2*a] )
00528 handle_mirror( c, i, a, -1, x2 );
00529 else {
00530 handle_collision( particle, 1+2*a, c, x2 );
00531 return( false );
00532 }
00533 } else if( i[a] >= (_pidata._geom->size(a)-1) ) {
00534 if( _mirror[2*a+1] )
00535 handle_mirror( c, i, a, +1, x2 );
00536 else {
00537 handle_collision( particle, 2+2*a, c, x2 );
00538 return( false );
00539 }
00540 }
00541 }
00542
00543 return( true );
00544 }
00545
00551 bool limit_trajectory_advance( const PP &x1, PP &x2 ) {
00552
00553 bool touched = false;
00554
00555 for( size_t a = 0; a < PP::dim(); a++ ) {
00556
00557 double lim1 = _pidata._geom->origo(a) -
00558 (_pidata._geom->size(a)-1)*_pidata._geom->h();
00559 double lim2 = _pidata._geom->origo(a) +
00560 2*(_pidata._geom->size(a)-1)*_pidata._geom->h();
00561
00562 if( x2[2*a+1] < lim1 ) {
00563
00564 double K = (lim1 - x1[2*a+1]) / (x2[2*a+1] - x1[2*a+1]);
00565 x2 = x1 + K*(x2-x1);
00566 touched = true;
00567 #ifdef DEBUG_PARTICLE_ITERATOR
00568 std::cout << "Limiting step to:\n";
00569 std::cout << " x2: " << x2 << "\n";
00570 #endif
00571 } else if(x2[2*a+1] > lim2 ) {
00572
00573 double K = (lim2 - x1[2*a+1]) / (x2[2*a+1] - x1[2*a+1]);
00574 x2 = x1 + K*(x2-x1);
00575 touched = true;
00576 #ifdef DEBUG_PARTICLE_ITERATOR
00577 std::cout << "Limiting step to:\n";
00578 std::cout << " x2: " << x2 << "\n";
00579 #endif
00580 }
00581 }
00582
00583 return( touched );
00584 }
00585
00590 void build_coldata( bool force_linear, const PP &x1, const PP &x2 ) {
00591
00592 try {
00593 if( _polyint && !force_linear )
00594 ColData<PP>::build_coldata_poly( _coldata, *_pidata._geom, x1, x2 );
00595 else
00596 ColData<PP>::build_coldata_linear( _coldata, *_pidata._geom, x1, x2 );
00597 } catch( std::bad_alloc ) {
00598 throw( ErrorNoMem( ERROR_LOCATION, "Out of memory building ColData" ) );
00599 }
00600
00601 }
00602
00622 bool handle_trajectory( Particle<PP> &particle, const PP &x1, PP &x2,
00623 bool force_linear, bool first_step ) {
00624
00625 #ifdef DEBUG_PARTICLE_ITERATOR
00626 std::cout << "Handle trajectory from x1 to x2:\n";
00627 std::cout << " x1: " << x1 << "\n";
00628 std::cout << " x2: " << x2 << "\n";
00629 #endif
00630
00631
00632
00633 if( limit_trajectory_advance( x1, x2 ) )
00634 force_linear = true;
00635
00636 build_coldata( force_linear, x1, x2 );
00637
00638
00639
00640
00641
00642
00643 if( _coldata.size() == 0 ) {
00644 #ifdef DEBUG_PARTICLE_ITERATOR
00645 std::cout << "No coldata\n";
00646 #endif
00647 return( true );
00648 }
00649
00650
00651 int i[3] = {0, 0, 0};
00652 for( size_t cdir = 0; cdir < PP::dim(); cdir++ )
00653 i[cdir] = (int)floor( (x1[2*cdir+1]-_pidata._geom->origo(cdir))/_pidata._geom->h() );
00654
00655
00656 #ifdef DEBUG_PARTICLE_ITERATOR
00657 std::cout << "Process coldata points:\n";
00658 #endif
00659 for( size_t a = 0; a < _coldata.size(); a++ ) {
00660
00661 #ifdef DEBUG_PARTICLE_ITERATOR
00662 std::cout << " Coldata " << std::setw(4) << a << ": "
00663 << _coldata[a]._x << ", "
00664 << std::setw(3) << i[0] << " "
00665 << std::setw(3) << i[1] << " "
00666 << std::setw(3) << i[2] << " "
00667 << std::setw(3) << _coldata[a]._dir << "\n";
00668 #endif
00669
00670
00671
00672 handle_trajectory_advance( particle, a, i, x2 );
00673
00674
00675 if( _pidata._scharge )
00676 scharge_add_from_trajectory( *_pidata._scharge, _scharge_mutex, particle.IQ(),
00677 _xi, _coldata[a]._x );
00678
00679
00680 if( _thand_cb )
00681 (*_thand_cb)( &particle, &_coldata[a]._x, &x2 );
00682
00683 #ifdef DEBUG_PARTICLE_ITERATOR
00684 if( particle.get_status() == PARTICLE_OUT ) {
00685 std::cout << " Particle out\n";
00686 std::cout << " x = " << x2 << "\n";
00687 } else if( particle.get_status() == PARTICLE_COLL ) {
00688 std::cout << " Particle collided\n";
00689 std::cout << " x = " << x2 << "\n";
00690 }
00691 #endif
00692
00693 if( particle.get_status() != PARTICLE_OK ) {
00694 save_trajectory_point( x2 );
00695 _coldata.clear();
00696 return( false );
00697 }
00698
00699
00700 _xi = _coldata[a]._x;
00701 }
00702
00703 #ifdef DEBUG_PARTICLE_ITERATOR
00704 std::cout << "Coldata done\n";
00705 #endif
00706 _coldata.clear();
00707 return( true );
00708 }
00709
00710
00713 bool axis_mirror_required( const PP &x2 ) {
00714 return( _pidata._geom->geom_mode() == MODE_CYL &&
00715 x2[4] < 0.0 &&
00716 x2[3] <= 0.01*_pidata._geom->h() &&
00717 x2[3]*fabs(x2[5]) <= 1.0e-9*fabs(x2[4]) );
00718
00719 }
00720
00721
00727 bool handle_axis_mirror_step( Particle<PP> &particle, const PP &x1, PP &x2 ) {
00728
00729
00730 double dxdt[5];
00731 PP::get_derivatives( x2[0], &x2[1], dxdt, (void *)&_pidata );
00732
00733
00734
00735 double dt = -x2[3]/x2[4];
00736 PP xc;
00737 xc[0] = x2[0]+dt;
00738 xc[1] = x2[1]+(x2[2]+0.5*dxdt[1]*dt)*dt;
00739 xc[2] = x2[2];
00740 xc[3] = x2[3]+x2[4]*dt;
00741 xc[4] = x2[4];
00742 xc[5] = x2[5];
00743
00744
00745 PP x3 = 2*xc - x2;
00746 x3[3] *= -1.0;
00747 x3[4] *= -1.0;
00748 x3[5] *= -1.0;
00749
00750 #ifdef DEBUG_PARTICLE_ITERATOR
00751 std::cout << "Particle mirror:\n";
00752 std::cout << " x1: " << x1 << "\n";
00753 std::cout << " x2: " << x2 << "\n";
00754 std::cout << " xc: " << xc << "\n";
00755 std::cout << " x3: " << x3 << "\n";
00756 #endif
00757
00758
00759 if( !handle_trajectory( particle, x2, x3, true, false ) )
00760 return( false );
00761
00762
00763 save_trajectory_point( x2 );
00764 save_trajectory_point( xc );
00765 xc[4] *= -1.0;
00766 xc[5] *= -1.0;
00767 save_trajectory_point( xc );
00768
00769
00770
00771 gsl_odeiv_step_reset( _step );
00772 gsl_odeiv_evolve_reset( _evolve );
00773
00774
00775 x2 = x3;
00776 return( true );
00777 }
00778
00792 bool check_particle_definition( PP &x ) {
00793
00794 #ifdef DEBUG_PARTICLE_ITERATOR
00795 std::cout << "Particle defined at:\n";
00796 std::cout << " x = " << x << "\n";
00797 #endif
00798
00799
00800 for( size_t a = 0; a < PP::size(); a++ ) {
00801 if( comp_isnan( x[a] ) )
00802 return( false );
00803 }
00804
00805
00806 if( _pidata._geom->inside( x.location() ) )
00807 return( false );
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841 return( true );
00842 }
00843
00844 double calculate_dt( const PP &x, const double *dxdt ) {
00845
00846 double spd = 0.0, acc = 0.0;
00847
00848 for( size_t a = 0; a < (PP::size()-1)/2; a++ ) {
00849
00850 spd += dxdt[2*a]*dxdt[2*a];
00851
00852 acc += dxdt[2*a+1]*dxdt[2*a+1];
00853 }
00854 if( _pidata._geom->geom_mode() == MODE_CYL ) {
00855
00856
00857 spd += x[3]*x[3]*x[5]*x[5];
00858
00859 acc += x[3]*x[3]*dxdt[4]*dxdt[4];
00860 }
00861
00862
00863 spd = _pidata._geom->h() / sqrt(spd);
00864 acc = sqrt( 2.0*_pidata._geom->h() / sqrt(acc) );
00865
00866 return( spd < acc ? spd : acc );
00867 }
00868
00869 public:
00870
00897 ParticleIterator( particle_iterator_type_e type, double epsabs, double epsrel,
00898 bool polyint, uint32_t maxsteps, double maxt,
00899 uint32_t trajdiv, bool mirror[6], ScalarField *scharge,
00900 pthread_mutex_t *scharge_mutex,
00901 const VectorField *efield, const VectorField *bfield,
00902 const Geometry *geom )
00903 : _type(type), _polyint(polyint), _epsabs(epsabs), _epsrel(epsrel), _maxsteps(maxsteps), _maxt(maxt),
00904 _trajdiv(trajdiv), _pidata(scharge,efield,bfield,geom),
00905 _thand_cb(0), _tend_cb(0), _bsup_cb(0), _pdb(0), _scharge_mutex(scharge_mutex),
00906 _stat(geom->number_of_boundaries()) {
00907
00908
00909 _mirror[0] = mirror[0];
00910 _mirror[1] = mirror[1];
00911 _mirror[2] = mirror[2];
00912 _mirror[3] = mirror[3];
00913 _mirror[4] = mirror[4];
00914 _mirror[5] = mirror[5];
00915
00916
00917 _system.jacobian = NULL;
00918 _system.params = (void *)&_pidata;
00919 _system.function = PP::get_derivatives;
00920 _system.dimension = PP::size()-1;
00921
00922
00923
00924
00925
00926 double scale_abs[PP::size()-1];
00927 for( uint32_t a = 0; a < (uint32_t)PP::size()-2; a+=2 ) {
00928 scale_abs[a+0] = 1.0;
00929 scale_abs[a+1] = 1.0e6;
00930 }
00931 if( _pidata._geom->geom_mode() == MODE_CYL )
00932 scale_abs[4] = 1.0;
00933
00934
00935 _step = gsl_odeiv_step_alloc( gsl_odeiv_step_rkck, _system.dimension );
00936
00937 _control = gsl_odeiv_control_scaled_new( _epsabs, _epsrel, 1.0, 1.0, scale_abs, PP::size()-1 );
00938 _evolve = gsl_odeiv_evolve_alloc( _system.dimension );
00939 }
00940
00941
00944 ~ParticleIterator() {
00945 gsl_odeiv_evolve_free( _evolve );
00946 gsl_odeiv_control_free( _control );
00947 gsl_odeiv_step_free( _step );
00948 }
00949
00950
00953 void set_trajectory_handler_callback( const TrajectoryHandlerCallback *thand_cb ) {
00954 _thand_cb = thand_cb;
00955 }
00956
00957
00960 void set_trajectory_end_callback( const TrajectoryEndCallback *tend_cb, ParticleDataBase *pdb ) {
00961 _tend_cb = tend_cb;
00962 _pdb = pdb;
00963 }
00964
00965
00968 void set_bfield_suppression_callback( const CallbackFunctorD_V *bsup_cb ) {
00969 _pidata.set_bfield_suppression_callback( bsup_cb );
00970 }
00971
00974 const ParticleStatistics &get_statistics( void ) const {
00975 return( _stat );
00976 }
00977
00978
00987 void operator()( Particle<PP> *particle, uint32_t pi ) {
00988
00989
00990 if( particle->get_status() != PARTICLE_OK )
00991 return;
00992
00993
00994 PP x = particle->x();
00995
00996
00997 if( !check_particle_definition( x ) ) {
00998 particle->set_status( PARTICLE_BADDEF );
00999 _stat.inc_end_baddef();
01000 return;
01001 }
01002 particle->x() = x;
01003
01004
01005 _traj.clear();
01006 save_trajectory_point( x );
01007 #ifdef DEBUG_PARTICLE_ITERATOR
01008 std::cout << x[0] << " "
01009 << x[1] << " "
01010 << x[2] << " "
01011 << x[3] << " "
01012 << x[4] << "\n";
01013 #endif
01014 _pidata._qm = particle->qm();
01015 _xi = x;
01016
01017
01018 gsl_odeiv_step_reset( _step );
01019 gsl_odeiv_evolve_reset( _evolve );
01020
01021
01022
01023 double dxdt[PP::size()-1];
01024 PP::get_derivatives( 0.0, &x[1], dxdt, (void *)&_pidata );
01025 double dt = calculate_dt( x, dxdt );
01026
01027 #ifdef DEBUG_PARTICLE_ITERATOR
01028 std::cout << "dxdt = ";
01029 for( size_t a = 0; a < PP::size()-1; a++ )
01030 std::cout << dxdt[a] << " ";
01031 std::cout << "\n";
01032 std::cout << "dt = " << dt << "\n";
01033 std::cout << "*** Starting iteration\n";
01034 #endif
01035
01036
01037
01038 PP x2;
01039 size_t nstp = 0;
01040 while( nstp < _maxsteps && x[0] < _maxt ) {
01041
01042 #ifdef DEBUG_PARTICLE_ITERATOR
01043 std::cout << "\n*** Step ***\n";
01044 std::cout << " x = " << x2 << "\n";
01045 std::cout << " dt = " << dt << " (proposed)\n";
01046 #endif
01047
01048
01049 x2 = x;
01050
01051 while( true ) {
01052 int retval = gsl_odeiv_evolve_apply( _evolve, _control, _step, &_system,
01053 &x2[0], _maxt, &dt, &x2[1] );
01054 if( retval == IBSIMU_DERIV_ERROR ) {
01055 #ifdef DEBUG_PARTICLE_ITERATOR
01056 std::cout << "Step rejected\n";
01057 std::cout << " x2 = " << x2 << "\n";
01058 std::cout << " dt = " << dt << "\n";
01059 #endif
01060 x2[0] = x[0];
01061
01062 dt *= 0.5;
01063 if( dt == 0.0 )
01064 throw( Error( ERROR_LOCATION, "too small step size" ) );
01065
01066 continue;
01067 } else if( retval == GSL_SUCCESS ) {
01068 break;
01069 } else {
01070 throw( Error( ERROR_LOCATION, "gsl_odeiv_evolve_apply failed" ) );
01071 }
01072 }
01073
01074
01075 if( nstp >= _maxsteps )
01076 break;
01077 if( x2[0] == x[0] )
01078 throw( Error( ERROR_LOCATION, "too small step size" ) );
01079
01080 #ifdef DEBUG_PARTICLE_ITERATOR
01081 std::cout << "Step accepted from x1 to x2:\n";
01082 std::cout << " dt = " << dt << " (taken)\n";
01083 std::cout << " x1 = " << x << "\n";
01084 std::cout << " x2 = " << x2 << "\n";
01085 #endif
01086
01087
01088 if( !handle_trajectory( *particle, x, x2, false, nstp == 0 ) ) {
01089 x = x2;
01090 break;
01091 }
01092
01093
01094
01095 if( axis_mirror_required( x2 ) ) {
01096 if( !handle_axis_mirror_step( *particle, x, x2 ) )
01097 break;
01098 }
01099
01100
01101 x = x2;
01102
01103
01104 save_trajectory_point( x2 );
01105
01106
01107 nstp++;
01108 }
01109
01110 #ifdef DEBUG_PARTICLE_ITERATOR
01111 std::cout << "\n*** Done stepping ***\n";
01112 #endif
01113
01114
01115 if( nstp == _maxsteps ) {
01116 particle->set_status( PARTICLE_NSTP );
01117 _stat.inc_end_step();
01118 } else if( x[0] >= _maxt ) {
01119 particle->set_status( PARTICLE_TIME );
01120 _stat.inc_end_time();
01121 }
01122
01123
01124 _stat.inc_sum_steps( nstp );
01125
01126
01127 if( _trajdiv != 0 && pi % _trajdiv == 0 )
01128 particle->copy_trajectory( _traj );
01129
01130
01131 particle->x() = x;
01132
01133
01134 if( _tend_cb )
01135 (*_tend_cb)( particle, _pdb );
01136 }
01137
01138 };
01139
01140
01141 #endif
01142