IBSimu 1.0.4
|
00001 00005 /* Copyright (c) 2005-2010 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 "particles.hpp" 00055 #include "efield.hpp" 00056 #include "scalarfield.hpp" 00057 #include "scharge.hpp" 00058 #include "scheduler.hpp" 00059 #include "polysolver.hpp" 00060 00061 00062 //#define DEBUG_PARTICLE_ITERATOR 1 00063 00064 00067 enum particle_iterator_type_e { 00068 PARTICLE_ITERATOR_ADAPTIVE = 0, 00069 PARTICLE_ITERATOR_FIXED_STEP_LEN 00070 }; 00071 00072 00078 template <class PP> class ParticleIterator { 00079 00082 struct ColData { 00083 PP _x; 00084 int _dir; 00087 ColData( PP x, int dir ) : _x(x), _dir(dir) {} 00088 00089 bool operator<( const ColData &cd ) const { 00090 return( _x[0] < cd._x[0] ); 00091 } 00092 }; 00093 00094 gsl_odeiv_system _system; 00095 gsl_odeiv_step *_step; 00096 gsl_odeiv_control *_control; 00097 gsl_odeiv_evolve *_evolve; 00099 particle_iterator_type_e _type; 00101 bool _polyint; 00102 double _epsabs; 00103 double _epsrel; 00104 uint32_t _maxsteps; 00105 double _maxt; 00106 uint32_t _trajdiv; 00108 bool _mirror[6]; 00110 Particle<PP> *_first; 00111 ParticleIteratorData _pidata; 00113 PP _xi; 00115 std::vector<PP> _traj; 00116 std::vector<ColData> _coldata; 00118 uint32_t _end_time; 00119 uint32_t _end_step; 00120 uint32_t _end_out; 00122 uint32_t _end_coll; 00124 uint32_t _end_baddef; 00125 uint32_t _sum_steps; 00136 bool check_collision( Particle<PP> &particle, const PP &x1, const PP &x2, PP &status_x ) { 00137 00138 size_t a; 00139 double K; 00140 Vec3D v1, v2, vc; 00141 00142 // Convert PP to Vec3D 00143 for( a = 0; a < (PP::size()-1)/2; a++ ) { 00144 v1[a] = x1[2*a+1]; 00145 v2[a] = x2[2*a+1]; 00146 } 00147 00148 // If inside solid, bracket for collision point 00149 if( (a = _pidata._g->inside( v2 )) >= 7 ) { 00150 K = _pidata._g->bracket_surface( a, v2, v1, vc ); 00151 } else { 00152 return( true ); // No collision happened. 00153 } 00154 00155 // Convert Vec3D to PP 00156 for( a = 0; a < PP::size(); a++ ) 00157 status_x[a] = x2[a] + K*(x1[a]-x2[a]); 00158 00159 // Remove all points from _traj after time status_x[0]. 00160 for( a = _traj.size()-1; a > 0; a-- ) { 00161 if( _traj[a][0] > status_x[0] ) 00162 _traj.pop_back(); 00163 else 00164 break; 00165 } 00166 00167 // Save last trajectory point and update status 00168 _traj.push_back( status_x ); 00169 particle.set_status( PARTICLE_COLL ); 00170 _end_coll++; 00171 00172 return( false ); // Collision happened. 00173 } 00174 00175 00183 void handle_mirror( size_t c, int i[3], size_t a, int border, PP &x2 ) { 00184 00185 #ifdef DEBUG_PARTICLE_ITERATOR 00186 std::cout << " handle_mirror( c = " << c 00187 << ", i = (" << i[0] << ", " << i[1] << ", " << i[2] 00188 << "), a = " << a << ", border = " << border 00189 << ")\n"; 00190 #endif 00191 00192 double xmirror; 00193 if( border < 0 ) { 00194 xmirror = _pidata._g->origo(a); 00195 i[a] = -i[a]-1; 00196 } else { 00197 xmirror = _pidata._g->max(a); 00198 i[a] = 2*_pidata._g->size(a)-i[a]-3; 00199 } 00200 00201 #ifdef DEBUG_PARTICLE_ITERATOR 00202 std::cout << " xmirror = " << xmirror << "\n"; 00203 std::cout << " i = (" << i[0] << ", " << i[1] << ", " << i[2] << ")\n"; 00204 std::cout << " xi = " << _xi << "\n"; 00205 #endif 00206 00207 // Check if found edge at first encounter 00208 bool caught_at_boundary = false; 00209 if( _coldata[c]._dir == border*((int)a+1) && 00210 ( i[a] == 0 || i[a] == (int)_pidata._g->size(a)-2 ) ) { 00211 caught_at_boundary = true; 00212 #ifdef DEBUG_PARTICLE_ITERATOR 00213 std::cout << " caught_at_boundary\n"; 00214 #endif 00215 } 00216 00217 // Mirror traj back to _xi 00218 if( caught_at_boundary ) { 00219 _traj.push_back( _coldata[c]._x ); 00220 } else { 00221 for( int b = _traj.size()-1; b > 0; b-- ) { 00222 if( _traj[b][0] >= _xi[0] ) { 00223 00224 #ifdef DEBUG_PARTICLE_ITERATOR 00225 std::cout << " mirroring traj[" << b << "] = " << _traj[b] << "\n"; 00226 #endif 00227 _traj[b][2*a+1] = 2.0*xmirror - _traj[b][2*a+1]; 00228 _traj[b][2*a+2] *= -1.0; 00229 } else 00230 break; 00231 } 00232 } 00233 00234 // Mirror rest of the coldata 00235 for( size_t b = c; b < _coldata.size(); b++ ) { 00236 if( (size_t)abs(_coldata[b]._dir) == a+1 ) 00237 _coldata[b]._dir *= -1; 00238 _coldata[b]._x[2*a+1] = 2.0*xmirror - _coldata[b]._x[2*a+1]; 00239 _coldata[b]._x[2*a+2] *= -1.0; 00240 } 00241 00242 if( caught_at_boundary ) 00243 _traj.push_back( _coldata[c]._x ); 00244 00245 // Mirror calculation point 00246 x2[2*a+1] = 2.0*xmirror - x2[2*a+1]; 00247 x2[2*a+2] *= -1.0; 00248 00249 // Coordinates changed, reset integrator 00250 gsl_odeiv_step_reset( _step ); 00251 gsl_odeiv_evolve_reset( _evolve ); 00252 } 00253 00254 00255 void handle_collision( Particle<PP> &particle, size_t c, PP &status_x ) { 00256 00257 #ifdef DEBUG_PARTICLE_ITERATOR 00258 std::cout << " handle_collision()\n"; 00259 #endif 00260 00261 _traj.push_back( _coldata[c]._x ); 00262 status_x = _coldata[c]._x; 00263 particle.set_status( PARTICLE_OUT ); 00264 _end_out++; 00265 } 00266 00267 00276 bool handle_trajectory_advance( Particle<PP> &particle, size_t c, int i[3], PP &x2 ) { 00277 00278 // Check for collisions with solids and advance coordinates i. 00279 if( PP::dim() == 2 ) { 00280 if( _coldata[c]._dir == -1 ) { 00281 if( ( abs(_pidata._g->mesh(i[0], i[1] )) >= 7 || 00282 abs(_pidata._g->mesh(i[0], i[1]+1)) >= 7 ) && 00283 !check_collision( particle, _xi, _coldata[c]._x, x2 ) ) 00284 return( false ); 00285 i[0]--; 00286 } else if( _coldata[c]._dir == +1 ) { 00287 if( ( abs(_pidata._g->mesh(i[0]+1,i[1] )) >= 7 || 00288 abs(_pidata._g->mesh(i[0]+1,i[1]+1)) >= 7 ) && 00289 !check_collision( particle, _xi, _coldata[c]._x, x2 ) ) 00290 return( false ); 00291 i[0]++; 00292 } else if( _coldata[c]._dir == -2 ) { 00293 if( ( abs(_pidata._g->mesh(i[0], i[1] )) >= 7 || 00294 abs(_pidata._g->mesh(i[0]+1,i[1] )) >= 7 ) && 00295 !check_collision( particle, _xi, _coldata[c]._x, x2 ) ) 00296 return( false ); 00297 i[1]--; 00298 } else { 00299 if( ( abs(_pidata._g->mesh(i[0], i[1]+1)) >= 7 || 00300 abs(_pidata._g->mesh(i[0]+1,i[1]+1)) >= 7 ) && 00301 !check_collision( particle, _xi, _coldata[c]._x, x2 ) ) 00302 return( false ); 00303 i[1]++; 00304 } 00305 } else if( PP::dim() == 3 ) { 00306 if( _coldata[c]._dir == -1 ) { 00307 if( ( abs(_pidata._g->mesh(i[0], i[1], i[2] )) >= 7 || 00308 abs(_pidata._g->mesh(i[0], i[1]+1,i[2] )) >= 7 || 00309 abs(_pidata._g->mesh(i[0], i[1], i[2]+1)) >= 7 || 00310 abs(_pidata._g->mesh(i[0], i[1]+1,i[2]+1)) >= 7 ) && 00311 !check_collision( particle, _xi, _coldata[c]._x, x2 ) ) 00312 return( false ); 00313 i[0]--; 00314 } else if( _coldata[c]._dir == +1 ) { 00315 if( ( abs(_pidata._g->mesh(i[0]+1,i[1], i[2] )) >= 7 || 00316 abs(_pidata._g->mesh(i[0]+1,i[1]+1,i[2] )) >= 7 || 00317 abs(_pidata._g->mesh(i[0]+1,i[1], i[2]+1)) >= 7 || 00318 abs(_pidata._g->mesh(i[0]+1,i[1]+1,i[2]+1)) >= 7 ) && 00319 !check_collision( particle, _xi, _coldata[c]._x, x2 ) ) 00320 return( false ); 00321 i[0]++; 00322 } else if( _coldata[c]._dir == -2 ) { 00323 if( ( abs(_pidata._g->mesh(i[0], i[1],i[2] )) >= 7 || 00324 abs(_pidata._g->mesh(i[0]+1,i[1],i[2] )) >= 7 || 00325 abs(_pidata._g->mesh(i[0], i[1],i[2]+1)) >= 7 || 00326 abs(_pidata._g->mesh(i[0]+1,i[1],i[2]+1)) >= 7 ) && 00327 !check_collision( particle, _xi, _coldata[c]._x, x2 ) ) 00328 return( false ); 00329 i[1]--; 00330 } else if( _coldata[c]._dir == +2 ) { 00331 if( ( abs(_pidata._g->mesh(i[0], i[1]+1,i[2] )) >= 7 || 00332 abs(_pidata._g->mesh(i[0]+1,i[1]+1,i[2] )) >= 7 || 00333 abs(_pidata._g->mesh(i[0], i[1]+1,i[2]+1)) >= 7 || 00334 abs(_pidata._g->mesh(i[0]+1,i[1]+1,i[2]+1)) >= 7 ) && 00335 !check_collision( particle, _xi, _coldata[c]._x, x2 ) ) 00336 return( false ); 00337 i[1]++; 00338 } else if( _coldata[c]._dir == -3 ) { 00339 if( ( abs(_pidata._g->mesh(i[0], i[1], i[2] )) >= 7 || 00340 abs(_pidata._g->mesh(i[0]+1,i[1], i[2] )) >= 7 || 00341 abs(_pidata._g->mesh(i[0], i[1]+1,i[2])) >= 7 || 00342 abs(_pidata._g->mesh(i[0]+1,i[1]+1,i[2])) >= 7 ) && 00343 !check_collision( particle, _xi, _coldata[c]._x, x2 ) ) 00344 return( false ); 00345 i[2]--; 00346 } else { 00347 if( ( abs(_pidata._g->mesh(i[0], i[1], i[2]+1)) >= 7 || 00348 abs(_pidata._g->mesh(i[0]+1,i[1], i[2]+1)) >= 7 || 00349 abs(_pidata._g->mesh(i[0], i[1]+1,i[2]+1)) >= 7 || 00350 abs(_pidata._g->mesh(i[0]+1,i[1]+1,i[2]+1)) >= 7 ) && 00351 !check_collision( particle, _xi, _coldata[c]._x, x2 ) ) 00352 return( false ); 00353 i[2]++; 00354 } 00355 } else { 00356 throw( Error( ERROR_LOCATION, "unsupported dimension number" ) ); 00357 } 00358 00359 // Check for collisions/mirroring with simulation boundary. Here 00360 // coordinates i are already advanced to next mesh. 00361 for( size_t a = 0; a < PP::dim(); a++ ) { 00362 00363 if( i[a] < 0 ) { 00364 if( _mirror[2*a] ) 00365 handle_mirror( c, i, a, -1, x2 ); 00366 else { 00367 handle_collision( particle, c, x2 ); 00368 return( false ); 00369 } 00370 } else if( i[a] >= (_pidata._g->size(a)-1) ) { 00371 if( _mirror[2*a+1] ) 00372 handle_mirror( c, i, a, +1, x2 ); 00373 else { 00374 handle_collision( particle, c, x2 ); 00375 return( false ); 00376 } 00377 } 00378 } 00379 00380 return( true ); 00381 } 00382 00391 void build_coldata_linear( Particle<PP> &particle, const PP &x1, const PP &x2 ) { 00392 00393 int a1, a2; 00394 00395 for( size_t a = 0; a < PP::dim(); a++ ) { 00396 00397 a1 = (int)floor( (x1[2*a+1]-_pidata._g->origo(a))/_pidata._g->h() ); 00398 a2 = (int)floor( (x2[2*a+1]-_pidata._g->origo(a))/_pidata._g->h() ); 00399 if( a1 > a2 ) { 00400 int a = a2; 00401 a2 = a1; 00402 a1 = a; 00403 } 00404 00405 for( int b = a1+1; b <= a2; b++ ) { 00406 00407 // Save intersection coordinates 00408 double K = (b*_pidata._g->h() + _pidata._g->origo(a) - x1[2*a+1]) / 00409 (x2[2*a+1] - x1[2*a+1]); 00410 if( K < 0.0 ) K = 0.0; 00411 else if( K > 1.0 ) K = 1.0; 00412 //std::cout << "Found valid root: " << K << "\n"; 00413 00414 if( x2[2*a+1] > x1[2*a+1] ) 00415 _coldata.push_back( ColData( x1 + (x2-x1)*K, a+1 ) ); 00416 else 00417 _coldata.push_back( ColData( x1 + (x2-x1)*K, -a-1 ) ); 00418 } 00419 } 00420 00421 #ifdef DEBUG_PARTICLE_ITERATOR 00422 std::cout << " Coldata linear built\n"; 00423 #endif 00424 } 00425 00434 void build_coldata_poly( Particle<PP> &particle, const PP &x1, const PP &x2 ) { 00435 00436 #ifdef DEBUG_PARTICLE_ITERATOR 00437 std::cout << "Building coldata using polynomial interpolation\n"; 00438 #endif 00439 00440 // Construct trajectory representation 00441 TrajectoryRep1D traj[PP::dim()]; 00442 for( size_t a = 0; a < PP::dim(); a++ ) { 00443 traj[a].construct( x2[0]-x1[0], 00444 x1[2*a+1], x1[2*a+2], 00445 x2[2*a+1], x2[2*a+2] ); 00446 } 00447 00448 // Solve trajectory intersections 00449 for( size_t a = 0; a < PP::dim(); a++ ) { 00450 00451 // Mesh number of x1 (start point) 00452 int i = (int)floor( (x1[2*a+1]-_pidata._g->origo(a))/_pidata._g->h() ); 00453 00454 // Search to negative(dj = -1) and positive (dj = +1) mesh directions 00455 for( int dj = -1; dj <= 1; dj += 2 ) { 00456 int j = i; 00457 if( dj == +1 ) 00458 j = i+1; 00459 int Kcount; // Solution counter 00460 double K[3]; // Solution array 00461 while( 1 ) { 00462 double val = _pidata._g->origo(a) + _pidata._g->h() * j; 00463 #ifdef DEBUG_PARTICLE_ITERATOR 00464 std::cout << " Searching intersections at coord(" << a << ") = " << val << "\n"; 00465 #endif 00466 Kcount = traj[a].solve( K, val ); 00467 if( Kcount == 0 ) 00468 break; // No valid roots 00469 00470 #ifdef DEBUG_PARTICLE_ITERATOR 00471 std::cout << " Found " << Kcount << " valid roots: "; 00472 for( int p = 0; p < Kcount; p++ ) 00473 std::cout << K[p] << " "; 00474 std::cout << "\n"; 00475 #endif 00476 00477 // Save roots to coldata 00478 for( int b = 0; b < Kcount; b++ ) { 00479 PP xcol; 00480 double x, v; 00481 xcol(0) = x1[0] + K[b]*(x2[0]-x1[0]); 00482 for( size_t c = 0; c < PP::dim(); c++ ) { 00483 traj[c].coord( x, v, K[b] ); 00484 if( a == c ) 00485 xcol[2*c+1] = val; // limit numerical inaccuracy 00486 else 00487 xcol[2*c+1] = x; 00488 xcol[2*c+2] = v; 00489 } 00490 if( _pidata._g->geom_mode() == MODE_CYL ) 00491 xcol[5] = x1[5] + K[b]*(x2[5]-x1[5]); 00492 if( xcol[2*a+2] >= 0.0 ) 00493 _coldata.push_back( ColData( xcol, a+1 ) ); 00494 else 00495 _coldata.push_back( ColData( xcol, -a-1 ) ); 00496 } 00497 00498 j += dj; 00499 } 00500 } 00501 } 00502 00503 #ifdef DEBUG_PARTICLE_ITERATOR 00504 std::cout << " Coldata built\n"; 00505 #endif 00506 } 00507 00513 bool limit_trajectory_advance( const PP &x1, PP &x2 ) { 00514 00515 bool touched = false; 00516 00517 for( size_t a = 0; a < PP::dim(); a++ ) { 00518 00519 double lim1 = _pidata._g->origo(a) - 00520 (_pidata._g->size(a)-1)*_pidata._g->h(); 00521 double lim2 = _pidata._g->origo(a) + 00522 2*(_pidata._g->size(a)-1)*_pidata._g->h(); 00523 00524 if( x2[2*a+1] < lim1 ) { 00525 00526 double K = (lim1 - x1[2*a+1]) / (x2[2*a+1] - x1[2*a+1]); 00527 x2 = x1 + K*(x2-x1); 00528 touched = true; 00529 #ifdef DEBUG_PARTICLE_ITERATOR 00530 std::cout << "Limiting step to:\n"; 00531 std::cout << " x2: " << x2 << "\n"; 00532 #endif 00533 } else if(x2[2*a+1] > lim2 ) { 00534 00535 double K = (lim2 - x1[2*a+1]) / (x2[2*a+1] - x1[2*a+1]); 00536 x2 = x1 + K*(x2-x1); 00537 touched = true; 00538 #ifdef DEBUG_PARTICLE_ITERATOR 00539 std::cout << "Limiting step to:\n"; 00540 std::cout << " x2: " << x2 << "\n"; 00541 #endif 00542 } 00543 } 00544 00545 return( touched ); 00546 } 00547 00563 bool handle_trajectory( Particle<PP> &particle, const PP &x1, PP &x2, 00564 bool force_linear=false ) { 00565 00566 #ifdef DEBUG_PARTICLE_ITERATOR 00567 std::cout << "Handle trajectory from x1 to x2:\n"; 00568 std::cout << " x1: " << x1 << "\n"; 00569 std::cout << " x2: " << x2 << "\n"; 00570 #endif 00571 00572 // Limit trajectory advance to double the simulation box 00573 // If limitation done, force to linear interpolation 00574 if( limit_trajectory_advance( x1, x2 ) ) 00575 force_linear = true; 00576 00577 // Make coldata 00578 if( _polyint && !force_linear ) 00579 build_coldata_poly( particle, x1, x2 ); 00580 else 00581 build_coldata_linear( particle, x1, x2 ); 00582 00583 // No intersections, nothing to do 00584 if( _coldata.size() == 0 ) { 00585 #ifdef DEBUG_PARTICLE_ITERATOR 00586 std::cout << "No coldata\n"; 00587 #endif 00588 return( true ); 00589 } 00590 00591 // Sort intersections in increasing time order 00592 sort( _coldata.begin(), _coldata.end() ); 00593 00594 // Starting mesh index 00595 int i[3] = {0, 0, 0}; 00596 for( size_t cdir = 0; cdir < PP::dim(); cdir++ ) 00597 i[cdir] = (int)floor( (x1[2*cdir+1]-_pidata._g->origo(cdir))/_pidata._g->h() ); 00598 00599 // Process intersection points 00600 #ifdef DEBUG_PARTICLE_ITERATOR 00601 std::cout << "Process coldata points:\n"; 00602 #endif 00603 for( size_t a = 0; a < _coldata.size(); a++ ) { 00604 00605 #ifdef DEBUG_PARTICLE_ITERATOR 00606 std::cout << " Coldata " << std::setw(4) << a << ": " 00607 << _coldata[a]._x << ", " 00608 << std::setw(3) << i[0] << " " 00609 << std::setw(3) << i[1] << " " 00610 << std::setw(3) << i[2] << " " 00611 << std::setw(3) << _coldata[a]._dir << "\n"; 00612 #endif 00613 00614 // Advance particle in mesh, check for possible collisions and 00615 // do mirroring. 00616 handle_trajectory_advance( particle, a, i, x2 ); 00617 00618 // Update space charge for one mesh. 00619 if( _pidata._scharge ) 00620 scharge_add_from_trajectory( *_pidata._scharge, particle.IQ(), 00621 _xi, _coldata[a]._x ); 00622 00623 #ifdef DEBUG_PARTICLE_ITERATOR 00624 if( particle.get_status() == PARTICLE_OUT ) { 00625 std::cout << " Particle out\n"; 00626 std::cout << " x = " << x2 << "\n"; 00627 } else if( particle.get_status() == PARTICLE_COLL ) { 00628 std::cout << " Particle collided\n"; 00629 std::cout << " x = " << x2 << "\n"; 00630 } 00631 #endif 00632 // Clear coldata and exit if particle collided. 00633 if( particle.get_status() != PARTICLE_OK ) { 00634 _coldata.clear(); 00635 return( false ); 00636 } 00637 00638 // Update last intersection point xi. 00639 _xi = _coldata[a]._x; 00640 } 00641 00642 #ifdef DEBUG_PARTICLE_ITERATOR 00643 std::cout << "Coldata done\n"; 00644 #endif 00645 _coldata.clear(); 00646 return( true ); 00647 } 00648 00649 00652 bool axis_mirror_required( const PP &x2 ) { 00653 return( _pidata._g->geom_mode() == MODE_CYL && 00654 x2[4] < 0.0 && 00655 x2[3] <= 0.01*_pidata._g->h() && 00656 x2[3]*fabs(x2[5]) <= 1.0e-9*fabs(x2[4]) ); 00657 00658 } 00659 00660 00666 bool handle_axis_mirror_step( Particle<PP> &particle, const PP &x1, PP &x2 ) { 00667 00668 // Get acceleration at x2 00669 double dxdt[5]; 00670 PP::get_derivatives( x2[0], &x2[1], dxdt, (void *)&_pidata ); 00671 00672 // Calculate crossover point assuming zero acceleration in 00673 // r-direction and constant acceleration in x-direction 00674 double dt = -x2[3]/x2[4]; 00675 PP xc; 00676 xc[0] = x2[0]+dt; 00677 xc[1] = x2[1]+(x2[2]+0.5*dxdt[1]*dt)*dt; 00678 xc[2] = x2[2]; 00679 xc[3] = x2[3]+x2[4]*dt; 00680 xc[4] = x2[4]; 00681 xc[5] = x2[5]; 00682 00683 // Mirror x2 to x3 00684 PP x3 = 2*xc - x2; 00685 x3[3] *= -1.0; 00686 x3[4] *= -1.0; 00687 x3[5] *= -1.0; 00688 00689 #ifdef DEBUG_PARTICLE_ITERATOR 00690 std::cout << "Particle mirror:\n"; 00691 std::cout << " x1: " << x1 << "\n"; 00692 std::cout << " x2: " << x2 << "\n"; 00693 std::cout << " xc: " << xc << "\n"; 00694 std::cout << " x3: " << x3 << "\n"; 00695 #endif 00696 00697 // Handle step with linear interpolation to avoid going to r<=0 00698 if( !handle_trajectory( particle, x2, x3, true ) ) 00699 return( false ); // Particle done 00700 00701 // Save trajectory calculation points 00702 _traj.push_back( x2 ); 00703 _traj.push_back( xc ); 00704 xc[4] *= -1.0; 00705 xc[5] *= -1.0; 00706 _traj.push_back( xc ); 00707 00708 // Next step not a continuation of previous one, reset 00709 // integrator 00710 gsl_odeiv_step_reset( _step ); 00711 gsl_odeiv_evolve_reset( _evolve ); 00712 00713 // Continue iteration at mirrored point 00714 x2 = x3; 00715 return( true ); 00716 } 00717 00728 bool check_particle_definition( PP &x ) { 00729 00730 #ifdef DEBUG_PARTICLE_ITERATOR 00731 std::cout << "Particle defined at:\n"; 00732 std::cout << " x = " << x << "\n"; 00733 #endif 00734 00735 // Check if inside solids of outside geometry. 00736 if( _pidata._g->inside( x.location() ) ) 00737 return( false ); 00738 00739 // Check if particle on simulation geometry border and directed outwards 00740 /* 00741 for( size_t a = 0; a < PP::dim(); a++ ) { 00742 if( x[2*a+1] == _pidata._g->origo(a) && x[2*a+2] < 0.0 ) { 00743 if( _mirror[2*a] ) { 00744 x[2*a+2] *= -1.0; 00745 #ifdef DEBUG_PARTICLE_ITERATOR 00746 std::cout << "Mirroring to:\n"; 00747 std::cout << " x = " << x << "\n"; 00748 #endif 00749 } else { 00750 return( false ); 00751 } 00752 00753 } else if( x[2*a+1] == _pidata._g->max(a) & x[2*a+2] > 0.0 ) { 00754 if( _mirror[2*a+1] ) { 00755 x[2*a+2] *= -1.0; 00756 #ifdef DEBUG_PARTICLE_ITERATOR 00757 std::cout << "Mirroring to:\n"; 00758 std::cout << " x = " << x << "\n"; 00759 #endif 00760 } else { 00761 return( false ); 00762 } 00763 } 00764 } 00765 00766 #ifdef DEBUG_PARTICLE_ITERATOR 00767 std::cout << "Definition ok\n"; 00768 #endif 00769 00770 */ 00771 return( true ); 00772 } 00773 00774 double calculate_dt( const PP &x, const double *dxdt ) { 00775 00776 double spd = 0.0, acc = 0.0; 00777 00778 for( size_t a = 0; a < (PP::size()-1)/2; a++ ) { 00779 //std::cout << "spd += " << dxdt[2*a]*dxdt[2*a] << "\n"; 00780 spd += dxdt[2*a]*dxdt[2*a]; 00781 //std::cout << "acc += " << dxdt[2*a+1]*dxdt[2*a+1] << "\n"; 00782 acc += dxdt[2*a+1]*dxdt[2*a+1]; 00783 } 00784 if( _pidata._g->geom_mode() == MODE_CYL ) { 00785 //std::cout << "MODE_CYL\n"; 00786 //std::cout << "spd += " << x[3]*x[3]*x[5]*x[5] << "\n"; 00787 spd += x[3]*x[3]*x[5]*x[5]; 00788 //std::cout << "acc += " << x[3]*x[3]*dxdt[4]*dxdt[4] << "\n"; 00789 acc += x[3]*x[3]*dxdt[4]*dxdt[4]; 00790 } 00791 //std::cout << "spd = " << sqrt(spd) << "\n"; 00792 //std::cout << "acc = " << sqrt(acc) << "\n"; 00793 spd = _pidata._g->h() / sqrt(spd); 00794 acc = sqrt( 2.0*_pidata._g->h() / sqrt(acc) ); 00795 00796 return( spd < acc ? spd : acc ); 00797 } 00798 00799 public: 00800 00827 ParticleIterator( particle_iterator_type_e type, double epsabs, double epsrel, 00828 bool polyint, uint32_t maxsteps, double maxt, 00829 uint32_t trajdiv, bool mirror[6], ScalarField *scharge, const Efield *efield, 00830 const VectorField *bfield, const Geometry *g, Particle<PP> *first ) 00831 : _type(type), _polyint(polyint), _epsabs(epsabs), _epsrel(epsrel), _maxsteps(maxsteps), _maxt(maxt), 00832 _trajdiv(trajdiv), _first(first), _pidata(scharge,efield,bfield,g), _end_time(0), 00833 _end_step(0), _end_out(0), _end_coll(0), _end_baddef(0), _sum_steps(0) { 00834 00835 // Initialize mirroring 00836 _mirror[0] = mirror[0]; 00837 _mirror[1] = mirror[1]; 00838 _mirror[2] = mirror[2]; 00839 _mirror[3] = mirror[3]; 00840 _mirror[4] = mirror[4]; 00841 _mirror[5] = mirror[5]; 00842 00843 // Initialize system of ordinary differential equations (ODE) 00844 _system.jacobian = NULL; 00845 _system.params = (void *)&_pidata; 00846 _system.function = PP::get_derivatives; 00847 _system.dimension = PP::size()-1; // Time is not part of differential equation dimensions 00848 00849 // Make scale 00850 // 2D: x vx y vy 00851 // Cyl: x vx r vr omega 00852 // 3D: x vx y vy z vz 00853 double scale_abs[PP::size()-1]; 00854 for( uint32_t a = 0; a < (uint32_t)PP::size()-2; a+=2 ) { 00855 scale_abs[a+0] = 1.0; 00856 scale_abs[a+1] = 1.0e6; 00857 } 00858 if( _pidata._g->geom_mode() == MODE_CYL ) 00859 scale_abs[4] = 1.0; 00860 00861 // Initialize ODE solver 00862 _step = gsl_odeiv_step_alloc( gsl_odeiv_step_rkck, _system.dimension ); 00863 //_control = gsl_odeiv_control_standard_new( _epsabs, _epsrel, 1.0, 1.0 ); 00864 _control = gsl_odeiv_control_scaled_new( _epsabs, _epsrel, 1.0, 1.0, scale_abs, PP::size()-1 ); 00865 _evolve = gsl_odeiv_evolve_alloc( _system.dimension ); 00866 } 00867 00868 00871 ~ParticleIterator() { 00872 gsl_odeiv_evolve_free( _evolve ); 00873 gsl_odeiv_control_free( _control ); 00874 gsl_odeiv_step_free( _step ); 00875 } 00876 00877 00880 void enable_nsimp_plasma_threshold( const ScalarField *epot, double phi_plasma ) { 00881 _pidata._epot = epot; 00882 _pidata._phi_plasma = phi_plasma; 00883 } 00884 00885 00893 void get_statistics( uint32_t &end_time, uint32_t &end_step, uint32_t &end_out, 00894 uint32_t &end_coll, uint32_t &end_baddef, uint32_t &sum_steps ) { 00895 end_time = _end_time; 00896 end_step = _end_step; 00897 end_out = _end_out; 00898 end_coll = _end_coll; 00899 end_baddef = _end_baddef; 00900 sum_steps = _sum_steps; 00901 } 00902 00911 void operator()( Particle<PP> *particle, 00912 Scheduler<ParticleIterator<PP>,Particle<PP>,Error> &scheduler ) { 00913 00914 // Copy starting point to x and 00915 PP x = particle->x(); 00916 00917 // Check particle definition 00918 if( !check_particle_definition( x ) ) { 00919 particle->set_status( PARTICLE_BADDEF ); 00920 _end_baddef++; 00921 return; 00922 } 00923 particle->x() = x; 00924 00925 // Reset trajectory and save first trajectory point. 00926 _traj.clear(); 00927 _traj.push_back( x ); 00928 #ifdef DEBUG_PARTICLE_ITERATOR 00929 std::cout << x[0] << " " 00930 << x[1] << " " 00931 << x[2] << " " 00932 << x[3] << " " 00933 << x[4] << "\n"; 00934 #endif 00935 _pidata._qm = particle->qm(); 00936 _xi = x; 00937 00938 // Reset integrator 00939 gsl_odeiv_step_reset( _step ); 00940 gsl_odeiv_evolve_reset( _evolve ); 00941 00942 // Make initial guess for step size 00943 //std::cout << "Guess dt ------------------------------------------------\n"; 00944 double dxdt[PP::size()-1]; 00945 PP::get_derivatives( 0.0, &x[1], dxdt, (void *)&_pidata ); 00946 double dt = calculate_dt( x, dxdt ); 00947 00948 #ifdef DEBUG_PARTICLE_ITERATOR 00949 std::cout << "dxdt = "; 00950 for( size_t a = 0; a < PP::size()-1; a++ ) 00951 std::cout << dxdt[a] << " "; 00952 std::cout << "\n"; 00953 std::cout << "dt = " << dt << "\n"; 00954 std::cout << "*** Starting iteration\n"; 00955 #endif 00956 00957 // Iterate ODEs until maximum steps are done, time is used 00958 // or particle collides. 00959 PP x2; 00960 size_t nstp = 0; 00961 while( nstp < _maxsteps && x[0] < _maxt ) { 00962 00963 #ifdef DEBUG_PARTICLE_ITERATOR 00964 std::cout << "\n*** Step ***\n"; 00965 std::cout << " x = " << x2 << "\n"; 00966 std::cout << " dt = " << dt << " (proposed)\n"; 00967 #endif 00968 00969 // Take a step. 00970 x2 = x; 00971 00972 while( nstp < _maxsteps ) { 00973 int retval = gsl_odeiv_evolve_apply( _evolve, _control, _step, &_system, 00974 &x2[0], _maxt, &dt, &x2[1] ); 00975 if( retval == IBSIMU_DERIV_ERROR ) { 00976 #ifdef DEBUG_PARTICLE_ITERATOR 00977 std::cout << "Step rejected\n"; 00978 std::cout << " x2 = " << x2 << "\n"; 00979 std::cout << " dt = " << dt << "\n"; 00980 #endif 00981 x2[0] = x[0]; // Reset time (this shouldn't be necessary - there 00982 // is a bug in GSL-1.12, report has been sent) 00983 dt *= 0.5; 00984 nstp++; 00985 continue; 00986 } else if( retval == GSL_SUCCESS ) { 00987 break; 00988 } else { 00989 throw( Error( ERROR_LOCATION, "gsl_odeiv_evolve_apply failed" ) ); 00990 } 00991 } 00992 00993 // Check step count number and step size validity 00994 if( nstp >= _maxsteps ) 00995 break; 00996 if( x2[0] == x[0] ) 00997 throw( Error( ERROR_LOCATION, "too small step size" ) ); 00998 00999 #ifdef DEBUG_PARTICLE_ITERATOR 01000 std::cout << "Step accepted from x1 to x2:\n"; 01001 std::cout << " dt = " << dt << " (taken)\n"; 01002 std::cout << " x1 = " << x << "\n"; 01003 std::cout << " x2 = " << x2 << "\n"; 01004 #endif 01005 01006 // Handle collisions and space charge of step. 01007 if( !handle_trajectory( *particle, x, x2 ) ) { 01008 x = x2; 01009 break; // Particle done 01010 } 01011 01012 // Check if particle mirroring is required to avoid 01013 // singularity at symmetry axis. 01014 if( axis_mirror_required( x2 ) ) { 01015 if( !handle_axis_mirror_step( *particle, x, x2 ) ) 01016 break; // Particle done 01017 } 01018 01019 // Propagate coordinates 01020 x = x2; 01021 01022 // Save trajectory point 01023 _traj.push_back( x2 ); 01024 01025 // Increase step count. 01026 nstp++; 01027 } 01028 01029 #ifdef DEBUG_PARTICLE_ITERATOR 01030 std::cout << "\n*** Done stepping ***\n"; 01031 #endif 01032 01033 // Check if step count or time limited 01034 if( nstp == _maxsteps ) { 01035 particle->set_status( PARTICLE_NSTP ); 01036 _end_step++; 01037 } else if( x[0] >= _maxt ) { 01038 particle->set_status( PARTICLE_TIME ); 01039 _end_time++; 01040 } 01041 01042 // Save step count 01043 _sum_steps += nstp; 01044 01045 // Save trajectory of current particle 01046 if( _trajdiv != 0 && (particle-_first) % _trajdiv == 0 ) 01047 particle->copy_trajectory( _traj ); 01048 01049 // Save last particle location 01050 particle->x() = x; 01051 } 01052 01053 /* 01054 std::cout << "Kala\n"; 01055 01056 std::cout << _coldata.capacity() << "\n"; 01057 std::cout << _traj.capacity() << "\n"; 01058 */ 01059 }; 01060 01061 01062 01063 01064 #endif 01065 01066 01067 01068 01069 01070 01071 01072 01073 01074 01075 01076 01077 01078 01079