OPAL (Object Oriented Parallel Accelerator Library) 2024.2
OPAL
ParallelCyclotronTracker.cpp
Go to the documentation of this file.
1//
2// Class ParallelCyclotronTracker
3// Tracker for OPAL-Cycl
4//
5// Copyright (c) 2007 - 2014, Jianjun Yang, Andreas Adelmann and Matthias Toggweiler,
6// Paul Scherrer Institut, Villigen PSI, Switzerland
7// Copyright (c) 2014, Daniel Winklehner, MIT, Cambridge, MA, USA
8// Copyright (c) 2012 - 2023, Paul Scherrer Institut, Villigen PSI, Switzerland
9// All rights reserved
10//
11// Implemented as part of the PhD thesis
12// "Beam dynamics in high intensity cyclotrons including neighboring bunch effects"
13// and the paper
14// "Beam dynamics in high intensity cyclotrons including neighboring bunch effects:
15// Model, implementation, and application"
16// (https://journals.aps.org/prab/pdf/10.1103/PhysRevSTAB.13.064201)
17//
18// This file is part of OPAL.
19//
20// OPAL is free software: you can redistribute it and/or modify
21// it under the terms of the GNU General Public License as published by
22// the Free Software Foundation, either version 3 of the License, or
23// (at your option) any later version.
24//
25// You should have received a copy of the GNU General Public License
26// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
27//
28#include "Algorithms/ParallelCyclotronTracker.h"
29
30#include "AbsBeamline/CCollimator.h"
31#include "AbsBeamline/Corrector.h"
32#include "AbsBeamline/Cyclotron.h"
33#include "AbsBeamline/Degrader.h"
34#include "AbsBeamline/Drift.h"
35#include "AbsBeamline/Marker.h"
36#include "AbsBeamline/Monitor.h"
37#include "AbsBeamline/Multipole.h"
38#include "AbsBeamline/MultipoleT.h"
39#include "AbsBeamline/MultipoleTBase.h"
40#include "AbsBeamline/Offset.h"
41#include "AbsBeamline/OutputPlane.h"
42#include "AbsBeamline/PluginElement.h"
43#include "AbsBeamline/Probe.h"
44#include "AbsBeamline/RBend.h"
45#include "AbsBeamline/RFCavity.h"
46#include "AbsBeamline/Ring.h"
47#include "AbsBeamline/SBend.h"
48#include "AbsBeamline/SBend3D.h"
49#include "AbsBeamline/ScalingFFAMagnet.h"
50#include "AbsBeamline/Septum.h"
51#include "AbsBeamline/Solenoid.h"
52#include "AbsBeamline/Stripper.h"
53#include "AbsBeamline/Vacuum.h"
54#include "AbsBeamline/VariableRFCavity.h"
55#include "AbsBeamline/VariableRFCavityFringeField.h"
56#include "AbsBeamline/VerticalFFAMagnet.h"
57
58#include "AbstractObjects/Element.h"
59#include "AbstractObjects/OpalData.h"
60
61#include "Algorithms/Ctunes.h"
62#include "Algorithms/PolynomialTimeDependence.h"
63
64#include "BasicActions/DumpEMFields.h"
65#include "BasicActions/DumpFields.h"
66
67#include "Beamlines/Beamline.h"
68#include "Beamlines/FlaggedBeamline.h"
69
70#include "Distribution/Distribution.h"
71
72#include "Elements/OpalBeamline.h"
73
74#include "Physics/Physics.h"
75#include "Physics/Units.h"
76
77#include "Structure/BoundaryGeometry.h"
78#include "Structure/DataSink.h"
79#include "Structure/LossDataSink.h"
80
81#include "Utilities/OpalException.h"
82#include "Utilities/Options.h"
83
84#include <cmath>
85#include <fstream>
86#include <iostream>
87#include <limits>
88#include <numeric>
89
90Vector_t const ParallelCyclotronTracker::xaxis = Vector_t(1.0, 0.0, 0.0);
91Vector_t const ParallelCyclotronTracker::yaxis = Vector_t(0.0, 1.0, 0.0);
92Vector_t const ParallelCyclotronTracker::zaxis = Vector_t(0.0, 0.0, 1.0);
93
94extern Inform *gmsg;
95extern Inform *gmsgALL;
96
97
110ParallelCyclotronTracker::ParallelCyclotronTracker(const Beamline& beamline,
111 PartBunchBase<double, 3>* bunch,
112 DataSink& ds,
113 const PartData& reference,
114 bool revBeam, bool revTrack,
115 int maxSTEPS,
116 Steppers::TimeIntegrator timeintegrator,
117 const int& numBunch,
118 const double& mbEta,
119 const double& mbPara,
120 const std::string& mbMode,
121 const std::string& mbBinning)
122 : Tracker(beamline, bunch, reference, revBeam, revTrack)
123 , bgf_m(nullptr)
124 , maxSteps_m(maxSTEPS)
125 , lastDumpedStep_m(0)
126 , myNode_m(Ippl::myNode())
127 , initialLocalNum_m(bunch->getLocalNum())
128 , initialTotalNum_m(bunch->getTotalNum())
129 , opalRing_m(nullptr)
130 , itsStepper_mp(nullptr)
131 , mode_m(TrackingMode::UNDEFINED)
132 , stepper_m(timeintegrator)
133 , rotation_m(3,3)
134{
135 itsBeamline = dynamic_cast<Beamline*>(beamline.clone());
136 itsDataSink = &ds;
137
138 if ( numBunch > 1 ) {
139 mbHandler_m = std::unique_ptr<MultiBunchHandler>(
140 new MultiBunchHandler(bunch, numBunch, mbEta,
141 mbPara, mbMode, mbBinning)
142 );
143 }
144
145 IntegrationTimer_m = IpplTimings::getTimer("Integration");
146 TransformTimer_m = IpplTimings::getTimer("Frametransform");
147 DumpTimer_m = IpplTimings::getTimer("Dump");
148 BinRepartTimer_m = IpplTimings::getTimer("Binaryrepart");
149 PluginElemTimer_m = IpplTimings::getTimer("PluginElements");
150 DelParticleTimer_m = IpplTimings::getTimer("DeleteParticles");
151
152 setTrackingMode();
153}
154
159ParallelCyclotronTracker::~ParallelCyclotronTracker() {
160 for (Component* component : myElements) {
161 delete(component);
162 }
163 for (auto fd : FieldDimensions) {
164 delete(fd);
165 }
166 delete itsBeamline;
167 // delete opalRing_m;
168}
169
170
171void ParallelCyclotronTracker::bgf_main_collision_test() {
172 if (!bgf_m) return;
173
174 Inform msg("bgf_main_collision_test ");
175
179
180 Vector_t intecoords = 0.0;
181
182 // This has to match the dT in the rk4 pusher
183 double dtime = itsBunch_m->getdT() * getHarmonicNumber();
184
185 int triId = 0;
186 for (size_t i = 0; i < itsBunch_m->getLocalNum(); i++) {
187 int res = bgf_m->partInside(itsBunch_m->R[i], itsBunch_m->P[i],
188 dtime, intecoords, triId);
189 if (res >= 0) {
190 lossDs_m->addParticle(OpalParticle(itsBunch_m->ID[i],
191 itsBunch_m->R[i], itsBunch_m->P[i],
192 itsBunch_m->getT(),
193 itsBunch_m->Q[i], itsBunch_m->M[i]),
194 std::make_pair(turnnumber_m, itsBunch_m->bunchNum[i]));
195 itsBunch_m->Bin[i] = -1;
196 *gmsgALL << level4 << "* Particle " << itsBunch_m->ID[i]
197 << " lost on boundary geometry" << endl;
198 }
199 }
200}
201
202// only used for dumping into stat file
203void ParallelCyclotronTracker::dumpAngle(const double& theta,
204 double& prevAzimuth,
205 double& azimuth,
206 const short& bunchNr) {
207 if ( prevAzimuth < 0.0 ) { // only at first occurrence
208 double plus = 0.0;
209 if ( OpalData::getInstance()->inRestartRun() ) {
210 plus = 360.0 * (turnnumber_m - bunchNr);
211 }
212 azimuth = theta + plus;
213 } else {
214 double dtheta = theta - prevAzimuth;
215 if ( dtheta < 0 ) {
216 dtheta += 360.0;
217 }
218 if ( dtheta > 180 ) { // rotating clockwise, reduce angle
219 dtheta -= 360;
220 }
221 azimuth += dtheta;
222 }
223 prevAzimuth = theta;
224}
225
226
227double ParallelCyclotronTracker::computeRadius(const Vector_t& meanR) const {
228 return Units::m2mm * std::sqrt(meanR(0) * meanR(0) + meanR(1) * meanR(1));
229}
230
231
232void ParallelCyclotronTracker::computePathLengthUpdate(std::vector<double>& dl,
233 const double& dt)
234{
235 // the last element in dotP is the dot-product over all particles
236 std::vector<double> dotP(dl.size());
237 if ( Options::psDumpFrame == DumpFrame::BUNCH_MEAN || isMultiBunch() ) {
238
239 for (unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) {
240 dotP[itsBunch_m->bunchNum[i]] += dot(itsBunch_m->P[i], itsBunch_m->P[i]);
241 }
242
243 allreduce(dotP.data(), dotP.size(), std::plus<double>());
244
245 // dot-product over all particles
246 double sum = std::accumulate(dotP.begin(), dotP.end() - 1, 0);
247 dotP.back() = sum / double(itsBunch_m->getTotalNum());
248
249 // bunch specific --> multi-bunches only
250 for (short b = 0; b < (short)dotP.size() - 1; ++b) {
251 dotP[b] /= double(itsBunch_m->getTotalNumPerBunch(b));
252 }
253
254 } else if ( itsBunch_m->getLocalNum() == 0 ) {
255 // here we are in DumpFrame::GLOBAL mode
256 dotP[0] = 0.0;
257 } else {
258 // here we are in DumpFrame::GLOBAL mode
259 dotP[0] = dot(itsBunch_m->P[0], itsBunch_m->P[0]);
260 }
261
262 for (size_t i = 0; i < dotP.size(); ++i) {
263 double const gamma = std::sqrt(1.0 + dotP[i]);
264 double const beta = std::sqrt(dotP[i]) / gamma;
265 dl[i] = Physics::c * dt * beta;
266 }
267}
268
269
275void ParallelCyclotronTracker::openFiles(size_t numFiles, std::string SfileName) {
276
277 for (unsigned int i=0; i<numFiles; i++) {
278 std::string SfileName2 = SfileName;
279 if (i<=2) {
280 SfileName2 += std::string("-Angle" + std::to_string(int(i)) + ".dat");
281 }
282 else {
283 // for single Particle Mode, output after each turn, to define matched initial phase ellipse.
284 SfileName2 += std::string("-afterEachTurn.dat");
285 }
286
287 outfTheta_m.emplace_back(new std::ofstream(SfileName2.c_str()));
288 outfTheta_m.back()->precision(8);
289 outfTheta_m.back()->setf(std::ios::scientific, std::ios::floatfield);
290 *outfTheta_m.back() << "# r [mm] beta_r*gamma "
291 << "theta [deg] beta_theta*gamma "
292 << "z [mm] beta_z*gamma" << std::endl;
293 }
294}
295
301void ParallelCyclotronTracker::closeFiles() {
302 for (auto & file : outfTheta_m) {
303 file->close();
304 }
305}
306
307
313void ParallelCyclotronTracker::visitCyclotron(const Cyclotron& cycl) {
314 *gmsg << "* ----------------------------- Cyclotron -------------------------------- *" << endl;
315
316 cycl_m = dynamic_cast<Cyclotron*>(cycl.clone());
317 myElements.push_back(cycl_m);
318
319 // Is this a Spiral Inflector Simulation? If yes, we'll give the user some
320 // useful information
321 spiral_flag = cycl_m->getSpiralFlag();
322 if (spiral_flag) {
323 *gmsg << endl << "* This is a Spiral Inflector Simulation! This means the following:" << endl;
324 *gmsg << "* 1.) It is up to the user to provide appropriate geometry, electric and magnetic fields!" << endl;
325 *gmsg << "* (Use BANDRF type cyclotron and use RFMAPFN to load both magnetic" << endl;
326 *gmsg << "* and electric fields, setting SUPERPOSE to an array of TRUE values.)" << endl;
327 *gmsg << "* 2.) For high currents it is strongly recommended to use the SAAMG fieldsolver," << endl;
328 *gmsg << "* FFT does not give the correct results (boundary conditions are missing)." << endl;
329 *gmsg << "* 3.) The whole geometry will be meshed and used for the fieldsolver." << endl;
330 *gmsg << "* There will be no transformations of the bunch into a local frame und consequently," << endl;
331 *gmsg << "* the problem will be treated non-relativistically!" << endl;
332 *gmsg << "* (This is not an issue for spiral inflectors as they are typically < 100 keV/amu.)" << endl;
333 *gmsg << endl << "* Note: For now, multi-bunch mode (MBM) needs to be de-activated for spiral inflector" << endl;
334 *gmsg << "* and space charge needs to be solved every time-step. numBunch_m and scSolveFreq are reset." << endl;
335 if (isMultiBunch()) {
336 mbHandler_m = nullptr;
337 }
338 }
339
340 // Fresh run (no restart):
341 if (!OpalData::getInstance()->inRestartRun()) {
342
343 // Get reference values from cyclotron element
344 referenceR = cycl_m->getRinit();
345 referenceTheta = cycl_m->getPHIinit();
346 referenceZ = cycl_m->getZinit();
347 referencePr = cycl_m->getPRinit();
348 referencePz = cycl_m->getPZinit();
349 referencePtot = itsReference.getGamma() * itsReference.getBeta();
350
351 // Calculate reference azimuthal (tangential) momentum from total-, z- and radial momentum:
352 float insqrt = referencePtot * referencePtot - \
353 referencePr * referencePr - referencePz * referencePz;
354
355 if (insqrt < 0) {
356 if (insqrt > -1.0e-10) {
357 referencePt = 0.0;
358 } else {
359 throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
360 "Pt imaginary");
361 }
362 } else {
363 referencePt = std::sqrt(insqrt);
364 }
365
366 if (referencePtot < 0.0) {
367 referencePt *= -1.0;
368 }
369 // End calculate referencePt
370
371 // Restart a run:
372 } else {
373
374 // If the user wants to save the restarted run in local frame,
375 // make sure the previous h5 file was local too
376 if (Options::psDumpFrame != DumpFrame::GLOBAL) {
377 if (!previousH5Local) {
378 throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
379 "You are trying a local restart from a global h5 file");
380 }
381 // Else, if the user wants to save the restarted run in global frame,
382 // make sure the previous h5 file was global too
383 } else {
384 if (previousH5Local) {
385 throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
386 "You are trying a global restart from a local h5 file");
387 }
388 }
389
390 // Adjust some of the reference variables from the h5 file
391 referenceR *= Units::mm2m;
392 referenceZ *= Units::mm2m;
393 referenceTheta *= Units::deg2rad;
394 referencePhi *= Units::deg2rad;
395 referencePsi *= Units::deg2rad;
396 referencePtot = bega;
397 }
398
399 cycl_m->checkInitialReferenceParticle(referenceR, referenceTheta, referenceZ);
400
401 sinRefTheta_m = std::sin(referenceTheta);
402 cosRefTheta_m = std::cos(referenceTheta);
403
404 *gmsg << endl;
405 *gmsg << "* Bunch global starting position:" << endl;
406 *gmsg << "* RINIT = " << referenceR << " [m]" << endl;
407 *gmsg << "* PHIINIT = " << referenceTheta * Units::rad2deg << " [deg]" << endl;
408 *gmsg << "* ZINIT = " << referenceZ << " [m]" << endl;
409 *gmsg << endl;
410 *gmsg << "* Bunch global starting momenta:" << endl;
411 *gmsg << "* Initial gamma = " << itsReference.getGamma() << endl;
412 *gmsg << "* Initial beta = " << itsReference.getBeta() << endl;
413 *gmsg << "* Reference total momentum = " << referencePtot << " [beta gamma]" << endl;
414 *gmsg << "* Reference azimuthal momentum (Pt) = " << referencePt << " [beta gamma]" << endl;
415 *gmsg << "* Reference radial momentum (Pr) = " << referencePr << " [beta gamma]" << endl;
416 *gmsg << "* Reference axial momentum (Pz) = " << referencePz << " [beta gamma]" << endl;
417 *gmsg << endl;
418
419 double sym = cycl_m->getSymmetry();
420 *gmsg << "* " << sym << "-fold field symmetry " << endl;
421
422 // ckr: this just returned the default value as defined in Component.h
423 // double rff = cycl_m->getRfFrequ(0);
424 // *gmsg << "* Rf frequency= " << rff << " [MHz]" << endl;
425
426 std::string fmfn = cycl_m->getFieldMapFN();
427 *gmsg << "* Field map file = '" << fmfn << "'" << endl;
428
429 std::string type = cycl_m->getCyclotronType();
430 *gmsg << "* Type of cyclotron = " << type << " " << endl;
431
432 double rmin = cycl_m->getMinR();
433 double rmax = cycl_m->getMaxR();
434 *gmsg << "* Radial aperture = " << rmin << " ... " << rmax << " [m]" << endl;
435
436 double zmin = cycl_m->getMinZ();
437 double zmax = cycl_m->getMaxZ();
438 *gmsg << "* Vertical aperture = " << zmin << " ... " << zmax << " [m]" << endl;
439
440 double h = cycl_m->getCyclHarm();
441 *gmsg << "* Number of trimcoils = " << cycl_m->getNumberOfTrimcoils() << endl;
442 *gmsg << "* Harmonic number h = " << h << endl;
443
444 std::vector<double> rffrequ = cycl_m->getRfFrequ();
445 *gmsg << "* RF frequency = " << Util::doubleVectorToString(rffrequ) << " [MHz]" << endl;
446
447 if (cycl_m->getBFieldType() == Cyclotron::BFieldType::BANDRF) {
448 std::vector<double> rfphi = cycl_m->getRfPhi();
449 for (size_t i = 0; i < rfphi.size(); ++i) {
450 rfphi[i] = rfphi[i] * Units::rad2deg;
451 }
452 *gmsg << "* RF inital phase = " << Util::doubleVectorToString(rfphi) << " [deg]" << endl;
453
454 std::vector<double> escale = cycl_m->getEScale();
455 *gmsg << "* E-field scale factor = " << Util::doubleVectorToString(escale) << endl;
456
457 std::vector<bool> superpose = cycl_m->getSuperpose();
458 *gmsg << "* Superpose electric field maps -> " << Util::boolVectorToUpperString(superpose) << endl;
459 }
460
461 // Read in cyclotron field maps
462 cycl_m->initialise(itsBunch_m, cycl_m->getBScale());
463
464 double BcParameter[8] = {};
465 BcParameter[0] = cycl_m->getRmin();
466 BcParameter[1] = cycl_m->getRmax();
467
468 // Store inner radius and outer radius of cyclotron field map in the list
469 buildupFieldList(BcParameter, ElementType::CYCLOTRON, cycl_m);
470}
471
477void ParallelCyclotronTracker::visitCCollimator(const CCollimator& coll) {
478 *gmsg << "* ----------------------------- Collimator ------------------------------- *" << endl;
479
480 CCollimator* elptr = dynamic_cast<CCollimator*>(coll.clone());
481 myElements.push_back(elptr);
482
483 *gmsg << "* Name = " << elptr->getName() << endl;
484
485 double xstart = elptr->getXStart();
486 *gmsg << "* XStart = " << xstart << " [m]" << endl;
487
488 double xend = elptr->getXEnd();
489 *gmsg << "* XEnd = " << xend << " [m]" << endl;
490
491 double ystart = elptr->getYStart();
492 *gmsg << "* YStart = " << ystart << " [m]" << endl;
493
494 double yend = elptr->getYEnd();
495 *gmsg << "* YEnd = " << yend << " [m]" << endl;
496
497 double zstart = elptr->getZStart();
498 *gmsg << "* ZStart = " << zstart << " [m]" << endl;
499
500 double zend = elptr->getZEnd();
501 *gmsg << "* ZEnd = " << zend << " [m]" << endl;
502
503 double width = elptr->getWidth();
504 *gmsg << "* Width = " << width << " [m]" << endl;
505
506 elptr->initialise(itsBunch_m);
507
508 double BcParameter[8] = {};
509
510 BcParameter[0] = xstart;
511 BcParameter[1] = xend;
512 BcParameter[2] = ystart;
513 BcParameter[3] = yend;
514 BcParameter[4] = width;
515
516 buildupFieldList(BcParameter, ElementType::CCOLLIMATOR, elptr);
517}
518
524void ParallelCyclotronTracker::visitCorrector(const Corrector& corr) {
525 *gmsg << "In Corrector; L= " << corr.getElementLength() << endl;
526 myElements.push_back(dynamic_cast<Corrector*>(corr.clone()));
527}
528
534void ParallelCyclotronTracker::visitDegrader(const Degrader& deg) {
535 *gmsg << "In Degrader; L= " << deg.getElementLength() << endl;
536 myElements.push_back(dynamic_cast<Degrader*>(deg.clone()));
537
538}
539
545void ParallelCyclotronTracker::visitDrift(const Drift& drift) {
546 *gmsg << "In drift L= " << drift.getElementLength() << endl;
547 myElements.push_back(dynamic_cast<Drift*>(drift.clone()));
548}
549
558
564void ParallelCyclotronTracker::visitOffset(const Offset& off) {
565 if (opalRing_m == nullptr)
566 throw OpalException(
567 "ParallelCylcotronTracker::visitOffset",
568 "Attempt to place an offset when Ring not defined");
569 opalRing_m->appendElement(off);
570}
571
577void ParallelCyclotronTracker::visitMarker(const Marker& marker) {
578 // *gmsg << "In Marker; L= " << marker.getElementLength() << endl;
579 myElements.push_back(dynamic_cast<Marker*>(marker.clone()));
580 // Do nothing.
581}
582
588void ParallelCyclotronTracker::visitMonitor(const Monitor& corr) {
589 // *gmsg << "In Monitor; L= " << corr.getElementLength() << endl;
590 myElements.push_back(dynamic_cast<Monitor*>(corr.clone()));
591 // applyDrift(flip_s * corr.getElementLength());
592}
593
599void ParallelCyclotronTracker::visitMultipole(const Multipole& mult) {
600 *gmsg << "In Multipole; L= " << mult.getElementLength() << " however the element is missing " << endl;
601 myElements.push_back(dynamic_cast<Multipole*>(mult.clone()));
602}
603
609void ParallelCyclotronTracker::visitMultipoleT(const MultipoleT& multT) {
610 *gmsg << "Adding MultipoleT" << endl;
611 if (opalRing_m != nullptr) {
612 opalRing_m->appendElement(multT);
613 } else {
614 throw OpalException("ParallelCyclotronTracker::visitMultipoleT",
615 "Need to define a RINGDEFINITION to use MultipoleT element");
616 }
617 myElements.push_back(dynamic_cast<MultipoleT*>(multT.clone()));
618}
619
625void ParallelCyclotronTracker::visitOutputPlane(const OutputPlane& plane) {
626 if (opalRing_m == nullptr) {
627 throw OpalException("ParallelCylcotronTracker::visitOutputPlane",
628 "Attempt to place an OutputPlane when Ring not defined");
629 }
630
631 OutputPlane* elptr = dynamic_cast<OutputPlane*>(plane.clone());
632 elptr->setGlobalFieldMap(opalRing_m);
633 myElements.push_back(elptr);
634 elptr->initialise(itsBunch_m);
635
636 double BcParameter[8] = {};
637
638 Vector_t centre = elptr->getCentre();
639 Vector_t norm = elptr->getNormal();
640 double width = elptr->getHorizontalExtent();
641 BcParameter[0] = centre[0]-width*norm[1];
642 BcParameter[1] = centre[0]+width*norm[1];
643 BcParameter[2] = centre[1]-width*norm[0];
644 BcParameter[3] = centre[1]+width*norm[0];
645 BcParameter[4] = Units::mm2m; // thickness, not used
646 buildupFieldList(BcParameter, ElementType::OUTPUTPLANE, elptr);
647}
648
649
655void ParallelCyclotronTracker::visitProbe(const Probe& prob) {
656 *gmsg << "* ----------------------------- Probe ------------------------------------ *" << endl;
657
658 Probe* elptr = dynamic_cast<Probe*>(prob.clone());
659 myElements.push_back(elptr);
660
661 *gmsg << "* Name = " << elptr->getName() << endl;
662
663 double xstart = elptr->getXStart();
664 *gmsg << "* XStart = " << xstart << " [m]" << endl;
665
666 double xend = elptr->getXEnd();
667 *gmsg << "* XEnd = " << xend << " [m]" << endl;
668
669 double ystart = elptr->getYStart();
670 *gmsg << "* YStart = " << ystart << " [m]" << endl;
671
672 double yend = elptr->getYEnd();
673 *gmsg << "* YEnd = " << yend << " [m]" << endl;
674
675 // initialise, do nothing
676 elptr->initialise(itsBunch_m);
677
678 double BcParameter[8] = {};
679
680 BcParameter[0] = xstart;
681 BcParameter[1] = xend;
682 BcParameter[2] = ystart;
683 BcParameter[3] = yend;
684 BcParameter[4] = 1 ; // width
685
686 buildupFieldList(BcParameter, ElementType::PROBE, elptr);
687}
688
694void ParallelCyclotronTracker::visitRBend(const RBend& bend) {
695 *gmsg << "In RBend; L= " << bend.getElementLength() << " however the element is missing " << endl;
696 myElements.push_back(dynamic_cast<RBend*>(bend.clone()));
697}
698
704void ParallelCyclotronTracker::visitRFCavity(const RFCavity& as) {
705 *gmsg << "* ----------------------------- RFCavity --------------------------------- * " << endl;
706
707 RFCavity* elptr = dynamic_cast<RFCavity*>(as.clone());
708 myElements.push_back(elptr);
709
710 if ( elptr->getCavityType() != CavityType::SGSW ) {
711 throw OpalException("ParallelCyclotronTracker::visitRFCavity",
712 "\"" + elptr->getCavityTypeString() + "\" is not valid \"TYPE\" for RFCavity.\n"
713 "The ParallelCyclotronTracker can only play with cyclotron type RF system currently...");
714 }
715
716 *gmsg << "* Name = " << elptr->getName() << endl;
717
718 std::string fmfn = elptr->getFieldMapFN();
719 *gmsg << "* RF Field map file = '" << fmfn << "'" << endl;
720
721 double rmin = elptr->getRmin();
722 *gmsg << "* Minimal radius of cavity = " << rmin << " [m]" << endl;
723
724 double rmax = elptr->getRmax();
725 *gmsg << "* Maximal radius of cavity = " << rmax << " [m]" << endl;
726
727 double rff = elptr->getCycFrequency();
728 *gmsg << "* RF frequency (2*pi*f) = " << rff << " [rad/s]" << endl;
729
730 double angle = elptr->getAzimuth();
731 *gmsg << "* Cavity azimuth position = " << angle * Units::rad2deg << " [deg] " << endl;
732
733 double gap = elptr->getGapWidth();
734 *gmsg << "* Cavity gap width = " << gap << " [m] " << endl;
735
736 double pdis = elptr->getPerpenDistance();
737 *gmsg << "* Cavity Shift distance = " << pdis << " [m] " << endl;
738
739 double phi0 = elptr->getPhi0();
740 *gmsg << "* Initial RF phase (t=0) = " << phi0 * Units::rad2deg << " [deg] " << endl;
741
742 /*
743 Setup time dependence and in case of no
744 timedependence use a polynom with a_0 = 1 and a_k = 0, k = 1,2,3.
745 */
746
747 std::shared_ptr<AbstractTimeDependence> freq_atd = nullptr;
748 std::shared_ptr<AbstractTimeDependence> ampl_atd = nullptr;
749 std::shared_ptr<AbstractTimeDependence> phase_atd = nullptr;
750
751 dvector_t unityVec;
752 unityVec.push_back(1.);
753 unityVec.push_back(0.);
754 unityVec.push_back(0.);
755 unityVec.push_back(0.);
756
757 std::string frequencyModelName = elptr->getFrequencyModelName();
758 std::string amplitudeModelName = elptr->getAmplitudeModelName();
759 std::string phaseModelName = elptr->getPhaseModelName();
760
761 if (!frequencyModelName.empty()) {
762 freq_atd = AbstractTimeDependence::getTimeDependence(frequencyModelName);
763 *gmsg << "* Variable frequency RF Model name " << frequencyModelName << endl;
764 } else {
765 freq_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
766 }
767
768 if (!amplitudeModelName.empty()) {
769 ampl_atd = AbstractTimeDependence::getTimeDependence(amplitudeModelName);
770 *gmsg << "* Variable amplitude RF Model name " << amplitudeModelName << endl;
771 } else {
772 ampl_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
773 }
774
775 if (!phaseModelName.empty()) {
776 phase_atd = AbstractTimeDependence::getTimeDependence(phaseModelName);
777 *gmsg << "* Variable phase RF Model name " << phaseModelName << endl;
778 } else {
779 phase_atd = std::shared_ptr<AbstractTimeDependence>(new PolynomialTimeDependence(unityVec));
780 }
781 // read cavity voltage profile data from file.
782 elptr->initialise(itsBunch_m, freq_atd, ampl_atd, phase_atd);
783
784 double BcParameter[8] = {};
785
786 BcParameter[0] = rmin;
787 BcParameter[1] = rmax;
788 BcParameter[2] = pdis;
789 BcParameter[3] = angle;
790
791 buildupFieldList(BcParameter, ElementType::RFCAVITY, elptr);
792}
793
798void ParallelCyclotronTracker::visitRing(const Ring& ring) {
799 *gmsg << "* ----------------------------- Ring ------------------------------------- *" << endl;
800
801 delete opalRing_m;
802
803 opalRing_m = dynamic_cast<Ring*>(ring.clone());
804
805 myElements.push_back(opalRing_m);
806
807 opalRing_m->initialise(itsBunch_m);
808
809 referenceR = opalRing_m->getBeamRInit();
810 referencePr = opalRing_m->getBeamPRInit();
811 referenceTheta = opalRing_m->getBeamPhiInit();
812 beamInitialRotation = opalRing_m->getBeamThetaInit();
813
814 if (!Util::angleBetweenAngles(referenceTheta, -Physics::pi, Physics::pi)) {
815 throw OpalException("Error in ParallelCyclotronTracker::visitRing",
816 "BEAM_PHIINIT is out of [-180, 180)");
817 }
818
819 referenceZ = 0.0;
820 referencePz = 0.0;
821
822 referencePtot = itsReference.getGamma() * itsReference.getBeta();
823 referencePt = std::sqrt(referencePtot * referencePtot - referencePr * referencePr);
824
825 if (referencePtot < 0.0)
826 referencePt *= -1.0;
827
828 sinRefTheta_m = std::sin(referenceTheta);
829 cosRefTheta_m = std::cos(referenceTheta);
830
831 double BcParameter[8] = {}; // zero initialise array
832
833 buildupFieldList(BcParameter, ElementType::RING, opalRing_m);
834
835 // Finally print some diagnostic
836 *gmsg << "* Initial beam radius = " << referenceR << " [m] " << endl;
837 *gmsg << "* Initial gamma = " << itsReference.getGamma() << endl;
838 *gmsg << "* Initial beta = " << itsReference.getBeta() << endl;
839 *gmsg << "* Total reference momentum = " << referencePtot << " [beta gamma]" << endl;
840 *gmsg << "* Reference azimuthal momentum = " << referencePt << " [beta gamma]" << endl;
841 *gmsg << "* Reference radial momentum = " << referencePr << " [beta gamma]" << endl;
842 *gmsg << "* " << opalRing_m->getSymmetry() << " fold field symmetry " << endl;
843 *gmsg << "* Harmonic number h = " << opalRing_m->getHarmonicNumber() << " " << endl;
844}
845
851void ParallelCyclotronTracker::visitSBend(const SBend& bend) {
852 *gmsg << "In SBend; L = " << bend.getElementLength() << " however the element is missing " << endl;
853 myElements.push_back(dynamic_cast<SBend*>(bend.clone()));
854}
855
856void ParallelCyclotronTracker::visitSBend3D(const SBend3D& bend) {
857 *gmsg << "Adding SBend3D" << endl;
858 if (opalRing_m != nullptr)
859 opalRing_m->appendElement(bend);
860 else
861 throw OpalException("ParallelCyclotronTracker::visitSBend3D",
862 "Need to define a RINGDEFINITION to use SBend3D element");
863}
864
865void ParallelCyclotronTracker::visitScalingFFAMagnet(const ScalingFFAMagnet& bend) {
866 *gmsg << "Adding ScalingFFAMagnet" << endl;
867 if (opalRing_m != nullptr) {
868 ScalingFFAMagnet* newBend = bend.clone(); // setup the end field, if required
869 newBend->setupEndField();
870 opalRing_m->appendElement(*newBend);
871 } else {
872 throw OpalException("ParallelCyclotronTracker::visitScalingFFAMagnet",
873 "Need to define a RINGDEFINITION to use ScalingFFAMagnet element");
874 }
875}
876
882void ParallelCyclotronTracker::visitSeptum(const Septum& sept) {
883 *gmsg << "* ----------------------------- Septum ----------------------------------- *" << endl;
884
885 Septum* elptr = dynamic_cast<Septum*>(sept.clone());
886 myElements.push_back(elptr);
887
888 *gmsg << "* Name = " << elptr->getName() << endl;
889
890 double xstart = elptr->getXStart();
891 *gmsg << "* XStart = " << xstart << " [m]" << endl;
892
893 double xend = elptr->getXEnd();
894 *gmsg << "* XEnd = " << xend << " [m]" << endl;
895
896 double ystart = elptr->getYStart();
897 *gmsg << "* YStart = " << ystart << " [m]" << endl;
898
899 double yend = elptr->getYEnd();
900 *gmsg << "* YEnd = " << yend << " [m]" << endl;
901
902 double width = elptr->getWidth();
903 *gmsg << "* Width = " << width << " [m]" << endl;
904
905 // initialise, do nothing
906 elptr->initialise(itsBunch_m);
907
908 double BcParameter[8] = {};
909
910 BcParameter[0] = xstart;
911 BcParameter[1] = xend;
912 BcParameter[2] = ystart;
913 BcParameter[3] = yend;
914 BcParameter[4] = width;
915
916 buildupFieldList(BcParameter, ElementType::SEPTUM, elptr);
917}
918
924void ParallelCyclotronTracker::visitSolenoid(const Solenoid& solenoid) {
925 myElements.push_back(dynamic_cast<Solenoid*>(solenoid.clone()));
926 Component* elptr = *(--myElements.end());
927 if (!elptr->hasAttribute("ELEMEDGE")) {
928 *gmsg << "Solenoid: no position of the element given" << endl;
929 return;
930 }
931}
932
938void ParallelCyclotronTracker::visitStripper(const Stripper& stripper) {
939 *gmsg << "* ----------------------------- Stripper --------------------------------- *" << endl;
940
941 Stripper* elptr = dynamic_cast<Stripper*>(stripper.clone());
942 myElements.push_back(elptr);
943
944 *gmsg << "* Name = " << elptr->getName() << endl;
945
946 double xstart = elptr->getXStart();
947 *gmsg << "* XStart = " << xstart << " [m]" << endl;
948
949 double xend = elptr->getXEnd();
950 *gmsg << "* XEnd = " << xend << " [m]" << endl;
951
952 double ystart = elptr->getYStart();
953 *gmsg << "* YStart = " << ystart << " [m]" << endl;
954
955 double yend = elptr->getYEnd();
956 *gmsg << "* YEnd = " << yend << " [m]" << endl;
957
958 double opcharge = elptr->getOPCharge();
959 *gmsg << "* Charge of outcoming particle = +e * " << opcharge << endl;
960
961 double opmass = elptr->getOPMass();
962 *gmsg << "* Mass of the outcoming particle = " << opmass << " [GeV/c^2]" << endl;
963
964 bool stop = elptr->getStop();
965 *gmsg << std::boolalpha
966 << "* Particles stripped will be deleted after interaction -> "
967 << stop << endl;
968
969 elptr->initialise(itsBunch_m);
970
971 double BcParameter[8] = {};
972
973 BcParameter[0] = xstart;
974 BcParameter[1] = xend;
975 BcParameter[2] = ystart;
976 BcParameter[3] = yend;
977 BcParameter[4] = 1; // width
978 BcParameter[5] = opcharge;
979 BcParameter[6] = opmass;
980
981 buildupFieldList(BcParameter, ElementType::STRIPPER, elptr);
982}
983
989void ParallelCyclotronTracker::visitVacuum(const Vacuum& vac) {
990 *gmsg << "* ----------------------------- Vacuum ----------------------------------- *" << endl;
991
992 Vacuum* elptr = dynamic_cast<Vacuum*>(vac.clone());
993 myElements.push_back(elptr);
994
995 double BcParameter[8] = {};
996
997 std::string gas = elptr->getResidualGasName();
998 *gmsg << "* Residual gas = " << gas << endl;
999
1000 double pressure = elptr->getPressure();
1001 if ( std::filesystem::exists(elptr->getPressureMapFN()) ) {
1002 std::string pmfn = elptr->getPressureMapFN();
1003 *gmsg << "* Pressure field map file = '" << pmfn << "'" << endl;
1004 *gmsg << "* Default pressure = " << pressure << " [mbar]" << endl;
1005 } else {
1006 *gmsg << "* Pressure = " << pressure << " [mbar]" << endl;
1007 }
1008 double pscale = elptr->getPScale();
1009
1010 double temperature = elptr->getTemperature();
1011 *gmsg << "* Temperature = " << temperature << " [K]" << endl;
1012
1013 bool stop = elptr->getStop();
1014 *gmsg << std::boolalpha
1015 << "* Particles stripped will be deleted after interaction -> "
1016 << stop << endl;
1017
1018 elptr->initialise(itsBunch_m);
1019
1020 BcParameter[0] = pressure;
1021 BcParameter[1] = pscale;
1022 BcParameter[2] = temperature;
1023
1024 buildupFieldList(BcParameter, ElementType::VACUUM, elptr);
1025}
1026
1032void ParallelCyclotronTracker::visitVariableRFCavity(const VariableRFCavity& cav) {
1033 *gmsg << "Adding Variable RF Cavity" << endl;
1034 if (opalRing_m != nullptr)
1035 opalRing_m->appendElement(cav);
1036 else
1037 throw OpalException("ParallelCyclotronTracker::visitVariableRFCavity",
1038 "Need to define a RINGDEFINITION to use VariableRFCavity element");
1039}
1040
1046void ParallelCyclotronTracker::visitVariableRFCavityFringeField
1047 (const VariableRFCavityFringeField& cav) {
1048 *gmsg << "Adding Variable RF Cavity with Fringe Field" << endl;
1049 if (opalRing_m != nullptr)
1050 opalRing_m->appendElement(cav);
1051 else
1052 throw OpalException("ParallelCyclotronTracker::visitVariableRFCavityFringeField",
1053 "Need to define a RINGDEFINITION to use VariableRFCavity element");
1054}
1055
1061void ParallelCyclotronTracker::visitVerticalFFAMagnet(const VerticalFFAMagnet& mag) {
1062 *gmsg << "Adding Vertical FFA Magnet" << endl;
1063 if (opalRing_m != nullptr)
1064 opalRing_m->appendElement(mag);
1065 else
1066 throw OpalException("ParallelCyclotronTracker::visitVerticalFFAMagnet",
1067 "Need to define a RINGDEFINITION to use VerticalFFAMagnet element");
1068}
1069
1070
1078void ParallelCyclotronTracker::buildupFieldList(double BcParameter[], ElementType elementType, Component *elptr) {
1079 beamline_list::iterator sindex;
1080
1081 type_pair *localpair = new type_pair();
1082 localpair->first = elementType;
1083
1084 for (int i = 0; i < 8; i++)
1085 *(((localpair->second).first) + i) = *(BcParameter + i);
1086
1087 (localpair->second).second = elptr;
1088
1089 // always put cyclotron as the first element in the list.
1090 if (elementType == ElementType::RING || elementType == ElementType::CYCLOTRON) {
1091 sindex = FieldDimensions.begin();
1092 } else {
1093 sindex = FieldDimensions.end();
1094 }
1095 FieldDimensions.insert(sindex, localpair);
1096}
1097
1103void ParallelCyclotronTracker::visitBeamline(const Beamline& bl) {
1104 const FlaggedBeamline* fbl = static_cast<const FlaggedBeamline*>(&bl);
1105 fbl->iterate(*this, false);//*dynamic_cast<BeamlineVisitor *>(this), false);
1106}
1107
1108void ParallelCyclotronTracker::checkNumPart(std::string s) {
1109 int nlp = itsBunch_m->getLocalNum();
1110 int minnlp = 0;
1111 int maxnlp = 111111;
1112 reduce(nlp, minnlp, OpMinAssign());
1113 reduce(nlp, maxnlp, OpMaxAssign());
1114 *gmsg << s << "min local particle number: "<< minnlp << endl;
1115 *gmsg << "* max local particle number: " << maxnlp << endl;
1116}
1117
1118
1119void ParallelCyclotronTracker::execute() {
1120 /*
1121 Initialize common variables and structures
1122 for the integrators
1123 */
1124 if (OpalData::getInstance()->inRestartRun()) {
1125 OpalData::getInstance()->setOpenMode(OpalData::OpenMode::APPEND);
1126 }
1127
1128 step_m = 0;
1129 restartStep0_m = 0;
1130 turnnumber_m = 1;
1131 azimuth_m = -1.0;
1132 prevAzimuth_m = -1.0;
1133
1134 // Record how many bunches have already been injected. ONLY FOR MPM
1135 if (isMultiBunch())
1136 mbHandler_m->setNumBunch(itsBunch_m->getNumBunch());
1137
1138 itsBeamline->accept(*this);
1139 if (opalRing_m != nullptr)
1140 opalRing_m->lockRing();
1141
1142 // Display the selected elements
1143 *gmsg << "* ------------------------------------------------------------------------ *" << endl;
1144 *gmsg << "* The selected Beam line elements are :" << endl;
1145 for (auto fd : FieldDimensions) {
1146 ElementType type = fd->first;
1147 *gmsg << "* -> " << ElementBase::getTypeString(type) << endl;
1148 if (type == ElementType::RFCAVITY) {
1149 RFCavity* cav = static_cast<RFCavity*>((fd->second).second);
1150 CavityCrossData ccd = {cav, cav->getSinAzimuth(), cav->getCosAzimuth(), cav->getPerpenDistance()};
1151 cavCrossDatas_m.push_back(ccd);
1152 } else if ( type == ElementType::CCOLLIMATOR ||
1153 type == ElementType::OUTPUTPLANE ||
1154 type == ElementType::PROBE ||
1155 type == ElementType::SEPTUM ||
1156 type == ElementType::STRIPPER) {
1157 PluginElement* element = static_cast<PluginElement*>((fd->second).second);
1158 pluginElements_m.push_back(element);
1159 }
1160 }
1161 *gmsg << "* ------------------------------------------------------------------------ *" << endl << endl;
1162
1163 // Get BoundaryGeometry that is already initialized
1164 bgf_m = OpalData::getInstance()->getGlobalGeometry();
1165 if (bgf_m)
1166 lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(bgf_m->getOpalName(),!Options::asciidump));
1167
1168 // External field arrays for dumping
1169 for (int k = 0; k < 2; k++) {
1170 FDext_m[k] = Vector_t(0.0, 0.0, 0.0);
1171 }
1172 extE_m = Vector_t(0.0, 0.0, 0.0);
1173 extB_m = Vector_t(0.0, 0.0, 0.0);
1174 DumpFields::writeFields((((*FieldDimensions.begin())->second).second));
1175 DumpEMFields::writeFields((((*FieldDimensions.begin())->second).second));
1176
1177 function_t func = std::bind(&ParallelCyclotronTracker::getFieldsAtPoint,
1178 this,
1179 std::placeholders::_1,
1180 std::placeholders::_2,
1181 std::placeholders::_3,
1182 std::placeholders::_4);
1183
1184 switch ( stepper_m ) {
1185 case Steppers::TimeIntegrator::LF2: {
1186 *gmsg << "* 2nd order Leap-Frog integrator" << endl;
1187 itsStepper_mp.reset(new LF2<function_t>(func));
1188 break;
1189 }
1190 case Steppers::TimeIntegrator::MTS: {
1191 *gmsg << "* Multiple time stepping (MTS) integrator" << endl;
1192 break;
1193 }
1194 case Steppers::TimeIntegrator::RK4:
1195 default: {
1196 *gmsg << "* 4th order Runge-Kutta integrator" << endl;
1197 itsStepper_mp.reset(new RK4<function_t>(func));
1198 break;
1199 }
1200 }
1201
1202 if ( stepper_m == Steppers::TimeIntegrator::MTS) {
1203 MtsTracker();
1204 } else {
1205 GenericTracker();
1206 }
1207
1208 *gmsg << "* ------------------------------------------------------------------------ *" << endl;
1209 *gmsg << "* Finalizing i.e. write data and close files :" << endl;
1210 for (auto fd : FieldDimensions) {
1211 ((fd->second).second)->finalise();
1212 }
1213 if (bgf_m && lossDs_m) {
1214 lossDs_m->save();
1215 }
1216 *gmsg << "* ------------------------------------------------------------------------ *" << endl;
1217}
1218
1219
1220void ParallelCyclotronTracker::MtsTracker() {
1221 /*
1222 * variable unit meaning
1223 * ------------------------------------------------
1224 * t [s] time
1225 * dt [s] time step
1226 * oldReferenceTheta [rad] azimuth angle
1227 * itsBunch_m->R [m] particle position
1228 *
1229 */
1230
1231 double t = 0, dt = 0, oldReferenceTheta = 0;
1232 std::tie(t, dt, oldReferenceTheta) = initializeTracking_m();
1233
1234 int const numSubsteps = std::max(Options::mtsSubsteps, 1);
1235 *gmsg << "* MTS: Number of substeps per step is " << numSubsteps << endl;
1236
1237 double const dt_inner = dt / double(numSubsteps);
1238 *gmsg << "* MTS: The inner time step is therefore " << dt_inner * Units::s2ns << " [ns]" << endl;
1239
1240// int SteptoLastInj = itsBunch_m->getSteptoLastInj();
1241
1242 bool flagTransition = false; // flag to determine when to transit from single-bunch to multi-bunches mode
1243
1244 *gmsg << "* ---------------------- Start tracking ---------------------------------- *" << endl;
1245
1246 if (itsBunch_m->hasFieldSolver()) {
1247 computeSpaceChargeFields_m();
1248 }
1249
1250 for (; (step_m < maxSteps_m) && (itsBunch_m->getTotalNum()>0); step_m++) {
1251
1252 bool finishedTurn = false;
1253
1254 if (step_m % Options::sptDumpFreq == 0) {
1255 singleParticleDump();
1256 }
1257
1258 Ippl::Comm->barrier();
1259
1260 // First half kick from space charge force
1261 if (itsBunch_m->hasFieldSolver()) {
1262 kick(0.5 * dt);
1263 }
1264
1265 // Substeps for external field integration
1266 for (int n = 0; n < numSubsteps; ++n) {
1267 borisExternalFields(dt_inner);
1268 }
1269 // bunch injection
1270 // TODO: Where is correct location for this piece of code? Beginning/end of step? Before field solve?
1271 injectBunch(flagTransition);
1272
1273 if ( itsBunch_m->hasFieldSolver() ) {
1274 computeSpaceChargeFields_m();
1275 } else {
1276 // if field solver is not available , only update bunch, to transfer particles between nodes if needed,
1277 // reset parameters such as LocalNum, initialTotalNum_m.
1278 // INFOMSG("No space charge Effects are included"<<endl;);
1279 if ((step_m % Options::repartFreq * 100) == 0) { //TODO: why * 100?
1280 Vector_t const meanP = calcMeanP();
1281 double const phi = calculateAngle(meanP(0), meanP(1)) - 0.5 * Physics::pi;
1282 Vector_t const meanR = calcMeanR();
1283 globalToLocal(itsBunch_m->R, phi, meanR);
1284 itsBunch_m->updateNumTotal();
1285 repartition();
1286 localToGlobal(itsBunch_m->R, phi, meanR);
1287 }
1288 }
1289
1290 // Second half kick from space charge force
1291 if (itsBunch_m->hasFieldSolver()) {
1292 kick(0.5 * dt);
1293 }
1294 // recalculate bingamma and reset the BinID for each particles according to its current gamma
1295 if (isMultiBunch() && (step_m % Options::rebinFreq == 0)) {
1296 mbHandler_m->updateParticleBins(itsBunch_m);
1297 }
1298
1299 // dump some data after one push in single particle tracking
1300 if ( mode_m == TrackingMode::SINGLE ) {
1301 unsigned int i = 0;
1302
1303 double temp_meanTheta = calculateAngle2(itsBunch_m->R[i](0),
1304 itsBunch_m->R[i](1)); // [-pi, pi]
1305
1306 dumpThetaEachTurn_m(t, itsBunch_m->R[i], itsBunch_m->P[i],
1307 temp_meanTheta, finishedTurn);
1308
1309 dumpAzimuthAngles_m(t, itsBunch_m->R[i], itsBunch_m->P[i],
1310 oldReferenceTheta, temp_meanTheta);
1311
1312 oldReferenceTheta = temp_meanTheta;
1313 } else if ( mode_m == TrackingMode::BUNCH ) {
1314 // both for single bunch and multi-bunch
1315 // avoid dump at the first step
1316 // finishedTurn has not been changed in first push
1317 if ( isTurnDone() ) {
1318 ++turnnumber_m;
1319 finishedTurn = true;
1320
1321 *gmsg << endl;
1322 *gmsg << "*** Finished turn " << turnnumber_m - 1
1323 << ", Total number of live particles: "
1324 << itsBunch_m->getTotalNum() << endl;
1325 }
1326
1327 // Recalculate bingamma and reset the BinID for each particles according to its current gamma
1328 if (isMultiBunch() && (step_m % Options::rebinFreq == 0)) {
1329 mbHandler_m->updateParticleBins(itsBunch_m);
1330 }
1331 }
1332
1333 // printing + updating bunch parameters + updating t
1334 update_m(t, dt, finishedTurn);
1335 }
1336
1337 // Some post-integration stuff
1338 *gmsg << endl;
1339 *gmsg << "* ---------------------- DONE TRACKING PARTICLES ------------------------- *" << endl;
1340
1341 //FIXME
1342 dvector_t Ttime = dvector_t();
1343 dvector_t Tdeltr = dvector_t();
1344 dvector_t Tdeltz = dvector_t();
1345 ivector_t TturnNumber = ivector_t();
1346
1347 finalizeTracking_m(Ttime, Tdeltr, Tdeltz, TturnNumber);
1348}
1349
1350
1351void ParallelCyclotronTracker::GenericTracker() {
1352 /*
1353 * variable unit meaning
1354 * ------------------------------------------------
1355 * t [s] time
1356 * dt [s] time step
1357 * oldReferenceTheta [rad] azimuth angle
1358 * itsBunch_m->R [m] particle position
1359 *
1360 */
1361 // Generic Tracker that has three modes defined by timeIntegrator_m:
1362 // 0 --- RK-4 (default)
1363 // 1 --- LF-2
1364 // (2 --- MTS ... not yet implemented in generic tracker)
1365 // mbHandler_m->getNumBunch() determines the number of bunches in multibunch mode (MBM, 1 for OFF)
1366 // Total number of particles determines single particle mode (SPM, 1 particle) or
1367 // Static Equilibrium Orbit Mode (SEO, 2 particles)
1368
1369 double t = 0, dt = 0, oldReferenceTheta = 0;
1370 std::tie(t, dt, oldReferenceTheta) = initializeTracking_m();
1371
1372 // If initialTotalNum_m = 2, trigger SEO mode and prepare for transverse tuning calculation
1373 // Where is the IF here? -DW
1374 dvector_t Ttime, Tdeltr, Tdeltz;
1375 ivector_t TturnNumber;
1376
1377 // Apply the plugin elements: probe, collimator, stripper, septum once before first step
1378 bool flagNeedUpdate = applyPluginElements(dt);
1379
1380 // Destroy particles if they are marked as Bin = -1 in the plugin elements
1381 // or out of global aperture
1382 deleteParticle(flagNeedUpdate);
1383
1384 /********************************
1385 * Main integration loop *
1386 ********************************/
1387 *gmsg << endl;
1388 *gmsg << "* ---------------------- Start tracking ---------------------------------- *" << endl;
1389
1390 for (; (step_m < maxSteps_m) && (itsBunch_m->getTotalNum()>0); step_m++) {
1391
1392 bool finishedTurn = false;
1393
1394 switch (mode_m) {
1395 case TrackingMode::SEO: {
1396 // initialTotalNum_m == 2
1397 seoMode_m(t, dt, finishedTurn, Ttime, Tdeltr, Tdeltz, TturnNumber);
1398 break;
1399 }
1400 case TrackingMode::SINGLE: {
1401 // initialTotalNum_m == 1
1402 singleMode_m(t, dt, finishedTurn, oldReferenceTheta);
1403 break;
1404 }
1405 case TrackingMode::BUNCH: {
1406 // initialTotalNum_m > 2
1407 // Start Tracking for number of particles > 2 (i.e. not single and not SEO mode)
1408 bunchMode_m(t, dt, finishedTurn);
1409 break;
1410 }
1411 case TrackingMode::UNDEFINED: {
1412 default:
1413 throw OpalException("ParallelCyclotronTracker::GenericTracker()",
1414 "No such tracking mode.");
1415 }
1416 }
1417 // Update bunch and some parameters and output some info
1418 update_m(t, dt, finishedTurn);
1419
1420 } // end for: the integration is DONE after maxSteps_m steps or if all particles are lost!
1421
1422 // Some post-integration stuff
1423 *gmsg << endl;
1424 *gmsg << "* ---------------------- DONE TRACKING PARTICLES ------------------------- *" << endl;
1425
1426 finalizeTracking_m(Ttime, Tdeltr, Tdeltz, TturnNumber);
1427}
1428
1429bool ParallelCyclotronTracker::getFieldsAtPoint(const double& t, const size_t& Pindex, Vector_t& Efield, Vector_t& Bfield) {
1430
1431 bool outOfBound = this->computeExternalFields_m(Pindex, t, Efield, Bfield);
1432
1433 // For runs without space charge effects, override this step to save time
1434 if (itsBunch_m->hasFieldSolver()) {
1435
1436 // Don't do for reference particle
1437 if (itsBunch_m->ID[Pindex] != 0) {
1438
1439 // add external Field and self space charge field
1440 Efield += itsBunch_m->Ef[Pindex];
1441 Bfield += itsBunch_m->Bf[Pindex];
1442 }
1443 }
1444
1445 return outOfBound;
1446}
1447
1448
1459bool ParallelCyclotronTracker::checkGapCross(Vector_t Rold, Vector_t Rnew,
1460 RFCavity* rfcavity, double& Dold) {
1461 bool flag = false;
1462 double sinx = rfcavity->getSinAzimuth();
1463 double cosx = rfcavity->getCosAzimuth();
1464
1465 double PerpenDistance = rfcavity->getPerpenDistance();
1466 double distNew = (Rnew[0] * sinx - Rnew[1] * cosx) - PerpenDistance;
1467 double distOld = (Rold[0] * sinx - Rold[1] * cosx) - PerpenDistance;
1468 if (distOld > 0.0 && distNew <= 0.0) flag = true;
1469 // This parameter is used correct cavity phase
1470 Dold = distOld;
1471 return flag;
1472}
1473
1474bool ParallelCyclotronTracker::RFkick(RFCavity* rfcavity, const double t, const double dt, const int Pindex) {
1475
1476 double radius = std::sqrt(std::pow(itsBunch_m->R[Pindex](0), 2.0) + std::pow(itsBunch_m->R[Pindex](1), 2.0)
1477 - std::pow(rfcavity->getPerpenDistance() , 2.0));
1478 double rmin = rfcavity->getRmin();
1479 double rmax = rfcavity->getRmax();
1480 double nomalRadius = (radius - rmin) / (rmax - rmin);
1481 double tempP[3];
1482 if (nomalRadius <= 1.0 && nomalRadius >= 0.0) {
1483
1484 for (int j = 0; j < 3; j++)
1485 tempP[j] = itsBunch_m->P[Pindex](j); //[px,py,pz] units: dimensionless
1486
1487 // here evaluate voltage and conduct momenta kicking;
1488 rfcavity->getMomentaKick(nomalRadius, tempP, t, dt, itsBunch_m->ID[Pindex], itsBunch_m->getM(), itsBunch_m->getQ()); // t : ns
1489
1490 for (int j = 0; j < 3; j++)
1491 itsBunch_m->P[Pindex](j) = tempP[j];
1492 return true;
1493 }
1494 return false;
1495}
1496
1497
1498struct adder {
1499 adder() : sum(0) {}
1500 double sum;
1501 void operator()(double x) { sum += x; }
1502};
1503
1516bool ParallelCyclotronTracker::getTunes(dvector_t& t, dvector_t& r, dvector_t& z,
1517 int lastTurn, double& /*nur*/, double& /*nuz*/) {
1518 TUNE_class *tune;
1519
1520 int Ndat = t.size();
1521
1522 /*
1523 remove mean
1524 */
1525 double rsum = for_each(r.begin(), r.end(), adder()).sum;
1526
1527 for (int i = 0; i < Ndat; i++)
1528 r[i] -= rsum;
1529
1530 double zsum = for_each(z.begin(), z.end(), adder()).sum;
1531
1532 for (int i = 0; i < Ndat; i++)
1533 z[i] -= zsum;
1534 double ti = *(t.begin());
1535 double tf = t[t.size()-1];
1536 double T = (tf - ti);
1537
1538 t.clear();
1539 for (int i = 0; i < Ndat; i++) {
1540 t.push_back(i);
1541 }
1542
1543 T = t[Ndat-1];
1544
1545 *gmsg << endl;
1546 *gmsg << "* ************************************* nuR ******************************************* *" << endl;
1547 *gmsg << endl << "* ===> " << Ndat << " data points Ti=" << ti << " Tf= " << tf << " -> T= " << T << endl;
1548
1549 int nhis_lomb = 10;
1550 int stat = 0;
1551 // book tune class
1552 tune = new TUNE_class();
1553 stat = tune->lombAnalysis(t, r, nhis_lomb, T / lastTurn);
1554 if (stat != 0)
1555 *gmsg << "* TUNE: Lomb analysis failed" << endl;
1556 *gmsg << "* ************************************************************************************* *" << endl;
1557
1558 delete tune;
1559 tune = nullptr;
1560 // FIXME: FixMe: need to come from the inputfile
1561 nhis_lomb = 100;
1562
1563 if (zsum != 0.0) {
1564
1565 *gmsg << endl;
1566 *gmsg << "* ************************************* nuZ ******************************************* *" << endl;
1567 *gmsg << endl << "* ===> " << Ndat << " data points Ti=" << ti << " Tf= " << tf << " -> T= " << T << endl;
1568
1569 // book tune class
1570 tune = new TUNE_class();
1571 stat = tune->lombAnalysis(t, z, nhis_lomb, T / lastTurn);
1572 if (stat != 0)
1573 *gmsg << "* TUNE: Lomb analysis failed" << endl;
1574 *gmsg << "* ************************************************************************************* *" << endl;
1575
1576 delete tune;
1577 tune = nullptr;
1578 }
1579 return true;
1580}
1581
1582
1583double ParallelCyclotronTracker::getHarmonicNumber() const {
1584 if (opalRing_m != nullptr)
1585 return opalRing_m->getHarmonicNumber();
1586 Cyclotron* elcycl = dynamic_cast<Cyclotron*>(((*FieldDimensions.begin())->second).second);
1587 if (elcycl != nullptr)
1588 return elcycl->getCyclHarm();
1589 throw OpalException("ParallelCyclotronTracker::getHarmonicNumber()",
1590 std::string("The first item in the FieldDimensions list does not ")
1591 +std::string("seem to be a Ring or a Cyclotron element"));
1592}
1593
1594
1595Vector_t ParallelCyclotronTracker::calcMeanR(short bunchNr) const {
1596 Vector_t meanR(0.0, 0.0, 0.0);
1597
1598 for (unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) {
1599 // take all particles if bunchNr <= -1
1600 if ( bunchNr > -1 && itsBunch_m->bunchNum[i] != bunchNr)
1601 continue;
1602
1603 for (int d = 0; d < 3; ++d) {
1604 meanR(d) += itsBunch_m->R[i](d);
1605 }
1606 }
1607
1608 reduce(meanR, meanR, OpAddAssign());
1609
1610 size_t n = itsBunch_m->getTotalNum();
1611
1612 if ( bunchNr > -1 )
1613 n = itsBunch_m->getTotalNumPerBunch(bunchNr);
1614
1615 return meanR / Vector_t(n);
1616}
1617
1618Vector_t ParallelCyclotronTracker::calcMeanP() const {
1619 Vector_t meanP(0.0, 0.0, 0.0);
1620
1621 for (unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) {
1622 for (int d = 0; d < 3; ++d) {
1623 meanP(d) += itsBunch_m->P[i](d);
1624 }
1625 }
1626
1627 reduce(meanP, meanP, OpAddAssign());
1628 return meanP / Vector_t(itsBunch_m->getTotalNum());
1629}
1630
1639
1640void ParallelCyclotronTracker::globalToLocal(ParticleAttrib<Vector_t>& particleVectors,
1641 double phi, Vector_t const translationToGlobal) {
1642 IpplTimings::startTimer(TransformTimer_m);
1643 particleVectors -= translationToGlobal;
1644
1645 rotation_m(0,0) = std::cos(phi);
1646 rotation_m(0,1) = std::sin(phi);
1647 rotation_m(0,2) = 0;
1648 rotation_m(1,0) = -std::sin(phi);
1649 rotation_m(1,1) = std::cos(phi);
1650 rotation_m(1,2) = 0;
1651 rotation_m(2,0) = 0;
1652 rotation_m(2,1) = 0;
1653 rotation_m(2,2) = 1;
1654
1655 for (unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) {
1656 particleVectors[i] = prod_boost_vector(rotation_m, particleVectors[i]);
1657 }
1658 IpplTimings::stopTimer(TransformTimer_m);
1659}
1660
1661void ParallelCyclotronTracker::localToGlobal(ParticleAttrib<Vector_t>& particleVectors,
1662 double phi, Vector_t const translationToGlobal) {
1663 IpplTimings::startTimer(TransformTimer_m);
1664
1665 rotation_m(0,0) = std::cos(phi);
1666 rotation_m(0,1) = -std::sin(phi);
1667 rotation_m(0,2) = 0;
1668 rotation_m(1,0) = std::sin(phi);
1669 rotation_m(1,1) = std::cos(phi);
1670 rotation_m(1,2) = 0;
1671 rotation_m(2,0) = 0;
1672 rotation_m(2,1) = 0;
1673 rotation_m(2,2) = 1;
1674
1675 for (unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) {
1676 particleVectors[i] = prod_boost_vector(rotation_m, particleVectors[i]);
1677 }
1678
1679 particleVectors += translationToGlobal;
1680 IpplTimings::stopTimer(TransformTimer_m);
1681}
1682
1683
1684inline void ParallelCyclotronTracker::globalToLocal(ParticleAttrib<Vector_t>& particleVectors,
1685 Quaternion_t const quaternion,
1686 Vector_t const meanR) {
1687 IpplTimings::startTimer(TransformTimer_m);
1688
1689 // Translation from global to local
1690 particleVectors -= meanR;
1691
1692 // Rotation to align P_mean with x-axis
1693 rotateWithQuaternion(particleVectors, quaternion);
1694 IpplTimings::stopTimer(TransformTimer_m);
1695}
1696
1697inline void ParallelCyclotronTracker::localToGlobal(ParticleAttrib<Vector_t>& particleVectors,
1698 Quaternion_t const quaternion,
1699 Vector_t const meanR) {
1700 IpplTimings::startTimer(TransformTimer_m);
1701 // Reverse the quaternion by multiplying the axis components (x,y,z) with -1
1702 Quaternion_t reverseQuaternion = quaternion * -1.0;
1703 reverseQuaternion(0) *= -1.0;
1704
1705 // Rotation back to original P_mean direction
1706 rotateWithQuaternion(particleVectors, reverseQuaternion);
1707
1708 // Translation from local to global
1709 particleVectors += meanR;
1710 IpplTimings::stopTimer(TransformTimer_m);
1711}
1712
1713inline void ParallelCyclotronTracker::globalToLocal(ParticleAttrib<Vector_t>& particleVectors,
1714 double const phi,
1715 double const psi,
1716 Vector_t const meanR) {
1717
1718 IpplTimings::startTimer(TransformTimer_m);
1719 //double const tolerance = 1.0e-4; // TODO What is a good angle threshold? -DW
1720
1721 // Translation from global to local
1722 particleVectors -= meanR;
1723
1724 // Rotation to align P_mean with x-axis
1725 rotateAroundZ(particleVectors, phi);
1726
1727 // If theta is large enough (i.e. there is significant momentum in z direction)
1728 // rotate around x-axis next
1729 //if (fabs(psi) > tolerance)
1730 rotateAroundX(particleVectors, psi);
1731 IpplTimings::stopTimer(TransformTimer_m);
1732}
1733
1734inline void ParallelCyclotronTracker::globalToLocal(Vector_t& myVector,
1735 double const phi,
1736 double const psi,
1737 Vector_t const meanR) {
1738
1739 IpplTimings::startTimer(TransformTimer_m);
1740 //double const tolerance = 1.0e-4; // TODO What is a good angle threshold? -DW
1741
1742 // Translation from global to local
1743 myVector -= meanR;
1744
1745 // Rotation to align P_mean with x-axis
1746 rotateAroundZ(myVector, phi);
1747
1748 // If theta is large enough (i.e. there is significant momentum in z direction)
1749 // rotate around x-axis next
1750 //if (fabs(psi) > tolerance)
1751 rotateAroundX(myVector, psi);
1752 IpplTimings::stopTimer(TransformTimer_m);
1753}
1754
1755inline void ParallelCyclotronTracker::localToGlobal(ParticleAttrib<Vector_t>& particleVectors,
1756 double const phi,
1757 double const psi,
1758 Vector_t const meanR) {
1759
1760 IpplTimings::startTimer(TransformTimer_m);
1761 //double const tolerance = 1.0e-4; // TODO What is a good angle threshold? -DW
1762
1763 // If theta is large enough (i.e. there is significant momentum in z direction)
1764 // rotate back around x-axis next
1765 //if (fabs(psi) > tolerance)
1766 rotateAroundX(particleVectors, -psi);
1767
1768 // Rotation to align P_mean with x-axis
1769 rotateAroundZ(particleVectors, -phi);
1770
1771 // Translation from local to global
1772 particleVectors += meanR;
1773 IpplTimings::stopTimer(TransformTimer_m);
1774}
1775
1776inline void ParallelCyclotronTracker::localToGlobal(Vector_t& myVector,
1777 double const phi,
1778 double const psi,
1779 Vector_t const meanR) {
1780
1781 IpplTimings::startTimer(TransformTimer_m);
1782 //double const tolerance = 1.0e-4; // TODO What is a good angle threshold? -DW
1783
1784 // If theta is large enough (i.e. there is significant momentum in z direction)
1785 // rotate back around x-axis next
1786 //if (fabs(psi) > tolerance)
1787 rotateAroundX(myVector, -psi);
1788
1789 // Rotation to align P_mean with x-axis
1790 rotateAroundZ(myVector, -phi);
1791
1792 // Translation from local to global
1793 myVector += meanR;
1794 IpplTimings::stopTimer(TransformTimer_m);
1795}
1796
1797inline void ParallelCyclotronTracker::rotateWithQuaternion(ParticleAttrib<Vector_t>& particleVectors, Quaternion_t const quaternion) {
1798
1799 Vector_t const quaternionVectorComponent = Vector_t(quaternion(1), quaternion(2), quaternion(3));
1800 double const quaternionScalarComponent = quaternion(0);
1801
1802 for (unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) {
1803
1804 particleVectors[i] = 2.0f * dot(quaternionVectorComponent, particleVectors[i]) * quaternionVectorComponent +
1805 (quaternionScalarComponent * quaternionScalarComponent -
1806 dot(quaternionVectorComponent, quaternionVectorComponent)) * particleVectors[i] + 2.0f *
1807 quaternionScalarComponent * cross(quaternionVectorComponent, particleVectors[i]);
1808 }
1809}
1810
1811inline void ParallelCyclotronTracker::normalizeQuaternion(Quaternion_t& quaternion){
1812
1813 double tolerance = 1.0e-10;
1814 double length2 = dot(quaternion, quaternion);
1815
1816 if (std::abs(length2) > tolerance && std::abs(length2 - 1.0f) > tolerance) {
1817
1818 double length = std::sqrt(length2);
1819 quaternion /= length;
1820 }
1821}
1822
1823inline void ParallelCyclotronTracker::normalizeVector(Vector_t& vector) {
1824
1825 double tolerance = 1.0e-10;
1826 double length2 = dot(vector, vector);
1827
1828 if (std::abs(length2) > tolerance && std::abs(length2 - 1.0f) > tolerance) {
1829
1830 double length = std::sqrt(length2);
1831 vector /= length;
1832 }
1833}
1834
1835inline void ParallelCyclotronTracker::rotateAroundZ(ParticleAttrib<Vector_t>& particleVectors, double const phi) {
1836 // Clockwise rotation of particles 'particleVectors' by 'phi' around Z axis
1837
1838 rotation_m(0,0) = std::cos(phi);
1839 rotation_m(0,1) = std::sin(phi);
1840 rotation_m(0,2) = 0;
1841 rotation_m(1,0) = -std::sin(phi);
1842 rotation_m(1,1) = std::cos(phi);
1843 rotation_m(1,2) = 0;
1844 rotation_m(2,0) = 0;
1845 rotation_m(2,1) = 0;
1846 rotation_m(2,2) = 1;
1847
1848 for (unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) {
1849 particleVectors[i] = prod_boost_vector(rotation_m, particleVectors[i]);
1850 }
1851}
1852
1853inline void ParallelCyclotronTracker::rotateAroundZ(Vector_t& myVector, double const phi) {
1854 // Clockwise rotation of single vector 'myVector' by 'phi' around Z axis
1855
1856 rotation_m(0,0) = std::cos(phi);
1857 rotation_m(0,1) = std::sin(phi);
1858 rotation_m(0,2) = 0;
1859 rotation_m(1,0) = -std::sin(phi);
1860 rotation_m(1,1) = std::cos(phi);
1861 rotation_m(1,2) = 0;
1862 rotation_m(2,0) = 0;
1863 rotation_m(2,1) = 0;
1864 rotation_m(2,2) = 1;
1865
1866 myVector = prod_boost_vector(rotation_m, myVector);
1867}
1868
1869inline void ParallelCyclotronTracker::rotateAroundX(ParticleAttrib<Vector_t>& particleVectors, double const psi) {
1870 // Clockwise rotation of particles 'particleVectors' by 'psi' around X axis
1871
1872 rotation_m(0,0) = 1;
1873 rotation_m(0,1) = 0;
1874 rotation_m(0,2) = 0;
1875 rotation_m(1,0) = 0;
1876 rotation_m(1,1) = std::cos(psi);
1877 rotation_m(1,2) = std::sin(psi);
1878 rotation_m(2,0) = 0;
1879 rotation_m(2,1) = -std::sin(psi);
1880 rotation_m(2,2) = std::cos(psi);
1881
1882 for (unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) {
1883 particleVectors[i] = prod_boost_vector(rotation_m, particleVectors[i]);
1884 }
1885}
1886
1887inline void ParallelCyclotronTracker::rotateAroundX(Vector_t& myVector, double const psi) {
1888 // Clockwise rotation of single vector 'myVector' by 'psi' around X axis
1889
1890 rotation_m(0,0) = 1;
1891 rotation_m(0,1) = 0;
1892 rotation_m(0,2) = 0;
1893 rotation_m(1,0) = 0;
1894 rotation_m(1,1) = std::cos(psi);
1895 rotation_m(1,2) = std::sin(psi);
1896 rotation_m(2,0) = 0;
1897 rotation_m(2,1) = -std::sin(psi);
1898 rotation_m(2,2) = std::cos(psi);
1899
1900 myVector = prod_boost_vector(rotation_m, myVector);
1901}
1902
1903inline void ParallelCyclotronTracker::getQuaternionTwoVectors(Vector_t u, Vector_t v, Quaternion_t& quaternion) {
1904 // four vector (w,x,y,z) of the quaternion of P_mean with the positive x-axis
1905
1906 normalizeVector(u);
1907 normalizeVector(v);
1908
1909 double k_cos_theta = dot(u, v);
1910 double k = std::sqrt(dot(u, u) * dot(v, v));
1911 double tolerance1 = 1.0e-5;
1912 double tolerance2 = 1.0e-8;
1913 Vector_t resultVectorComponent;
1914
1915 if (std::abs(k_cos_theta / k + 1.0) < tolerance1) {
1916 // u and v are almost exactly antiparallel so we need to do
1917 // 180 degree rotation around any vector orthogonal to u
1918
1919 resultVectorComponent = cross(u, xaxis);
1920
1921 // If by chance u is parallel to xaxis, use zaxis instead
1922 if (dot(resultVectorComponent, resultVectorComponent) < tolerance2) {
1923 resultVectorComponent = cross(u, zaxis);
1924 }
1925
1926 double halfAngle = 0.5 * Physics::pi;
1927 double sinHalfAngle = std::sin(halfAngle);
1928
1929 resultVectorComponent *= sinHalfAngle;
1930
1931 k = 0.0;
1932 k_cos_theta = std::cos(halfAngle);
1933
1934 } else {
1935 resultVectorComponent = cross(u, v);
1936 }
1937
1938 quaternion(0) = k_cos_theta + k;
1939 quaternion(1) = resultVectorComponent(0);
1940 quaternion(2) = resultVectorComponent(1);
1941 quaternion(3) = resultVectorComponent(2);
1942
1943 normalizeQuaternion(quaternion);
1944}
1945
1946
1947bool ParallelCyclotronTracker::push(double h) {
1948
1949 IpplTimings::startTimer(IntegrationTimer_m);
1950
1951 bool flagNeedUpdate = false;
1952 for (unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) {
1953
1954 Vector_t const oldR = itsBunch_m->R[i];
1955 double const gamma = Util::getGamma(itsBunch_m->P[i]);
1956 double const c_gamma = Physics::c / gamma;
1957 Vector_t const v = itsBunch_m->P[i] * c_gamma;
1958
1959 itsBunch_m->R[i] += h * v;
1960
1961 for (const auto & ccd : cavCrossDatas_m) {
1962 double const distNew = (itsBunch_m->R[i][0] * ccd.sinAzimuth - itsBunch_m->R[i][1] * ccd.cosAzimuth) - ccd.perpenDistance;
1963 bool tagCrossing = false;
1964 double distOld;
1965 if (distNew <= 0.0) {
1966 distOld = (oldR[0] * ccd.sinAzimuth - oldR[1] * ccd.cosAzimuth) - ccd.perpenDistance;
1967 if (distOld > 0.0) tagCrossing = true;
1968 }
1969 if (tagCrossing) {
1970 double const dt1 = distOld / euclidean_norm(v);
1971 double const dt2 = h - dt1;
1972
1973 // Re-track particle from the old position to cavity gap point
1974 itsBunch_m->R[i] = oldR + dt1 * v;
1975
1976 // Momentum kick
1977 RFkick(ccd.cavity, itsBunch_m->getT(), dt1, i);
1978
1979 itsBunch_m->R[i] += dt2 * itsBunch_m->P[i] * c_gamma;
1980 }
1981 }
1982 flagNeedUpdate |= (itsBunch_m->Bin[i] < 0);
1983 }
1984
1985 IpplTimings::stopTimer(IntegrationTimer_m);
1986 return flagNeedUpdate;
1987}
1988
1989
1990bool ParallelCyclotronTracker::kick(double h) {
1991 IpplTimings::startTimer(IntegrationTimer_m);
1992
1993 bool flagNeedUpdate = false;
1994 BorisPusher pusher;
1995 double const q = itsBunch_m->Q[0] / Physics::q_e; // For now all particles have the same charge
1996 double const M = itsBunch_m->M[0] * Units::GeV2eV; // For now all particles have the same rest energy
1997
1998 for (unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) {
1999
2000 pusher.kick(itsBunch_m->R[i], itsBunch_m->P[i],
2001 itsBunch_m->Ef[i], itsBunch_m->Bf[i],
2002 h, M, q);
2003 flagNeedUpdate |= (itsBunch_m->Bin[i] < 0);
2004 }
2005 IpplTimings::stopTimer(IntegrationTimer_m);
2006 return flagNeedUpdate;
2007}
2008
2009
2010void ParallelCyclotronTracker::borisExternalFields(double h) {
2011
2012 // push particles for first half step
2013 bool flagNeedUpdate = push(0.5 * h);
2014
2015 // Evaluate external fields
2016 IpplTimings::startTimer(IntegrationTimer_m);
2017 for (unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) {
2018
2019 itsBunch_m->Ef[i] = Vector_t(0.0, 0.0, 0.0);
2020 itsBunch_m->Bf[i] = Vector_t(0.0, 0.0, 0.0);
2021
2022 computeExternalFields_m(i, itsBunch_m->getT(),
2023 itsBunch_m->Ef[i], itsBunch_m->Bf[i]);
2024 }
2025 IpplTimings::stopTimer(IntegrationTimer_m);
2026
2027 // Kick particles for full step
2028 flagNeedUpdate |= kick(h);
2029
2030 // push particles for second half step
2031 flagNeedUpdate |= push(0.5 * h);
2032
2033 // apply the plugin elements: probe, collimator, stripper, septum
2034 flagNeedUpdate |= applyPluginElements(h);
2035 // destroy particles if they are marked as Bin=-1 in the plugin elements or out of global aperture
2036 deleteParticle(flagNeedUpdate);
2037}
2038
2039
2040bool ParallelCyclotronTracker::applyPluginElements(const double dt) {
2041 IpplTimings::startTimer(PluginElemTimer_m);
2042
2043 for (beamline_list::iterator sindex = ++(FieldDimensions.begin()); sindex != FieldDimensions.end(); ++sindex) {
2044 if (((*sindex)->first) == ElementType::VACUUM) {
2045 Vacuum* vac = static_cast<Vacuum*>(((*sindex)->second).second);
2046 vac->checkVacuum(itsBunch_m, cycl_m);
2047 }
2048 }
2049
2050 bool flag = false;
2051 for (PluginElement* element : pluginElements_m) {
2052 bool tmp = element->check(itsBunch_m,
2053 turnnumber_m,
2054 itsBunch_m->getT(),
2055 dt);
2056 flag |= tmp;
2057
2058 if ( tmp ) {
2059 itsBunch_m->updateNumTotal();
2060 *gmsg << "* Total number of particles after PluginElement= "
2061 << itsBunch_m->getTotalNum() << endl;
2062 }
2063 }
2064
2065 IpplTimings::stopTimer(PluginElemTimer_m);
2066 return flag;
2067}
2068
2069bool ParallelCyclotronTracker::deleteParticle(bool flagNeedUpdate) {
2070 IpplTimings::startTimer(DelParticleTimer_m);
2071 // Update immediately if any particles are lost during this step
2072
2073 if ((step_m + 1) % Options::delPartFreq != 0) {
2074 IpplTimings::stopTimer(DelParticleTimer_m);
2075 return false;
2076 }
2077
2078 allreduce(flagNeedUpdate, 1, std::logical_or<bool>());
2079
2080 if (flagNeedUpdate) {
2081 short bunchCount = itsBunch_m->getNumBunch();
2082 std::vector<size_t> locLostParticleNum(bunchCount, 0);
2083
2084 const int leb = itsBunch_m->getLastemittedBin();
2085 std::unique_ptr<size_t[]> localBinCount;
2086
2087 if ( isMultiBunch() ) {
2088 localBinCount = std::unique_ptr<size_t[]>(new size_t[leb]);
2089 for (int i = 0; i < leb; ++i)
2090 localBinCount[i] = 0;
2091 }
2092
2093 for (unsigned int i = 0; i < itsBunch_m->getLocalNum(); i++) {
2094 if (itsBunch_m->Bin[i] < 0) {
2095 ++locLostParticleNum[itsBunch_m->bunchNum[i]];
2096 itsBunch_m->destroy(1, i);
2097 } else if ( isMultiBunch() ) {
2098 /* we need to count the local number of particles
2099 * per energy bin
2100 */
2101 ++localBinCount[itsBunch_m->Bin[i]];
2102 }
2103 }
2104
2105 if ( isMultiBunch() ) {
2106 // set the local bin count
2107 for (int i = 0; i < leb; ++i) {
2108 itsBunch_m->setLocalBinCount(localBinCount[i], i);
2109 }
2110 }
2111
2112 std::vector<size_t> localnum(bunchCount + 1);
2113 for (size_t i = 0; i < localnum.size() - 1; ++i) {
2114 localnum[i] = itsBunch_m->getLocalNumPerBunch(i) - locLostParticleNum[i];
2115 itsBunch_m->setLocalNumPerBunch(localnum[i], i);
2116 }
2117
2118 /* We need to destroy the particles now
2119 * before we compute the means. We also
2120 * have to update the total number of particles
2121 * otherwise the statistics are wrong.
2122 */
2123 itsBunch_m->performDestroy(true);
2124
2125 /* total number of particles of individual bunches
2126 * last index of vector contains total number over all
2127 * bunches, used as a check
2128 */
2129 std::vector<size_t> totalnum(bunchCount + 1);
2130 localnum[bunchCount] = itsBunch_m->getLocalNum();
2131
2132 allreduce(localnum.data(), totalnum.data(), localnum.size(), std::plus<size_t>());
2133 itsBunch_m->setTotalNum(totalnum[bunchCount]);
2134
2135 for (short i = 0; i < bunchCount; ++i) {
2136 itsBunch_m->setTotalNumPerBunch(totalnum[i], i);
2137 }
2138
2139 size_t sum = std::accumulate(totalnum.begin(),
2140 totalnum.end() - 1, 0);
2141
2142 if ( sum != totalnum[bunchCount] ) {
2143 throw OpalException("ParallelCyclotronTracker::deleteParticle()",
2144 "Total number of particles " + std::to_string(totalnum[bunchCount]) +
2145 " != " + std::to_string(sum) + " (sum over all bunches)");
2146 }
2147
2148 size_t globLostParticleNum = 0;
2149 size_t locNumLost = std::accumulate(locLostParticleNum.begin(),
2150 locLostParticleNum.end(), 0);
2151
2152 reduce(locNumLost, globLostParticleNum, 1, std::plus<size_t>());
2153
2154 if ( globLostParticleNum > 0 ) {
2155 *gmsg << level3 << "At step " << step_m
2156 << ", lost " << globLostParticleNum << " particles" << endl;
2157 }
2158
2159 if (totalnum[bunchCount] == 0) {
2160 IpplTimings::stopTimer(DelParticleTimer_m);
2161 return flagNeedUpdate;
2162 }
2163
2164 Vector_t const meanR = calcMeanR();
2165 Vector_t const meanP = calcMeanP();
2166
2167 // Bunch (local) azimuth at meanR w.r.t. y-axis
2168 double const phi = calculateAngle(meanP(0), meanP(1)) - 0.5 * Physics::pi;
2169
2170 // Bunch (local) elevation at meanR w.r.t. xy plane
2171 double const psi = 0.5 * Physics::pi - std::acos(meanP(2) / euclidean_norm(meanP));
2172
2173 // For statistics data, the bunch is transformed into a local coordinate system
2174 // with meanP in direction of y-axis -DW
2175 globalToLocal(itsBunch_m->R, phi, psi, meanR);
2176 globalToLocal(itsBunch_m->P, phi, psi, Vector_t(0.0)); // P should be rotated, but not shifted
2177
2178 // Now destroy particles and update pertinent parameters in local frame
2179 // Note that update() will be called within boundp() -DW
2180 itsBunch_m->boundp();
2181 //itsBunch_m->update();
2182
2183 itsBunch_m->calcBeamParameters();
2184
2185 localToGlobal(itsBunch_m->R, phi, psi, meanR);
2186 localToGlobal(itsBunch_m->P, phi, psi, Vector_t(0.0));
2187
2188 // If particles were deleted, recalculate bingamma and reset BinID for remaining particles
2189 if ( isMultiBunch() )
2190 mbHandler_m->updateParticleBins(itsBunch_m);
2191 }
2192
2193 IpplTimings::stopTimer(DelParticleTimer_m);
2194 return flagNeedUpdate;
2195}
2196
2197void ParallelCyclotronTracker::initTrackOrbitFile() {
2198
2199 std::string f = OpalData::getInstance()->getInputBasename() + std::string("-trackOrbit.dat");
2200
2201 outfTrackOrbit_m.setf(std::ios::scientific, std::ios::floatfield);
2202 outfTrackOrbit_m.precision(8);
2203
2204 if (myNode_m == 0) {
2205 if (OpalData::getInstance()->inRestartRun()) {
2206 outfTrackOrbit_m.open(f.c_str(), std::ios::app);
2207 outfTrackOrbit_m << "# Restart at integration step " << itsBunch_m->getLocalTrackStep() << std::endl;
2208 } else {
2209 outfTrackOrbit_m.open(f.c_str());
2210 outfTrackOrbit_m << "# The six-dimensional phase space data in the global Cartesian coordinates" << std::endl;
2211 outfTrackOrbit_m << "# Part. ID x [m] beta_x*gamma y [m] beta_y*gamma z [m] beta_z*gamma" << std::endl;
2212 }
2213 }
2214}
2215
2216void ParallelCyclotronTracker::initDistInGlobalFrame() {
2217
2218 if (!OpalData::getInstance()->inRestartRun()) {
2219 // Start a new run (no restart)
2220
2221 double const initialReferenceTheta = referenceTheta;
2222
2223 // TODO: Replace with TracerParticle
2224 // Force the initial phase space values of the particle with ID = 0 to zero,
2225 // to set it as a reference particle.
2226 if (initialTotalNum_m > 2) {
2227 for (size_t i = 0; i < initialLocalNum_m; ++i) {
2228 if (itsBunch_m->ID[i] == 0) {
2229 itsBunch_m->R[i] = Vector_t(0.0);
2230 itsBunch_m->P[i] = Vector_t(0.0);
2231 }
2232 }
2233 }
2234
2235 // Initialize global R
2236 Vector_t const initMeanR = Vector_t(referenceR * cosRefTheta_m,
2237 referenceR * sinRefTheta_m,
2238 referenceZ);
2239
2240 localToGlobal(itsBunch_m->R, initialReferenceTheta, initMeanR);
2241
2242 // Initialize global P (Cartesian, but input P_ref is in Pr, Ptheta, Pz,
2243 // so translation has to be done before the rotation this once)
2244 // Cave: In the local frame, the positive y-axis is the direction of movement -DW
2245 for (size_t i = 0; i < initialLocalNum_m; ++i) {
2246 itsBunch_m->P[i](0) += referencePr;
2247 itsBunch_m->P[i](1) += referencePt;
2248 itsBunch_m->P[i](2) += referencePz;
2249 }
2250
2251 // Out of the three coordinates of meanR (R, Theta, Z) only the angle
2252 // changes the momentum vector...
2253 double angle = initialReferenceTheta + beamInitialRotation;
2254 localToGlobal(itsBunch_m->P, angle);
2255
2256 DistributionType distType = itsBunch_m->getDistType();
2257 if (distType == DistributionType::FROMFILE) {
2258 checkFileMomentum();
2259 }
2260
2261 // Initialize the bin number of the first bunch to 0
2262 for (size_t i = 0; i < initialLocalNum_m; ++i) {
2263 itsBunch_m->Bin[i] = 0;
2264 }
2265
2266 // Set time step per particle
2267 setTimeStep();
2268
2269 // Backup initial distribution if multi bunch mode
2270 if ((initialTotalNum_m > 2) && isMultiBunch() && mbHandler_m->isForceMode()) {
2271 mbHandler_m->saveBunch(itsBunch_m);
2272 }
2273
2274 // Else: Restart from the distribution in the h5 file
2275 } else {
2276
2277 // Do a local frame restart (we have already checked that the old h5 file was saved in local
2278 // frame as well).
2279 if ((Options::psDumpFrame != DumpFrame::GLOBAL)) {
2280
2281 *gmsg << "* Restart in the local frame" << endl;
2282
2283 Vector_t const initMeanR = Vector_t(referenceR * cosRefTheta_m,
2284 referenceR * sinRefTheta_m,
2285 referenceZ);
2286
2287 // Do the transformations
2288 localToGlobal(itsBunch_m->R, referencePhi, referencePsi, initMeanR);
2289 localToGlobal(itsBunch_m->P, referencePhi, referencePsi);
2290
2291 // Initialize the bin number of the first bunch to 0
2292 for (size_t i = 0; i < initialLocalNum_m; ++i) {
2293 itsBunch_m->Bin[i] = 0;
2294 }
2295
2296 // Or do a global frame restart (no transformations necessary)
2297 } else {
2298 *gmsg << "* Restart in the global frame" << endl;
2299
2300 pathLength_m = itsBunch_m->get_sPos();
2301 }
2302 }
2303
2304 // set the number of particles per bunch
2305 itsBunch_m->countTotalNumPerBunch();
2306
2307 // ------------ Get some Values ---------------------------------------------------------- //
2308 Vector_t const meanR = calcMeanR();
2309 Vector_t const meanP = calcMeanP();
2310
2311 // Bunch (local) azimuth at meanR w.r.t. y-axis
2312 double const phi = calculateAngle(meanP(0), meanP(1)) - 0.5 * Physics::pi;
2313
2314 // Bunch (local) elevation at meanR w.r.t. xy plane
2315 double const psi = 0.5 * Physics::pi - std::acos(meanP(2) / euclidean_norm(meanP));
2316
2317 double radius = std::sqrt(meanR[0] * meanR[0] + meanR[1] * meanR[1]);
2318
2319 if ( isMultiBunch() ) {
2320 mbHandler_m->setRadiusTurns(radius);
2321 }
2322
2323 // Do boundp and repartition in the local frame at beginning of this run
2324 globalToLocal(itsBunch_m->R, phi, psi, meanR);
2325 globalToLocal(itsBunch_m->P, phi, psi); // P should be rotated, but not shifted
2326
2327 itsBunch_m->boundp();
2328
2329 checkNumPart(std::string("\n* Before repartition: "));
2330 repartition();
2331 checkNumPart(std::string("\n* After repartition: "));
2332
2333 itsBunch_m->calcBeamParameters();
2334
2335 *gmsg << endl << "* *********************** Bunch information in local frame: ************************";
2336 *gmsg << *itsBunch_m << endl;
2337
2338 localToGlobal(itsBunch_m->R, phi, psi, meanR);
2339 localToGlobal(itsBunch_m->P, phi, psi);
2340
2341 // Save initial distribution if not a restart
2342 if (!OpalData::getInstance()->inRestartRun()) {
2343 step_m -= 1;
2344
2345 bunchDumpPhaseSpaceData();
2346 bunchDumpStatData();
2347
2348 step_m += 1;
2349 }
2350
2351 // Print out the Bunch information at beginning of the run.
2352 itsBunch_m->calcBeamParameters();
2353
2354 // multi-bunch simulation only
2355 saveInjectValues();
2356
2357 *gmsg << endl << "* *********************** Bunch information in global frame: ***********************";
2358 *gmsg << *itsBunch_m << endl;
2359}
2360
2361void ParallelCyclotronTracker::setTimeStep() {
2362 for (size_t i = 0; i < initialLocalNum_m; ++i) {
2363 itsBunch_m->dt[i] = itsBunch_m->getdT();
2364 }
2365}
2366
2367void ParallelCyclotronTracker::checkFileMomentum() {
2368
2369 double pTotalMean = 0.0;
2370 for (size_t i = 0; i < initialLocalNum_m; ++i) {
2371 pTotalMean += euclidean_norm(itsBunch_m->P[i]);
2372 }
2373
2374 allreduce(pTotalMean, 1, std::plus<double>());
2375
2376 pTotalMean /= initialTotalNum_m;
2377
2378 double tolerance = itsBunch_m->getMomentumTolerance();
2379 if (tolerance > 0 &&
2380 std::abs(pTotalMean - referencePtot) / pTotalMean > tolerance) {
2381 throw OpalException("ParallelCyclotronTracker::checkFileMomentum",
2382 "The total momentum of the particle distribution\n"
2383 "in the global reference frame: " +
2384 std::to_string(pTotalMean) + ",\n"
2385 "is different from the momentum given\n"
2386 "in the \"BEAM\" command: " +
2387 std::to_string(referencePtot) + ".\n"
2388 "In Opal-cycl the initial distribution\n"
2389 "is specified in the local reference frame.\n"
2390 "When using a \"FROMFILE\" type distribution, the momentum \n"
2391 "must be the same as the specified in the \"BEAM\" command,\n"
2392 "which is in global reference frame.");
2393 }
2394}
2395
2396
2397//TODO: This can be completely rewritten with TracerParticle -DW
2398void ParallelCyclotronTracker::singleParticleDump() {
2399 IpplTimings::startTimer(DumpTimer_m);
2400
2401 if (Ippl::getNodes() > 1 ) {
2402
2403 double x;
2404 int id;
2405 dvector_t tmpr;
2406 ivector_t tmpi;
2407
2408 int tag = Ippl::Comm->next_tag(IPPL_APP_TAG4, IPPL_APP_CYCLE);
2409
2410 // for all nodes, find the location of particle with ID = 0 & 1 in bunch containers
2411 int found[2] = {-1, -1};
2412 int counter = 0;
2413
2414 for (size_t i = 0; i < itsBunch_m->getLocalNum(); ++i) {
2415 if (itsBunch_m->ID[i] == 0) {
2416 found[counter] = i;
2417 counter++;
2418 }
2419 if (itsBunch_m->ID[i] == 1) {
2420 found[counter] = i;
2421 counter++;
2422 }
2423 }
2424
2425 if (myNode_m == 0) {
2426 int notReceived = Ippl::getNodes() - 1;
2427 int numberOfPart = 0;
2428 // receive other nodes
2429 while(notReceived > 0) {
2430
2431 int node = COMM_ANY_NODE;
2432 Message *rmsg = Ippl::Comm->receive_block(node, tag);
2433
2434 if (rmsg == nullptr)
2435 ERRORMSG("Could not receive from client nodes in main." << endl);
2436
2437 --notReceived;
2438
2439 rmsg->get(&numberOfPart);
2440
2441 for (int i = 0; i < numberOfPart; ++i) {
2442 rmsg->get(&id);
2443 tmpi.push_back(id);
2444 for (int ii = 0; ii < 6; ii++) {
2445 rmsg->get(&x);
2446 tmpr.push_back(x);
2447 }
2448 }
2449 delete rmsg;
2450 }
2451 // own node
2452 for (int i = 0; i < counter; ++i) {
2453
2454 tmpi.push_back(itsBunch_m->ID[found[i]]);
2455
2456 for (int j = 0; j < 3; ++j) {
2457 tmpr.push_back(itsBunch_m->R[found[i]](j));
2458 tmpr.push_back(itsBunch_m->P[found[i]](j));
2459 }
2460 }
2461 // store
2462 dvector_t::iterator itParameter = tmpr.begin();
2463
2464 for (auto tmpid : tmpi) {
2465
2466 outfTrackOrbit_m << "ID" << tmpid;
2467
2468 if (tmpid == 0) { // for stat file
2469 itsBunch_m->RefPartR_m[0] = *itParameter;
2470 itsBunch_m->RefPartR_m[1] = *(itParameter + 2);
2471 itsBunch_m->RefPartR_m[2] = *(itParameter + 4);
2472 itsBunch_m->RefPartP_m[0] = *(itParameter + 1);
2473 itsBunch_m->RefPartP_m[1] = *(itParameter + 3);
2474 itsBunch_m->RefPartP_m[2] = *(itParameter + 5);
2475 }
2476 for (int ii = 0; ii < 6; ii++) {
2477 outfTrackOrbit_m << " " << *itParameter;
2478 ++itParameter;
2479 }
2480
2481 outfTrackOrbit_m << std::endl;
2482 }
2483 } else {
2484 // for other nodes
2485 Message *smsg = new Message();
2486 smsg->put(counter);
2487
2488 for (int i = 0; i < counter; i++) {
2489
2490 smsg->put(itsBunch_m->ID[found[i]]);
2491
2492 for (int j = 0; j < 3; j++) {
2493 smsg->put(itsBunch_m->R[found[i]](j));
2494 smsg->put(itsBunch_m->P[found[i]](j));
2495 }
2496 }
2497
2498 if (!Ippl::Comm->send(smsg, 0, tag)) {
2499 ERRORMSG("Ippl::Comm->send(smsg, 0, tag) failed " << endl);
2500 }
2501 }
2502
2503 Ippl::Comm->barrier();
2504
2505 } else {
2506
2507 for (size_t i = 0; i < itsBunch_m->getLocalNum(); i++) {
2508 if (itsBunch_m->ID[i] == 0 || itsBunch_m->ID[i] == 1) {
2509
2510 outfTrackOrbit_m << "ID" << itsBunch_m->ID[i] << " ";
2511 outfTrackOrbit_m << itsBunch_m->R[i](0) << " " << itsBunch_m->P[i](0) << " ";
2512 outfTrackOrbit_m << itsBunch_m->R[i](1) << " " << itsBunch_m->P[i](1) << " ";
2513 outfTrackOrbit_m << itsBunch_m->R[i](2) << " " << itsBunch_m->P[i](2) << std::endl;
2514
2515 if (itsBunch_m->ID[i] == 0) { // for stat file
2516 itsBunch_m->RefPartR_m = itsBunch_m->R[i];
2517 itsBunch_m->RefPartP_m = itsBunch_m->P[i];
2518 }
2519 }
2520 }
2521 }
2522
2523 IpplTimings::stopTimer(DumpTimer_m);
2524}
2525
2526void ParallelCyclotronTracker::bunchDumpStatData(){
2527
2528 IpplTimings::startTimer(DumpTimer_m);
2529
2530 // dump stat file per bunch in case of multi-bunch mode
2531 if (isMultiBunch()) {
2532 double phi = 0.0, psi = 0.0;
2533 Vector_t meanR = calcMeanR();
2534
2535 // Bunch (global) angle w.r.t. x-axis (cylinder coordinates)
2536 double theta = calculateAngle(meanR(0), meanR(1)) * Units::rad2deg;
2537
2538 // fix azimuth_m --> make monotonically increasing
2539 dumpAngle(theta, prevAzimuth_m, azimuth_m);
2540
2541 updateAzimuthAndRadius();
2542
2543 if (Options::psDumpFrame != DumpFrame::GLOBAL) {
2544 Vector_t meanP = calcMeanP();
2545
2546 // Bunch (local) azimuth at meanR w.r.t. y-axis
2547 phi = calculateAngle(meanP(0), meanP(1)) - 0.5 * Physics::pi;
2548
2549 // Bunch (local) elevation at meanR w.r.t. xy plane
2550 psi = 0.5 * Physics::pi - std::acos(meanP(2) / euclidean_norm(meanP));
2551
2552 // Rotate so Pmean is in positive y direction. No shift, so that normalized emittance and
2553 // unnormalized emittance as well as centroids are calculated correctly in
2554 // PartBunch::calcBeamParameters()
2555 globalToLocal(itsBunch_m->R, phi, psi, meanR);
2556 globalToLocal(itsBunch_m->P, phi, psi);
2557 }
2558
2559 itsDataSink->writeMultiBunchStatistics(itsBunch_m, mbHandler_m.get());
2560
2561 if (Options::psDumpFrame != DumpFrame::GLOBAL) {
2562 localToGlobal(itsBunch_m->R, phi, psi, meanR);
2563 localToGlobal(itsBunch_m->P, phi, psi);
2564 }
2565
2566 IpplTimings::stopTimer(DumpTimer_m);
2567 return;
2568 }
2569
2570 // --------------------------------- Get some Values ---------------------------------------- //
2571 double const temp_t = itsBunch_m->getT();
2572 Vector_t meanR;
2573 Vector_t meanP;
2574 if (Options::psDumpFrame == DumpFrame::BUNCH_MEAN) {
2575 meanR = calcMeanR();
2576 meanP = calcMeanP();
2577 } else if (itsBunch_m->getLocalNum() > 0) {
2578 meanR = itsBunch_m->R[0];
2579 meanP = itsBunch_m->P[0];
2580 }
2581 double phi = 0;
2582 double psi = 0;
2583
2584 // Bunch (global) angle w.r.t. x-axis (cylinder coordinates)
2585 double azimuth = calculateAngle(meanR(0), meanR(1)) * Units::rad2deg;
2586
2587 // fix azimuth_m --> make monotonically increasing
2588 dumpAngle(azimuth, prevAzimuth_m, azimuth_m);
2589
2590 // -------------- Calculate the external fields at the center of the bunch ----------------- //
2591 beamline_list::iterator DumpSindex = FieldDimensions.begin();
2592
2593 extE_m = Vector_t(0.0, 0.0, 0.0);
2594 extB_m = Vector_t(0.0, 0.0, 0.0);
2595
2596 (((*DumpSindex)->second).second)->apply(meanR, meanP, temp_t, extE_m, extB_m);
2597
2598 // If we are saving in local frame, bunch and fields at the bunch center have to be rotated
2599 // TODO: Make decision if we maybe want to always save statistics data in local frame? -DW
2600 if (Options::psDumpFrame != DumpFrame::GLOBAL) {
2601 // -------------------- ----------- Do Transformations ---------------------------------- //
2602 // Bunch (local) azimuth at meanR w.r.t. y-axis
2603 phi = calculateAngle(meanP(0), meanP(1)) - 0.5 * Physics::pi;
2604
2605 // Bunch (local) elevation at meanR w.r.t. xy plane
2606 psi = 0.5 * Physics::pi - std::acos(meanP(2) / euclidean_norm(meanP));
2607
2608 // Rotate so Pmean is in positive y direction. No shift, so that normalized emittance and
2609 // unnormalized emittance as well as centroids are calculated correctly in
2610 // PartBunch::calcBeamParameters()
2611 globalToLocal(extB_m, phi, psi);
2612 globalToLocal(extE_m, phi, psi);
2613 globalToLocal(itsBunch_m->R, phi, psi);
2614 globalToLocal(itsBunch_m->P, phi, psi);
2615 }
2616
2617 FDext_m[0] = extB_m * Units::kG2T;
2618 FDext_m[1] = extE_m; // kV/mm? -DW
2619
2620 // Save the stat file
2621 itsDataSink->dumpSDDS(itsBunch_m, FDext_m, azimuth_m);
2622
2623 // If we are in local mode, transform back after saving
2624 if (Options::psDumpFrame != DumpFrame::GLOBAL) {
2625 localToGlobal(itsBunch_m->R, phi, psi);
2626 localToGlobal(itsBunch_m->P, phi, psi);
2627 }
2628
2629 IpplTimings::stopTimer(DumpTimer_m);
2630}
2631
2632
2633void ParallelCyclotronTracker::bunchDumpPhaseSpaceData() {
2634 // --------------------------------- Particle dumping --------------------------------------- //
2635 // Note: Don't dump when
2636 // 1. after one turn
2637 // in order to sychronize the dump step for multi-bunch and single bunch for comparison
2638 // with each other during post-process phase.
2639 // ------------------------------------------------------------------------------------------ //
2640 IpplTimings::startTimer(DumpTimer_m);
2641
2642 // --------------------------------- Get some Values ---------------------------------------- //
2643 double const temp_t = itsBunch_m->getT();
2644
2645 Vector_t meanR;
2646 Vector_t meanP;
2647
2648 // in case of multi-bunch mode take always bunch mean (although it takes all bunches)
2649 if (Options::psDumpFrame == DumpFrame::BUNCH_MEAN || isMultiBunch()) {
2650 meanR = calcMeanR();
2651 meanP = calcMeanP();
2652 } else if (itsBunch_m->getLocalNum() > 0) {
2653 meanR = itsBunch_m->R[0];
2654 meanP = itsBunch_m->P[0];
2655 }
2656
2657 double const betagamma_temp = euclidean_norm(meanP);
2658
2659 double const E = itsBunch_m->get_meanKineticEnergy();
2660
2661 // Bunch (global) angle w.r.t. x-axis (cylinder coordinates)
2662 double const theta = std::atan2(meanR(1), meanR(0));
2663
2664 // Bunch (local) azimuth at meanR w.r.t. y-axis
2665 double const phi = calculateAngle(meanP(0), meanP(1)) - 0.5 * Physics::pi;
2666
2667 // Bunch (local) elevation at meanR w.r.t. xy plane
2668 double const psi = 0.5 * Physics::pi - std::acos(meanP(2) / betagamma_temp);
2669
2670 // ---------------- Re-calculate reference values in format of input values ----------------- //
2671 // Position:
2672 // New OPAL 2.0: Init in m (back to mm just for output) -DW
2673 referenceR = computeRadius(meanR); // includes m --> mm conversion
2674 referenceTheta = Units::rad2deg * theta;
2675 referenceZ = Units::m2mm * meanR(2);
2676
2677 // Momentum in Theta-hat, R-hat, Z-hat coordinates at position meanR:
2678 referencePtot = betagamma_temp;
2679 referencePz = meanP(2);
2680 referencePr = meanP(0) * std::cos(theta) + meanP(1) * std::sin(theta);
2681 referencePt = std::sqrt(referencePtot * referencePtot - \
2682 referencePz * referencePz - referencePr * referencePr);
2683
2684 // ----- Calculate the external fields at the center of the bunch (Cave: Global Frame) ----- //
2685 beamline_list::iterator DumpSindex = FieldDimensions.begin();
2686
2687 extE_m = Vector_t(0.0, 0.0, 0.0);
2688 extB_m = Vector_t(0.0, 0.0, 0.0);
2689
2690 (((*DumpSindex)->second).second)->apply(meanR, meanP, temp_t, extE_m, extB_m);
2691 FDext_m[0] = extB_m * Units::kG2T;
2692 FDext_m[1] = extE_m;
2693
2694 // -------------- If flag psDumpFrame is global, dump bunch in global frame ------------- //
2695 if (Options::psDumpFreq < std::numeric_limits<int>::max() ){
2696 bool dumpLocalFrame = true;
2697 std::string dumpString = "local";
2698 if (Options::psDumpFrame == DumpFrame::GLOBAL) {
2699 dumpLocalFrame = false;
2700 dumpString = "global";
2701 }
2702
2703 if (dumpLocalFrame == true) {
2704 // ---------------- If flag psDumpFrame is local, dump bunch in local frame ---------------- //
2705
2706 // The bunch is transformed into a local coordinate system with meanP in direction of y-axis
2707 globalToLocal(itsBunch_m->R, phi, psi, meanR);
2708 globalToLocal(itsBunch_m->P, phi, psi); // P should only be rotated
2709
2710 globalToLocal(extB_m, phi, psi);
2711 globalToLocal(extE_m, phi, psi);
2712 }
2713
2714 FDext_m[0] = extB_m * Units::kG2T;
2715 FDext_m[1] = extE_m;
2716
2717 lastDumpedStep_m = itsDataSink->dumpH5(itsBunch_m, // Local and in m
2718 FDext_m, E,
2719 referencePr,
2720 referencePt,
2721 referencePz,
2722 referenceR,
2723 referenceTheta,
2724 referenceZ,
2725 phi * Units::rad2deg, // P_mean azimuth
2726 // at ref. R/Th/Z
2727 psi * Units::rad2deg, // P_mean elevation
2728 // at ref. R/Th/Z
2729 dumpLocalFrame); // Flag localFrame
2730
2731 if (dumpLocalFrame == true) {
2732 // Return to global frame
2733 localToGlobal(itsBunch_m->R, phi, psi, meanR);
2734 localToGlobal(itsBunch_m->P, phi, psi);
2735 }
2736
2737 // Tell user in which mode we are dumping
2738 // New: no longer dumping for num part < 3, omit phase space dump number info
2739 if (lastDumpedStep_m == -1) {
2740 *gmsg << endl << "* Integration step " << step_m + 1
2741 << " (no phase space dump for <= 2 particles)" << endl;
2742 } else {
2743 *gmsg << endl << "* Phase space dump " << lastDumpedStep_m
2744 << " (" << dumpString << " frame) at integration step " << step_m + 1 << endl;
2745 }
2746 }
2747
2748 // Print dump information on screen
2749 *gmsg << "* T = " << temp_t << " s"
2750 << ", Live Particles: " << itsBunch_m->getTotalNum() << endl;
2751 *gmsg << "* E = " << E << " MeV"
2752 << ", beta * gamma = " << betagamma_temp << endl;
2753 *gmsg << "* Bunch position: R = " << referenceR << " mm"
2754 << ", Theta = " << referenceTheta << " Deg"
2755 << ", Z = " << referenceZ << " mm" << endl;
2756 *gmsg << "* Local Azimuth = " << phi * Units::rad2deg << " Deg"
2757 << ", Local Elevation = " << psi * Units::rad2deg << " Deg" << endl;
2758
2759 IpplTimings::stopTimer(DumpTimer_m);
2760}
2761
2762bool ParallelCyclotronTracker::isTurnDone() {
2763 return (step_m > 10) && (((step_m + 1) %setup_m.stepsPerTurn) == 0);
2764}
2765
2766void ParallelCyclotronTracker::update_m(double& t, const double& dt,
2767 const bool& finishedTurn)
2768{
2769 // Reference parameters
2770 t += dt;
2771
2772 updateTime(dt);
2773
2774 itsBunch_m->setLocalTrackStep((step_m + 1));
2775 if (!(step_m + 1 % 1000)) {
2776 *gmsg << "Step " << step_m + 1 << endl;
2777 }
2778
2779 updatePathLength(dt);
2780
2781 // Here is global frame, don't do itsBunch_m->boundp();
2782
2783 if (itsBunch_m->getTotalNum()>0) {
2784 // Only dump last step if we have particles left.
2785 // Check separately for phase space (ps) and statistics (stat) data dump frequency
2786 if ( mode_m != TrackingMode::SEO && ( ((step_m + 1) % Options::psDumpFreq == 0) ||
2787 (Options::psDumpEachTurn && finishedTurn)))
2788 {
2789 // Write phase space data to h5 file
2790 bunchDumpPhaseSpaceData();
2791 }
2792
2793 if ( mode_m != TrackingMode::SEO && ( ((step_m + 1) % Options::statDumpFreq == 0) ||
2794 (Options::psDumpEachTurn && finishedTurn)))
2795 {
2796 // Write statistics data to stat file
2797 bunchDumpStatData();
2798 }
2799 }
2800
2801 if (Options::psDumpEachTurn && finishedTurn) {
2802 for (PluginElement* element : pluginElements_m) {
2803 element->save();
2804 }
2805 }
2806}
2807
2808
2809std::tuple<double, double, double> ParallelCyclotronTracker::initializeTracking_m() {
2810 // Read in some control parameters
2811 setup_m.scSolveFreq = (spiral_flag) ? 1 : Options::scSolveFreq;
2812 setup_m.stepsPerTurn = itsBunch_m->getStepsPerTurn();
2813
2814 // Define 3 special azimuthal angles where dump particle's six parameters
2815 // at each turn into 3 ASCII files. only important in single particle tracking
2816 azimuth_angle_m.resize(3);
2817 azimuth_angle_m[0] = 0.0;
2818 azimuth_angle_m[1] = 22.5 * Units::deg2rad;
2819 azimuth_angle_m[2] = 45.0 * Units::deg2rad;
2820
2821 double harm = getHarmonicNumber();
2822 double dt = itsBunch_m->getdT() * harm;
2823 double t = itsBunch_m->getT();
2824
2825 double oldReferenceTheta = referenceTheta; // init here, reset each step
2826 setup_m.deltaTheta = Physics::pi / (setup_m.stepsPerTurn); // half of the average angle per step
2827
2828 //int stepToLastInj = itsBunch_m->getSteptoLastInj(); // TODO: Do we need this? -DW
2829
2830 // Record how many bunches have already been injected. ONLY FOR MBM
2831 if (isMultiBunch()) {
2832 mbHandler_m->setNumBunch(itsBunch_m->getNumBunch());
2833 }
2834
2835 initTrackOrbitFile();
2836
2837 // Get data from h5 file for restart run and reset current step
2838 // to last step of previous simulation
2839 if (OpalData::getInstance()->inRestartRun()) {
2840
2841 restartStep0_m = itsBunch_m->getLocalTrackStep();
2842 step_m = restartStep0_m;
2843 turnnumber_m = step_m / setup_m.stepsPerTurn + 1;
2844
2845 *gmsg << "* Restart at integration step " << restartStep0_m
2846 << " at turn " << turnnumber_m - 1 << endl;
2847
2848 initPathLength();
2849 }
2850
2851 setup_m.stepsNextCheck = step_m + setup_m.stepsPerTurn; // Steps to next check for transition
2852
2853 initDistInGlobalFrame();
2854
2855 if (isMultiBunch()) {
2856 mbHandler_m->updateParticleBins(itsBunch_m);
2857 }
2858
2859 // --- Output to user --- //
2860 *gmsg << "* Beginning of this run is at t = " << t * Units::s2ns << " [ns]" << endl;
2861 *gmsg << "* The time step is set to dt = " << dt * Units::s2ns << " [ns]" << endl;
2862
2863 if ( isMultiBunch() ) {
2864 *gmsg << "* MBM: Time interval between neighbour bunches is set to "
2865 << setup_m.stepsPerTurn * dt * Units::s2ns << "[ns]" << endl;
2866 *gmsg << "* MBM: The particles energy bin reset frequency is set to "
2867 << Options::rebinFreq << endl;
2868 }
2869
2870 *gmsg << "* Single particle trajectory dump frequency is set to " << Options::sptDumpFreq << endl;
2871 *gmsg << "* The frequency to solve space charge fields is set to " << setup_m.scSolveFreq << endl;
2872 *gmsg << "* The repartition frequency is set to " << Options::repartFreq << endl;
2873
2874 switch ( mode_m ) {
2875 case TrackingMode::SEO: {
2876 *gmsg << endl;
2877 *gmsg << "* ------------------------- STATIC EQUILIBRIUM ORBIT MODE ----------------------------- *" << endl
2878 << "* Instruction: When the total particle number is equal to 2, SEO mode is triggered *" << endl
2879 << "* automatically. This mode does NOT include any RF cavities. The initial distribution *" << endl
2880 << "* file must be specified. In the file the first line is for reference particle and the *" << endl
2881 << "* second line is for off-center particle. The tune is calculated by FFT routines based *" << endl
2882 << "* on these two particles. *" << endl
2883 << "* ---------------- NOTE: SEO MODE ONLY WORKS SERIALLY ON SINGLE NODE ------------------ *" << endl;
2884
2885 if (Ippl::getNodes() != 1)
2886 throw OpalException("Error in ParallelCyclotronTracker::initializeTracking_m",
2887 "SEO MODE ONLY WORKS SERIALLY ON SINGLE NODE");
2888 break;
2889 }
2890 case TrackingMode::SINGLE: {
2891 *gmsg << endl;
2892 *gmsg << "* ------------------------------ SINGLE PARTICLE MODE --------------------------------- *" << endl
2893 << "* Instruction: When the total particle number is equal to 1, single particle mode is *" << endl
2894 << "* triggered automatically. The initial distribution file must be specified which should *" << endl
2895 << "* contain only one line for the single particle *" << endl
2896 << "* ---------NOTE: SINGLE PARTICLE MODE ONLY WORKS SERIALLY ON A SINGLE NODE ------------ *" << endl;
2897
2898 if (Ippl::getNodes() != 1)
2899 throw OpalException("Error in ParallelCyclotronTracker::initializeTracking_m",
2900 "SINGLE PARTICLE MODE ONLY WORKS SERIALLY ON A SINGLE NODE");
2901
2902 // For single particle mode open output files
2903 openFiles(azimuth_angle_m.size() + 1, OpalData::getInstance()->getInputBasename());
2904 break;
2905 }
2906 case TrackingMode::BUNCH: {
2907 break;
2908 }
2909 case TrackingMode::UNDEFINED: {}
2910 default: {
2911 throw OpalException("ParallelCyclotronTracker::initializeTracking_m()",
2912 "No such tracking mode.");
2913 }
2914 }
2915
2916 return std::make_tuple(t, dt, oldReferenceTheta);
2917}
2918
2919
2920void ParallelCyclotronTracker::finalizeTracking_m(dvector_t& Ttime,
2921 dvector_t& Tdeltr,
2922 dvector_t& Tdeltz, ivector_t& TturnNumber) {
2923
2924 for (size_t ii = 0; ii < (itsBunch_m->getLocalNum()); ii++) {
2925 if (itsBunch_m->ID[ii] == 0) {
2926 double FinalMomentum2 = std::pow(itsBunch_m->P[ii](0), 2.0) + std::pow(itsBunch_m->P[ii](1), 2.0) + std::pow(itsBunch_m->P[ii](2), 2.0);
2927 double FinalEnergy = (std::sqrt(1.0 + FinalMomentum2) - 1.0) * itsBunch_m->getM() * Units::eV2MeV;
2928 *gmsg << "* Final energy of reference particle = " << FinalEnergy << " [MeV]" << endl;
2929 *gmsg << "* Total phase space dump number(includes the initial distribution) = " << lastDumpedStep_m + 1 << endl;
2930 *gmsg << "* One can restart simulation from the last dump step (--restart " << lastDumpedStep_m << ")" << endl;
2931 }
2932 }
2933
2934 Ippl::Comm->barrier();
2935
2936 switch ( mode_m ) {
2937 case TrackingMode::SEO: {
2938 // Calculate tunes after tracking.
2939 *gmsg << endl;
2940 *gmsg << "* **************** The result for tune calculation (NO space charge) ******************* *" << endl
2941 << "* Number of tracked turns: " << TturnNumber.back() << endl;
2942 double nur, nuz;
2943 getTunes(Ttime, Tdeltr, Tdeltz, TturnNumber.back(), nur, nuz);
2944 break;
2945 }
2946 case TrackingMode::SINGLE:
2947 closeFiles();
2948 // fall through
2949 case TrackingMode::BUNCH: // we do nothing
2950 case TrackingMode::UNDEFINED:
2951 default: {
2952 // not for multibunch
2953 if (!isMultiBunch()) {
2954 *gmsg << "*" << endl;
2955 *gmsg << "* Finished during turn " << turnnumber_m << " (" << turnnumber_m - 1 << " turns completed)" << endl;
2956 *gmsg << "* Cave: Turn number is not correct for restart mode"<< endl;
2957 }
2958 }
2959 }
2960
2961 Ippl::Comm->barrier();
2962
2963 if (myNode_m == 0) {
2964 outfTrackOrbit_m.close();
2965 }
2966
2967 *gmsg << endl << "* *********************** Bunch information in global frame: ***********************";
2968
2969 if (itsBunch_m->getTotalNum() > 0) {
2970 // Print out the Bunch information at end of the run.
2971 itsBunch_m->calcBeamParameters();
2972 *gmsg << *itsBunch_m << endl;
2973 } else {
2974 *gmsg << endl << "* No Particles left in bunch" << endl;
2975 *gmsg << "* **********************************************************************************" << endl;
2976 }
2977}
2978
2979void ParallelCyclotronTracker::setTrackingMode() {
2980 if ( initialTotalNum_m == 1 ) {
2981 mode_m = TrackingMode::SINGLE;
2982 } else if ( initialTotalNum_m == 2 ) {
2983 mode_m = TrackingMode::SEO;
2984 } else if ( initialTotalNum_m > 2 ) {
2985 mode_m = TrackingMode::BUNCH;
2986 } else {
2987 mode_m = TrackingMode::UNDEFINED;
2988 }
2989}
2990
2991void ParallelCyclotronTracker::seoMode_m(double& t, const double dt, bool& /*finishedTurn*/,
2992 dvector_t& Ttime, dvector_t& Tdeltr,
2993 dvector_t& Tdeltz, ivector_t& TturnNumber) {
2994
2995 // 2 particles: Trigger SEO mode
2996 // (Switch off cavity and calculate betatron oscillation tuning)
2997 double r_tuning[2], z_tuning[2];
2998
2999 IpplTimings::startTimer(IntegrationTimer_m);
3000 for (size_t i = 0; i < (itsBunch_m->getLocalNum()); i++) {
3001
3002 if ((step_m % Options::sptDumpFreq == 0)) {
3003 outfTrackOrbit_m << "ID" << (itsBunch_m->ID[i]);
3004 outfTrackOrbit_m << " " << itsBunch_m->R[i](0)
3005 << " " << itsBunch_m->P[i](0)
3006 << " " << itsBunch_m->R[i](1)
3007 << " " << itsBunch_m->P[i](1)
3008 << " " << itsBunch_m->R[i](2)
3009 << " " << itsBunch_m->P[i](2)
3010 << std::endl;
3011 }
3012
3013 double OldTheta = calculateAngle(itsBunch_m->R[i](0), itsBunch_m->R[i](1));
3014 r_tuning[i] = itsBunch_m->R[i](0) * std::cos(OldTheta) +
3015 itsBunch_m->R[i](1) * std::sin(OldTheta);
3016
3017 z_tuning[i] = itsBunch_m->R[i](2);
3018
3019 // Integrate for one step in the lab Cartesian frame (absolute value).
3020 itsStepper_mp->advance(itsBunch_m, i, t, dt);
3021
3022 if ( (i == 0) && isTurnDone() ) {
3023 turnnumber_m++;
3024 }
3025
3026 } // end for: finish one step tracking for all particles for initialTotalNum_m == 2 mode
3027 IpplTimings::stopTimer(IntegrationTimer_m);
3028
3029 // store dx and dz for future tune calculation if higher precision needed, reduce freqSample.
3030 if (step_m % Options::sptDumpFreq == 0) {
3031 Ttime.push_back(t);
3032 Tdeltz.push_back(z_tuning[1]);
3033 Tdeltr.push_back(r_tuning[1] - r_tuning[0]);
3034 TturnNumber.push_back(turnnumber_m);
3035 }
3036}
3037
3038
3039void ParallelCyclotronTracker::singleMode_m(double& t, const double dt,
3040 bool& finishedTurn, double& oldReferenceTheta) {
3041 // 1 particle: Trigger single particle mode
3042
3043 // ********************************************************************************** //
3044 // * This was moved here b/c collision should be tested before the actual * //
3045 // * timestep (bgf_main_collision_test() predicts the next step automatically) * //
3046
3047 // apply the plugin elements: probe, collimator, stripper, septum
3048 bool flagNeedUpdate = applyPluginElements(dt);
3049
3050 // check if we lose particles at the boundary
3051 bgf_main_collision_test();
3052 // ********************************************************************************** //
3053
3054 if (itsBunch_m->getLocalNum() == 0) return; // might happen if particle is in collimator
3055
3056 IpplTimings::startTimer(IntegrationTimer_m);
3057
3058 unsigned int i = 0; // we only have a single particle
3059 if ( step_m % Options::sptDumpFreq == 0 ) {
3060 outfTrackOrbit_m << "ID" <<itsBunch_m->ID[i]
3061 << " " << itsBunch_m->R[i](0)
3062 << " " << itsBunch_m->P[i](0)
3063 << " " << itsBunch_m->R[i](1)
3064 << " " << itsBunch_m->P[i](1)
3065 << " " << itsBunch_m->R[i](2)
3066 << " " << itsBunch_m->P[i](2)
3067 << std::endl;
3068 }
3069
3070 double temp_meanTheta = calculateAngle2(itsBunch_m->R[i](0),
3071 itsBunch_m->R[i](1)); // [-pi, pi]
3072
3073 dumpThetaEachTurn_m(t, itsBunch_m->R[i], itsBunch_m->P[i],
3074 temp_meanTheta, finishedTurn);
3075
3076 dumpAzimuthAngles_m(t, itsBunch_m->R[i], itsBunch_m->P[i],
3077 oldReferenceTheta, temp_meanTheta);
3078
3079 oldReferenceTheta = temp_meanTheta;
3080
3081 // used for gap crossing checking
3082 Vector_t Rold = itsBunch_m->R[i]; // [x,y,z] (m)
3083 Vector_t Pold = itsBunch_m->P[i]; // [px,py,pz] (beta*gamma)
3084
3085 // integrate for one step in the lab Cartesian frame (absolute value).
3086 itsStepper_mp->advance(itsBunch_m, i, t, dt);
3087
3088 // If gap crossing happens, do momenta kicking (not if gap crossing just happened)
3089 if (itsBunch_m->cavityGapCrossed[i] == true) {
3090 itsBunch_m->cavityGapCrossed[i] = false;
3091 } else {
3092 gapCrossKick_m(i, t, dt, Rold, Pold);
3093 }
3094 IpplTimings::stopTimer(IntegrationTimer_m);
3095
3096 // Destroy particles if they are marked as Bin = -1 in the plugin elements
3097 // or out of the global aperture
3098 flagNeedUpdate |= (itsBunch_m->Bin[i] < 0);
3099 deleteParticle(flagNeedUpdate);
3100}
3101
3102
3103void ParallelCyclotronTracker::bunchMode_m(double& t, const double dt, bool& finishedTurn) {
3104
3105 // Flag for transition from single-bunch to multi-bunches mode
3106 static bool flagTransition = false;
3107
3108 // single particle dumping
3109 if (step_m % Options::sptDumpFreq == 0)
3110 singleParticleDump();
3111
3112 // Find out if we need to do bunch injection
3113 injectBunch(flagTransition);
3114
3115// oldReferenceTheta = calculateAngle2(meanR(0), meanR(1));
3116
3117 // Calculate SC field before each time step and keep constant during integration.
3118 // Space Charge effects are included only when total macroparticles number is NOT LESS THAN 1000.
3119 if (itsBunch_m->hasFieldSolver()) {
3120
3121 if (step_m % setup_m.scSolveFreq == 0) {
3122 computeSpaceChargeFields_m();
3123 } else {
3124 // If we are not solving for the space charge fields at this time step
3125 // we will apply the fields from the previous step and have to rotate them
3126 // accordingly. For this we find the quaternion between the previous mean momentum (PreviousMeanP)
3127 // and the current mean momentum (meanP) and rotate the fields with this quaternion.
3128
3129 Vector_t const meanP = calcMeanP();
3130
3131 Quaternion_t quaternionToNewMeanP;
3132
3133 getQuaternionTwoVectors(PreviousMeanP, meanP, quaternionToNewMeanP);
3134
3135 // Reset PreviousMeanP. Cave: This HAS to be after the quaternion is calculated!
3136 PreviousMeanP = calcMeanP();
3137
3138 // Rotate the fields
3139 globalToLocal(itsBunch_m->Ef, quaternionToNewMeanP);
3140 globalToLocal(itsBunch_m->Bf, quaternionToNewMeanP);
3141 }
3142 }
3143
3144 // *** This was moved here b/c collision should be tested before the **********************
3145 // *** actual timestep (bgf_main_collision_test() predicts the next step automatically) -DW
3146 // Apply the plugin elements: probe, collimator, stripper, septum
3147 bool flagNeedUpdate = applyPluginElements(dt);
3148
3149 // check if we lose particles at the boundary
3150 bgf_main_collision_test();
3151
3152 IpplTimings::startTimer(IntegrationTimer_m);
3153
3154 for (size_t i = 0; i < itsBunch_m->getLocalNum(); i++) {
3155
3156 // used for gap crossing checking
3157 Vector_t Rold = itsBunch_m->R[i]; // [x,y,z] (m)
3158 Vector_t Pold = itsBunch_m->P[i]; // [px,py,pz] (beta*gamma)
3159
3160 // Integrate for one step in the lab Cartesian frame (absolute value).
3161 itsStepper_mp->advance(itsBunch_m, i, t, dt);
3162
3163 // If gap crossing happens, do momenta kicking (not if gap crossing just happened)
3164 if (itsBunch_m->cavityGapCrossed[i] == true) {
3165 itsBunch_m->cavityGapCrossed[i] = false;
3166 } else {
3167 gapCrossKick_m(i, t, dt, Rold, Pold);
3168 }
3169 flagNeedUpdate |= (itsBunch_m->Bin[i] < 0);
3170 }
3171
3172 IpplTimings::stopTimer(IntegrationTimer_m);
3173
3174 // Destroy particles if they are marked as Bin = -1 in the plugin elements
3175 // or out of global aperture
3176 deleteParticle(flagNeedUpdate);
3177
3178 // Recalculate bingamma and reset the BinID for each particles according to its current gamma
3179 if (isMultiBunch() && step_m % Options::rebinFreq == 0) {
3180 mbHandler_m->updateParticleBins(itsBunch_m);
3181 }
3182
3183 // Some status output for user after each turn
3184 if ( isTurnDone() ) {
3185 turnnumber_m++;
3186 finishedTurn = true;
3187
3188 *gmsg << endl;
3189 *gmsg << "*** Finished turn " << turnnumber_m - 1
3190 << ", Total number of live particles: "
3191 << itsBunch_m->getTotalNum() << endl;
3192 }
3193
3194 Ippl::Comm->barrier();
3195}
3196
3197
3198void ParallelCyclotronTracker::gapCrossKick_m(size_t i, double t,
3199 double dt,
3200 const Vector_t& Rold,
3201 const Vector_t& Pold) {
3202
3203 for (beamline_list::iterator sindex = ++(FieldDimensions.begin());
3204 sindex != FieldDimensions.end(); ++sindex)
3205 {
3206 bool tag_crossing = false;
3207 double DistOld = 0.0;
3208 RFCavity* rfcav;
3209
3210 if (((*sindex)->first) == ElementType::RFCAVITY) {
3211 // here check gap cross in the list, if do, set tag_crossing to TRUE
3212 rfcav = static_cast<RFCavity*>(((*sindex)->second).second);
3213 tag_crossing = checkGapCross(Rold, itsBunch_m->R[i], rfcav, DistOld);
3214 }
3215
3216 if ( tag_crossing ) {
3217 itsBunch_m->cavityGapCrossed[i] = true;
3218
3219 double oldBeta = euclidean_norm(Util::getBeta(Pold));
3220 double dt1 = DistOld / (oldBeta * Physics::c);
3221 double dt2 = dt - dt1;
3222
3223 // retrack particle from the old postion to cavity gap point
3224 // restore the old coordinates and momenta
3225 itsBunch_m->R[i] = Rold;
3226 itsBunch_m->P[i] = Pold;
3227
3228 if (dt / dt1 < 1.0e9) {
3229 itsStepper_mp->advance(itsBunch_m, i, t, dt1);
3230 }
3231 // Momentum kick
3232 RFkick(rfcav, t, dt1, i);
3233
3234 /* Retrack particle from cavity gap point for
3235 * the left time to finish the entire timestep
3236 */
3237 if (dt / dt2 < 1.0e9) {
3238 itsStepper_mp->advance(itsBunch_m, i, t, dt2);
3239 }
3240 }
3241 }
3242}
3243
3244
3245void ParallelCyclotronTracker::dumpAzimuthAngles_m(const double& t,
3246 const Vector_t& R,
3247 const Vector_t& P,
3248 const double& oldReferenceTheta,
3249 const double& temp_meanTheta) {
3250
3251 for (unsigned int i=0; i<=2; i++) {
3252 if ((oldReferenceTheta < azimuth_angle_m[i] - setup_m.deltaTheta) &&
3253 ( temp_meanTheta >= azimuth_angle_m[i] - setup_m.deltaTheta))
3254 {
3255 *(outfTheta_m[i]) << "#Turn number = " << turnnumber_m
3256 << ", Time = " << t * Units::s2ns << " [ns]" << std::endl
3257 << " " << std::hypot(R(0), R(1))
3258 << " " << P(0) * std::cos(temp_meanTheta) + P(1) * std::sin(temp_meanTheta)
3259 << " " << temp_meanTheta * Units::rad2deg
3260 << " " << - P(0) * std::sin(temp_meanTheta) + P(1) * std::cos(temp_meanTheta)
3261 << " " << R(2)
3262 << " " << P(2) << std::endl;
3263 }
3264 }
3265}
3266
3267void ParallelCyclotronTracker::dumpThetaEachTurn_m(const double& t,
3268 const Vector_t& R,
3269 const Vector_t& P,
3270 const double& temp_meanTheta,
3271 bool& finishedTurn) {
3272 if ( isTurnDone() ) {
3273 ++turnnumber_m;
3274 finishedTurn = true;
3275
3276 *gmsg << "* SPT: Finished turn " << turnnumber_m - 1 << endl;
3277
3278 *(outfTheta_m[3]) << "#Turn number = " << turnnumber_m
3279 << ", Time = " << t * Units::s2ns << " [ns]" << std::endl
3280 << " " << std::sqrt(R(0) * R(0) + R(1) * R(1))
3281 << " " << P(0) * std::cos(temp_meanTheta) +
3282 P(1) * std::sin(temp_meanTheta)
3283 << " " << temp_meanTheta / Units::deg2rad
3284 << " " << - P(0) * std::sin(temp_meanTheta) +
3285 P(1) * std::cos(temp_meanTheta)
3286 << " " << R(2)
3287 << " " << P(2) << std::endl;
3288 }
3289}
3290
3291
3292void ParallelCyclotronTracker::computeSpaceChargeFields_m() {
3293 // Firstly reset E and B to zero before fill new space charge field data for each track step
3294 itsBunch_m->Bf = Vector_t(0.0);
3295 itsBunch_m->Ef = Vector_t(0.0);
3296
3297 if (spiral_flag && itsBunch_m->getFieldSolverType() == FieldSolverType::SAAMG) {
3298 // --- Single bunch mode with spiral inflector --- //
3299
3300 // If we are doing a spiral inflector simulation and are using the SAAMG solver
3301 // we don't rotate or translate the bunch and gamma is 1.0 (non-relativistic).
3302
3303 itsBunch_m->setGlobalMeanR(Vector_t(0.0, 0.0, 0.0));
3304 itsBunch_m->setGlobalToLocalQuaternion(Quaternion_t(1.0, 0.0, 0.0, 0.0));
3305
3306 itsBunch_m->computeSelfFields_cycl(1.0);
3307
3308 } else {
3309
3310 Vector_t const meanR = calcMeanR();
3311
3312 PreviousMeanP = calcMeanP();
3313
3314 // Since calcMeanP takes into account all particles of all bins (TODO: Check this! -DW)
3315 // Using the quaternion method with PreviousMeanP and yaxis should give the correct result
3316 Quaternion_t quaternionToYAxis;
3317
3318 getQuaternionTwoVectors(PreviousMeanP, yaxis, quaternionToYAxis);
3319
3320 globalToLocal(itsBunch_m->R, quaternionToYAxis, meanR);
3321
3322 if ((step_m + 1) % Options::boundpDestroyFreq == 0) {
3323 itsBunch_m->boundp_destroyCycl();
3324 } else {
3325 itsBunch_m->boundp();
3326 }
3327
3328 if (hasMultiBunch()) {
3329 // --- Multibunch mode --- //
3330
3331 // Calculate gamma for each energy bin
3332 itsBunch_m->calcGammas_cycl();
3333
3334 repartition();
3335
3336 // Calculate space charge field for each energy bin
3337 for (int b = 0; b < itsBunch_m->getLastemittedBin(); b++) {
3338 itsBunch_m->setBinCharge(b, itsBunch_m->getChargePerParticle());
3339 itsBunch_m->setGlobalMeanR(meanR);
3340 itsBunch_m->setGlobalToLocalQuaternion(quaternionToYAxis);
3341 itsBunch_m->computeSelfFields_cycl(b);
3342 }
3343
3344 itsBunch_m->Q = itsBunch_m->getChargePerParticle();
3345
3346 } else {
3347 // --- Single bunch mode --- //
3348 double temp_meangamma = Util::getGamma(PreviousMeanP);
3349
3350 repartition();
3351
3352 itsBunch_m->setGlobalMeanR(meanR);
3353 itsBunch_m->setGlobalToLocalQuaternion(quaternionToYAxis);
3354 itsBunch_m->computeSelfFields_cycl(temp_meangamma);
3355 }
3356
3357 // Transform coordinates back to global
3358 localToGlobal(itsBunch_m->R, quaternionToYAxis, meanR);
3359
3360 // Transform self field back to global frame (rotate only)
3361 localToGlobal(itsBunch_m->Ef, quaternionToYAxis);
3362 localToGlobal(itsBunch_m->Bf, quaternionToYAxis);
3363 }
3364}
3365
3366
3367bool ParallelCyclotronTracker::computeExternalFields_m(const size_t& i, const double& t,
3368 Vector_t& Efield, Vector_t& Bfield) {
3369
3370 beamline_list::iterator sindex = FieldDimensions.begin();
3371
3372 // Flag whether a particle is out of field
3373 bool outOfBound = (((*sindex)->second).second)->apply(i, t, Efield, Bfield);
3374
3375 Bfield *= Units::kG2T;
3376 Efield *= Units::kV2V / Units::mm2m;
3377
3378 return outOfBound;
3379}
3380
3381bool ParallelCyclotronTracker::computeExternalFields_m(const Vector_t& R, const Vector_t& P, const double& t,
3382 Vector_t& Efield, Vector_t& Bfield) {
3383
3384 beamline_list::iterator sindex = FieldDimensions.begin();
3385 // Flag whether a particle is out of field
3386 bool outOfBound = (((*sindex)->second).second)->apply(R, P, t, Efield, Bfield);
3387
3388 Bfield *= Units::kG2T;
3389 Efield *= Units::kV2V / Units::mm2m;
3390
3391 return outOfBound;
3392}
3393
3394
3395
3396void ParallelCyclotronTracker::injectBunch(bool& flagTransition) {
3397 if (!isMultiBunch() || step_m != setup_m.stepsNextCheck) {
3398 return;
3399 }
3400
3401 const short result = mbHandler_m->injectBunch(itsBunch_m, itsReference,
3402 flagTransition);
3403
3404 switch ( result ) {
3405 case 0: {
3406 // nothing happened
3407 break;
3408 }
3409 case 1: {
3410 // bunch got saved
3411 saveInjectValues();
3412 setup_m.stepsNextCheck += setup_m.stepsPerTurn;
3413 if ( flagTransition ) {
3414 *gmsg << "* MBM: Saving beam distribution at turn " << turnnumber_m << endl;
3415 *gmsg << "* MBM: After one revolution, Multi-Bunch Mode will be invoked" << endl;
3416 }
3417 break;
3418 }
3419 case 2: {
3420 // bunch got injected
3421 setup_m.stepsNextCheck += setup_m.stepsPerTurn;
3422 break;
3423 }
3424 default: {
3425 throw OpalException("ParallelCyclotronTracker::injectBunch()",
3426 "Unknown return value " + std::to_string(result));
3427 }
3428 }
3429}
3430
3431
3432void ParallelCyclotronTracker::saveInjectValues() {
3433 if ( !isMultiBunch() ) {
3434 return;
3435 }
3436
3437 Vector_t meanR = calcMeanR();
3438
3439 // Bunch (global) angle w.r.t. x-axis (cylinder coordinates)
3440 double theta = calculateAngle(meanR(0), meanR(1)) * Units::rad2deg;
3441
3442 // fix azimuth_m --> make monotonically increasing
3443 dumpAngle(theta, prevAzimuth_m, azimuth_m);
3444
3445 const double radius = computeRadius(meanR);
3446
3447 MultiBunchHandler::injection_t& inj = mbHandler_m->getInjectionValues();
3448
3449 inj.time = itsBunch_m->getT();
3450 inj.pathlength = itsBunch_m->get_sPos();
3451 inj.azimuth = azimuth_m;
3452 inj.radius = radius;
3453}
3454
3455
3456void ParallelCyclotronTracker::updatePathLength(const double& dt) {
3457 /* the last element includes all particles,
3458 * all other elements are bunch number specific
3459 */
3460 std::vector<double> lpaths(1);
3461
3462 if ( isMultiBunch() ) {
3463 lpaths.resize(mbHandler_m->getNumBunch() + 1);
3464 }
3465
3466 computePathLengthUpdate(lpaths, dt);
3467
3468 pathLength_m += lpaths.back();
3469 itsBunch_m->set_sPos(pathLength_m);
3470
3471 if ( isMultiBunch() ) {
3472 mbHandler_m->updatePathLength(lpaths);
3473 }
3474}
3475
3476
3477void ParallelCyclotronTracker::updateTime(const double& dt) {
3478 double t = itsBunch_m->getT();
3479
3480 itsBunch_m->setT(t + dt);
3481
3482 if ( isMultiBunch() ) {
3483 mbHandler_m->updateTime(dt);
3484 }
3485}
3486
3487
3488void ParallelCyclotronTracker::updateAzimuthAndRadius() {
3489 if (!isMultiBunch()) {
3490 return;
3491 }
3492
3493 for (short b = 0; b < mbHandler_m->getNumBunch(); ++b) {
3494 Vector_t meanR = calcMeanR(b);
3495 MultiBunchHandler::beaminfo_t& binfo = mbHandler_m->getBunchInfo(b);
3496
3497 binfo.radius = computeRadius(meanR);
3498 double azimuth = calculateAngle(meanR(0), meanR(1)) * Units::rad2deg;
3499 dumpAngle(azimuth, binfo.prevAzimuth, binfo.azimuth, b);
3500 }
3501}
3502
3503
3504void ParallelCyclotronTracker::initPathLength() {
3505 if ( isMultiBunch() ) {
3506 // we need to reset the path length of each bunch
3507 itsDataSink->setMultiBunchInitialPathLengh(mbHandler_m.get());
3508 }
3509}
cross
Vector3D cross(const Vector3D &lhs, const Vector3D &rhs)
Vector cross product.
Definition Vector3D.cpp:111
dot
double dot(const Vector3D &lhs, const Vector3D &rhs)
Vector dot product.
Definition Vector3D.cpp:118
ElementType
ElementType
Definition ElementBase.h:88
Quaternion_t
Quaternion Quaternion_t
Definition Quaternion.h:42
DumpFrame::BUNCH_MEAN
@ BUNCH_MEAN
Definition Options.h:28
FlaggedBeamline
TBeamline< FlaggedElmPtr > FlaggedBeamline
A beam line with flagged elements.
Definition FlaggedBeamline.h:34
euclidean_norm
T euclidean_norm(const Vector< T > &)
Euclidean norm.
Definition Vector.h:243
gmsgALL
Inform * gmsgALL
Definition Main.cpp:71
prod_boost_vector
T prod_boost_vector(const boost::numeric::ublas::matrix< double > &rotation, const T &vector)
Definition BoostMatrix.h:26
DistributionType
DistributionType
Definition Distribution.h:50
gmsg
Inform * gmsg
Definition Main.cpp:70
allreduce
void allreduce(const T *input, T *output, int count, Op op)
Definition GlobalComm.hpp:510
reduce
bool reduce(Communicate &, InputIterator, InputIterator, OutputIterator, const ReduceOp &, bool *IncludeVal=0)
Definition GlobalComm.hpp:55
COMM_ANY_NODE
const int COMM_ANY_NODE
Definition Communicate.h:40
IPPL_APP_CYCLE
#define IPPL_APP_CYCLE
Definition Tags.h:103
IPPL_APP_TAG4
#define IPPL_APP_TAG4
Definition Tags.h:97
for_each
bool for_each(const BareFieldIterator< T, D > &p, SameFieldID s, C)
Definition AssignDefs.h:30
sum
T::PETE_Expr_t::PETE_Return_t sum(const PETE_Expr< T > &expr)
Definition PETE.h:1111
level4
Inform & level4(Inform &inf)
Definition Inform.cpp:48
endl
Inform & endl(Inform &inf)
Definition Inform.cpp:42
level3
Inform & level3(Inform &inf)
Definition Inform.cpp:47
ERRORMSG
#define ERRORMSG(msg)
Definition IpplInfo.h:350
Ippl
IpplInfo Ippl
Definition IpplInfo.h:353
Physics::q_e
constexpr double q_e
The elementary charge in As.
Definition Physics.h:69
Physics::c
constexpr double c
The velocity of light in m/s.
Definition Physics.h:45
Physics::pi
constexpr double pi
The value of.
Definition Physics.h:30
Units::mm2m
constexpr double mm2m
Definition Units.h:29
Units::m2mm
constexpr double m2mm
Definition Units.h:26
Units::eV2MeV
constexpr double eV2MeV
Definition Units.h:77
Units::kV2V
constexpr double kV2V
Definition Units.h:62
Units::GeV2eV
constexpr double GeV2eV
Definition Units.h:68
Units::kG2T
constexpr double kG2T
Definition Units.h:59
Units::rad2deg
constexpr double rad2deg
Definition Units.h:146
Units::deg2rad
constexpr double deg2rad
Definition Units.h:143
Units::s2ns
constexpr double s2ns
Definition Units.h:44
Options::mtsSubsteps
int mtsSubsteps
Definition Options.cpp:59
Options::psDumpFreq
int psDumpFreq
The frequency to dump the phase space, i.e.dump data when steppsDumpFreq==0.
Definition Options.cpp:39
Options::boundpDestroyFreq
int boundpDestroyFreq
Definition Options.cpp:87
Options::psDumpEachTurn
bool psDumpEachTurn
phase space dump flag for OPAL-cycl
Definition Options.cpp:43
Options::asciidump
bool asciidump
Definition Options.cpp:85
Options::sptDumpFreq
int sptDumpFreq
The frequency to dump single particle trajectory of particles with ID = 0 & 1.
Definition Options.cpp:47
Options::delPartFreq
unsigned int delPartFreq
The frequency to delete particles (currently: OPAL-cycl only).
Definition Options.cpp:109
Options::scSolveFreq
int scSolveFreq
The frequency to solve space charge fields.
Definition Options.cpp:57
Options::repartFreq
int repartFreq
The frequency to do particles repartition for better load balance between nodes.
Definition Options.cpp:49
Options::rebinFreq
int rebinFreq
The frequency to reset energy bin ID for all particles.
Definition Options.cpp:55
Options::psDumpFrame
DumpFrame psDumpFrame
flag to decide in which coordinate frame the phase space will be dumped for OPAL-cycl
Definition Options.cpp:45
Options::statDumpFreq
int statDumpFreq
The frequency to dump statistical values, e.e. dump data when stepstatDumpFreq==0.
Definition Options.cpp:41
Util::angleBetweenAngles
bool angleBetweenAngles(const double angle, const double min, const double max)
check if angle (in rad and in range [0,2pi]) is within [min, max]
Definition Util.h:206
Util::doubleVectorToString
std::string doubleVectorToString(const std::vector< double > &v)
Definition Util.cpp:176
Util::boolVectorToUpperString
std::string boolVectorToUpperString(const std::vector< bool > &b)
Definition Util.cpp:161
Util::getBeta
Vector_t getBeta(Vector_t p)
Definition Util.h:51
Util::getGamma
double getGamma(Vector_t p)
Definition Util.h:46
Steppers::TimeIntegrator
TimeIntegrator
Definition Steppers.h:25
OpalData::getInputBasename
std::string getInputBasename()
get input file name without extension
Definition OpalData.cpp:674
OpalData::getInstance
static OpalData * getInstance()
Definition OpalData.cpp:196
OpalData::getGlobalGeometry
BoundaryGeometry * getGlobalGeometry()
Definition OpalData.cpp:461
OpalData::setOpenMode
void setOpenMode(OpenMode openMode)
Definition OpalData.cpp:349
TUNE_class::lombAnalysis
int lombAnalysis(double *x, double *y, int Ndat, int nhis)
Definition Ctunes.cpp:140
adder::operator()
void operator()(double x)
Definition ParallelCyclotronTracker.cpp:1501
ParallelCyclotronTracker::computePathLengthUpdate
void computePathLengthUpdate(std::vector< double > &dl, const double &dt)
Definition ParallelCyclotronTracker.cpp:232
ParallelCyclotronTracker::pluginElements_m
std::vector< PluginElement * > pluginElements_m
Definition ParallelCyclotronTracker.h:217
ParallelCyclotronTracker::dumpAzimuthAngles_m
void dumpAzimuthAngles_m(const double &t, const Vector_t &R, const Vector_t &P, const double &oldReferenceTheta, const double &temp_meanTheta)
Definition ParallelCyclotronTracker.cpp:3245
ParallelCyclotronTracker::IntegrationTimer_m
IpplTimings::TimerRef IntegrationTimer_m
Definition ParallelCyclotronTracker.h:329
ParallelCyclotronTracker::PluginElemTimer_m
IpplTimings::TimerRef PluginElemTimer_m
Definition ParallelCyclotronTracker.h:333
ParallelCyclotronTracker::visitDegrader
virtual void visitDegrader(const Degrader &)
Apply the algorithm to a degrader.
Definition ParallelCyclotronTracker.cpp:534
ParallelCyclotronTracker::visitRBend
virtual void visitRBend(const RBend &)
Apply the algorithm to a rectangular bend.
Definition ParallelCyclotronTracker.cpp:694
ParallelCyclotronTracker::xaxis
static Vector_t const xaxis
The positive axes unit vectors.
Definition ParallelCyclotronTracker.h:230
ParallelCyclotronTracker::visitProbe
virtual void visitProbe(const Probe &)
Apply the algorithm to a probe.
Definition ParallelCyclotronTracker.cpp:655
ParallelCyclotronTracker::visitVacuum
virtual void visitVacuum(const Vacuum &)
Apply the algorithm to a vacuum space.
Definition ParallelCyclotronTracker.cpp:989
ParallelCyclotronTracker::function_t
std::function< bool(const double &, const size_t &, Vector_t &, Vector_t &)> function_t
Definition ParallelCyclotronTracker.h:455
ParallelCyclotronTracker::computeExternalFields_m
bool computeExternalFields_m(const Vector_t &R, const Vector_t &P, const double &t, Vector_t &Efield, Vector_t &Bfield)
Calculate the field map at an arbitrary point.
Definition ParallelCyclotronTracker.cpp:3381
ParallelCyclotronTracker::azimuth_angle_m
std::vector< double > azimuth_angle_m
the different azimuthal angles for the outfTheta_m output files
Definition ParallelCyclotronTracker.h:309
ParallelCyclotronTracker::mbHandler_m
std::unique_ptr< MultiBunchHandler > mbHandler_m
Definition ParallelCyclotronTracker.h:259
ParallelCyclotronTracker::injectBunch
void injectBunch(bool &flagTransition)
Definition ParallelCyclotronTracker.cpp:3396
ParallelCyclotronTracker::rotateAroundZ
void rotateAroundZ(ParticleAttrib< Vector_t > &particleVectors, double const phi)
Definition ParallelCyclotronTracker.cpp:1835
ParallelCyclotronTracker::visitStripper
virtual void visitStripper(const Stripper &)
Apply the algorithm to a particle stripper.
Definition ParallelCyclotronTracker.cpp:938
ParallelCyclotronTracker::normalizeQuaternion
void normalizeQuaternion(Quaternion_t &quaternion)
Definition ParallelCyclotronTracker.cpp:1811
ParallelCyclotronTracker::visitCCollimator
virtual void visitCCollimator(const CCollimator &)
Apply the algorithm to a collimator.
Definition ParallelCyclotronTracker.cpp:477
ParallelCyclotronTracker::singleMode_m
void singleMode_m(double &t, const double dt, bool &finishedTurn, double &oldReferenceTheta)
Definition ParallelCyclotronTracker.cpp:3039
ParallelCyclotronTracker::update_m
void update_m(double &t, const double &dt, const bool &finishedTurn)
Update time and path length, write to output files.
Definition ParallelCyclotronTracker.cpp:2766
ParallelCyclotronTracker::outfTheta_m
std::vector< std::unique_ptr< std::ofstream > > outfTheta_m
output coordinates at different azimuthal angles and one after every turn
Definition ParallelCyclotronTracker.h:307
ParallelCyclotronTracker::rotateAroundX
void rotateAroundX(ParticleAttrib< Vector_t > &particleVectors, double const psi)
Definition ParallelCyclotronTracker.cpp:1869
ParallelCyclotronTracker::calculateAngle
double calculateAngle(double x, double y)
Definition ParallelCyclotronTracker.h:559
ParallelCyclotronTracker::isTurnDone
bool isTurnDone()
Check if turn done.
Definition ParallelCyclotronTracker.cpp:2762
ParallelCyclotronTracker::visitBeamline
virtual void visitBeamline(const Beamline &)
Apply the algorithm to a beam line.
Definition ParallelCyclotronTracker.cpp:1103
ParallelCyclotronTracker::checkGapCross
bool checkGapCross(Vector_t Rold, Vector_t Rnew, RFCavity *rfcavity, double &DistOld)
Definition ParallelCyclotronTracker.cpp:1459
ParallelCyclotronTracker::visitOffset
virtual void visitOffset(const Offset &)
Apply the algorithm to a offset (placement).
Definition ParallelCyclotronTracker.cpp:564
ParallelCyclotronTracker::DelParticleTimer_m
IpplTimings::TimerRef DelParticleTimer_m
Definition ParallelCyclotronTracker.h:334
ParallelCyclotronTracker::RFkick
bool RFkick(RFCavity *rfcavity, const double t, const double dt, const int Pindex)
Definition ParallelCyclotronTracker.cpp:1474
ParallelCyclotronTracker::visitScalingFFAMagnet
virtual void visitScalingFFAMagnet(const ScalingFFAMagnet &bend)
Apply the algorithm to a scaling FFA magnet.
Definition ParallelCyclotronTracker.cpp:865
ParallelCyclotronTracker::computeRadius
double computeRadius(const Vector_t &meanR) const
Definition ParallelCyclotronTracker.cpp:227
ParallelCyclotronTracker::updatePathLength
void updatePathLength(const double &dt)
Definition ParallelCyclotronTracker.cpp:3456
ParallelCyclotronTracker::openFiles
void openFiles(size_t numFiles, std::string fn)
@ open / close output coordinate files
Definition ParallelCyclotronTracker.cpp:275
ParallelCyclotronTracker::visitSBend3D
virtual void visitSBend3D(const SBend3D &)
Apply the algorithm to a sector bend with 3D field map.
Definition ParallelCyclotronTracker.cpp:856
ParallelCyclotronTracker::bega
double bega
The reference variables.
Definition ParallelCyclotronTracker.h:235
ParallelCyclotronTracker::dumpAngle
void dumpAngle(const double &theta, double &prevAzimuth, double &azimuth, const short &bunchNr=0)
Definition ParallelCyclotronTracker.cpp:203
ParallelCyclotronTracker::applyPluginElements
bool applyPluginElements(const double dt)
Definition ParallelCyclotronTracker.cpp:2040
ParallelCyclotronTracker::visitCorrector
virtual void visitCorrector(const Corrector &)
Apply the algorithm to a closed orbit corrector.
Definition ParallelCyclotronTracker.cpp:524
ParallelCyclotronTracker::initializeTracking_m
std::tuple< double, double, double > initializeTracking_m()
Definition ParallelCyclotronTracker.cpp:2809
ParallelCyclotronTracker::type_pair
std::pair< ElementType, element_pair > type_pair
Definition ParallelCyclotronTracker.h:61
ParallelCyclotronTracker::visitRing
virtual void visitRing(const Ring &ring)
Apply the algorithm to a ring.
Definition ParallelCyclotronTracker.cpp:798
ParallelCyclotronTracker::BinRepartTimer_m
IpplTimings::TimerRef BinRepartTimer_m
Definition ParallelCyclotronTracker.h:332
ParallelCyclotronTracker::globalToLocal
void globalToLocal(ParticleAttrib< Vector_t > &vectorArray, double phi, Vector_t const translationToGlobal=0)
Definition ParallelCyclotronTracker.cpp:1640
ParallelCyclotronTracker::getFieldsAtPoint
bool getFieldsAtPoint(const double &t, const size_t &Pindex, Vector_t &Efield, Vector_t &Bfield)
Definition ParallelCyclotronTracker.cpp:1429
ParallelCyclotronTracker::cavCrossDatas_m
std::vector< CavityCrossData > cavCrossDatas_m
Definition ParallelCyclotronTracker.h:218
ParallelCyclotronTracker::visitDrift
virtual void visitDrift(const Drift &)
Apply the algorithm to a drift space.
Definition ParallelCyclotronTracker.cpp:545
ParallelCyclotronTracker::itsStepper_mp
std::unique_ptr< Stepper< function_t > > itsStepper_mp
Definition ParallelCyclotronTracker.h:457
ParallelCyclotronTracker::setup_m
struct ParallelCyclotronTracker::settings setup_m
ParallelCyclotronTracker::finalizeTracking_m
void finalizeTracking_m(dvector_t &Ttime, dvector_t &Tdeltr, dvector_t &Tdeltz, ivector_t &TturnNumber)
Definition ParallelCyclotronTracker.cpp:2920
ParallelCyclotronTracker::visitVariableRFCavityFringeField
virtual void visitVariableRFCavityFringeField(const VariableRFCavityFringeField &cav)
Apply the algorithm to a variable RF cavity with Fringe Field.
Definition ParallelCyclotronTracker.cpp:1047
ParallelCyclotronTracker::gapCrossKick_m
void gapCrossKick_m(size_t i, double t, double dt, const Vector_t &Rold, const Vector_t &Pold)
Definition ParallelCyclotronTracker.cpp:3198
ParallelCyclotronTracker::getQuaternionTwoVectors
void getQuaternionTwoVectors(Vector_t u, Vector_t v, Quaternion_t &quaternion)
Definition ParallelCyclotronTracker.cpp:1903
ParallelCyclotronTracker::visitSBend
virtual void visitSBend(const SBend &)
Apply the algorithm to a sector bend.
Definition ParallelCyclotronTracker.cpp:851
ParallelCyclotronTracker::visitMarker
virtual void visitMarker(const Marker &)
Apply the algorithm to a marker.
Definition ParallelCyclotronTracker.cpp:577
ParallelCyclotronTracker::execute
virtual void execute()
Apply the algorithm to the top-level beamline.
Definition ParallelCyclotronTracker.cpp:1119
ParallelCyclotronTracker::dumpThetaEachTurn_m
void dumpThetaEachTurn_m(const double &t, const Vector_t &R, const Vector_t &P, const double &temp_meanTheta, bool &finishedTurn)
Definition ParallelCyclotronTracker.cpp:3267
ParallelCyclotronTracker::bunchMode_m
void bunchMode_m(double &t, const double dt, bool &finishedTurn)
Definition ParallelCyclotronTracker.cpp:3103
ParallelCyclotronTracker::myElements
std::list< Component * > myElements
Definition ParallelCyclotronTracker.h:215
ParallelCyclotronTracker::lossDs_m
std::unique_ptr< LossDataSink > lossDs_m
Definition ParallelCyclotronTracker.h:446
ParallelCyclotronTracker::visitMultipoleT
virtual void visitMultipoleT(const MultipoleT &)
Apply the algorithm to an arbitrary multipole.
Definition ParallelCyclotronTracker.cpp:609
ParallelCyclotronTracker::visitVariableRFCavity
virtual void visitVariableRFCavity(const VariableRFCavity &cav)
Apply the algorithm to a variabel RF cavity.
Definition ParallelCyclotronTracker.cpp:1032
ParallelCyclotronTracker::buildupFieldList
void buildupFieldList(double BcParameter[], ElementType elementType, Component *elptr)
Definition ParallelCyclotronTracker.cpp:1078
ParallelCyclotronTracker::rotateWithQuaternion
void rotateWithQuaternion(ParticleAttrib< Vector_t > &vectorArray, Quaternion_t const quaternion)
Definition ParallelCyclotronTracker.cpp:1797
ParallelCyclotronTracker::visitMultipole
virtual void visitMultipole(const Multipole &)
Apply the algorithm to a multipole.
Definition ParallelCyclotronTracker.cpp:599
ParallelCyclotronTracker::visitSolenoid
virtual void visitSolenoid(const Solenoid &)
Apply the algorithm to a solenoid.
Definition ParallelCyclotronTracker.cpp:924
ParallelCyclotronTracker::stepper_m
Steppers::TimeIntegrator stepper_m
Definition ParallelCyclotronTracker.h:468
ParallelCyclotronTracker::TransformTimer_m
IpplTimings::TimerRef TransformTimer_m
Definition ParallelCyclotronTracker.h:331
ParallelCyclotronTracker::calcMeanR
Vector_t calcMeanR(short bunchNr=-1) const
Definition ParallelCyclotronTracker.cpp:1595
ParallelCyclotronTracker::calculateAngle2
double calculateAngle2(double x, double y)
Definition ParallelCyclotronTracker.h:573
ParallelCyclotronTracker::visitCyclotron
virtual void visitCyclotron(const Cyclotron &cycl)
Apply the algorithm to a cyclotron.
Definition ParallelCyclotronTracker.cpp:313
ParallelCyclotronTracker::localToGlobal
void localToGlobal(ParticleAttrib< Vector_t > &vectorArray, double phi, Vector_t const translationToGlobal=0)
Definition ParallelCyclotronTracker.cpp:1661
ParallelCyclotronTracker::visitRFCavity
virtual void visitRFCavity(const RFCavity &)
Apply the algorithm to a RF cavity.
Definition ParallelCyclotronTracker.cpp:704
ParallelCyclotronTracker::getTunes
bool getTunes(dvector_t &t, dvector_t &r, dvector_t &z, int lastTurn, double &nur, double &nuz)
Definition ParallelCyclotronTracker.cpp:1516
ParallelCyclotronTracker::visitFlexibleCollimator
virtual void visitFlexibleCollimator(const FlexibleCollimator &)
Apply the algorithm to a flexible collimator.
Definition ParallelCyclotronTracker.cpp:555
ParallelCyclotronTracker::visitVerticalFFAMagnet
virtual void visitVerticalFFAMagnet(const VerticalFFAMagnet &bend)
Apply the algorithm to a vertical FFA magnet.
Definition ParallelCyclotronTracker.cpp:1061
ParallelCyclotronTracker::visitSeptum
virtual void visitSeptum(const Septum &)
Apply the algorithm to a septum.
Definition ParallelCyclotronTracker.cpp:882
ParallelCyclotronTracker::seoMode_m
void seoMode_m(double &t, const double dt, bool &finishedTurn, dvector_t &Ttime, dvector_t &Tdeltr, dvector_t &Tdeltz, ivector_t &TturnNumber)
Definition ParallelCyclotronTracker.cpp:2991
ParallelCyclotronTracker::maxSteps_m
int maxSteps_m
The maximal number of steps the system is integrated.
Definition ParallelCyclotronTracker.h:227
ParallelCyclotronTracker::visitOutputPlane
virtual void visitOutputPlane(const OutputPlane &)
Apply the algorithm to a outputplane.
Definition ParallelCyclotronTracker.cpp:625
ParallelCyclotronTracker::visitMonitor
virtual void visitMonitor(const Monitor &)
Apply the algorithm to a beam position monitor.
Definition ParallelCyclotronTracker.cpp:588
DumpEMFields::writeFields
static void writeFields(Component *field)
Definition DumpEMFields.cpp:209
DumpFields::writeFields
static void writeFields(Component *field)
Definition DumpFields.cpp:128
CCollimator::getWidth
double getWidth()
Definition CCollimator.h:100
CCollimator::getZStart
double getZStart()
Member variable access.
Definition CCollimator.h:90
CCollimator::getZEnd
double getZEnd()
Definition CCollimator.h:95
Component
Interface for a single beam element.
Definition Component.h:50
Corrector
Interface for general corrector.
Definition Corrector.h:35
Cyclotron::getCyclHarm
virtual double getCyclHarm() const
Definition Cyclotron.cpp:300
Drift
Interface for drift space.
Definition Drift.h:33
ElementBase::getName
virtual const std::string & getName() const
Get element name.
Definition ElementBase.cpp:162
ElementBase::getElementLength
virtual double getElementLength() const
Get design length.
Definition ElementBase.h:415
ElementBase::clone
virtual ElementBase * clone() const =0
Return clone.
ElementBase::hasAttribute
virtual bool hasAttribute(const std::string &aKey) const
Test for existence of an attribute.
Definition ElementBase.cpp:200
ElementBase::getTypeString
std::string getTypeString() const
Definition ElementBase.h:579
Marker
Interface for a marker.
Definition Marker.h:32
Multipole
Interface for general multipole.
Definition Multipole.h:47
MultipoleT::clone
ElementBase * clone() const override
Definition MultipoleT.cpp:69
OutputPlane::getCentre
Vector_t getCentre() const
Definition OutputPlane.h:280
OutputPlane::clone
ElementBase * clone() const override
Definition OutputPlane.cpp:58
OutputPlane::getHorizontalExtent
double getHorizontalExtent() const
Definition OutputPlane.h:304
OutputPlane::setGlobalFieldMap
void setGlobalFieldMap(Component *field)
Definition OutputPlane.h:267
OutputPlane::getNormal
Vector_t getNormal() const
Definition OutputPlane.h:271
PluginElement::getYStart
double getYStart() const
Definition PluginElement.cpp:184
PluginElement::initialise
virtual void initialise(PartBunchBase< double, 3 > *bunch, double &startField, double &endField) override
Pure virtual implementation of Component.
Definition PluginElement.cpp:48
PluginElement::getXEnd
double getXEnd() const
Definition PluginElement.cpp:180
PluginElement::getXStart
double getXStart() const
Member variable access.
Definition PluginElement.cpp:176
PluginElement::getYEnd
double getYEnd() const
Definition PluginElement.cpp:188
Probe
Definition Probe.h:28
RBend
Definition RBend.h:58
RFCavity::getPhaseModelName
std::string getPhaseModelName()
Definition RFCavity.h:464
RFCavity::getRmax
virtual double getRmax() const
Definition RFCavity.cpp:299
RFCavity::getMomentaKick
void getMomentaKick(const double normalRadius, double momentum[], const double t, const double dtCorrt, const int PID, const double restMass, const int chargenumber)
used in OPAL-cycl
Definition RFCavity.cpp:373
RFCavity::getAzimuth
virtual double getAzimuth() const
Definition RFCavity.cpp:308
RFCavity::getAmplitudeModelName
std::string getAmplitudeModelName()
Definition RFCavity.h:449
RFCavity::getCosAzimuth
virtual double getCosAzimuth() const
Definition RFCavity.cpp:316
RFCavity::initialise
virtual void initialise(PartBunchBase< double, 3 > *bunch, double &startField, double &endField) override
Definition RFCavity.cpp:167
RFCavity::getSinAzimuth
virtual double getSinAzimuth() const
Definition RFCavity.cpp:312
RFCavity::getFieldMapFN
virtual std::string getFieldMapFN() const
Definition RFCavity.cpp:345
RFCavity::getCycFrequency
virtual double getCycFrequency() const
Definition RFCavity.cpp:359
RFCavity::getGapWidth
virtual double getGapWidth() const
Definition RFCavity.cpp:324
RFCavity::getFrequencyModelName
std::string getFrequencyModelName()
Definition RFCavity.h:479
RFCavity::getPerpenDistance
virtual double getPerpenDistance() const
Definition RFCavity.cpp:320
RFCavity::getCavityType
CavityType getCavityType() const
Definition RFCavity.h:414
RFCavity::getPhi0
virtual double getPhi0() const
Definition RFCavity.cpp:328
RFCavity::getCavityTypeString
std::string getCavityTypeString() const
Definition RFCavity.cpp:341
RFCavity::getRmin
virtual double getRmin() const
Definition RFCavity.cpp:290
Ring
Ring describes a ring type geometry for tracking.
Definition Ring.h:53
Ring::clone
virtual ElementBase * clone() const override
Definition Ring.h:145
SBend
Definition SBend.h:68
ScalingFFAMagnet::clone
ScalingFFAMagnet * clone() const override
Definition ScalingFFAMagnet.cpp:52
Septum::initialise
virtual void initialise(PartBunchBase< double, 3 > *bunch, double &startField, double &endField) override
Pure virtual implementation of Component.
Definition Septum.cpp:47
Septum::getWidth
double getWidth() const
Definition Septum.cpp:57
Solenoid
Interface for solenoids.
Definition Solenoid.h:36
Stripper::getOPMass
double getOPMass() const
Definition Stripper.cpp:86
Stripper::getOPCharge
double getOPCharge() const
Definition Stripper.cpp:82
Stripper::getStop
bool getStop() const
Definition Stripper.cpp:94
Vacuum::getPressure
double getPressure() const
Definition Vacuum.cpp:128
Vacuum::getPressureMapFN
std::string getPressureMapFN() const
Definition Vacuum.cpp:141
Vacuum::getPScale
double getPScale() const
Definition Vacuum.cpp:149
Vacuum::initialise
virtual void initialise(PartBunchBase< double, 3 > *bunch, double &startField, double &endField) override
Definition Vacuum.cpp:233
Vacuum::checkVacuum
virtual bool checkVacuum(PartBunchBase< double, 3 > *bunch, Cyclotron *cycl)
Definition Vacuum.cpp:201
Vacuum::getTemperature
double getTemperature() const
Definition Vacuum.cpp:162
Vacuum::getResidualGasName
std::string getResidualGasName()
Definition Vacuum.cpp:115
Vacuum::getStop
bool getStop() const
Definition Vacuum.cpp:175
AbstractTimeDependence::getTimeDependence
static std::shared_ptr< AbstractTimeDependence > getTimeDependence(std::string name)
Definition AbstractTimeDependence.cpp:36
AbstractTracker::itsReference
const PartData itsReference
The reference information.
Definition AbstractTracker.h:54
Tracker::Tracker
Tracker(const Beamline &, const PartData &, bool backBeam, bool backTrack)
Constructor.
Definition Tracker.cpp:82
Tracker::itsBunch_m
PartBunchBase< double, 3 > * itsBunch_m
The bunch of particles to be tracked.
Definition Tracker.h:151
Beamline
An abstract sequence of beam line components.
Definition Beamline.h:34
TBeamline::iterate
virtual void iterate(BeamlineVisitor &, bool r2l) const
Apply visitor to all elements of the line.
Definition TBeamline.h:205
BorisPusher::kick
void kick(const Vector_t &R, Vector_t &P, const Vector_t &Ef, const Vector_t &Bf, const double &dt) const
Definition BorisPusher.h:65
LF2
Leap-Frog 2nd order.
Definition LF2.h:26
RK4
4-th order Runnge-Kutta stepper
Definition RK4.h:27
OpalException
The base class for all OPAL exceptions.
Definition OpalException.h:28
Message::put
Message & put(const T &val)
Definition Message.h:406
Message::get
Message & get(const T &cval)
Definition Message.h:476
IpplInfo::getNodes
static int getNodes()
Definition IpplInfo.cpp:670
IpplInfo::Comm
static Communicate * Comm
Definition IpplInfo.h:84
IpplTimings::getTimer
static TimerRef getTimer(const char *nm)
Definition IpplTimings.h:182
IpplTimings::stopTimer
static void stopTimer(TimerRef t)
Definition IpplTimings.h:192
IpplTimings::startTimer
static void startTimer(TimerRef t)
Definition IpplTimings.h:187
Vector_t
Vektor< double, 3 > Vector_t
Definition src/Classic/Algorithms/Vektor.h:6