particleiterator.hpp
Go to the documentation of this file.
1 
5 /* Copyright (c) 2005-2011 Taneli Kalvas. All rights reserved.
6  *
7  * You can redistribute this software and/or modify it under the terms
8  * of the GNU General Public License as published by the Free Software
9  * Foundation; either version 2 of the License, or (at your option)
10  * any later version.
11  *
12  * This library is distributed in the hope that it will be useful, but
13  * WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15  * General Public License for more details.
16  *
17  * You should have received a copy of the GNU General Public License
18  * along with this library (file "COPYING" included in the package);
19  * if not, write to the Free Software Foundation, Inc., 51 Franklin
20  * Street, Fifth Floor, Boston, MA 02110-1301 USA
21  *
22  * If you have questions about your rights to use or distribute this
23  * software, please contact Berkeley Lab's Technology Transfer
24  * Department at TTD@lbl.gov. Other questions, comments and bug
25  * reports should be sent directly to the author via email at
26  * taneli.kalvas@jyu.fi.
27  *
28  * NOTICE. This software was developed under partial funding from the
29  * U.S. Department of Energy. As such, the U.S. Government has been
30  * granted for itself and others acting on its behalf a paid-up,
31  * nonexclusive, irrevocable, worldwide license in the Software to
32  * reproduce, prepare derivative works, and perform publicly and
33  * display publicly. Beginning five (5) years after the date
34  * permission to assert copyright is obtained from the U.S. Department
35  * of Energy, and subject to any subsequent five (5) year renewals,
36  * the U.S. Government is granted for itself and others acting on its
37  * behalf a paid-up, nonexclusive, irrevocable, worldwide license in
38  * the Software to reproduce, prepare derivative works, distribute
39  * copies to the public, perform publicly and display publicly, and to
40  * permit others to do so.
41  */
42 
43 #ifndef PARTICLEITERATOR_HPP
44 #define PARTICLEITERATOR_HPP 1
45 
46 
47 #include <vector>
48 #include <iostream>
49 #include <algorithm>
50 #include <iomanip>
51 #include <gsl/gsl_odeiv.h>
52 #include <gsl/gsl_poly.h>
53 #include "geometry.hpp"
54 #include "compmath.hpp"
55 #include "trajectory.hpp"
56 #include "particles.hpp"
57 #include "vectorfield.hpp"
58 #include "scalarfield.hpp"
59 #include "scharge.hpp"
60 #include "scheduler.hpp"
61 #include "polysolver.hpp"
62 #include "particledatabase.hpp"
63 
64 
65 //#define DEBUG_PARTICLE_ITERATOR 1
66 
67 
73 };
74 
75 
84 template <class PP> class ColData {
85 public:
86  PP _x;
87  int _dir;
92  ColData( PP x, int dir ) : _x(x), _dir(dir) {}
93 
98  bool operator<( const ColData &cd ) const {
99  return( _x[0] < cd._x[0] );
100  }
101 
110  static void build_coldata_linear( std::vector<ColData> &coldata, const Mesh &mesh,
111  const PP &x1, const PP &x2 ) {
112 
113  coldata.clear();
114 
115  for( size_t a = 0; a < PP::dim(); a++ ) {
116 
117  int a1 = (int)floor( (x1[2*a+1]-mesh.origo(a))/mesh.h() );
118  int a2 = (int)floor( (x2[2*a+1]-mesh.origo(a))/mesh.h() );
119  if( a1 > a2 ) {
120  int a = a2;
121  a2 = a1;
122  a1 = a;
123  }
124 
125  for( int b = a1+1; b <= a2; b++ ) {
126 
127  // Save intersection coordinates
128  double K = (b*mesh.h() + mesh.origo(a) - x1[2*a+1]) /
129  (x2[2*a+1] - x1[2*a+1]);
130  if( K < 0.0 ) K = 0.0;
131  else if( K > 1.0 ) K = 1.0;
132  //std::cout << "Found valid root: " << K << "\n";
133 
134  if( x2[2*a+1] > x1[2*a+1] )
135  coldata.push_back( ColData( x1 + (x2-x1)*K, a+1 ) );
136  else
137  coldata.push_back( ColData( x1 + (x2-x1)*K, -a-1 ) );
138  }
139  }
140 
141  // Sort intersections in increasing time order
142  sort( coldata.begin(), coldata.end() );
143  }
144 
153  static void build_coldata_poly( std::vector<ColData> &coldata, const Mesh &mesh,
154  const PP &x1, const PP &x2 ) {
155 
156 #ifdef DEBUG_PARTICLE_ITERATOR
157  std::cout << "Building coldata using polynomial interpolation\n";
158 #endif
159 
160  coldata.clear();
161 
162  // Construct trajectory representation
163  TrajectoryRep1D traj[PP::dim()];
164  for( size_t a = 0; a < PP::dim(); a++ ) {
165  traj[a].construct( x2[0]-x1[0],
166  x1[2*a+1], x1[2*a+2],
167  x2[2*a+1], x2[2*a+2] );
168  }
169 
170  // Solve trajectory intersections
171  for( size_t a = 0; a < PP::dim(); a++ ) {
172 
173  // Mesh number of x1 (start point)
174  int i = (int)floor( (x1[2*a+1]-mesh.origo(a))/mesh.h() );
175 
176  // Search to negative (dj = -1) and positive (dj = +1) mesh directions
177  for( int dj = -1; dj <= 1; dj += 2 ) {
178  int j = i;
179  if( dj == +1 )
180  j = i+1;
181  int Kcount; // Solution counter
182  double K[3]; // Solution array
183  while( 1 ) {
184 
185  // Intersection point
186  double val = mesh.origo(a) + mesh.h() * j;
187  if( val < mesh.origo(a) )
188  break;
189  else if( val > mesh.max(a) )
190  break;
191 
192 #ifdef DEBUG_PARTICLE_ITERATOR
193  std::cout << " Searching intersections at coord(" << a << ") = " << val << "\n";
194 #endif
195  Kcount = traj[a].solve( K, val );
196  if( Kcount == 0 )
197  break; // No valid roots
198 
199 #ifdef DEBUG_PARTICLE_ITERATOR
200  std::cout << " Found " << Kcount << " valid roots: ";
201  for( int p = 0; p < Kcount; p++ )
202  std::cout << K[p] << " ";
203  std::cout << "\n";
204 #endif
205 
206  // Save roots to coldata
207  for( int b = 0; b < Kcount; b++ ) {
208  PP xcol;
209  double x, v;
210  xcol(0) = x1[0] + K[b]*(x2[0]-x1[0]);
211  for( size_t c = 0; c < PP::dim(); c++ ) {
212  traj[c].coord( x, v, K[b] );
213  if( a == c )
214  xcol[2*c+1] = val; // limit numerical inaccuracy
215  else
216  xcol[2*c+1] = x;
217  xcol[2*c+2] = v;
218  }
219  if( mesh.geom_mode() == MODE_CYL )
220  xcol[5] = x1[5] + K[b]*(x2[5]-x1[5]);
221  if( xcol[2*a+2] >= 0.0 )
222  coldata.push_back( ColData( xcol, a+1 ) );
223  else
224  coldata.push_back( ColData( xcol, -a-1 ) );
225  }
226 
227  j += dj;
228  }
229  }
230  }
231 
232  // Sort intersections in increasing time order
233  sort( coldata.begin(), coldata.end() );
234 
235 #ifdef DEBUG_PARTICLE_ITERATOR
236  std::cout << " Coldata built\n";
237 #endif
238  }
239 
240 };
241 
242 
250 template <class PP> class ParticleIterator {
251 
252  gsl_odeiv_system _system;
253  gsl_odeiv_step *_step;
254  gsl_odeiv_control *_control;
255  gsl_odeiv_evolve *_evolve;
259  bool _polyint;
260  double _epsabs;
261  double _epsrel;
262  uint32_t _maxsteps;
263  double _maxt;
264  uint32_t _trajdiv;
266  bool _mirror[6];
268  ParticleIteratorData _pidata;
269  const TrajectoryHandlerCallback *_thand_cb;
270  const TrajectoryEndCallback *_tend_cb;
271  const TrajectoryEndCallback *_bsup_cb;
272  ParticleDataBase *_pdb;
273  pthread_mutex_t *_scharge_mutex;
275  PP _xi;
277  std::vector<PP> _traj;
278  std::vector<ColData<PP> > _coldata;
280  ParticleStatistics _stat;
288  void save_trajectory_point( PP x ) {
289 
290  try {
291  _traj.push_back( x );
292  } catch( std::bad_alloc ) {
293  throw( ErrorNoMem( ERROR_LOCATION, "Out of memory saving trajectory" ) );
294  }
295  }
296 
305  bool check_collision( Particle<PP> &particle, const PP &x1, const PP &x2, PP &status_x ) {
306 
307  // If inside solid, bracket for collision point
308  Vec3D v2 = x2.location();
309  int32_t bound = _pidata._geom->inside( v2 );
310  if( bound < 7 )
311  return( true ); // No collision happened.
312  Vec3D vc;
313  Vec3D v1 = x1.location();
314  double K = _pidata._geom->bracket_surface( bound, v2, v1, vc );
315 
316  // Calculate new PP
317  for( size_t a = 0; a < PP::size(); a++ )
318  status_x[a] = x2[a] + K*(x1[a]-x2[a]);
319 
320  // Remove all points from trajectory after time status_x[0].
321  for( size_t a = _traj.size()-1; a > 0; a-- ) {
322  if( _traj[a][0] > status_x[0] )
323  _traj.pop_back();
324  else
325  break;
326  }
327 
328  // Save last trajectory point and update status
329  //save_trajectory_point( status_x );
330  particle.set_status( PARTICLE_COLL );
331 
332  // Update collision statistics for boundary
333  _stat.add_bound_collision( bound, particle.IQ() );
334 
335  return( false ); // Collision happened.
336  }
337 
338 
346  void handle_mirror( size_t c, int i[3], size_t a, int border, PP &x2 ) {
347 
348 #ifdef DEBUG_PARTICLE_ITERATOR
349  std::cout << " handle_mirror( c = " << c
350  << ", i = (" << i[0] << ", " << i[1] << ", " << i[2]
351  << "), a = " << a << ", border = " << border
352  << ")\n";
353 #endif
354 
355  double xmirror;
356  if( border < 0 ) {
357  xmirror = _pidata._geom->origo(a);
358  i[a] = -i[a]-1;
359  } else {
360  xmirror = _pidata._geom->max(a);
361  i[a] = 2*_pidata._geom->size(a)-i[a]-3;
362  }
363 
364 #ifdef DEBUG_PARTICLE_ITERATOR
365  std::cout << " xmirror = " << xmirror << "\n";
366  std::cout << " i = (" << i[0] << ", " << i[1] << ", " << i[2] << ")\n";
367  std::cout << " xi = " << _xi << "\n";
368 #endif
369 
370  // Check if found edge at first encounter
371  bool caught_at_boundary = false;
372  if( _coldata[c]._dir == border*((int)a+1) &&
373  ( i[a] == 0 || i[a] == (int)_pidata._geom->size(a)-2 ) ) {
374  caught_at_boundary = true;
375 #ifdef DEBUG_PARTICLE_ITERATOR
376  std::cout << " caught_at_boundary\n";
377 #endif
378  }
379 
380  // Mirror traj back to _xi
381  if( caught_at_boundary ) {
382  save_trajectory_point( _coldata[c]._x );
383  } else {
384  for( int b = _traj.size()-1; b > 0; b-- ) {
385  if( _traj[b][0] >= _xi[0] ) {
386 
387 #ifdef DEBUG_PARTICLE_ITERATOR
388  std::cout << " mirroring traj[" << b << "] = " << _traj[b] << "\n";
389 #endif
390  _traj[b][2*a+1] = 2.0*xmirror - _traj[b][2*a+1];
391  _traj[b][2*a+2] *= -1.0;
392  } else
393  break;
394  }
395  }
396 
397  // Mirror rest of the coldata
398  for( size_t b = c; b < _coldata.size(); b++ ) {
399  if( (size_t)abs(_coldata[b]._dir) == a+1 )
400  _coldata[b]._dir *= -1;
401  _coldata[b]._x[2*a+1] = 2.0*xmirror - _coldata[b]._x[2*a+1];
402  _coldata[b]._x[2*a+2] *= -1.0;
403  }
404 
405  if( caught_at_boundary )
406  save_trajectory_point( _coldata[c]._x );
407 
408  // Mirror calculation point
409  x2[2*a+1] = 2.0*xmirror - x2[2*a+1];
410  x2[2*a+2] *= -1.0;
411 
412  // Coordinates changed, reset integrator
413  gsl_odeiv_step_reset( _step );
414  gsl_odeiv_evolve_reset( _evolve );
415  }
416 
417 
418  void handle_collision( Particle<PP> &particle, uint32_t bound, size_t c, PP &status_x ) {
419 
420 #ifdef DEBUG_PARTICLE_ITERATOR
421  std::cout << " handle_collision()\n";
422 #endif
423 
424  //save_trajectory_point( _coldata[c]._x );
425  status_x = _coldata[c]._x;
426  particle.set_status( PARTICLE_OUT );
427  _stat.add_bound_collision( bound, particle.IQ() );
428  }
429 
430 
439  bool handle_trajectory_advance( Particle<PP> &particle, size_t c, int i[3], PP &x2 ) {
440 
441  // Check for collisions with solids and advance coordinates i.
442  if( PP::dim() == 2 ) {
443  if( _coldata[c]._dir == -1 ) {
444  if( ( abs(_pidata._geom->mesh(i[0], i[1] )) >= 7 ||
445  abs(_pidata._geom->mesh(i[0], i[1]+1)) >= 7 ) &&
446  !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
447  return( false );
448  i[0]--;
449  } else if( _coldata[c]._dir == +1 ) {
450  if( ( abs(_pidata._geom->mesh(i[0]+1,i[1] )) >= 7 ||
451  abs(_pidata._geom->mesh(i[0]+1,i[1]+1)) >= 7 ) &&
452  !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
453  return( false );
454  i[0]++;
455  } else if( _coldata[c]._dir == -2 ) {
456  if( ( abs(_pidata._geom->mesh(i[0], i[1] )) >= 7 ||
457  abs(_pidata._geom->mesh(i[0]+1,i[1] )) >= 7 ) &&
458  !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
459  return( false );
460  i[1]--;
461  } else {
462  if( ( abs(_pidata._geom->mesh(i[0], i[1]+1)) >= 7 ||
463  abs(_pidata._geom->mesh(i[0]+1,i[1]+1)) >= 7 ) &&
464  !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
465  return( false );
466  i[1]++;
467  }
468  } else if( PP::dim() == 3 ) {
469  if( _coldata[c]._dir == -1 ) {
470  if( ( abs(_pidata._geom->mesh(i[0], i[1], i[2] )) >= 7 ||
471  abs(_pidata._geom->mesh(i[0], i[1]+1,i[2] )) >= 7 ||
472  abs(_pidata._geom->mesh(i[0], i[1], i[2]+1)) >= 7 ||
473  abs(_pidata._geom->mesh(i[0], i[1]+1,i[2]+1)) >= 7 ) &&
474  !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
475  return( false );
476  i[0]--;
477  } else if( _coldata[c]._dir == +1 ) {
478  if( ( abs(_pidata._geom->mesh(i[0]+1,i[1], i[2] )) >= 7 ||
479  abs(_pidata._geom->mesh(i[0]+1,i[1]+1,i[2] )) >= 7 ||
480  abs(_pidata._geom->mesh(i[0]+1,i[1], i[2]+1)) >= 7 ||
481  abs(_pidata._geom->mesh(i[0]+1,i[1]+1,i[2]+1)) >= 7 ) &&
482  !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
483  return( false );
484  i[0]++;
485  } else if( _coldata[c]._dir == -2 ) {
486  if( ( abs(_pidata._geom->mesh(i[0], i[1],i[2] )) >= 7 ||
487  abs(_pidata._geom->mesh(i[0]+1,i[1],i[2] )) >= 7 ||
488  abs(_pidata._geom->mesh(i[0], i[1],i[2]+1)) >= 7 ||
489  abs(_pidata._geom->mesh(i[0]+1,i[1],i[2]+1)) >= 7 ) &&
490  !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
491  return( false );
492  i[1]--;
493  } else if( _coldata[c]._dir == +2 ) {
494  if( ( abs(_pidata._geom->mesh(i[0], i[1]+1,i[2] )) >= 7 ||
495  abs(_pidata._geom->mesh(i[0]+1,i[1]+1,i[2] )) >= 7 ||
496  abs(_pidata._geom->mesh(i[0], i[1]+1,i[2]+1)) >= 7 ||
497  abs(_pidata._geom->mesh(i[0]+1,i[1]+1,i[2]+1)) >= 7 ) &&
498  !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
499  return( false );
500  i[1]++;
501  } else if( _coldata[c]._dir == -3 ) {
502  if( ( abs(_pidata._geom->mesh(i[0], i[1], i[2] )) >= 7 ||
503  abs(_pidata._geom->mesh(i[0]+1,i[1], i[2] )) >= 7 ||
504  abs(_pidata._geom->mesh(i[0], i[1]+1,i[2])) >= 7 ||
505  abs(_pidata._geom->mesh(i[0]+1,i[1]+1,i[2])) >= 7 ) &&
506  !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
507  return( false );
508  i[2]--;
509  } else {
510  if( ( abs(_pidata._geom->mesh(i[0], i[1], i[2]+1)) >= 7 ||
511  abs(_pidata._geom->mesh(i[0]+1,i[1], i[2]+1)) >= 7 ||
512  abs(_pidata._geom->mesh(i[0], i[1]+1,i[2]+1)) >= 7 ||
513  abs(_pidata._geom->mesh(i[0]+1,i[1]+1,i[2]+1)) >= 7 ) &&
514  !check_collision( particle, _xi, _coldata[c]._x, x2 ) )
515  return( false );
516  i[2]++;
517  }
518  } else {
519  throw( Error( ERROR_LOCATION, "unsupported dimension number" ) );
520  }
521 
522  // Check for collisions/mirroring with simulation boundary. Here
523  // coordinates i are already advanced to next mesh.
524  for( size_t a = 0; a < PP::dim(); a++ ) {
525 
526  if( i[a] < 0 ) {
527  if( _mirror[2*a] )
528  handle_mirror( c, i, a, -1, x2 );
529  else {
530  handle_collision( particle, 1+2*a, c, x2 );
531  return( false );
532  }
533  } else if( i[a] >= (_pidata._geom->size(a)-1) ) {
534  if( _mirror[2*a+1] )
535  handle_mirror( c, i, a, +1, x2 );
536  else {
537  handle_collision( particle, 2+2*a, c, x2 );
538  return( false );
539  }
540  }
541  }
542 
543  return( true );
544  }
545 
551  bool limit_trajectory_advance( const PP &x1, PP &x2 ) {
552 
553  bool touched = false;
554 
555  for( size_t a = 0; a < PP::dim(); a++ ) {
556 
557  double lim1 = _pidata._geom->origo(a) -
558  (_pidata._geom->size(a)-1)*_pidata._geom->h();
559  double lim2 = _pidata._geom->origo(a) +
560  2*(_pidata._geom->size(a)-1)*_pidata._geom->h();
561 
562  if( x2[2*a+1] < lim1 ) {
563 
564  double K = (lim1 - x1[2*a+1]) / (x2[2*a+1] - x1[2*a+1]);
565  x2 = x1 + K*(x2-x1);
566  touched = true;
567 #ifdef DEBUG_PARTICLE_ITERATOR
568  std::cout << "Limiting step to:\n";
569  std::cout << " x2: " << x2 << "\n";
570 #endif
571  } else if(x2[2*a+1] > lim2 ) {
572 
573  double K = (lim2 - x1[2*a+1]) / (x2[2*a+1] - x1[2*a+1]);
574  x2 = x1 + K*(x2-x1);
575  touched = true;
576 #ifdef DEBUG_PARTICLE_ITERATOR
577  std::cout << "Limiting step to:\n";
578  std::cout << " x2: " << x2 << "\n";
579 #endif
580  }
581  }
582 
583  return( touched );
584  }
585 
590  void build_coldata( bool force_linear, const PP &x1, const PP &x2 ) {
591 
592  try {
593  if( _polyint && !force_linear )
594  ColData<PP>::build_coldata_poly( _coldata, *_pidata._geom, x1, x2 );
595  else
596  ColData<PP>::build_coldata_linear( _coldata, *_pidata._geom, x1, x2 );
597  } catch( std::bad_alloc ) {
598  throw( ErrorNoMem( ERROR_LOCATION, "Out of memory building ColData" ) );
599  }
600 
601  }
602 
622  bool handle_trajectory( Particle<PP> &particle, const PP &x1, PP &x2,
623  bool force_linear, bool first_step ) {
624 
625 #ifdef DEBUG_PARTICLE_ITERATOR
626  std::cout << "Handle trajectory from x1 to x2:\n";
627  std::cout << " x1: " << x1 << "\n";
628  std::cout << " x2: " << x2 << "\n";
629 #endif
630 
631  // Limit trajectory advance to double the simulation box
632  // If limitation done, force to linear interpolation
633  if( limit_trajectory_advance( x1, x2 ) )
634  force_linear = true;
635 
636  build_coldata( force_linear, x1, x2 );
637 
638  // TODO
639  // Remove entrance to geometry if coming from outside or make
640  // code to skip the collision detection for these particles
641 
642  // No intersections, nothing to do
643  if( _coldata.size() == 0 ) {
644 #ifdef DEBUG_PARTICLE_ITERATOR
645  std::cout << "No coldata\n";
646 #endif
647  return( true );
648  }
649 
650  // Starting mesh index
651  int i[3] = {0, 0, 0};
652  for( size_t cdir = 0; cdir < PP::dim(); cdir++ )
653  i[cdir] = (int)floor( (x1[2*cdir+1]-_pidata._geom->origo(cdir))/_pidata._geom->h() );
654 
655  // Process intersection points
656 #ifdef DEBUG_PARTICLE_ITERATOR
657  std::cout << "Process coldata points:\n";
658 #endif
659  for( size_t a = 0; a < _coldata.size(); a++ ) {
660 
661 #ifdef DEBUG_PARTICLE_ITERATOR
662  std::cout << " Coldata " << std::setw(4) << a << ": "
663  << _coldata[a]._x << ", "
664  << std::setw(3) << i[0] << " "
665  << std::setw(3) << i[1] << " "
666  << std::setw(3) << i[2] << " "
667  << std::setw(3) << _coldata[a]._dir << "\n";
668 #endif
669 
670  // Advance particle in mesh, check for possible collisions and
671  // do mirroring.
672  handle_trajectory_advance( particle, a, i, x2 );
673 
674  // Update space charge for one mesh.
675  if( _pidata._scharge )
676  scharge_add_from_trajectory( *_pidata._scharge, _scharge_mutex, particle.IQ(),
677  _xi, _coldata[a]._x );
678 
679  // Call trajectory handler callback
680  if( _thand_cb )
681  (*_thand_cb)( &particle, &_coldata[a]._x, &x2 );
682 
683 #ifdef DEBUG_PARTICLE_ITERATOR
684  if( particle.get_status() == PARTICLE_OUT ) {
685  std::cout << " Particle out\n";
686  std::cout << " x = " << x2 << "\n";
687  } else if( particle.get_status() == PARTICLE_COLL ) {
688  std::cout << " Particle collided\n";
689  std::cout << " x = " << x2 << "\n";
690  }
691 #endif
692  // Clear coldata and exit if particle collided.
693  if( particle.get_status() != PARTICLE_OK ) {
694  save_trajectory_point( x2 );
695  _coldata.clear();
696  return( false );
697  }
698 
699  // Update last intersection point xi.
700  _xi = _coldata[a]._x;
701  }
702 
703 #ifdef DEBUG_PARTICLE_ITERATOR
704  std::cout << "Coldata done\n";
705 #endif
706  _coldata.clear();
707  return( true );
708  }
709 
710 
713  bool axis_mirror_required( const PP &x2 ) {
714  return( _pidata._geom->geom_mode() == MODE_CYL &&
715  x2[4] < 0.0 &&
716  x2[3] <= 0.01*_pidata._geom->h() &&
717  x2[3]*fabs(x2[5]) <= 1.0e-9*fabs(x2[4]) );
718 
719  }
720 
721 
727  bool handle_axis_mirror_step( Particle<PP> &particle, const PP &x1, PP &x2 ) {
728 
729  // Get acceleration at x2
730  double dxdt[5];
731  PP::get_derivatives( x2[0], &x2[1], dxdt, (void *)&_pidata );
732 
733  // Calculate crossover point assuming zero acceleration in
734  // r-direction and constant acceleration in x-direction
735  double dt = -x2[3]/x2[4];
736  PP xc;
737  xc[0] = x2[0]+dt;
738  xc[1] = x2[1]+(x2[2]+0.5*dxdt[1]*dt)*dt;
739  xc[2] = x2[2];
740  xc[3] = x2[3]+x2[4]*dt;
741  xc[4] = x2[4];
742  xc[5] = x2[5];
743 
744  // Mirror x2 to x3
745  PP x3 = 2*xc - x2;
746  x3[3] *= -1.0;
747  x3[4] *= -1.0;
748  x3[5] *= -1.0;
749 
750 #ifdef DEBUG_PARTICLE_ITERATOR
751  std::cout << "Particle mirror:\n";
752  std::cout << " x1: " << x1 << "\n";
753  std::cout << " x2: " << x2 << "\n";
754  std::cout << " xc: " << xc << "\n";
755  std::cout << " x3: " << x3 << "\n";
756 #endif
757 
758  // Handle step with linear interpolation to avoid going to r<=0
759  if( !handle_trajectory( particle, x2, x3, true, false ) )
760  return( false ); // Particle done
761 
762  // Save trajectory calculation points
763  save_trajectory_point( x2 );
764  save_trajectory_point( xc );
765  xc[4] *= -1.0;
766  xc[5] *= -1.0;
767  save_trajectory_point( xc );
768 
769  // Next step not a continuation of previous one, reset
770  // integrator
771  gsl_odeiv_step_reset( _step );
772  gsl_odeiv_evolve_reset( _evolve );
773 
774  // Continue iteration at mirrored point
775  x2 = x3;
776  return( true );
777  }
778 
792  bool check_particle_definition( PP &x ) {
793 
794 #ifdef DEBUG_PARTICLE_ITERATOR
795  std::cout << "Particle defined at:\n";
796  std::cout << " x = " << x << "\n";
797 #endif
798 
799  // Check for NaN
800  for( size_t a = 0; a < PP::size(); a++ ) {
801  if( comp_isnan( x[a] ) )
802  return( false );
803  }
804 
805  // Check if inside solids or outside simulation geometry box.
806  if( _pidata._geom->inside( x.location() ) )
807  return( false );
808 
809  // Check if particle on simulation geometry border and directed outwards
810  /*
811  for( size_t a = 0; a < PP::dim(); a++ ) {
812  if( x[2*a+1] == _pidata._geom->origo(a) && x[2*a+2] < 0.0 ) {
813  if( _mirror[2*a] ) {
814  x[2*a+2] *= -1.0;
815 #ifdef DEBUG_PARTICLE_ITERATOR
816  std::cout << "Mirroring to:\n";
817  std::cout << " x = " << x << "\n";
818 #endif
819  } else {
820  return( false );
821  }
822 
823  } else if( x[2*a+1] == _pidata._geom->max(a) & x[2*a+2] > 0.0 ) {
824  if( _mirror[2*a+1] ) {
825  x[2*a+2] *= -1.0;
826 #ifdef DEBUG_PARTICLE_ITERATOR
827  std::cout << "Mirroring to:\n";
828  std::cout << " x = " << x << "\n";
829 #endif
830  } else {
831  return( false );
832  }
833  }
834  }
835 
836 #ifdef DEBUG_PARTICLE_ITERATOR
837  std::cout << "Definition ok\n";
838 #endif
839 
840  */
841  return( true );
842  }
843 
844  double calculate_dt( const PP &x, const double *dxdt ) {
845 
846  double spd = 0.0, acc = 0.0;
847 
848  for( size_t a = 0; a < (PP::size()-1)/2; a++ ) {
849  //std::cout << "spd += " << dxdt[2*a]*dxdt[2*a] << "\n";
850  spd += dxdt[2*a]*dxdt[2*a];
851  //std::cout << "acc += " << dxdt[2*a+1]*dxdt[2*a+1] << "\n";
852  acc += dxdt[2*a+1]*dxdt[2*a+1];
853  }
854  if( _pidata._geom->geom_mode() == MODE_CYL ) {
855  //std::cout << "MODE_CYL\n";
856  //std::cout << "spd += " << x[3]*x[3]*x[5]*x[5] << "\n";
857  spd += x[3]*x[3]*x[5]*x[5];
858  //std::cout << "acc += " << x[3]*x[3]*dxdt[4]*dxdt[4] << "\n";
859  acc += x[3]*x[3]*dxdt[4]*dxdt[4];
860  }
861  //std::cout << "spd = " << sqrt(spd) << "\n";
862  //std::cout << "acc = " << sqrt(acc) << "\n";
863  spd = _pidata._geom->h() / sqrt(spd);
864  acc = sqrt( 2.0*_pidata._geom->h() / sqrt(acc) );
865 
866  return( spd < acc ? spd : acc );
867  }
868 
869 public:
870 
897  ParticleIterator( particle_iterator_type_e type, double epsabs, double epsrel,
898  bool polyint, uint32_t maxsteps, double maxt,
899  uint32_t trajdiv, bool mirror[6], ScalarField *scharge,
900  pthread_mutex_t *scharge_mutex,
901  const VectorField *efield, const VectorField *bfield,
902  const Geometry *geom )
903  : _type(type), _polyint(polyint), _epsabs(epsabs), _epsrel(epsrel), _maxsteps(maxsteps), _maxt(maxt),
904  _trajdiv(trajdiv), _pidata(scharge,efield,bfield,geom),
905  _thand_cb(0), _tend_cb(0), _bsup_cb(0), _pdb(0), _scharge_mutex(scharge_mutex),
906  _stat(geom->number_of_boundaries()) {
907 
908  // Initialize mirroring
909  _mirror[0] = mirror[0];
910  _mirror[1] = mirror[1];
911  _mirror[2] = mirror[2];
912  _mirror[3] = mirror[3];
913  _mirror[4] = mirror[4];
914  _mirror[5] = mirror[5];
915 
916  // Initialize system of ordinary differential equations (ODE)
917  _system.jacobian = NULL;
918  _system.params = (void *)&_pidata;
919  _system.function = PP::get_derivatives;
920  _system.dimension = PP::size()-1; // Time is not part of differential equation dimensions
921 
922  // Make scale
923  // 2D: x vx y vy
924  // Cyl: x vx r vr omega
925  // 3D: x vx y vy z vz
926  double scale_abs[PP::size()-1];
927  for( uint32_t a = 0; a < (uint32_t)PP::size()-2; a+=2 ) {
928  scale_abs[a+0] = 1.0;
929  scale_abs[a+1] = 1.0e6;
930  }
931  if( _pidata._geom->geom_mode() == MODE_CYL )
932  scale_abs[4] = 1.0;
933 
934  // Initialize ODE solver
935  _step = gsl_odeiv_step_alloc( gsl_odeiv_step_rkck, _system.dimension );
936  //_control = gsl_odeiv_control_standard_new( _epsabs, _epsrel, 1.0, 1.0 );
937  _control = gsl_odeiv_control_scaled_new( _epsabs, _epsrel, 1.0, 1.0, scale_abs, PP::size()-1 );
938  _evolve = gsl_odeiv_evolve_alloc( _system.dimension );
939  }
940 
941 
945  gsl_odeiv_evolve_free( _evolve );
946  gsl_odeiv_control_free( _control );
947  gsl_odeiv_step_free( _step );
948  }
949 
950 
954  _thand_cb = thand_cb;
955  }
956 
957 
961  _tend_cb = tend_cb;
962  _pdb = pdb;
963  }
964 
965 
969  _pidata.set_bfield_suppression_callback( bsup_cb );
970  }
971 
974  const ParticleStatistics &get_statistics( void ) const {
975  return( _stat );
976  }
977 
978 
987  void operator()( Particle<PP> *particle, uint32_t pi ) {
988 
989  // Check particle status
990  if( particle->get_status() != PARTICLE_OK )
991  return;
992 
993  // Copy starting point to x and
994  PP x = particle->x();
995 
996  // Check particle definition
997  if( !check_particle_definition( x ) ) {
998  particle->set_status( PARTICLE_BADDEF );
999  _stat.inc_end_baddef();
1000  return;
1001  }
1002  particle->x() = x;
1003 
1004  // Reset trajectory and save first trajectory point.
1005  _traj.clear();
1006  save_trajectory_point( x );
1007 #ifdef DEBUG_PARTICLE_ITERATOR
1008  std::cout << x[0] << " "
1009  << x[1] << " "
1010  << x[2] << " "
1011  << x[3] << " "
1012  << x[4] << "\n";
1013 #endif
1014  _pidata._qm = particle->qm();
1015  _xi = x;
1016 
1017  // Reset integrator
1018  gsl_odeiv_step_reset( _step );
1019  gsl_odeiv_evolve_reset( _evolve );
1020 
1021  // Make initial guess for step size
1022  //std::cout << "Guess dt ------------------------------------------------\n";
1023  double dxdt[PP::size()-1];
1024  PP::get_derivatives( 0.0, &x[1], dxdt, (void *)&_pidata );
1025  double dt = calculate_dt( x, dxdt );
1026 
1027 #ifdef DEBUG_PARTICLE_ITERATOR
1028  std::cout << "dxdt = ";
1029  for( size_t a = 0; a < PP::size()-1; a++ )
1030  std::cout << dxdt[a] << " ";
1031  std::cout << "\n";
1032  std::cout << "dt = " << dt << "\n";
1033  std::cout << "*** Starting iteration\n";
1034 #endif
1035 
1036  // Iterate ODEs until maximum steps are done, time is used
1037  // or particle collides.
1038  PP x2;
1039  size_t nstp = 0; // Steps taken
1040  while( nstp < _maxsteps && x[0] < _maxt ) {
1041 
1042 #ifdef DEBUG_PARTICLE_ITERATOR
1043  std::cout << "\n*** Step ***\n";
1044  std::cout << " x = " << x2 << "\n";
1045  std::cout << " dt = " << dt << " (proposed)\n";
1046 #endif
1047 
1048  // Take a step.
1049  x2 = x;
1050 
1051  while( true ) {
1052  int retval = gsl_odeiv_evolve_apply( _evolve, _control, _step, &_system,
1053  &x2[0], _maxt, &dt, &x2[1] );
1054  if( retval == IBSIMU_DERIV_ERROR ) {
1055 #ifdef DEBUG_PARTICLE_ITERATOR
1056  std::cout << "Step rejected\n";
1057  std::cout << " x2 = " << x2 << "\n";
1058  std::cout << " dt = " << dt << "\n";
1059 #endif
1060  x2[0] = x[0]; // Reset time (this shouldn't be necessary - there
1061  // is a bug in GSL-1.12, report has been sent)
1062  dt *= 0.5;
1063  if( dt == 0.0 )
1064  throw( Error( ERROR_LOCATION, "too small step size" ) );
1065  //nstp++;
1066  continue;
1067  } else if( retval == GSL_SUCCESS ) {
1068  break;
1069  } else {
1070  throw( Error( ERROR_LOCATION, "gsl_odeiv_evolve_apply failed" ) );
1071  }
1072  }
1073 
1074  // Check step count number and step size validity
1075  if( nstp >= _maxsteps )
1076  break;
1077  if( x2[0] == x[0] )
1078  throw( Error( ERROR_LOCATION, "too small step size" ) );
1079 
1080 #ifdef DEBUG_PARTICLE_ITERATOR
1081  std::cout << "Step accepted from x1 to x2:\n";
1082  std::cout << " dt = " << dt << " (taken)\n";
1083  std::cout << " x1 = " << x << "\n";
1084  std::cout << " x2 = " << x2 << "\n";
1085 #endif
1086 
1087  // Handle collisions and space charge of step.
1088  if( !handle_trajectory( *particle, x, x2, false, nstp == 0 ) ) {
1089  x = x2;
1090  break; // Particle done
1091  }
1092 
1093  // Check if particle mirroring is required to avoid
1094  // singularity at symmetry axis.
1095  if( axis_mirror_required( x2 ) ) {
1096  if( !handle_axis_mirror_step( *particle, x, x2 ) )
1097  break; // Particle done
1098  }
1099 
1100  // Propagate coordinates
1101  x = x2;
1102 
1103  // Save trajectory point
1104  save_trajectory_point( x2 );
1105 
1106  // Increase step count.
1107  nstp++;
1108  }
1109 
1110 #ifdef DEBUG_PARTICLE_ITERATOR
1111  std::cout << "\n*** Done stepping ***\n";
1112 #endif
1113 
1114  // Check if step count or time limited
1115  if( nstp == _maxsteps ) {
1116  particle->set_status( PARTICLE_NSTP );
1117  _stat.inc_end_step();
1118  } else if( x[0] >= _maxt ) {
1119  particle->set_status( PARTICLE_TIME );
1120  _stat.inc_end_time();
1121  }
1122 
1123  // Save step count
1124  _stat.inc_sum_steps( nstp );
1125 
1126  // Save trajectory of current particle
1127  if( _trajdiv != 0 && pi % _trajdiv == 0 )
1128  particle->copy_trajectory( _traj );
1129 
1130  // Save last particle location
1131  particle->x() = x;
1132 
1133  // Call trajectory end callback
1134  if( _tend_cb )
1135  (*_tend_cb)( particle, _pdb );
1136  }
1137 
1138 };
1139 
1140 
1141 #endif
1142 
#define IBSIMU_DERIV_ERROR
Definition: particles.hpp:58
uint32_t inside(const Vec3D &x) const
Returns 0 if point x is vacuum or the number of solid of x is inside a defined solid. Returns a number from 1 to 6 if point x is outside the defined geometry. If the point is inside several defined solids, the solid with the highest solid number is returned.
void add_bound_collision(uint32_t bound, double IQ)
Geometry definition
static void build_coldata_poly(std::vector< ColData > &coldata, const Mesh &mesh, const PP &x1, const PP &x2)
Find mesh intersections of polynomially interpolated particle trajectory segment. ...
Definition: particleiterator.hpp:153
Basic error class.
Definition: error.hpp:142
Trajectory interpolation solver.
Scalar fields.
void inc_sum_steps(void)
Definition: particlestatistics.hpp:99
Abstract base class for vector field.
Definition: vectorfield.hpp:53
void inc_end_step(void)
Definition: particlestatistics.hpp:97
Definition: particles.hpp:70
double bracket_surface(uint32_t n, const Vec3D &xin, const Vec3D &xout, Vec3D &xsurf) const
Find solid n surface location by bracketing.
Particle and particle point objects
Definition: particles.hpp:72
Temporary data bundle for particle iterators.
Definition: particles.hpp:900
Mesh geometry definion.
Definition: mesh.hpp:67
Vec3D origo(void) const
Returns origo vector of geometry.
Definition: mesh.hpp:128
void set_bfield_suppression_callback(const CallbackFunctorD_V *bsup_cb)
Set B-field potential dependent suppression callback.
Definition: particles.hpp:915
Trajectory representation between two calculated points in 1d.
Definition: trajectory.hpp:61
void inc_end_baddef(void)
Definition: particlestatistics.hpp:98
Trajectory end callback.
Definition: particledatabase.hpp:71
Particle databases
Polynomial solver.
Definition: particles.hpp:71
Particle iterator class for continuous Vlasov-type iteration.
Definition: particleiterator.hpp:250
Compatibility math.
Definition: particleiterator.hpp:71
Trajectory handler callback.
Definition: particledatabase.hpp:57
static void build_coldata_linear(std::vector< ColData > &coldata, const Mesh &mesh, const PP &x1, const PP &x2)
Find mesh intersections of linearly interpolated particle trajectory segment.
Definition: particleiterator.hpp:110
Geometry defining class.
Definition: geometry.hpp:131
PP _x
Mesh intersection coordinates.
Definition: particleiterator.hpp:86
void set_trajectory_handler_callback(const TrajectoryHandlerCallback *thand_cb)
Set trajectory handler callback.
Definition: particleiterator.hpp:953
#define ERROR_LOCATION
Macro for setting error location when throwing errors.
Definition: error.hpp:72
Job scheduler for parallel processing.
double h(void) const
Returns mesh cell size.
Definition: mesh.hpp:146
PP & x()
Return reference to coordinate data.
Definition: particles.hpp:784
Space charge deposition functions.
particle_status_e get_status()
Return particle status.
Definition: particles.hpp:680
Definition: callback.hpp:61
void set_trajectory_end_callback(const TrajectoryEndCallback *tend_cb, ParticleDataBase *pdb)
Set trajectory end callback.
Definition: particleiterator.hpp:960
~ParticleIterator()
Destructor.
Definition: particleiterator.hpp:944
double qm() const
Return charge per mass ratio (q/m) [C/kg].
Definition: particles.hpp:703
Definition: particles.hpp:69
void set_bfield_suppression_callback(const CallbackFunctorD_V *bsup_cb)
Set B-field potential dependent suppression callback.
Definition: particleiterator.hpp:968
void operator()(Particle< PP > *particle, uint32_t pi)
Iterate a particle from start to end.
Definition: particleiterator.hpp:987
Particle class in some geometry.
Definition: particles.hpp:724
ColData(PP x, int dir)
Constructor for collision at x into direction dir.
Definition: particleiterator.hpp:92
Definition: particles.hpp:74
Definition: particles.hpp:73
void set_status(particle_status_e status)
Set particle status.
Definition: particles.hpp:684
ScalarField * _scharge
Space charge field or NULL.
Definition: particles.hpp:901
Particle iteration statistics.
Definition: particlestatistics.hpp:57
Vec3D max(void) const
Returns vector pointing to the last mesh point opposite of origo.
Definition: mesh.hpp:137
Mesh intersection (collision) coordinate data
Definition: particleiterator.hpp:84
Error class for memory allocation errors.
Definition: error.hpp:185
void scharge_add_from_trajectory(ScalarField &scharge, pthread_mutex_t *mutex, double IQ, const ParticleP2D &x1, const ParticleP2D &x2)
Function for adding charge to space charge density map from particle trajectory in 2d simulation...
const ParticleStatistics & get_statistics(void) const
Get particle iterator statistics.
Definition: particleiterator.hpp:974
particle_iterator_type_e
Particle iterator type.
Definition: particleiterator.hpp:70
void coord(double &x, double &v, double K)
Calculate location x and velocity v at parametric time K.
Definition: particleiterator.hpp:72
double _qm
Precalculated q/m.
Definition: particles.hpp:905
bool operator<(const ColData &cd) const
Compare coldata entry times.
Definition: particleiterator.hpp:98
double IQ() const
Return current or charge carried by trajectory or particle cloud [A/C].
Definition: particles.hpp:691
const signed char & mesh(int32_t i) const
Returns a const reference to solid mesh.
Definition: geometry.hpp:282
int solve(double K[3], double x, int extrapolate=0)
Solves for trajectory intersection with location.
void construct(double dt, double x1, double v1, double x2, double v2, trajectory_rep_e force=TRAJ_EMPTY)
Construct representation of trajectory from (x1,v1) to (x2,v2) in time dt.
Int3D size(void) const
Returns size array of geometry.
Definition: mesh.hpp:116
const Geometry * _geom
Geometry.
Definition: particles.hpp:904
Cylindrically symmetric geometry.
Definition: types.hpp:62
Three dimensional vector.
Definition: vec3d.hpp:58
void copy_trajectory(const std::vector< PP > &traj)
Define trajectory by copying.
Definition: particles.hpp:808
ParticleIterator(particle_iterator_type_e type, double epsabs, double epsrel, bool polyint, uint32_t maxsteps, double maxt, uint32_t trajdiv, bool mirror[6], ScalarField *scharge, pthread_mutex_t *scharge_mutex, const VectorField *efield, const VectorField *bfield, const Geometry *geom)
Constructor for new particle iterator.
Definition: particleiterator.hpp:897
Scalar field class.
Definition: scalarfield.hpp:70
int _dir
Direction of particle at intersection. i: -1/+1, j: -2/+2, k: -3:/+3.
Definition: particleiterator.hpp:87
geom_mode_e geom_mode(void) const
Returns geometry mode.
Definition: mesh.hpp:108
Particle database base class.
Definition: particledatabase.hpp:167
void inc_end_time(void)
Definition: particlestatistics.hpp:96
int comp_isnan(double x)
Vector field base.