particleiterator.hpp
Go to the documentation of this file.
00001 
00005 /* Copyright (c) 2005-2011 Taneli Kalvas. All rights reserved.
00006  *
00007  * You can redistribute this software and/or modify it under the terms
00008  * of the GNU General Public License as published by the Free Software
00009  * Foundation; either version 2 of the License, or (at your option)
00010  * any later version.
00011  * 
00012  * This library is distributed in the hope that it will be useful, but
00013  * WITHOUT ANY WARRANTY; without even the implied warranty of
00014  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
00015  * General Public License for more details.
00016  * 
00017  * You should have received a copy of the GNU General Public License
00018  * along with this library (file "COPYING" included in the package);
00019  * if not, write to the Free Software Foundation, Inc., 51 Franklin
00020  * Street, Fifth Floor, Boston, MA 02110-1301 USA
00021  * 
00022  * If you have questions about your rights to use or distribute this
00023  * software, please contact Berkeley Lab's Technology Transfer
00024  * Department at TTD@lbl.gov. Other questions, comments and bug
00025  * reports should be sent directly to the author via email at
00026  * taneli.kalvas@jyu.fi.
00027  * 
00028  * NOTICE. This software was developed under partial funding from the
00029  * U.S.  Department of Energy.  As such, the U.S. Government has been
00030  * granted for itself and others acting on its behalf a paid-up,
00031  * nonexclusive, irrevocable, worldwide license in the Software to
00032  * reproduce, prepare derivative works, and perform publicly and
00033  * display publicly.  Beginning five (5) years after the date
00034  * permission to assert copyright is obtained from the U.S. Department
00035  * of Energy, and subject to any subsequent five (5) year renewals,
00036  * the U.S. Government is granted for itself and others acting on its
00037  * behalf a paid-up, nonexclusive, irrevocable, worldwide license in
00038  * the Software to reproduce, prepare derivative works, distribute
00039  * copies to the public, perform publicly and display publicly, and to
00040  * permit others to do so.
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 //#define DEBUG_PARTICLE_ITERATOR 1
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                 // Save intersection coordinates
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                 //std::cout << "Found valid root: " << K << "\n";
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         // Sort intersections in increasing time order
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         // Construct trajectory representation
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         // Solve trajectory intersections
00171         for( size_t a = 0; a < PP::dim(); a++ ) {
00172 
00173             // Mesh number of x1 (start point)
00174             int i = (int)floor( (x1[2*a+1]-mesh.origo(a))/mesh.h() );
00175             
00176             // Search to negative (dj = -1) and positive (dj = +1) mesh directions
00177             for( int dj = -1; dj <= 1; dj += 2 ) {
00178                 int j = i;
00179                 if( dj == +1 )
00180                     j = i+1;
00181                 int Kcount;  // Solution counter
00182                 double K[3]; // Solution array
00183                 while( 1 ) {
00184 
00185                     // Intersection point
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; // No valid roots
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                     // Save roots to coldata
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; // limit numerical inaccuracy
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         // Sort intersections in increasing time order
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         // If inside solid, bracket for collision point
00308         Vec3D v2 = x2.location();
00309         int32_t bound = _pidata._geom->inside( v2 );
00310         if( bound < 7 )
00311             return( true ); // No collision happened.
00312         Vec3D vc;
00313         Vec3D v1 = x1.location();
00314         double K = _pidata._geom->bracket_surface( bound, v2, v1, vc );
00315 
00316         // Calculate new PP
00317         for( size_t a = 0; a < PP::size(); a++ )
00318             status_x[a] = x2[a] + K*(x1[a]-x2[a]);
00319 
00320         // Remove all points from trajectory after time status_x[0].
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         // Save last trajectory point and update status
00329         //save_trajectory_point( status_x );
00330         particle.set_status( PARTICLE_COLL );
00331 
00332         // Update collision statistics for boundary
00333         _stat.add_bound_collision( bound, particle.IQ() );
00334 
00335         return( false ); // Collision happened.
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         // Check if found edge at first encounter
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         // Mirror traj back to _xi
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         // Mirror rest of the coldata
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         // Mirror calculation point
00409         x2[2*a+1] = 2.0*xmirror - x2[2*a+1];
00410         x2[2*a+2] *= -1.0;
00411         
00412         // Coordinates changed, reset integrator
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         //save_trajectory_point( _coldata[c]._x );
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         // Check for collisions with solids and advance coordinates i.
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         // Check for collisions/mirroring with simulation boundary. Here
00523         // coordinates i are already advanced to next mesh.
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         // Limit trajectory advance to double the simulation box
00632         // If limitation done, force to linear interpolation
00633         if( limit_trajectory_advance( x1, x2 ) )
00634             force_linear = true;
00635 
00636         build_coldata( force_linear, x1, x2 );
00637 
00638         // TODO
00639         // Remove entrance to geometry if coming from outside or make
00640         // code to skip the collision detection for these particles
00641 
00642         // No intersections, nothing to do
00643         if( _coldata.size() == 0 ) {
00644 #ifdef DEBUG_PARTICLE_ITERATOR
00645             std::cout << "No coldata\n";
00646 #endif
00647             return( true );
00648         }
00649 
00650         // Starting mesh index
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         // Process intersection points
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             // Advance particle in mesh, check for possible collisions and
00671             // do mirroring.
00672             handle_trajectory_advance( particle, a, i, x2 );
00673 
00674             // Update space charge for one mesh.
00675             if( _pidata._scharge )
00676                 scharge_add_from_trajectory( *_pidata._scharge, _scharge_mutex, particle.IQ(), 
00677                                              _xi, _coldata[a]._x );
00678 
00679             // Call trajectory handler callback
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             // Clear coldata and exit if particle collided.
00693             if( particle.get_status() != PARTICLE_OK ) {
00694                 save_trajectory_point( x2 );
00695                 _coldata.clear();
00696                 return( false );
00697             }
00698 
00699             // Update last intersection point xi.
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         // Get acceleration at x2
00730         double dxdt[5];
00731         PP::get_derivatives( x2[0], &x2[1], dxdt, (void *)&_pidata );
00732 
00733         // Calculate crossover point assuming zero acceleration in
00734         // r-direction and constant acceleration in x-direction
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         // Mirror x2 to x3
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         // Handle step with linear interpolation to avoid going to r<=0
00759         if( !handle_trajectory( particle, x2, x3, true, false ) )
00760             return( false ); // Particle done
00761 
00762         // Save trajectory calculation points
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         // Next step not a continuation of previous one, reset
00770         // integrator
00771         gsl_odeiv_step_reset( _step );
00772         gsl_odeiv_evolve_reset( _evolve );
00773 
00774         // Continue iteration at mirrored point
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         // Check for NaN
00800         for( size_t a = 0; a < PP::size(); a++ ) {
00801             if( comp_isnan( x[a] ) )
00802                 return( false );
00803         }
00804 
00805         // Check if inside solids or outside simulation geometry box.
00806         if( _pidata._geom->inside( x.location() ) )
00807             return( false );
00808 
00809         // Check if particle on simulation geometry border and directed outwards
00810         /*
00811         for( size_t a = 0; a < PP::dim(); a++ ) {
00812             if( x[2*a+1] == _pidata._geom->origo(a) && x[2*a+2] < 0.0 ) {
00813                 if( _mirror[2*a] ) {
00814                     x[2*a+2] *= -1.0;
00815 #ifdef DEBUG_PARTICLE_ITERATOR
00816         std::cout << "Mirroring to:\n";
00817         std::cout << "  x = " << x << "\n";
00818 #endif
00819                 } else {
00820                     return( false );
00821                 }
00822 
00823             } else if( x[2*a+1] == _pidata._geom->max(a) & x[2*a+2] > 0.0 ) {
00824                 if( _mirror[2*a+1] ) {
00825                     x[2*a+2] *= -1.0;
00826 #ifdef DEBUG_PARTICLE_ITERATOR
00827         std::cout << "Mirroring to:\n";
00828         std::cout << "  x = " << x << "\n";
00829 #endif
00830                 } else {
00831                     return( false );
00832                 }
00833             }
00834         }
00835 
00836 #ifdef DEBUG_PARTICLE_ITERATOR
00837         std::cout << "Definition ok\n";
00838 #endif
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             //std::cout << "spd += " << dxdt[2*a]*dxdt[2*a] << "\n";
00850             spd += dxdt[2*a]*dxdt[2*a];
00851             //std::cout << "acc += " << dxdt[2*a+1]*dxdt[2*a+1] << "\n";
00852             acc += dxdt[2*a+1]*dxdt[2*a+1];
00853         }
00854         if( _pidata._geom->geom_mode() == MODE_CYL ) {
00855             //std::cout << "MODE_CYL\n";
00856             //std::cout << "spd += " << x[3]*x[3]*x[5]*x[5] << "\n";
00857             spd += x[3]*x[3]*x[5]*x[5];
00858             //std::cout << "acc += " << x[3]*x[3]*dxdt[4]*dxdt[4] << "\n";
00859             acc += x[3]*x[3]*dxdt[4]*dxdt[4];
00860         }
00861         //std::cout << "spd = " << sqrt(spd) << "\n";
00862         //std::cout << "acc = " << sqrt(acc) << "\n";
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         // Initialize mirroring
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         // Initialize system of ordinary differential equations (ODE)
00917         _system.jacobian  = NULL;
00918         _system.params    = (void *)&_pidata;
00919         _system.function  = PP::get_derivatives;
00920         _system.dimension = PP::size()-1; // Time is not part of differential equation dimensions
00921 
00922         // Make scale
00923         // 2D:  x vx y vy
00924         // Cyl: x vx r vr omega
00925         // 3D:  x vx y vy z vz
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         // Initialize ODE solver
00935         _step    = gsl_odeiv_step_alloc( gsl_odeiv_step_rkck, _system.dimension );
00936         //_control = gsl_odeiv_control_standard_new( _epsabs, _epsrel, 1.0, 1.0 );
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         // Check particle status
00990         if( particle->get_status() != PARTICLE_OK )
00991             return;
00992 
00993         // Copy starting point to x and 
00994         PP x = particle->x();
00995 
00996         // Check particle definition
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         // Reset trajectory and save first trajectory point.
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         // Reset integrator
01018         gsl_odeiv_step_reset( _step );
01019         gsl_odeiv_evolve_reset( _evolve );
01020         
01021         // Make initial guess for step size
01022         //std::cout << "Guess dt ------------------------------------------------\n";
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         // Iterate ODEs until maximum steps are done, time is used 
01037         // or particle collides.
01038         PP x2;
01039         size_t nstp = 0; // Steps taken
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             // Take a step.
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]; // Reset time (this shouldn't be necessary - there 
01061                                   // is a bug in GSL-1.12, report has been sent)
01062                     dt *= 0.5;
01063                     if( dt == 0.0 )
01064                         throw( Error( ERROR_LOCATION, "too small step size" ) );
01065                     //nstp++;
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             // Check step count number and step size validity
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             // Handle collisions and space charge of step.
01088             if( !handle_trajectory( *particle, x, x2, false, nstp == 0 ) ) {
01089                 x = x2;
01090                 break; // Particle done
01091             }
01092 
01093             // Check if particle mirroring is required to avoid 
01094             // singularity at symmetry axis.
01095             if( axis_mirror_required( x2 ) ) {
01096                 if( !handle_axis_mirror_step( *particle, x, x2 ) )
01097                     break; // Particle done
01098             }
01099 
01100             // Propagate coordinates
01101             x = x2;
01102 
01103             // Save trajectory point
01104             save_trajectory_point( x2 );
01105             
01106             // Increase step count.
01107             nstp++;
01108         }
01109 
01110 #ifdef DEBUG_PARTICLE_ITERATOR
01111         std::cout << "\n*** Done stepping ***\n";
01112 #endif
01113 
01114         // Check if step count or time limited 
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         // Save step count
01124         _stat.inc_sum_steps( nstp );
01125 
01126         // Save trajectory of current particle
01127         if( _trajdiv != 0 && pi % _trajdiv == 0 )
01128             particle->copy_trajectory( _traj );
01129 
01130         // Save last particle location
01131         particle->x() = x;
01132 
01133         // Call trajectory end callback
01134         if( _tend_cb )
01135             (*_tend_cb)( particle, _pdb );
01136     }
01137 
01138 };
01139 
01140 
01141 #endif
01142