00001
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043 #ifndef PARTICLEDATABASE_HPP
00044 #define PARTICLEDATABASE_HPP 1
00045
00046
00047 #include <vector>
00048 #include "timer.hpp"
00049 #include "ibsimu.hpp"
00050 #include "trajectory.hpp"
00051 #include "particles.hpp"
00052 #include "particleiterator.hpp"
00053 #include "trajectorydiagnostics.hpp"
00054
00055
00056
00057
00058
00059
00060
00061
00069 class ParticleDataBase {
00070
00071 protected:
00072
00073 uint32_t _threadcount;
00074 double _epsabs;
00075 double _epsrel;
00076 bool _polyint;
00077 uint32_t _maxsteps;
00078 double _maxt;
00079 uint32_t _trajdiv;
00081 bool _mirror[6];
00083 double _rhosum;
00085 uint32_t _end_time;
00086 uint32_t _end_step;
00087 uint32_t _end_out;
00089 uint32_t _end_coll;
00091 uint32_t _end_baddef;
00092 uint32_t _sum_steps;
00094 int _iteration;
00096 bool _nsimp;
00097 const ScalarField *_epot;
00098 double _phi_plasma;
00099
00102 ParticleDataBase()
00103 : _threadcount(1), _epsabs(1e-6), _epsrel(1e-6), _polyint(true), _maxsteps(1000),
00104 _maxt(1e-3), _trajdiv(1), _rhosum(0.0), _end_time(0), _end_step(0), _end_out(0),
00105 _end_coll(0), _end_baddef(0), _sum_steps(0), _iteration(-1), _nsimp(false),
00106 _epot(NULL), _phi_plasma(0.0) {
00107 _mirror[0] = false;
00108 _mirror[1] = false;
00109 _mirror[2] = false;
00110 _mirror[3] = false;
00111 _mirror[4] = false;
00112 _mirror[5] = false;
00113 }
00114
00115 public:
00116
00117
00118
00119
00120
00123 virtual ~ParticleDataBase() {}
00124
00125
00126
00127
00128
00131 void set_thread_count( uint32_t threadcount ) {
00132 if( threadcount <= 0 )
00133 throw( Error( ERROR_LOCATION, "invalid parameter" ) );
00134 _threadcount = threadcount;
00135 }
00136
00142 void set_accuracy( double epsabs, double epsrel ) {
00143 _epsabs = epsabs;
00144 _epsrel = epsrel;
00145 }
00146
00147 void enable_nsimp_plasma_threshold( const ScalarField *epot, double phi_plasma ) {
00148 _nsimp = true;
00149 _epot = epot;
00150 _phi_plasma = phi_plasma;
00151 }
00152
00157 void set_polyint( bool polyint ) {
00158 _polyint = polyint;
00159 }
00160
00166 bool get_polyint( void ) const {
00167 return( _polyint );
00168 }
00169
00174 void set_max_steps( uint32_t maxsteps ) {
00175 if( maxsteps <= 0 )
00176 throw( Error( ERROR_LOCATION, "invalid parameter" ) );
00177 _maxsteps = maxsteps;
00178 }
00179
00184 void set_max_time( double maxt ) {
00185 if( maxt <= 0.0 )
00186 throw( Error( ERROR_LOCATION, "invalid parameter" ) );
00187 _maxt = maxt;
00188 }
00189
00196 void set_save_trajectories( uint32_t div ) {
00197 _trajdiv = div;
00198 }
00199
00204 void set_mirror( const bool mirror[6] ) {
00205 _mirror[0] = mirror[0];
00206 _mirror[1] = mirror[1];
00207 _mirror[2] = mirror[2];
00208 _mirror[3] = mirror[3];
00209 _mirror[4] = mirror[4];
00210 _mirror[5] = mirror[5];
00211 }
00212
00217 void get_mirror( bool mirror[6] ) const {
00218 mirror[0] = _mirror[0];
00219 mirror[1] = _mirror[1];
00220 mirror[2] = _mirror[2];
00221 mirror[3] = _mirror[3];
00222 mirror[4] = _mirror[4];
00223 mirror[5] = _mirror[5];
00224 }
00225
00226 int get_iteration_number( void ) const {
00227 return( _iteration );
00228 }
00229
00243 double get_rhosum( void ) {
00244 return( _rhosum );
00245 }
00246
00247
00248
00249
00250
00253 virtual size_t size( void ) const = 0;
00254
00257 virtual ParticleBase &particle( uint32_t i ) = 0;
00258
00261 virtual const ParticleBase &particle( uint32_t i ) const = 0;
00262
00265 virtual size_t traj_size( uint32_t i ) const = 0;
00266
00269 virtual void trajectory_point( double &t, Vec3D &loc, Vec3D &vel, uint32_t i, uint32_t j ) const = 0;
00270
00274 virtual void trajectories_at_plane( TrajectoryDiagnosticData &tdata,
00275 coordinate_axis_e axis,
00276 double val,
00277 const std::vector<trajectory_diagnostic_e> &diagnostics ) const = 0;
00278
00279
00280
00281
00282
00288 virtual void clear( void ) = 0;
00289
00295 virtual void clear_trajectories( void ) = 0;
00296
00297
00298
00299
00300
00303 virtual void reserve( size_t size ) = 0;
00304
00305
00306
00307
00308
00309 virtual void debug_print( void ) const = 0;
00310
00311 };
00312
00313
00328 template<class PP> class ParticleDataBasePP : public ParticleDataBase {
00329
00332 static void add_diagnostics( TrajectoryDiagnosticData &tdata, const PP &x,
00333 const Particle<PP> &p, int crd ) {
00334
00335 for( size_t a = 0; a < tdata.diag_size(); a++ ) {
00336
00337
00338 double data = 0.0;
00339 switch( tdata.diagnostic( a ) ) {
00340 case DIAG_NONE:
00341 data = 0.0;
00342 break;
00343 case DIAG_T:
00344 data = x[0];
00345 break;
00346 case DIAG_X:
00347 data = x[1];
00348 break;
00349 case DIAG_VX:
00350 data = x[2];
00351 break;
00352 case DIAG_Y:
00353 case DIAG_R:
00354 data = x[3];
00355 break;
00356 case DIAG_VY:
00357 case DIAG_VR:
00358 data = x[4];
00359 break;
00360 case DIAG_Z:
00361 data = x[5];
00362 break;
00363 case DIAG_VZ:
00364 data = x[6];
00365 break;
00366 case DIAG_W:
00367 data = x[5];
00368 break;
00369 case DIAG_VTHETA:
00370 data = x[5]*x[3];
00371 break;
00372 case DIAG_XP:
00373 data = x[2]/x[2*crd+2];
00374 break;
00375 case DIAG_YP:
00376 case DIAG_RP:
00377 data = x[4]/x[2*crd+2];
00378 break;
00379 case DIAG_AP:
00380 data = x[3]*x[5]/x[2*crd+2];
00381 break;
00382 case DIAG_ZP:
00383 data = x[6]/x[2*crd+2];
00384 break;
00385 case DIAG_CURR:
00386 data = p.IQ();
00387 break;
00388 case DIAG_QM:
00389 data = p.qm();
00390 break;
00391 case DIAG_EK:
00392 Vec3D velocity = x.velocity();
00393 data = velocity.norm2();
00394 break;
00395 }
00396
00397 tdata.add_data( a, data );
00398 }
00399 }
00400
00401 protected:
00402
00403 std::vector<Particle<PP> > _particles;
00407 ParticleDataBasePP() {}
00408
00409 public:
00410
00411
00412
00413
00414
00417 ~ParticleDataBasePP() {}
00418
00419
00420
00421
00422
00425 virtual size_t size( void ) const { return( _particles.size() ); }
00426
00429 virtual Particle<PP> &particle( uint32_t i ) { return( _particles[i] ); }
00430
00433 virtual const Particle<PP> &particle( uint32_t i ) const { return( _particles[i] ); }
00434
00437 virtual size_t traj_size( uint32_t i ) const { return( _particles[i].traj_size() ); }
00438
00441 virtual void trajectory_point( double &t, Vec3D &loc, Vec3D &vel, uint32_t i, uint32_t j ) const {
00442 PP x = _particles[i].traj(j);
00443 t = x[0];
00444 loc = x.location();
00445 vel = x.velocity();
00446 }
00447
00451 virtual void trajectories_at_plane( TrajectoryDiagnosticData &tdata,
00452 coordinate_axis_e axis,
00453 double val,
00454 const std::vector<trajectory_diagnostic_e> &diagnostics ) const {
00455
00456 if( ibsimu.get_verbose_output() )
00457 std::cout << "Making trajectory diagnostics at "
00458 << coordinate_axis_string[axis] << " = " << val << "\n";
00459
00460
00461 switch( PP::geom_mode() ) {
00462 case MODE_1D:
00463 throw( Error( ERROR_LOCATION, "unsupported dimension number" ) );
00464 break;
00465 case MODE_2D:
00466 if( axis == AXIS_R || axis == AXIS_Z )
00467 throw( Error( ERROR_LOCATION, "nonexistent axis" ) );
00468 break;
00469 case MODE_CYL:
00470 if( axis == AXIS_Y || axis == AXIS_Z )
00471 throw( Error( ERROR_LOCATION, "nonexistent axis" ) );
00472 break;
00473 case MODE_3D:
00474 if( axis == AXIS_R )
00475 throw( Error( ERROR_LOCATION, "nonexistent axis" ) );
00476 break;
00477 default:
00478 throw( Error( ERROR_LOCATION, "unsupported dimension number" ) );
00479 }
00480
00481
00482 for( size_t a = 0; a < diagnostics.size(); a++ ) {
00483 if( diagnostics[a] == DIAG_NONE )
00484 throw( Error( ERROR_LOCATION, "invalid diagnostics query \'DIAG_NONE\'" ) );
00485 else if( PP::geom_mode() != MODE_CYL && (diagnostics[a] == DIAG_R ||
00486 diagnostics[a] == DIAG_VR ||
00487 diagnostics[a] == DIAG_RP ||
00488 diagnostics[a] == DIAG_W ||
00489 diagnostics[a] == DIAG_VTHETA ||
00490 diagnostics[a] == DIAG_AP) )
00491 throw( Error( ERROR_LOCATION, "invalid diagnostics query for geometry type" ) );
00492 }
00493
00494
00495 tdata.clear();
00496 for( size_t a = 0; a < diagnostics.size(); a++ ) {
00497 tdata.add_data_column( diagnostics[a] );
00498 }
00499
00500
00501 int crd;
00502 switch( axis ) {
00503 case AXIS_X:
00504 crd = 0;
00505 break;
00506 case AXIS_Y:
00507 case AXIS_R:
00508 crd = 1;
00509 break;
00510 case AXIS_Z:
00511 crd = 2;
00512 break;
00513 default:
00514 throw( Error( ERROR_LOCATION, "unsupported axis" ) );
00515 }
00516
00517
00518 double Isum = 0.0;
00519 std::vector<PP> intsc;
00520 for( size_t a = 0; a < _particles.size(); a++ ) {
00521 size_t N = _particles[a].traj_size();
00522 if( N < 2 )
00523 continue;
00524 PP x1 = _particles[a].traj(0);
00525 size_t nintsum = 0;
00526 for( size_t b = 1; b < N; b++ ) {
00527 PP x2 = _particles[a].traj(b);
00528 intsc.clear();
00529 size_t nintsc = PP::trajectory_intersections_at_plane( intsc, crd, val, x1, x2 );
00530 nintsum += nintsc;
00531 for( size_t c = 0; c < nintsc; c++ ) {
00532 Isum += _particles[a].IQ();
00533 add_diagnostics( tdata, intsc[c], _particles[a], crd );
00534 }
00535
00536 x1 = x2;
00537 }
00538 }
00539
00540 if( ibsimu.get_verbose_output() ) {
00541 std::cout << " number of trajectories = " << tdata.traj_size() << "\n";
00542 if( PP::geom_mode() == MODE_2D )
00543 std::cout << " total current = " << Isum << " A/m\n";
00544 else
00545 std::cout << " total current = " << Isum << " A\n";
00546 }
00547 }
00548
00549
00550
00551
00552
00555 virtual void clear( void ) {
00556 _particles.clear();
00557 _rhosum = 0.0;
00558 }
00559
00565 virtual void clear_trajectories( void ) {
00566 for( uint32_t a = 0; a < _particles.size(); a++ )
00567 _particles[a].clear_trajectory();
00568 }
00569
00570
00571
00572
00573
00576 virtual void reserve( size_t size ) { _particles.reserve( size ); }
00577
00589 void add_particle( double IQ, double q, double m, const PP &x ) {
00590 _particles.push_back( Particle<PP>( IQ, CHARGE_E*q, MASS_U*m, x ) );
00591 }
00592
00597 void add_particle( const Particle<PP> &pp ) {
00598 _particles.push_back( pp );
00599 }
00600
00601
00602
00603
00604
00612 virtual void iterate_trajectories( ScalarField &scharge, const Efield &efield,
00613 const VectorField &bfield, const Geometry &g ) {
00614
00615 ScalarField *schmap[_threadcount];
00616 std::vector<ParticleIterator<PP> *> iterators;
00617
00618 Timer t;
00619 if( ibsimu.get_verbose_output() )
00620 std::cout << "Calculating particle trajectories\n";
00621 _iteration++;
00622
00623
00624 if( g.geom_mode() != PP::geom_mode() )
00625 throw( Error( ERROR_LOCATION, "Differing geometry modes" ) );
00626
00627
00628 scharge.clear();
00629
00630
00631 _end_time = 0;
00632 _end_step = 0;
00633 _end_out = 0;
00634 _end_coll = 0;
00635 _end_baddef = 0;
00636 _sum_steps = 0;
00637
00638
00639 if( _particles.size() == 0 ) {
00640 std::cout << " no particles to calculate\n";
00641 return;
00642 }
00643
00644
00645 for( uint32_t a = 0; a < _threadcount; a++ ) {
00646 if( a == 0 ) schmap[a] = &scharge;
00647 else schmap[a] = new ScalarField( g );
00648 iterators.push_back( new ParticleIterator<PP>( PARTICLE_ITERATOR_ADAPTIVE, _epsabs, _epsrel,
00649 _polyint, _maxsteps,
00650 _maxt, _trajdiv, _mirror, schmap[a],
00651 &efield, &bfield, &g, &_particles[0] ) );
00652 if( _nsimp )
00653 iterators[a]->enable_nsimp_plasma_threshold( _epot, _phi_plasma );
00654 }
00655
00656
00657 Scheduler<ParticleIterator<PP>,Particle<PP>,Error> scheduler( iterators );
00658
00659
00660 for( size_t a = 0; a < _particles.size(); a++ )
00661 scheduler.add_problem( &_particles[a] );
00662
00663
00664 scheduler.run();
00665 scheduler.finish();
00666
00667 if( scheduler.is_error() ) {
00668
00669 std::vector<Error> err;
00670 std::vector<Particle<PP> *> part;
00671 scheduler.get_errors( err, part );
00672 throw( err[0] );
00673 }
00674
00675
00676
00677 for( uint32_t a = 0; a < _threadcount; a++ ) {
00678 if( a != 0 ) {
00679 scharge += *schmap[a];
00680 delete schmap[a];
00681 }
00682 uint32_t end_time, end_step, end_out, end_coll, end_baddef, sum_steps;
00683 iterators[a]->get_statistics( end_time, end_step, end_out,
00684 end_coll, end_baddef, sum_steps );
00685 _end_time += end_time;
00686 _end_step += end_step;
00687 _end_out += end_out;
00688 _end_coll += end_coll;
00689 _end_baddef += end_baddef;
00690 _sum_steps += sum_steps;
00691 delete iterators[a];
00692 }
00693
00694 scharge_finalize( scharge );
00695
00696 t.stop();
00697 if( ibsimu.get_verbose_output() ) {
00698 std::cout << " Particle histories (" << _particles.size() << " total):\n";
00699 std::cout << " time limited = " << _end_time << "\n";
00700 std::cout << " step count limited = " << _end_step << "\n";
00701 std::cout << " out of geometry = " << _end_out << "\n";
00702 std::cout << " collided = " << _end_coll << "\n";
00703 std::cout << " bad definitions = " << _end_baddef << "\n";
00704 std::cout << " total steps = " << _sum_steps << "\n";
00705 std::cout << " steps per particle (ave) = " << _sum_steps/(double)_particles.size() << "\n";
00706 std::cout << " time used = " << t << "\n";
00707 }
00708 }
00709
00716 virtual void step_particles( const Efield &efield, const Geometry &g, double dt ) {
00717
00718 }
00719
00720
00721
00722
00723
00726 virtual void debug_print( void ) const {
00727 std::cout << "threadcount = " << _threadcount << "\n";
00728 std::cout << "epsabs = " << _epsabs << "\n";
00729 std::cout << "epsrel = " << _epsrel << "\n";
00730 std::cout << "maxsteps = " << _maxsteps << "\n";
00731 std::cout << "maxt = " << _maxt << "\n";
00732 std::cout << "trajdiv = " << _trajdiv << "\n";
00733 std::cout << "mirror = (";
00734 for( uint32_t a = 0; a < 5; a++ )
00735 std::cout << _mirror[a] << ", ";
00736 std::cout << _mirror[5] << ")\n";
00737
00738 for( uint32_t a = 0; a < _particles.size(); a++ ) {
00739 std::cout << "Particle " << a << ":\n";
00740 _particles[a].debug_print();
00741 }
00742 }
00743
00744 };
00745
00746
00747
00759 class ParticleDataBase2D : public ParticleDataBasePP<ParticleP2D> {
00760
00761
00762 public:
00763
00764
00765
00766
00767
00770 ParticleDataBase2D() {}
00771
00774 ~ParticleDataBase2D() {}
00775
00776
00777
00778
00779
00799 void add_2d_beam_with_energy( uint32_t N, double J, double q, double m,
00800 double E, double Tp, double Tt,
00801 double x1, double y1, double x2, double y2 );
00802
00818 void add_2d_beam_with_velocity( uint32_t N, double J, double q, double m,
00819 double v, double dvp, double dvt,
00820 double x1, double y1, double x2, double y2 );
00821
00838 void add_2d_KV_beam_with_emittance( uint32_t N, double I, double q, double m,
00839 double a, double b, double e,
00840 double Ex, double x0, double y0 );
00841
00858 void add_2d_gaussian_beam_with_emittance( uint32_t N, double I, double q, double m,
00859 double a, double b, double e,
00860 double Ex, double x0, double y0 );
00861 };
00862
00863
00864
00865
00877 class ParticleDataBaseCyl : public ParticleDataBasePP<ParticlePCyl> {
00878
00879
00880 public:
00881
00882
00883
00884
00885
00888 ParticleDataBaseCyl() {}
00889
00892 ~ParticleDataBaseCyl () {}
00893
00894
00895
00896
00897
00917 void add_2d_beam_with_energy( uint32_t N, double J, double q, double m,
00918 double E, double Tp, double Tt,
00919 double x1, double y1, double x2, double y2 );
00920
00936 void add_2d_beam_with_velocity( uint32_t N, double J, double q, double m,
00937 double v, double dvp, double dvt,
00938 double x1, double y1, double x2, double y2 );
00939
00942 void add_2d_full_gaussian_beam( uint32_t N, double I, double q, double m,
00943 double Ex, double Tp, double Tt,
00944 double x0, double dr );
00945
00964 void add_2d_gaussian_beam_with_emittance( uint32_t N, double I, double q, double m,
00965 double a, double b, double e,
00966 double Ex, double x0 );
00967 };
00968
00969
00970
00982 class ParticleDataBase3D : public ParticleDataBasePP<ParticleP3D> {
00983
00984
00985 public:
00986
00987
00988
00989
00990
00993 ParticleDataBase3D() {}
00994
00997 ~ParticleDataBase3D() {}
00998
00999
01000
01001
01002
01026 void add_cylindrical_beam_with_energy( uint32_t N, double J, double q, double m,
01027 double E, double Tp, double Tt, Vec3D c,
01028 Vec3D dir1, Vec3D dir2, double r );
01029
01052 void add_cylindrical_beam_with_velocity( uint32_t N, double J, double q, double m,
01053 double v, double dvp, double dvt, Vec3D c,
01054 Vec3D dir1, Vec3D dir2, double r );
01055
01064 void add_rectangular_beam_with_energy( uint32_t N, double J, double q, double m,
01065 double E, double Tp, double Tt, Vec3D c,
01066 Vec3D dir1, Vec3D dir2, double size1, double size2 );
01067
01076 void add_rectangular_beam_with_velocity( uint32_t N, double J, double q, double m,
01077 double v, double dvp, double dvt, Vec3D c,
01078 Vec3D dir1, Vec3D dir2, double size1, double size2 );
01079
01080
01098 void add_3d_KV_beam_with_emittance( uint32_t N, double I, double q, double m,
01099 double ay, double by, double ey,
01100 double az, double bz, double ez,
01101 double Ex, double x0, double y0, double z0 );
01102
01120 void add_3d_gaussian_beam_with_emittance( uint32_t N, double I, double q, double m,
01121 double ay, double by, double ey,
01122 double az, double bz, double ez,
01123 double Ex, double x0, double y0, double z0 );
01124 };
01125
01126
01127
01128
01129
01130
01131 #endif
01132
01133
01134
01135
01136
01137
01138
01139
01140
01141
01142
01143
01144
01145
01146
01147
01148
01149