// ------------------------------------------------------------------------ // $RCSfile: ParallelCyclotronTracker.cpp,v $ // ------------------------------------------------------------------------ // $Revision: 1.1 $initialLocalNum_m // ------------------------------------------------------------------------ // Copyright: see Copyright.readme // ------------------------------------------------------------------------ // // Class: ParallelCyclotronTracker // The class for tracking particles with 3D space charge in Cyclotrons and FFAs // // ------------------------------------------------------------------------ // // $Date: 2007/10/17 04:00:08 $ // $Author: adelmann, yang, winklehner $ // // ------------------------------------------------------------------------ #include "Algorithms/ParallelCyclotronTracker.h" #include #include #include #include #include #include "AbstractObjects/Element.h" #include "AbstractObjects/OpalData.h" #include "AbsBeamline/BeamStripping.h" #include "AbsBeamline/CCollimator.h" #include "AbsBeamline/Corrector.h" #include "AbsBeamline/Cyclotron.h" #include "AbsBeamline/CyclotronValley.h" #include "AbsBeamline/Degrader.h" #include "AbsBeamline/Diagnostic.h" #include "AbsBeamline/Drift.h" #include "AbsBeamline/ElementBase.h" #include "AbsBeamline/Lambertson.h" #include "AbsBeamline/Offset.h" #include "AbsBeamline/Marker.h" #include "AbsBeamline/Monitor.h" #include "AbsBeamline/Multipole.h" #include "AbsBeamline/MultipoleT.h" #include "AbsBeamline/MultipoleTBase.h" #include "AbsBeamline/MultipoleTStraight.h" #include "AbsBeamline/MultipoleTCurvedConstRadius.h" #include "AbsBeamline/MultipoleTCurvedVarRadius.h" #include "AbsBeamline/PluginElement.h" #include "AbsBeamline/Probe.h" #include "AbsBeamline/RBend.h" #include "AbsBeamline/Ring.h" #include "AbsBeamline/RFCavity.h" #include "AbsBeamline/RFQuadrupole.h" #include "AbsBeamline/SBend.h" #include "AbsBeamline/SBend3D.h" #include "AbsBeamline/ScalingFFAMagnet.h" #include "AbsBeamline/Separator.h" #include "AbsBeamline/Septum.h" #include "AbsBeamline/Solenoid.h" #include "AbsBeamline/Stripper.h" #include "AbsBeamline/VariableRFCavity.h" #include "AbsBeamline/VariableRFCavityFringeField.h" #include "Algorithms/Ctunes.h" #include "Algorithms/PolynomialTimeDependence.h" #include "Beamlines/Beamline.h" #include "Beamlines/FlaggedBeamline.h" #include "Elements/OpalBeamline.h" #include "Physics/Physics.h" #include "Utilities/OpalException.h" #include "Utilities/Options.h" #include "BasicActions/DumpFields.h" #include "BasicActions/DumpEMFields.h" #include "Structure/BoundaryGeometry.h" #include "Structure/DataSink.h" #include "Structure/LossDataSink.h" using Physics::pi; using Physics::q_e; constexpr double c_mmtns = Physics::c * 1.0e-6; // m/s --> mm/ns Vector_t const ParallelCyclotronTracker::xaxis = Vector_t(1.0, 0.0, 0.0); Vector_t const ParallelCyclotronTracker::yaxis = Vector_t(0.0, 1.0, 0.0); Vector_t const ParallelCyclotronTracker::zaxis = Vector_t(0.0, 0.0, 1.0); extern Inform *gmsg; /** * Constructor ParallelCyclotronTracker * * @param beamline * @param bunch * @param ds * @param reference * @param revBeam * @param revTrack * @param maxSTEPS * @param timeIntegrator */ ParallelCyclotronTracker::ParallelCyclotronTracker(const Beamline &beamline, PartBunchBase *bunch, DataSink &ds, const PartData &reference, bool revBeam, bool revTrack, int maxSTEPS, int timeIntegrator, const int& numBunch, const double& mbEta, const double& mbPara, const std::string& mbMode, const std::string& mbBinning) : Tracker(beamline, bunch, reference, revBeam, revTrack) , bgf_m(nullptr) , maxSteps_m(maxSTEPS) , lastDumpedStep_m(0) , myNode_m(Ippl::myNode()) , initialLocalNum_m(bunch->getLocalNum()) , initialTotalNum_m(bunch->getTotalNum()) , opalRing_m(nullptr) , itsStepper_mp(nullptr) { itsBeamline = dynamic_cast(beamline.clone()); itsDataSink = &ds; if ( numBunch > 1 ) { mbHandler_m = std::unique_ptr( new MultiBunchHandler(bunch, numBunch, mbEta, mbPara, mbMode, mbBinning) ); } IntegrationTimer_m = IpplTimings::getTimer("Integration"); TransformTimer_m = IpplTimings::getTimer("Frametransform"); DumpTimer_m = IpplTimings::getTimer("Dump"); BinRepartTimer_m = IpplTimings::getTimer("Binaryrepart"); PluginElemTimer_m = IpplTimings::getTimer("PluginElements"); DelParticleTimer_m = IpplTimings::getTimer("DeleteParticles"); // FIXME Change track command if ( initialTotalNum_m == 1 ) { mode_m = MODE::SINGLE; } else if ( initialTotalNum_m == 2 ) { mode_m = MODE::SEO; } else if ( initialTotalNum_m > 2 ) { mode_m = MODE::BUNCH; } else mode_m = MODE::UNDEFINED; if ( timeIntegrator == 0 ) { stepper_m = stepper::INTEGRATOR::RK4; } else if ( timeIntegrator == 1) { stepper_m = stepper::INTEGRATOR::LF2; } else if ( timeIntegrator == 2) { stepper_m = stepper::INTEGRATOR::MTS; } else stepper_m = stepper::INTEGRATOR::UNDEFINED; } /** * Destructor ParallelCyclotronTracker * */ ParallelCyclotronTracker::~ParallelCyclotronTracker() { if(bgf_m) lossDs_m->save(); for(Component* component : myElements) { delete(component); } for(auto fd : FieldDimensions) { delete(fd); } delete itsBeamline; // delete opalRing_m; } void ParallelCyclotronTracker::bgf_main_collision_test() { if(!bgf_m) return; Inform msg("bgf_main_collision_test "); /** *Here we check if a particle is outside the domain, flag it for deletion */ Vector_t intecoords = 0.0; // This has to match the dT in the rk4 pusher double dtime = itsBunch_m->getdT() * getHarmonicNumber(); int triId = 0; for(size_t i = 0; i < itsBunch_m->getLocalNum(); i++) { int res = bgf_m->partInside(itsBunch_m->R[i], itsBunch_m->P[i], dtime, itsBunch_m->PType[i], itsBunch_m->Q[i], intecoords, triId); //int res = bgf_m->partInside(itsBunch_m->R[i]*1.0e-3, itsBunch_m->P[i], dtime, itsBunch_m->PType[i], itsBunch_m->Q[i], intecoords, triId); if(res >= 0) { lossDs_m->addParticle(itsBunch_m->R[i], itsBunch_m->P[i], itsBunch_m->ID[i], itsBunch_m->getT()*1e9, turnnumber_m, itsBunch_m->bunchNum[i]); itsBunch_m->Bin[i] = -1; Inform gmsgALL("OPAL", INFORM_ALL_NODES); gmsgALL << level4 << "* Particle " << itsBunch_m->ID[i] << " lost on boundary geometry" << endl; } } } // only used for dumping into stat file void ParallelCyclotronTracker::dumpAngle(const double& theta, double& prevAzimuth, double& azimuth, const short& bunchNr) { if ( prevAzimuth < 0.0 ) { // only at first occurrence double plus = 0.0; if ( OpalData::getInstance()->inRestartRun() ) { plus = 360.0 * (turnnumber_m - bunchNr); } azimuth = theta + plus; } else { double dtheta = theta - prevAzimuth; if ( dtheta < 0 ) { dtheta += 360.0; } if ( dtheta > 180 ) { // rotating clockwise, reduce angle dtheta -= 360; } azimuth += dtheta; } prevAzimuth = theta; } double ParallelCyclotronTracker::computeRadius(const Vector_t& meanR) const { // New OPAL 2.0: m --> mm return 1000.0 * std::sqrt(meanR(0) * meanR(0) + meanR(1) * meanR(1)); } void ParallelCyclotronTracker::computePathLengthUpdate(std::vector& dl, const double& dt) { // the last element in dotP is the dot-product over all particles std::vector dotP(dl.size()); if ( Options::psDumpFrame == Options::BUNCH_MEAN || isMultiBunch()) { for(unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) { dotP[itsBunch_m->bunchNum[i]] += dot(itsBunch_m->P[i], itsBunch_m->P[i]); } allreduce(dotP.data(), dotP.size(), std::plus()); // dot-product over all particles double sum = std::accumulate(dotP.begin(), dotP.end() - 1, 0); dotP.back() = sum / double(itsBunch_m->getTotalNum()); // bunch specific --> multi-bunches only for (short b = 0; b < (short)dotP.size() - 1; ++b) { dotP[b] /= double(itsBunch_m->getTotalNumPerBunch(b)); } } else if ( itsBunch_m->getLocalNum() == 0 ) { // here we are in Options::GLOBAL mode dotP[0] = 0.0; } else { // here we are in Options::GLOBAL mode dotP[0] = dot(itsBunch_m->P[0], itsBunch_m->P[0]); } for (size_t i = 0; i < dotP.size(); ++i) { double const gamma = std::sqrt(1.0 + dotP[i]); double const beta = std::sqrt(dotP[i]) / gamma; dl[i] = c_mmtns * dt * 1.0e-3 * beta; // unit: m } } /** * * * @param fn Base file name */ void ParallelCyclotronTracker::openFiles(size_t numFiles, std::string SfileName) { for (unsigned int i=0; iprecision(8); outfTheta_m.back()->setf(std::ios::scientific, std::ios::floatfield); *outfTheta_m.back() << "# r [mm] beta_r*gamma " << "theta [deg] beta_theta*gamma " << "z [mm] beta_z*gamma" << std::endl; } } /** * Close all files related to * special output in the Cyclotron * mode. */ void ParallelCyclotronTracker::closeFiles() { for (auto & file : outfTheta_m) { file->close(); } } /** * * @param ring */ void ParallelCyclotronTracker::visitRing(const Ring &ring) { *gmsg << "* ----------------------------- Adding Ring ------------------------------ *" << endl; delete opalRing_m; opalRing_m = dynamic_cast(ring.clone()); myElements.push_back(opalRing_m); opalRing_m->initialise(itsBunch_m); referenceR = opalRing_m->getBeamRInit(); referencePr = opalRing_m->getBeamPRInit(); referenceTheta = opalRing_m->getBeamPhiInit(); if(referenceTheta <= -180.0 || referenceTheta > 180.0) { throw OpalException("Error in ParallelCyclotronTracker::visitRing", "PHIINIT is out of [-180, 180)!"); } referenceZ = 0.0; referencePz = 0.0; referencePtot = itsReference.getGamma() * itsReference.getBeta(); referencePt = sqrt(referencePtot * referencePtot - referencePr * referencePr); if(referencePtot < 0.0) referencePt *= -1.0; sinRefTheta_m = sin(referenceTheta * Physics::deg2rad); cosRefTheta_m = cos(referenceTheta * Physics::deg2rad); double BcParameter[8] = {}; // zero initialise array buildupFieldList(BcParameter, ElementBase::RING, opalRing_m); // Finally print some diagnostic *gmsg << "* Initial beam radius = " << referenceR << " [mm] " << endl; *gmsg << "* Initial gamma = " << itsReference.getGamma() << endl; *gmsg << "* Initial beta = " << itsReference.getBeta() << endl; *gmsg << "* Total reference momentum = " << referencePtot * 1000.0 << " [MCU]" << endl; *gmsg << "* Reference azimuthal momentum = " << referencePt * 1000.0 << " [MCU]" << endl; *gmsg << "* Reference radial momentum = " << referencePr * 1000.0 << " [MCU]" << endl; *gmsg << "* " << opalRing_m->getSymmetry() << " fold field symmetry " << endl; *gmsg << "* Harmonic number h= " << opalRing_m->getHarmonicNumber() << " " << endl; } /** * * * @param cycl */ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) { *gmsg << "* -------------------------- Adding Cyclotron ---------------------------- *" << endl; cycl_m = dynamic_cast(cycl.clone()); myElements.push_back(cycl_m); // Is this a Spiral Inflector Simulation? If yes, we'll give the user some // useful information spiral_flag = cycl_m->getSpiralFlag(); if(spiral_flag) { *gmsg << endl << "* This is a Spiral Inflector Simulation! This means the following:" << endl; *gmsg << "* 1.) It is up to the user to provide appropriate geometry, electric and magnetic fields!" << endl; *gmsg << "* (Use BANDRF type cyclotron and use RFMAPFN to load both magnetic" << endl; *gmsg << "* and electric fields, setting SUPERPOSE to an array of TRUE values.)" << endl; *gmsg << "* 2.) For high currents it is strongly recommended to use the SAAMG fieldsolver," << endl; *gmsg << "* FFT does not give the correct results (boundary conditions are missing)." << endl; *gmsg << "* 3.) The whole geometry will be meshed and used for the fieldsolver." << endl; *gmsg << "* There will be no transformations of the bunch into a local frame und consequently," << endl; *gmsg << "* the problem will be treated non-relativistically!" << endl; *gmsg << "* (This is not an issue for spiral inflectors as they are typically < 100 keV/amu.)" << endl; *gmsg << endl << "* Note: For now, multi-bunch mode (MBM) needs to be de-activated for spiral inflector" << endl; *gmsg << "* and space charge needs to be solved every time-step. numBunch_m and scSolveFreq are reset." << endl; if (isMultiBunch()) { mbHandler_m = nullptr; } } // Fresh run (no restart): if(!OpalData::getInstance()->inRestartRun()) { // Get reference values from cyclotron element // For now, these are still stored in mm. should be the only ones. -DW referenceR = cycl_m->getRinit(); referenceTheta = cycl_m->getPHIinit(); referenceZ = cycl_m->getZinit(); referencePr = cycl_m->getPRinit(); referencePz = cycl_m->getPZinit(); if(referenceTheta <= -180.0 || referenceTheta > 180.0) { throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron", "PHIINIT is out of [-180, 180)!"); } referencePtot = itsReference.getGamma() * itsReference.getBeta(); // Calculate reference azimuthal (tangential) momentum from total-, z- and radial momentum: float insqrt = referencePtot * referencePtot - \ referencePr * referencePr - referencePz * referencePz; if(insqrt < 0) { if(insqrt > -1.0e-10) { referencePt = 0.0; } else { throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron", "Pt imaginary!"); } } else { referencePt = sqrt(insqrt); } if(referencePtot < 0.0) referencePt *= -1.0; // End calculate referencePt // Restart a run: } else { // If the user wants to save the restarted run in local frame, // make sure the previous h5 file was local too if (Options::psDumpFrame != Options::GLOBAL) { if (!previousH5Local) { throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron", "You are trying a local restart from a global h5 file!"); } // Else, if the user wants to save the restarted run in global frame, // make sure the previous h5 file was global too } else { if (previousH5Local) { throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron", "You are trying a global restart from a local h5 file!"); } } // Adjust some of the reference variables from the h5 file referencePhi *= Physics::deg2rad; referencePsi *= Physics::deg2rad; referencePtot = bega; if(referenceTheta <= -180.0 || referenceTheta > 180.0) { throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron", "PHIINIT is out of [-180, 180)!"); } } sinRefTheta_m = sin(referenceTheta * Physics::deg2rad); cosRefTheta_m = cos(referenceTheta * Physics::deg2rad); *gmsg << endl; *gmsg << "* Bunch global starting position:" << endl; *gmsg << "* RINIT = " << referenceR << " [mm]" << endl; *gmsg << "* PHIINIT = " << referenceTheta << " [deg]" << endl; *gmsg << "* ZINIT = " << referenceZ << " [mm]" << endl; *gmsg << endl; *gmsg << "* Bunch global starting momenta:" << endl; *gmsg << "* Initial gamma = " << itsReference.getGamma() << endl; *gmsg << "* Initial beta = " << itsReference.getBeta() << endl; *gmsg << "* Reference total momentum (beta * gamma) = " << referencePtot * 1000.0 << " [MCU]" << endl; *gmsg << "* Reference azimuthal momentum (Pt) = " << referencePt * 1000.0 << " [MCU]" << endl; *gmsg << "* Reference radial momentum (Pr) = " << referencePr * 1000.0 << " [MCU]" << endl; *gmsg << "* Reference axial momentum (Pz) = " << referencePz * 1000.0 << " [MCU]" << endl; *gmsg << endl; double sym = cycl_m->getSymmetry(); *gmsg << "* " << sym << "-fold field symmetry " << endl; // ckr: this just returned the default value as defined in Component.h // double rff = cycl_m->getRfFrequ(); // *gmsg << "* Rf frequency= " << rff << " [MHz]" << endl; std::string fmfn = cycl_m->getFieldMapFN(); *gmsg << "* Field map file name = " << fmfn << " " << endl; std::string type = cycl_m->getCyclotronType(); *gmsg << "* Type of cyclotron = " << type << " " << endl; double rmin = cycl_m->getMinR(); double rmax = cycl_m->getMaxR(); *gmsg << "* Radial aperture = " << rmin << " ... " << rmax<<" [m] "<< endl; double zmin = cycl_m->getMinZ(); double zmax = cycl_m->getMaxZ(); *gmsg << "* Vertical aperture = " << zmin << " ... " << zmax<<" [m]"<< endl; double h = cycl_m->getCyclHarm(); *gmsg << "* Number of trimcoils = " << cycl_m->getNumberOfTrimcoils() << endl; *gmsg << "* Harmonic number h = " << h << " " << endl; /** * To ease the initialise() function, set a integral parameter fieldflag internally. * Its value is by the option "TYPE" of the element "CYCLOTRON" * fieldflag = 1, readin PSI format measured field file (default) * fieldflag = 2, readin carbon cyclotron field file created by Jianjun Yang, TYPE=CARBONCYCL * fieldflag = 3, readin ANSYS format file for CYCIAE-100 created by Jianjun Yang, TYPE=CYCIAE * fieldflag = 4, readin AVFEQ format file for Riken cyclotrons * fieldflag = 5, readin FFA format file for MSU/FNAL FFA * fieldflag = 6, readin both median plane B field map and 3D E field map of RF cavity for compact cyclotron * fieldflag = 7, read in fields for Daniel's synchrocyclotron simulations */ int fieldflag = cycl_m->getFieldFlag(type); // Read in cyclotron field maps (midplane + 3D fields if desired). cycl_m->initialise(itsBunch_m, fieldflag, cycl_m->getBScale()); double BcParameter[8] = {}; BcParameter[0] = 0.001 * cycl_m->getRmin(); BcParameter[1] = 0.001 * cycl_m->getRmax(); // Store inner radius and outer radius of cyclotron field map in the list buildupFieldList(BcParameter, ElementBase::CYCLOTRON, cycl_m); } /** * Not implemented and most probable never used * */ void ParallelCyclotronTracker::visitBeamBeam(const BeamBeam &) { *gmsg << "In BeamBeam tracker is missing " << endl; } void ParallelCyclotronTracker::visitBeamStripping(const BeamStripping &bstp) { *gmsg << "* ------------------------------ Beam Stripping ------------------------------" << endl; BeamStripping* elptr = dynamic_cast(bstp.clone()); myElements.push_back(elptr); double BcParameter[8] = {}; if(elptr->getPressureMapFN() == "") { double pressure = elptr->getPressure(); *gmsg << "* Pressure = " << pressure << " [mbar]" << endl; BcParameter[0] = pressure; } double temperature = elptr->getTemperature(); *gmsg << "* Temperature = " << temperature << " [K]" << endl; std::string gas = elptr->getResidualGas(); *gmsg << "* Residual gas = " << gas << endl; bool stop = elptr->getStop(); *gmsg << std::boolalpha << "* Particles stripped will be deleted after interaction -> " << stop << endl; elptr->initialise(itsBunch_m, elptr->getPScale()); BcParameter[1] = temperature; BcParameter[2] = stop; buildupFieldList(BcParameter, ElementBase::BEAMSTRIPPING, elptr); } /** * * * @param coll */ void ParallelCyclotronTracker::visitCCollimator(const CCollimator &coll) { *gmsg << "* ------------------------------ Collimator ------------------------------" << endl; CCollimator* elptr = dynamic_cast(coll.clone()); myElements.push_back(elptr); double xstart = elptr->getXStart(); *gmsg << "* Xstart = " << xstart << " [mm]" << endl; double xend = elptr->getXEnd(); *gmsg << "* Xend = " << xend << " [mm]" << endl; double ystart = elptr->getYStart(); *gmsg << "* Ystart = " << ystart << " [mm]" << endl; double yend = elptr->getYEnd(); *gmsg << "* Yend = " << yend << " [mm]" << endl; double zstart = elptr->getZStart(); *gmsg << "* Zstart = " << zstart << " [mm]" << endl; double zend = elptr->getZEnd(); *gmsg << "* Zend = " << zend << " [mm]" << endl; double width = elptr->getWidth(); *gmsg << "* Width = " << width << " [mm]" << endl; elptr->initialise(itsBunch_m); double BcParameter[8] = {}; BcParameter[0] = 0.001 * xstart ; BcParameter[1] = 0.001 * xend; BcParameter[2] = 0.001 * ystart ; BcParameter[3] = 0.001 * yend; BcParameter[4] = 0.001 * width ; buildupFieldList(BcParameter, ElementBase::CCOLLIMATOR, elptr); } /** * * * @param corr */ void ParallelCyclotronTracker::visitCorrector(const Corrector &corr) { *gmsg << "In Corrector; L= " << corr.getElementLength() << endl; myElements.push_back(dynamic_cast(corr.clone())); } /** * * * @param degrader */ void ParallelCyclotronTracker::visitDegrader(const Degrader °) { *gmsg << "In Degrader; L= " << deg.getElementLength() << endl; myElements.push_back(dynamic_cast(deg.clone())); } /** * * * @param diag */ void ParallelCyclotronTracker::visitDiagnostic(const Diagnostic &diag) { *gmsg << "In Diagnostic; L= " << diag.getElementLength() << endl; myElements.push_back(dynamic_cast(diag.clone())); } /** * * * @param drift */ void ParallelCyclotronTracker::visitDrift(const Drift &drift) { *gmsg << "In drift L= " << drift.getElementLength() << endl; myElements.push_back(dynamic_cast(drift.clone())); } /** * * * @param */ void ParallelCyclotronTracker::visitFlexibleCollimator(const FlexibleCollimator &) { } /** * * * @param lamb */ void ParallelCyclotronTracker::visitLambertson(const Lambertson &lamb) { *gmsg << "In Lambertson; L= " << lamb.getElementLength() << endl; myElements.push_back(dynamic_cast(lamb.clone())); } void ParallelCyclotronTracker::visitOffset(const Offset & off) { if (opalRing_m == NULL) throw OpalException( "ParallelCylcotronTracker::visitOffset", "Attempt to place an offset when Ring not defined"); Offset* offNonConst = const_cast(&off); offNonConst->updateGeometry(opalRing_m->getNextPosition(), opalRing_m->getNextNormal()); opalRing_m->appendElement(off); } /** * * * @param marker */ void ParallelCyclotronTracker::visitMarker(const Marker &marker) { // *gmsg << "In Marker; L= " << marker.getElementLength() << endl; myElements.push_back(dynamic_cast(marker.clone())); // Do nothing. } /** * * * @param corr */ void ParallelCyclotronTracker::visitMonitor(const Monitor &corr) { // *gmsg << "In Monitor; L= " << corr.getElementLength() << endl; myElements.push_back(dynamic_cast(corr.clone())); // applyDrift(flip_s * corr.getElementLength()); } /** * * * @param mult */ void ParallelCyclotronTracker::visitMultipole(const Multipole &mult) { *gmsg << "In Multipole; L= " << mult.getElementLength() << " however the element is missing " << endl; myElements.push_back(dynamic_cast(mult.clone())); } /** * * * @param multT */ void ParallelCyclotronTracker::visitMultipoleT(const MultipoleT &multT) { *gmsg << "Adding MultipoleT" << endl; if (opalRing_m != NULL) { opalRing_m->appendElement(multT); } else { throw OpalException("ParallelCyclotronTracker::visitMultipoleT", "Need to define a RINGDEFINITION to use MultipoleT element"); } myElements.push_back(dynamic_cast(multT.clone())); } /** * * * @param multTstraight */ void ParallelCyclotronTracker::visitMultipoleTStraight(const MultipoleTStraight &multTstraight) { *gmsg << "Adding MultipoleTStraight" << endl; if (opalRing_m != NULL) { opalRing_m->appendElement(multTstraight); } else { throw OpalException("ParallelCyclotronTracker::visitMultipoleTStraight", "Need to define a RINGDEFINITION to use MultipoleTStraight element"); } myElements.push_back(dynamic_cast(multTstraight.clone())); } /** * * * @param multTccurv */ void ParallelCyclotronTracker::visitMultipoleTCurvedConstRadius(const MultipoleTCurvedConstRadius &multTccurv) { *gmsg << "Adding MultipoleTCurvedConstRadius" << endl; if (opalRing_m != NULL) { opalRing_m->appendElement(multTccurv); } else { throw OpalException("ParallelCyclotronTracker::visitMultipoleTCurvedConstRadius", "Need to define a RINGDEFINITION to use MultipoleTCurvedConstRadius element"); } myElements.push_back(dynamic_cast(multTccurv.clone())); } /** * * * @param multTvcurv */ void ParallelCyclotronTracker::visitMultipoleTCurvedVarRadius(const MultipoleTCurvedVarRadius &multTvcurv) { *gmsg << "Adding MultipoleTCurvedVarRadius" << endl; if (opalRing_m != NULL) { opalRing_m->appendElement(multTvcurv); } else { throw OpalException("ParallelCyclotronTracker::visitMultipoleTCurvedVarRadius", "Need to define a RINGDEFINITION to use MultipoleTCurvedVarRadius element"); } myElements.push_back(dynamic_cast(multTvcurv.clone())); } /** * * * @param prob */ void ParallelCyclotronTracker::visitProbe(const Probe &prob) { *gmsg << "* ------------------------------ Probe ------------------------------" << endl; Probe *elptr = dynamic_cast(prob.clone()); myElements.push_back(elptr); *gmsg << "* Name = " << elptr->getName() << endl; double xstart = elptr->getXStart(); *gmsg << "* XStart = " << xstart << " [mm]" << endl; double xend = elptr->getXEnd(); *gmsg << "* XEnd = " << xend << " [mm]" << endl; double ystart = elptr->getYStart(); *gmsg << "* YStart = " << ystart << " [mm]" << endl; double yend = elptr->getYEnd(); *gmsg << "* YEnd = " << yend << " [mm]" << endl; // initialise, do nothing elptr->initialise(itsBunch_m); double BcParameter[8] = {}; BcParameter[0] = 0.001 * xstart ; BcParameter[1] = 0.001 * xend; BcParameter[2] = 0.001 * ystart ; BcParameter[3] = 0.001 * yend; BcParameter[4] = 0.001 ; // width // store probe parameters in the list buildupFieldList(BcParameter, ElementBase::PROBE, elptr); } /** * * * @param bend */ void ParallelCyclotronTracker::visitRBend(const RBend &bend) { *gmsg << "In RBend; L= " << bend.getElementLength() << " however the element is missing " << endl; myElements.push_back(dynamic_cast(bend.clone())); } void ParallelCyclotronTracker::visitSBend3D(const SBend3D &bend) { *gmsg << "Adding SBend3D" << endl; if (opalRing_m != NULL) opalRing_m->appendElement(bend); else throw OpalException("ParallelCyclotronTracker::visitSBend3D", "Need to define a RINGDEFINITION to use SBend3D element"); } void ParallelCyclotronTracker::visitScalingFFAMagnet(const ScalingFFAMagnet &bend) { *gmsg << "Adding ScalingFFAMagnet" << endl; if (opalRing_m != NULL) { opalRing_m->appendElement(bend); } else { throw OpalException("ParallelCyclotronTracker::visitScalingFFAMagnet", "Need to define a RINGDEFINITION to use ScalingFFAMagnet element"); } } void ParallelCyclotronTracker::visitVariableRFCavity(const VariableRFCavity &cav) { *gmsg << "Adding Variable RF Cavity" << endl; if (opalRing_m != NULL) opalRing_m->appendElement(cav); else throw OpalException("ParallelCyclotronTracker::visitVariableRFCavity", "Need to define a RINGDEFINITION to use VariableRFCavity element"); } void ParallelCyclotronTracker::visitVariableRFCavityFringeField (const VariableRFCavityFringeField &cav) { *gmsg << "Adding Variable RF Cavity with Fringe Field" << endl; if (opalRing_m != NULL) opalRing_m->appendElement(cav); else throw OpalException("ParallelCyclotronTracker::visitVariableRFCavityFringeField", "Need to define a RINGDEFINITION to use VariableRFCavity element"); } /** * * * @param as */ void ParallelCyclotronTracker::visitRFCavity(const RFCavity &as) { *gmsg << "* ------------------------------ RFCavity ------------------------------" << endl; RFCavity *elptr = dynamic_cast(as.clone()); myElements.push_back(elptr); if((elptr->getComponentType() != "SINGLEGAP") && (elptr->getComponentType() != "DOUBLEGAP")) { *gmsg << (elptr->getComponentType()) << endl; throw OpalException("ParallelCyclotronTracker::visitRFCavity", "The ParallelCyclotronTracker can only play with cyclotron type RF system currently ..."); } double rmin = elptr->getRmin(); *gmsg << "* Minimal radius of cavity= " << rmin << " [mm]" << endl; double rmax = elptr->getRmax(); *gmsg << "* Maximal radius of cavity= " << rmax << " [mm]" << endl; double rff = elptr->getCycFrequency(); *gmsg << "* RF frequency (2*pi*f)= " << rff << " [rad/s]" << endl; std::string fmfn = elptr->getFieldMapFN(); *gmsg << "* RF Field map file name= " << fmfn << endl; double angle = elptr->getAzimuth(); *gmsg << "* Cavity azimuth position= " << angle << " [deg] " << endl; double gap = elptr->getGapWidth(); *gmsg << "* Cavity gap width= " << gap << " [mm] " << endl; double pdis = elptr->getPerpenDistance(); *gmsg << "* Cavity Shift distance= " << pdis << " [mm] " << endl; double phi0 = elptr->getPhi0(); *gmsg << "* Initial RF phase (t=0)= " << phi0 << " [deg] " << endl; /* Setup time dependence and in case of no timedependence use a polynom with a_0 = 1 and a_k = 0, k = 1,2,3. */ std::shared_ptr freq_atd = nullptr; std::shared_ptr ampl_atd = nullptr; std::shared_ptr phase_atd = nullptr; dvector_t unityVec; unityVec.push_back(1.); unityVec.push_back(0.); unityVec.push_back(0.); unityVec.push_back(0.); if (elptr->getFrequencyModelName() != "") { freq_atd = AbstractTimeDependence::getTimeDependence(elptr->getFrequencyModelName()); *gmsg << "* Variable frequency RF Model name " << elptr->getFrequencyModelName() << endl; } else freq_atd = std::shared_ptr(new PolynomialTimeDependence(unityVec)); if (elptr->getAmplitudeModelName() != "") { ampl_atd = AbstractTimeDependence::getTimeDependence(elptr->getAmplitudeModelName()); *gmsg << "* Variable amplitude RF Model name " << elptr->getAmplitudeModelName() << endl; } else ampl_atd = std::shared_ptr(new PolynomialTimeDependence(unityVec)); if (elptr->getPhaseModelName() != "") { phase_atd = AbstractTimeDependence::getTimeDependence(elptr->getPhaseModelName()); *gmsg << "* Variable phase RF Model name " << elptr->getPhaseModelName() << endl; } else phase_atd = std::shared_ptr(new PolynomialTimeDependence(unityVec)); // read cavity voltage profile data from file. elptr->initialise(itsBunch_m, freq_atd, ampl_atd, phase_atd); double BcParameter[8] = {}; BcParameter[0] = 0.001 * rmin; BcParameter[1] = 0.001 * rmax; BcParameter[2] = 0.001 * pdis; BcParameter[3] = angle; buildupFieldList(BcParameter, ElementBase::RFCAVITY, elptr); } /** * * * @param rfq */ void ParallelCyclotronTracker::visitRFQuadrupole(const RFQuadrupole &rfq) { *gmsg << "In RFQuadrupole; L = " << rfq.getElementLength() << " however the element is missing " << endl; myElements.push_back(dynamic_cast(rfq.clone())); } /** * * * @param bend */ void ParallelCyclotronTracker::visitSBend(const SBend &bend) { *gmsg << "In SBend; L = " << bend.getElementLength() << " however the element is missing " << endl; myElements.push_back(dynamic_cast(bend.clone())); } /** * * * @param sep */ void ParallelCyclotronTracker::visitSeparator(const Separator &sep) { *gmsg << "In Separator L= " << sep.getElementLength() << " however the element is missing " << endl; myElements.push_back(dynamic_cast(sep.clone())); } /** * * * @param sept */ void ParallelCyclotronTracker::visitSeptum(const Septum &sept) { *gmsg << endl << "* ------------------------------ Septum ------------------------------- *" << endl; Septum *elptr = dynamic_cast(sept.clone()); myElements.push_back(elptr); double xstart = elptr->getXStart(); *gmsg << "* XStart = " << xstart << " [mm]" << endl; double xend = elptr->getXEnd(); *gmsg << "* XEnd = " << xend << " [mm]" << endl; double ystart = elptr->getYStart(); *gmsg << "* YStart = " << ystart << " [mm]" << endl; double yend = elptr->getYEnd(); *gmsg << "* YEnd = " << yend << " [mm]" << endl; double width = elptr->getWidth(); *gmsg << "* Width = " << width << " [mm]" << endl; // initialise, do nothing elptr->initialise(itsBunch_m); double BcParameter[8] = {}; BcParameter[0] = 0.001 * xstart ; BcParameter[1] = 0.001 * xend; BcParameter[2] = 0.001 * ystart ; BcParameter[3] = 0.001 * yend; BcParameter[4] = 0.001 * width ; // store septum parameters in the list buildupFieldList(BcParameter, ElementBase::SEPTUM, elptr); } /** * * * @param solenoid */ void ParallelCyclotronTracker::visitSolenoid(const Solenoid &solenoid) { myElements.push_back(dynamic_cast(solenoid.clone())); Component *elptr = *(--myElements.end()); if(!elptr->hasAttribute("ELEMEDGE")) { *gmsg << "Solenoid: no position of the element given!" << endl; return; } } /** * * * @param pplate */ void ParallelCyclotronTracker::visitParallelPlate(const ParallelPlate &pplate) {//do nothing //*gmsg << "ParallelPlate: not in use in ParallelCyclotronTracker!" << endl; //buildupFieldList(startField, endField, elptr); } /** * * * @param cv */ void ParallelCyclotronTracker::visitCyclotronValley(const CyclotronValley &cv) { // Do nothing here. } /** * * * @param stripper */ void ParallelCyclotronTracker::visitStripper(const Stripper &stripper) { *gmsg << "* ------------------------------ Stripper ------------------------------" << endl; Stripper *elptr = dynamic_cast(stripper.clone()); myElements.push_back(elptr); *gmsg << "* Name = " << elptr->getName() << endl; double xstart = elptr->getXStart(); *gmsg << "* XStart = " << xstart << " [mm]" << endl; double xend = elptr->getXEnd(); *gmsg << "* XEnd = " << xend << " [mm]" << endl; double ystart = elptr->getYStart(); *gmsg << "* YStart = " << ystart << " [mm]" << endl; double yend = elptr->getYEnd(); *gmsg << "* YEnd = " << yend << " [mm]" << endl; double opcharge = elptr->getOPCharge(); *gmsg << "* Charge of outcoming particle = +e * " << opcharge << endl; double opmass = elptr->getOPMass(); *gmsg << "* Mass of the outcoming particle = " << opmass << " [GeV/c^2]" << endl; bool stop = elptr->getStop(); *gmsg << std::boolalpha << "* Particles stripped will be deleted after interaction -> " << stop << endl; elptr->initialise(itsBunch_m); double BcParameter[8] = {}; BcParameter[0] = 0.001 * xstart ; BcParameter[1] = 0.001 * xend; BcParameter[2] = 0.001 * ystart ; BcParameter[3] = 0.001 * yend; BcParameter[4] = 0.001; // width BcParameter[5] = opcharge; BcParameter[6] = opmass; buildupFieldList(BcParameter, ElementBase::STRIPPER, elptr); } /** * * * @param BcParameter * @param ElementType * @param elptr */ void ParallelCyclotronTracker::buildupFieldList(double BcParameter[], ElementBase::ElementType elementType, Component *elptr) { beamline_list::iterator sindex; type_pair *localpair = new type_pair(); localpair->first = elementType; for(int i = 0; i < 8; i++) *(((localpair->second).first) + i) = *(BcParameter + i); (localpair->second).second = elptr; // always put cyclotron as the first element in the list. if(elementType == ElementBase::RING || elementType == ElementBase::CYCLOTRON) { sindex = FieldDimensions.begin(); } else { sindex = FieldDimensions.end(); } FieldDimensions.insert(sindex, localpair); } /** * * * @param bl */ void ParallelCyclotronTracker::visitBeamline(const Beamline &bl) { const FlaggedBeamline* fbl = static_cast(&bl); fbl->iterate(*this, false);//*dynamic_cast(this), false); } void ParallelCyclotronTracker::checkNumPart(std::string s) { int nlp = itsBunch_m->getLocalNum(); int minnlp = 0; int maxnlp = 111111; reduce(nlp, minnlp, OpMinAssign()); reduce(nlp, maxnlp, OpMaxAssign()); *gmsg << s << " min local particle number " << minnlp << " max local particle number: " << maxnlp << endl; } /** * * */ void ParallelCyclotronTracker::execute() { /* Initialize common variables and structures for the integrators */ if (OpalData::getInstance()->inRestartRun()) { OpalData::getInstance()->setOpenMode(OpalData::OPENMODE::APPEND); } step_m = 0; restartStep0_m = 0; turnnumber_m = 1; azimuth_m = -1.0; prevAzimuth_m = -1.0; // Record how many bunches have already been injected. ONLY FOR MPM if (isMultiBunch()) mbHandler_m->setNumBunch(itsBunch_m->getNumBunch()); itsBeamline->accept(*this); if (opalRing_m != NULL) opalRing_m->lockRing(); // Display the selected elements *gmsg << "* ---------------------------------------------------" << endl; *gmsg << "* The selected Beam line elements are :" << endl; for(auto fd : FieldDimensions) { ElementBase::ElementType type = fd->first; *gmsg << "* -> " << ElementBase::getTypeString(type) << endl; if(type == ElementBase::RFCAVITY) { RFCavity* cav = static_cast((fd->second).second); CavityCrossData ccd = {cav, cav->getSinAzimuth(), cav->getCosAzimuth(), cav->getPerpenDistance() * 0.001}; cavCrossDatas_m.push_back(ccd); } else if( type == ElementBase::CCOLLIMATOR || type == ElementBase::PROBE || type == ElementBase::SEPTUM || type == ElementBase::STRIPPER) { PluginElement* element = static_cast((fd->second).second); pluginElements_m.push_back(element); } } *gmsg << "* ---------------------------------------------------" << endl; // Get BoundaryGeometry that is already initialized bgf_m = OpalData::getInstance()->getGlobalGeometry(); if (bgf_m) lossDs_m = std::unique_ptr(new LossDataSink("GEOM",!Options::asciidump)); // External field arrays for dumping for(int k = 0; k < 2; k++) FDext_m[k] = Vector_t(0.0, 0.0, 0.0); extE_m = Vector_t(0.0, 0.0, 0.0); extB_m = Vector_t(0.0, 0.0, 0.0); DumpFields::writeFields((((*FieldDimensions.begin())->second).second)); DumpEMFields::writeFields((((*FieldDimensions.begin())->second).second)); function_t func = std::bind(&ParallelCyclotronTracker::getFieldsAtPoint, this, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4); switch ( stepper_m ) { case stepper::INTEGRATOR::RK4: *gmsg << "* 4th order Runge-Kutta integrator" << endl; itsStepper_mp.reset(new RK4(func)); break; case stepper::INTEGRATOR::LF2: *gmsg << "* 2nd order Leap-Frog integrator" << endl; itsStepper_mp.reset(new LF2(func)); break; case stepper::INTEGRATOR::MTS: *gmsg << "* Multiple time stepping (MTS) integrator" << endl; break; case stepper::INTEGRATOR::UNDEFINED: default: itsStepper_mp.reset(nullptr); throw OpalException("ParallelCyclotronTracker::execute", "Invalid name of TIMEINTEGRATOR in Track command"); } if ( stepper_m == stepper::INTEGRATOR::MTS) MtsTracker(); else GenericTracker(); *gmsg << "* ----------------------------------------------- *" << endl; *gmsg << "* Finalizing i.e. write data and close files :" << endl; for(auto fd : FieldDimensions) { ((fd->second).second)->finalise(); } *gmsg << "* ----------------------------------------------- *" << endl; } void ParallelCyclotronTracker::MtsTracker() { /* * variable unit meaning * ------------------------------------------------ * t [ns] time * dt [ns] time step * oldReferenceTheta [rad] azimuth angle * itsBunch_m->R [m] particle position * */ double t = 0, dt = 0, oldReferenceTheta = 0; std::tie(t, dt, oldReferenceTheta) = initializeTracking_m(); int const numSubsteps = std::max(Options::mtsSubsteps, 1); *gmsg << "MTS: Number of substeps per step is " << numSubsteps << endl; double const dt_inner = dt / double(numSubsteps); *gmsg << "MTS: The inner time step is therefore " << dt_inner << " [ns]" << endl; // int SteptoLastInj = itsBunch_m->getSteptoLastInj(); bool flagTransition = false; // flag to determine when to transit from single-bunch to multi-bunches mode *gmsg << "* ---------------------------- Start tracking ----------------------------" << endl; if ( itsBunch_m->hasFieldSolver() ) computeSpaceChargeFields_m(); for(; (step_m < maxSteps_m) && (itsBunch_m->getTotalNum()>0); step_m++) { bool finishedTurn = false; if(step_m % Options::sptDumpFreq == 0) { singleParticleDump(); } Ippl::Comm->barrier(); // First half kick from space charge force if(itsBunch_m->hasFieldSolver()) { kick(0.5 * dt); } // Substeps for external field integration for(int n = 0; n < numSubsteps; ++n) borisExternalFields(dt_inner); // bunch injection // TODO: Where is correct location for this piece of code? Beginning/end of step? Before field solve? injectBunch(flagTransition); if ( itsBunch_m->hasFieldSolver() ) { computeSpaceChargeFields_m(); } else { // if field solver is not available , only update bunch, to transfer particles between nodes if needed, // reset parameters such as LocalNum, initialTotalNum_m. // INFOMSG("No space charge Effects are included!"<R, phi, meanR); itsBunch_m->updateNumTotal(); repartition(); localToGlobal(itsBunch_m->R, phi, meanR); } } // Second half kick from space charge force if(itsBunch_m->hasFieldSolver()) kick(0.5 * dt); // recalculate bingamma and reset the BinID for each particles according to its current gamma if (isMultiBunch() && (step_m % Options::rebinFreq == 0)) { mbHandler_m->updateParticleBins(itsBunch_m); } // dump some data after one push in single particle tracking if ( mode_m == MODE::SINGLE ) { unsigned int i = 0; double temp_meanTheta = calculateAngle2(itsBunch_m->R[i](0), itsBunch_m->R[i](1)); // [-pi, pi] dumpThetaEachTurn_m(t, itsBunch_m->R[i], itsBunch_m->P[i], temp_meanTheta, finishedTurn); dumpAzimuthAngles_m(t, itsBunch_m->R[i], itsBunch_m->P[i], oldReferenceTheta, temp_meanTheta); oldReferenceTheta = temp_meanTheta; } else if ( mode_m == MODE::BUNCH ) { // both for single bunch and multi-bunch // avoid dump at the first step // finishedTurn has not been changed in first push if( isTurnDone() ) { ++turnnumber_m; finishedTurn = true; *gmsg << endl; *gmsg << "*** Finished turn " << turnnumber_m - 1 << ", Total number of live particles: " << itsBunch_m->getTotalNum() << endl; } // Recalculate bingamma and reset the BinID for each particles according to its current gamma if (isMultiBunch() && (step_m % Options::rebinFreq == 0)) { mbHandler_m->updateParticleBins(itsBunch_m); } } // printing + updating bunch parameters + updating t update_m(t, dt, finishedTurn); } // Some post-integration stuff *gmsg << endl; *gmsg << "* ---------------------------- DONE TRACKING PARTICLES -------------------------------- * " << endl; //FIXME dvector_t Ttime = dvector_t(); dvector_t Tdeltr = dvector_t(); dvector_t Tdeltz = dvector_t(); ivector_t TturnNumber = ivector_t(); finalizeTracking_m(Ttime, Tdeltr, Tdeltz, TturnNumber); } /** * * */ void ParallelCyclotronTracker::GenericTracker() { /* * variable unit meaning * ------------------------------------------------ * t [ns] time * dt [ns] time step * oldReferenceTheta [rad] azimuth angle * itsBunch_m->R [m] particle position * */ // Generic Tracker that has three modes defined by timeIntegrator_m: // 0 --- RK-4 (default) // 1 --- LF-2 // (2 --- MTS ... not yet implemented in generic tracker) // mbHandler_m->getNumBunch() determines the number of bunches in multibunch mode (MBM, 1 for OFF) // Total number of particles determines single particle mode (SPM, 1 particle) or // Static Equilibrium Orbit Mode (SEO, 2 particles) double t = 0, dt = 0, oldReferenceTheta = 0; std::tie(t, dt, oldReferenceTheta) = initializeTracking_m(); // If initialTotalNum_m = 2, trigger SEO mode and prepare for transverse tuning calculation // Where is the IF here? -DW dvector_t Ttime, Tdeltr, Tdeltz; ivector_t TturnNumber; // Apply the plugin elements: probe, collimator, stripper, septum once before first step bool flagNeedUpdate = applyPluginElements(dt); // Destroy particles if they are marked as Bin = -1 in the plugin elements // or out of global aperture deleteParticle(flagNeedUpdate); /******************************** * Main integration loop * ********************************/ *gmsg << endl; *gmsg << "* --------------------------------- Start tracking ------------------------------------ *" << endl; for(; (step_m < maxSteps_m) && (itsBunch_m->getTotalNum()>0); step_m++) { bool finishedTurn = false; switch (mode_m) { case MODE::SEO: { // initialTotalNum_m == 2 seoMode_m(t, dt, finishedTurn, Ttime, Tdeltr, Tdeltz, TturnNumber); break; } case MODE::SINGLE: { // initialTotalNum_m == 1 singleMode_m(t, dt, finishedTurn, oldReferenceTheta); break; } case MODE::BUNCH: { // initialTotalNum_m > 2 // Start Tracking for number of particles > 2 (i.e. not single and not SEO mode) bunchMode_m(t, dt, finishedTurn); break; } case MODE::UNDEFINED: default: throw OpalException("ParallelCyclotronTracker::GenericTracker()", "No such tracking mode."); } // Update bunch and some parameters and output some info update_m(t, dt, finishedTurn); } // end for: the integration is DONE after maxSteps_m steps or if all particles are lost! // Some post-integration stuff *gmsg << endl; *gmsg << "* ---------------------------- DONE TRACKING PARTICLES -------------------------------- * " << endl; finalizeTracking_m(Ttime, Tdeltr, Tdeltz, TturnNumber); } bool ParallelCyclotronTracker::getFieldsAtPoint(const double &t, const size_t &Pindex, Vector_t &Efield, Vector_t &Bfield) { bool outOfBound = this->computeExternalFields_m(Pindex, t, Efield, Bfield); // For runs without space charge effects, override this step to save time if(itsBunch_m->hasFieldSolver()) { // Don't do for reference particle if(itsBunch_m->ID[Pindex] != 0) { // add external Field and self space charge field Efield += itsBunch_m->Ef[Pindex]; Bfield += itsBunch_m->Bf[Pindex]; } } return outOfBound; } /** * * * @param Rold * @param Rnew * @param elptr * @param Dold * * @return */ bool ParallelCyclotronTracker::checkGapCross(Vector_t Rold, Vector_t Rnew, RFCavity * rfcavity, double &Dold) { bool flag = false; double sinx = rfcavity->getSinAzimuth(); double cosx = rfcavity->getCosAzimuth(); // TODO: Presumably this is still in mm, so for now, change to m -DW double PerpenDistance = 0.001 * rfcavity->getPerpenDistance(); double distNew = (Rnew[0] * sinx - Rnew[1] * cosx) - PerpenDistance; double distOld = (Rold[0] * sinx - Rold[1] * cosx) - PerpenDistance; if(distOld > 0.0 && distNew <= 0.0) flag = true; // This parameter is used correct cavity phase Dold = 1.0e3 * distOld; // m --> mm return flag; } bool ParallelCyclotronTracker::RFkick(RFCavity * rfcavity, const double t, const double dt, const int Pindex){ // For OPAL 2.0: As long as the RFCavity is in mm, we have to change R to mm here -DW double radius = sqrt(pow(1000.0 * itsBunch_m->R[Pindex](0), 2.0) + pow(1000.0 * itsBunch_m->R[Pindex](1), 2.0) - pow(rfcavity->getPerpenDistance() , 2.0)); double rmin = rfcavity->getRmin(); double rmax = rfcavity->getRmax(); double nomalRadius = (radius - rmin) / (rmax - rmin); double tempP[3]; if(nomalRadius <= 1.0 && nomalRadius >= 0.0) { for(int j = 0; j < 3; j++) tempP[j] = itsBunch_m->P[Pindex](j); //[px,py,pz] units: dimensionless // here evaluate voltage and conduct momenta kicking; rfcavity->getMomentaKick(nomalRadius, tempP, t, dt, itsBunch_m->ID[Pindex], itsBunch_m->getM(), itsBunch_m->getQ()); // t : ns for(int j = 0; j < 3; j++) itsBunch_m->P[Pindex](j) = tempP[j]; return true; } return false; } struct adder : public std::unary_function { adder() : sum(0) {} double sum; void operator()(double x) { sum += x; } }; /** * * * @param t * @param r * @param z * @param lastTurn * @param nur * @param nuz * * @return */ bool ParallelCyclotronTracker::getTunes(dvector_t &t, dvector_t &r, dvector_t &z, int lastTurn, double &nur, double &nuz) { TUNE_class *tune; int Ndat = t.size(); /* remove mean */ double rsum = for_each(r.begin(), r.end(), adder()).sum; for(int i = 0; i < Ndat; i++) r[i] -= rsum; double zsum = for_each(z.begin(), z.end(), adder()).sum; for(int i = 0; i < Ndat; i++) z[i] -= zsum; double ti = *(t.begin()); double tf = t[t.size()-1]; double T = (tf - ti); t.clear(); for(int i = 0; i < Ndat; i++) { t.push_back(i); } T = t[Ndat-1]; *gmsg << endl; *gmsg << "* ************************************* nuR ******************************************* *" << endl; *gmsg << endl << "* ===> " << Ndat << " data points Ti=" << ti << " Tf= " << tf << " -> T= " << T << endl; int nhis_lomb = 10; int stat = 0; // book tune class tune = new TUNE_class(); stat = tune->lombAnalysis(t, r, nhis_lomb, T / lastTurn); if(stat != 0) *gmsg << "* TUNE: Lomb analysis failed" << endl; *gmsg << "* ************************************************************************************* *" << endl; delete tune; tune = NULL; // FIXME: FixMe: need to come from the inputfile nhis_lomb = 100; if(zsum != 0.0) { *gmsg << endl; *gmsg << "* ************************************* nuZ ******************************************* *" << endl; *gmsg << endl << "* ===> " << Ndat << " data points Ti=" << ti << " Tf= " << tf << " -> T= " << T << endl; // book tune class tune = new TUNE_class(); stat = tune->lombAnalysis(t, z, nhis_lomb, T / lastTurn); if(stat != 0) *gmsg << "* TUNE: Lomb analysis failed" << endl; *gmsg << "* ************************************************************************************* *" << endl; delete tune; tune = NULL; } return true; } double ParallelCyclotronTracker::getHarmonicNumber() const { if (opalRing_m != NULL) return opalRing_m->getHarmonicNumber(); Cyclotron* elcycl = dynamic_cast(((*FieldDimensions.begin())->second).second); if (elcycl != NULL) return elcycl->getCyclHarm(); throw OpalException("ParallelCyclotronTracker::getHarmonicNumber()", std::string("The first item in the FieldDimensions list does not ") +std::string("seem to be an Ring or a Cyclotron element")); } Vector_t ParallelCyclotronTracker::calcMeanR(short bunchNr) const { Vector_t meanR(0.0, 0.0, 0.0); for(unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) { // take all particles if bunchNr <= -1 if ( bunchNr > -1 && itsBunch_m->bunchNum[i] != bunchNr) continue; for(int d = 0; d < 3; ++d) { meanR(d) += itsBunch_m->R[i](d); } } reduce(meanR, meanR, OpAddAssign()); size_t n = itsBunch_m->getTotalNum(); if ( bunchNr > -1 ) n = itsBunch_m->getTotalNumPerBunch(bunchNr); return meanR / Vector_t(n); } Vector_t ParallelCyclotronTracker::calcMeanP() const { Vector_t meanP(0.0, 0.0, 0.0); for(unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) { for(int d = 0; d < 3; ++d) { meanP(d) += itsBunch_m->P[i](d); } } reduce(meanP, meanP, OpAddAssign()); return meanP / Vector_t(itsBunch_m->getTotalNum()); } void ParallelCyclotronTracker::repartition() { if((step_m % Options::repartFreq) == 0) { IpplTimings::startTimer(BinRepartTimer_m); itsBunch_m->do_binaryRepart(); Ippl::Comm->barrier(); IpplTimings::stopTimer(BinRepartTimer_m); } } void ParallelCyclotronTracker::globalToLocal(ParticleAttrib & particleVectors, double phi, Vector_t const translationToGlobal) { IpplTimings::startTimer(TransformTimer_m); particleVectors -= translationToGlobal; Tenzor const rotation( cos(phi), sin(phi), 0, -sin(phi), cos(phi), 0, 0, 0, 1); // clockwise rotation for(unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) { particleVectors[i] = dot(rotation, particleVectors[i]); } IpplTimings::stopTimer(TransformTimer_m); } void ParallelCyclotronTracker::localToGlobal(ParticleAttrib & particleVectors, double phi, Vector_t const translationToGlobal) { IpplTimings::startTimer(TransformTimer_m); Tenzor const rotation(cos(phi), -sin(phi), 0, sin(phi), cos(phi), 0, 0, 0, 1); // counter-clockwise rotation for(unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) { particleVectors[i] = dot(rotation, particleVectors[i]); } particleVectors += translationToGlobal; IpplTimings::stopTimer(TransformTimer_m); } inline void ParallelCyclotronTracker::globalToLocal(ParticleAttrib & particleVectors, Quaternion_t const quaternion, Vector_t const meanR) { IpplTimings::startTimer(TransformTimer_m); // Translation from global to local particleVectors -= meanR; // Rotation to align P_mean with x-axis rotateWithQuaternion(particleVectors, quaternion); IpplTimings::stopTimer(TransformTimer_m); } inline void ParallelCyclotronTracker::localToGlobal(ParticleAttrib & particleVectors, Quaternion_t const quaternion, Vector_t const meanR) { IpplTimings::startTimer(TransformTimer_m); // Reverse the quaternion by multiplying the axis components (x,y,z) with -1 Quaternion_t reverseQuaternion = quaternion * -1.0; reverseQuaternion(0) *= -1.0; // Rotation back to original P_mean direction rotateWithQuaternion(particleVectors, reverseQuaternion); // Translation from local to global particleVectors += meanR; IpplTimings::stopTimer(TransformTimer_m); } inline void ParallelCyclotronTracker::globalToLocal(ParticleAttrib & particleVectors, double const phi, double const psi, Vector_t const meanR) { IpplTimings::startTimer(TransformTimer_m); //double const tolerance = 1.0e-4; // TODO What is a good angle threshold? -DW // Translation from global to local particleVectors -= meanR; // Rotation to align P_mean with x-axis rotateAroundZ(particleVectors, phi); // If theta is large enough (i.e. there is significant momentum in z direction) // rotate around x-axis next //if (fabs(psi) > tolerance) rotateAroundX(particleVectors, psi); IpplTimings::stopTimer(TransformTimer_m); } inline void ParallelCyclotronTracker::globalToLocal(Vector_t & myVector, double const phi, double const psi, Vector_t const meanR) { IpplTimings::startTimer(TransformTimer_m); //double const tolerance = 1.0e-4; // TODO What is a good angle threshold? -DW // Translation from global to local myVector -= meanR; // Rotation to align P_mean with x-axis rotateAroundZ(myVector, phi); // If theta is large enough (i.e. there is significant momentum in z direction) // rotate around x-axis next //if (fabs(psi) > tolerance) rotateAroundX(myVector, psi); IpplTimings::stopTimer(TransformTimer_m); } inline void ParallelCyclotronTracker::localToGlobal(ParticleAttrib & particleVectors, double const phi, double const psi, Vector_t const meanR) { IpplTimings::startTimer(TransformTimer_m); //double const tolerance = 1.0e-4; // TODO What is a good angle threshold? -DW // If theta is large enough (i.e. there is significant momentum in z direction) // rotate back around x-axis next //if (fabs(psi) > tolerance) rotateAroundX(particleVectors, -psi); // Rotation to align P_mean with x-axis rotateAroundZ(particleVectors, -phi); // Translation from local to global particleVectors += meanR; IpplTimings::stopTimer(TransformTimer_m); } inline void ParallelCyclotronTracker::localToGlobal(Vector_t & myVector, double const phi, double const psi, Vector_t const meanR) { IpplTimings::startTimer(TransformTimer_m); //double const tolerance = 1.0e-4; // TODO What is a good angle threshold? -DW // If theta is large enough (i.e. there is significant momentum in z direction) // rotate back around x-axis next //if (fabs(psi) > tolerance) rotateAroundX(myVector, -psi); // Rotation to align P_mean with x-axis rotateAroundZ(myVector, -phi); // Translation from local to global myVector += meanR; IpplTimings::stopTimer(TransformTimer_m); } inline void ParallelCyclotronTracker::rotateWithQuaternion(ParticleAttrib & particleVectors, Quaternion_t const quaternion) { Vector_t const quaternionVectorComponent = Vector_t(quaternion(1), quaternion(2), quaternion(3)); double const quaternionScalarComponent = quaternion(0); for(unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) { particleVectors[i] = 2.0f * dot(quaternionVectorComponent, particleVectors[i]) * quaternionVectorComponent + (quaternionScalarComponent * quaternionScalarComponent - dot(quaternionVectorComponent, quaternionVectorComponent)) * particleVectors[i] + 2.0f * quaternionScalarComponent * cross(quaternionVectorComponent, particleVectors[i]); } } inline void ParallelCyclotronTracker::normalizeQuaternion(Quaternion_t & quaternion){ double tolerance = 1.0e-10; double length2 = dot(quaternion, quaternion); if (fabs(length2) > tolerance && fabs(length2 - 1.0f) > tolerance) { double length = sqrt(length2); quaternion /= length; } } inline void ParallelCyclotronTracker::normalizeVector(Vector_t & vector) { double tolerance = 1.0e-10; double length2 = dot(vector, vector); if (fabs(length2) > tolerance && fabs(length2 - 1.0f) > tolerance) { double length = sqrt(length2); vector /= length; } } inline void ParallelCyclotronTracker::rotateAroundZ(ParticleAttrib & particleVectors, double const phi) { // Clockwise rotation of particles 'particleVectors' by 'phi' around Z axis Tenzor const rotation( cos(phi), sin(phi), 0, -sin(phi), cos(phi), 0, 0, 0, 1); for(unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) { particleVectors[i] = dot(rotation, particleVectors[i]); } } inline void ParallelCyclotronTracker::rotateAroundZ(Vector_t & myVector, double const phi) { // Clockwise rotation of single vector 'myVector' by 'phi' around Z axis Tenzor const rotation( cos(phi), sin(phi), 0, -sin(phi), cos(phi), 0, 0, 0, 1); myVector = dot(rotation, myVector); } inline void ParallelCyclotronTracker::rotateAroundX(ParticleAttrib & particleVectors, double const psi) { // Clockwise rotation of particles 'particleVectors' by 'psi' around X axis Tenzor const rotation(1, 0, 0, 0, cos(psi), sin(psi), 0, -sin(psi), cos(psi)); for(unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) { particleVectors[i] = dot(rotation, particleVectors[i]); } } inline void ParallelCyclotronTracker::rotateAroundX(Vector_t & myVector, double const psi) { // Clockwise rotation of single vector 'myVector' by 'psi' around X axis Tenzor const rotation(1, 0, 0, 0, cos(psi), sin(psi), 0, -sin(psi), cos(psi)); myVector = dot(rotation, myVector); } inline void ParallelCyclotronTracker::getQuaternionTwoVectors(Vector_t u, Vector_t v, Quaternion_t & quaternion) { // four vector (w,x,y,z) of the quaternion of P_mean with the positive x-axis normalizeVector(u); normalizeVector(v); double k_cos_theta = dot(u, v); double k = sqrt(dot(u, u) * dot(v, v)); double tolerance1 = 1.0e-5; double tolerance2 = 1.0e-8; Vector_t resultVectorComponent; if (fabs(k_cos_theta / k + 1.0) < tolerance1) { // u and v are almost exactly antiparallel so we need to do // 180 degree rotation around any vector orthogonal to u resultVectorComponent = cross(u, xaxis); // If by chance u is parallel to xaxis, use zaxis instead if (dot(resultVectorComponent, resultVectorComponent) < tolerance2) { resultVectorComponent = cross(u, zaxis); } double halfAngle = 0.5 * pi; double sinHalfAngle = sin(halfAngle); resultVectorComponent *= sinHalfAngle; k = 0.0; k_cos_theta = cos(halfAngle); } else { resultVectorComponent = cross(u, v); } quaternion(0) = k_cos_theta + k; quaternion(1) = resultVectorComponent(0); quaternion(2) = resultVectorComponent(1); quaternion(3) = resultVectorComponent(2); normalizeQuaternion(quaternion); } bool ParallelCyclotronTracker::push(double h) { /* h [ns] * dt1 [ns] * dt2 [ns] */ IpplTimings::startTimer(IntegrationTimer_m); // h [ns] --> h [s] h *= 1.0e-9; bool flagNeedUpdate = false; for(unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) { Vector_t const oldR = itsBunch_m->R[i]; double const gamma = sqrt(1.0 + dot(itsBunch_m->P[i], itsBunch_m->P[i])); double const c_gamma = Physics::c / gamma; Vector_t const v = itsBunch_m->P[i] * c_gamma; itsBunch_m->R[i] += h * v; for(const auto & ccd : cavCrossDatas_m) { double const distNew = (itsBunch_m->R[i][0] * ccd.sinAzimuth - itsBunch_m->R[i][1] * ccd.cosAzimuth) - ccd.perpenDistance; bool tagCrossing = false; double distOld; if(distNew <= 0.0) { distOld = (oldR[0] * ccd.sinAzimuth - oldR[1] * ccd.cosAzimuth) - ccd.perpenDistance; if(distOld > 0.0) tagCrossing = true; } if(tagCrossing) { double const dt1 = distOld / sqrt(dot(v, v)); double const dt2 = h - dt1; // Retrack particle from the old postion to cavity gap point itsBunch_m->R[i] = oldR + dt1 * v; // Momentum kick //itsBunch_m->R[i] *= 1000.0; // RFkick uses [itsBunch_m->R] == mm RFkick(ccd.cavity, itsBunch_m->getT() * 1.0e9, dt1, i); //itsBunch_m->R[i] *= 0.001; itsBunch_m->R[i] += dt2 * itsBunch_m->P[i] * c_gamma; } } flagNeedUpdate |= (itsBunch_m->Bin[i] < 0); } IpplTimings::stopTimer(IntegrationTimer_m); return flagNeedUpdate; } bool ParallelCyclotronTracker::kick(double h) { IpplTimings::startTimer(IntegrationTimer_m); bool flagNeedUpdate = false; BorisPusher pusher; double const q = itsBunch_m->Q[0] / q_e; // For now all particles have the same charge double const M = itsBunch_m->M[0] * 1.0e9; // For now all particles have the same rest energy for(unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) { pusher.kick(itsBunch_m->R[i], itsBunch_m->P[i], itsBunch_m->Ef[i], itsBunch_m->Bf[i], h * 1.0e-9, M, q); flagNeedUpdate |= (itsBunch_m->Bin[i] < 0); } IpplTimings::stopTimer(IntegrationTimer_m); return flagNeedUpdate; } void ParallelCyclotronTracker::borisExternalFields(double h) { // h in [ns] // push particles for first half step bool flagNeedUpdate = push(0.5 * h); // Evaluate external fields IpplTimings::startTimer(IntegrationTimer_m); for(unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) { itsBunch_m->Ef[i] = Vector_t(0.0, 0.0, 0.0); itsBunch_m->Bf[i] = Vector_t(0.0, 0.0, 0.0); computeExternalFields_m(i, itsBunch_m->getT() * 1e9 /*ns*/, itsBunch_m->Ef[i], itsBunch_m->Bf[i]); } IpplTimings::stopTimer(IntegrationTimer_m); // Kick particles for full step flagNeedUpdate |= kick(h); // push particles for second half step flagNeedUpdate |= push(0.5 * h); // apply the plugin elements: probe, collimator, stripper, septum flagNeedUpdate |= applyPluginElements(h); // destroy particles if they are marked as Bin=-1 in the plugin elements or out of global apeture deleteParticle(flagNeedUpdate); } bool ParallelCyclotronTracker::applyPluginElements(const double dt) { IpplTimings::startTimer(PluginElemTimer_m); // Plugin Elements are all defined in mm, change beam to mm before applying itsBunch_m->R *= Vector_t(1000.0); for(beamline_list::iterator sindex = ++(FieldDimensions.begin()); sindex != FieldDimensions.end(); ++sindex) { if(((*sindex)->first) == ElementBase::BEAMSTRIPPING) { BeamStripping *bstp = static_cast(((*sindex)->second).second); bstp->checkBeamStripping(itsBunch_m, cycl_m, turnnumber_m, itsBunch_m->getT()*1e9 /*[ns]*/, dt); } } bool flag = false; for (PluginElement* element : pluginElements_m) { bool tmp = element->check(itsBunch_m, turnnumber_m, itsBunch_m->getT() * 1e9 /*[ns]*/, dt); flag |= tmp; if ( tmp ) { itsBunch_m->updateNumTotal(); *gmsg << "* Total number of particles after PluginElement= " << itsBunch_m->getTotalNum() << endl; } } itsBunch_m->R *= Vector_t(0.001); IpplTimings::stopTimer(PluginElemTimer_m); return flag; } bool ParallelCyclotronTracker::deleteParticle(bool flagNeedUpdate){ IpplTimings::startTimer(DelParticleTimer_m); // Update immediately if any particles are lost during this step allreduce(flagNeedUpdate, 1, std::logical_or()); if(flagNeedUpdate) { short bunchCount = itsBunch_m->getNumBunch(); std::vector locLostParticleNum(bunchCount, 0); const int leb = itsBunch_m->getLastemittedBin(); std::unique_ptr localBinCount; if ( isMultiBunch() ) { localBinCount = std::unique_ptr(new size_t[leb]); for (int i = 0; i < leb; ++i) localBinCount[i] = 0; } for(unsigned int i = 0; i < itsBunch_m->getLocalNum(); i++) { if(itsBunch_m->Bin[i] < 0) { ++locLostParticleNum[itsBunch_m->bunchNum[i]]; itsBunch_m->destroy(1, i); } else if ( isMultiBunch() ) { /* we need to count the local number of particles * per energy bin */ ++localBinCount[itsBunch_m->Bin[i]]; } } if ( isMultiBunch() ) { // set the local bin count for (int i = 0; i < leb; ++i) { itsBunch_m->setLocalBinCount(localBinCount[i], i); } } std::vector localnum(bunchCount + 1); for (size_t i = 0; i < localnum.size() - 1; ++i) { localnum[i] = itsBunch_m->getLocalNumPerBunch(i) - locLostParticleNum[i]; itsBunch_m->setLocalNumPerBunch(localnum[i], i); } /* We need to destroy the particles now * before we compute the means. We also * have to update the total number of particles * otherwise the statistics are wrong. */ itsBunch_m->performDestroy(true); /* total number of particles of individual bunches * last index of vector contains total number over all * bunches, used as a check */ std::vector totalnum(bunchCount + 1); localnum[bunchCount] = itsBunch_m->getLocalNum(); allreduce(localnum.data(), totalnum.data(), localnum.size(), std::plus()); itsBunch_m->setTotalNum(totalnum[bunchCount]); for (short i = 0; i < bunchCount; ++i) { itsBunch_m->setTotalNumPerBunch(totalnum[i], i); } size_t sum = std::accumulate(totalnum.begin(), totalnum.end() - 1, 0); if ( sum != totalnum[bunchCount] ) { throw OpalException("ParallelCyclotronTracker::deleteParticle()", "Total number of particles " + std::to_string(totalnum[bunchCount]) + " != " + std::to_string(sum) + " (sum over all bunches)"); } size_t globLostParticleNum = 0; size_t locNumLost = std::accumulate(locLostParticleNum.begin(), locLostParticleNum.end(), 0); reduce(locNumLost, globLostParticleNum, 1, std::plus()); if ( globLostParticleNum > 0 ) { *gmsg << "At step " << step_m << ", lost " << globLostParticleNum << " particles " << "on stripper, collimator, septum, geometry, " << "out of cyclotron aperture or beam stripping" << endl; } if (totalnum[bunchCount] == 0) { IpplTimings::stopTimer(DelParticleTimer_m); return flagNeedUpdate; } Vector_t const meanR = calcMeanR(); Vector_t const meanP = calcMeanP(); // Bunch (local) azimuth at meanR w.r.t. y-axis double const phi = calculateAngle(meanP(0), meanP(1)) - 0.5 * pi; // Bunch (local) elevation at meanR w.r.t. xy plane double const psi = 0.5 * pi - acos(meanP(2) / sqrt(dot(meanP, meanP))); // For statistics data, the bunch is transformed into a local coordinate system // with meanP in direction of y-axis -DW globalToLocal(itsBunch_m->R, phi, psi, meanR); globalToLocal(itsBunch_m->P, phi, psi, Vector_t(0.0)); // P should be rotated, but not shifted //itsBunch_m->R *= Vector_t(0.001); // mm --> m // Now destroy particles and update pertinent parameters in local frame // Note that update() will be called within boundp() -DW itsBunch_m->boundp(); //itsBunch_m->update(); itsBunch_m->calcBeamParameters(); //itsBunch_m->R *= Vector_t(1000.0); // m --> mm localToGlobal(itsBunch_m->R, phi, psi, meanR); localToGlobal(itsBunch_m->P, phi, psi, Vector_t(0.0)); // If particles were deleted, recalculate bingamma and reset BinID for remaining particles if ( isMultiBunch() ) mbHandler_m->updateParticleBins(itsBunch_m); } IpplTimings::stopTimer(DelParticleTimer_m); return flagNeedUpdate; } void ParallelCyclotronTracker::initTrackOrbitFile() { std::string f = OpalData::getInstance()->getInputBasename() + std::string("-trackOrbit.dat"); outfTrackOrbit_m.setf(std::ios::scientific, std::ios::floatfield); outfTrackOrbit_m.precision(8); if(myNode_m == 0) { if(OpalData::getInstance()->inRestartRun()) { outfTrackOrbit_m.open(f.c_str(), std::ios::app); outfTrackOrbit_m << "# Restart at integration step " << itsBunch_m->getLocalTrackStep() << std::endl; } else { outfTrackOrbit_m.open(f.c_str()); outfTrackOrbit_m << "# The six-dimensional phase space data in the global Cartesian coordinates" << std::endl; outfTrackOrbit_m << "# Part. ID x [m] beta_x*gamma y [m] beta_y*gamma z [m] beta_z*gamma" << std::endl; } } } void ParallelCyclotronTracker::initDistInGlobalFrame() { if(!OpalData::getInstance()->inRestartRun()) { // Start a new run (no restart) double const initialReferenceTheta = referenceTheta * Physics::deg2rad; // TODO: Replace with TracerParticle // Force the initial phase space values of the particle with ID = 0 to zero, // to set it as a reference particle. if(initialTotalNum_m > 2) { for(size_t i = 0; i < initialLocalNum_m; ++i) { if(itsBunch_m->ID[i] == 0) { itsBunch_m->R[i] = Vector_t(0.0); itsBunch_m->P[i] = Vector_t(0.0); } } } // Initialize global R //itsBunch_m->R *= Vector_t(1000.0); // m --> mm // NEW OPAL 2.0: Immediately change to m -DW Vector_t const initMeanR = Vector_t(0.001 * referenceR * cosRefTheta_m, // mm --> m 0.001 * referenceR * sinRefTheta_m, // mm --> m 0.001 * referenceZ); // mm --> m localToGlobal(itsBunch_m->R, initialReferenceTheta, initMeanR); // Initialize global P (Cartesian, but input P_ref is in Pr, Ptheta, Pz, // so translation has to be done before the rotation this once) // Cave: In the local frame, the positive y-axis is the direction of movement -DW for(size_t i = 0; i < initialLocalNum_m; ++i) { itsBunch_m->P[i](0) += referencePr; itsBunch_m->P[i](1) += referencePt; itsBunch_m->P[i](2) += referencePz; } // Out of the three coordinates of meanR (R, Theta, Z) only the angle // changes the momentum vector... localToGlobal(itsBunch_m->P, initialReferenceTheta); // Initialize the bin number of the first bunch to 0 for(size_t i = 0; i < initialLocalNum_m; ++i) { itsBunch_m->Bin[i] = 0; } // Backup initial distribution if multi bunch mode if ((initialTotalNum_m > 2) && isMultiBunch() && mbHandler_m->isForceMode()) { mbHandler_m->saveBunch(itsBunch_m); } // Else: Restart from the distribution in the h5 file } else { // Do a local frame restart (we have already checked that the old h5 file was saved in local // frame as well). if((Options::psDumpFrame != Options::GLOBAL)) { *gmsg << "* Restart in the local frame" << endl; //itsBunch_m->R *= Vector_t(1000.0); // m --> mm // referenceR and referenceZ are already in mm // New OPAL 2.0: Init in m -DW Vector_t const initMeanR = Vector_t(0.001 * referenceR * cosRefTheta_m, 0.001 * referenceR * sinRefTheta_m, 0.001 * referenceZ); // Do the tranformations localToGlobal(itsBunch_m->R, referencePhi, referencePsi, initMeanR); localToGlobal(itsBunch_m->P, referencePhi, referencePsi); // Initialize the bin number of the first bunch to 0 for(size_t i = 0; i < initialLocalNum_m; ++i) { itsBunch_m->Bin[i] = 0; } // Or do a global frame restart (no transformations necessary) } else { *gmsg << "* Restart in the global frame" << endl; pathLength_m = itsBunch_m->get_sPos(); //itsBunch_m->R *= Vector_t(1000.0); // m --> mm } } // set the number of particles per bunch itsBunch_m->countTotalNumPerBunch(); // ------------ Get some Values ---------------------------------------------------------- // Vector_t const meanR = calcMeanR(); Vector_t const meanP = calcMeanP(); // Bunch (local) azimuth at meanR w.r.t. y-axis double const phi = calculateAngle(meanP(0), meanP(1)) - 0.5 * pi; // Bunch (local) elevation at meanR w.r.t. xy plane double const psi = 0.5 * pi - acos(meanP(2) / sqrt(dot(meanP, meanP))); double radius = std::sqrt(meanR[0] * meanR[0] + meanR[1] * meanR[1]); // [m] if ( isMultiBunch() ) mbHandler_m->setRadiusTurns(radius); // Do boundp and repartition in the local frame at beginning of this run globalToLocal(itsBunch_m->R, phi, psi, meanR); globalToLocal(itsBunch_m->P, phi, psi); // P should be rotated, but not shifted //itsBunch_m->R *= Vector_t(0.001); // mm --> m itsBunch_m->boundp(); checkNumPart(std::string("* Before repartition: ")); repartition(); checkNumPart(std::string("* After repartition: ")); itsBunch_m->calcBeamParameters(); *gmsg << endl << "* *********************** Bunch information in local frame: ************************"; *gmsg << *itsBunch_m << endl; //itsBunch_m->R *= Vector_t(1000.0); // m --> mm localToGlobal(itsBunch_m->R, phi, psi, meanR); localToGlobal(itsBunch_m->P, phi, psi); // Save initial distribution if not a restart if(!OpalData::getInstance()->inRestartRun()) { step_m -= 1; bunchDumpPhaseSpaceData(); bunchDumpStatData(); step_m += 1; } // Print out the Bunch information at beginning of the run. Because the bunch information // displays in units of m we have to change back and forth one more time. //itsBunch_m->R *= Vector_t(0.001); // mm --> m itsBunch_m->calcBeamParameters(); // multi-bunch simulation only saveInjectValues(); *gmsg << endl << "* *********************** Bunch information in global frame: ***********************"; *gmsg << *itsBunch_m << endl; //itsBunch_m->R *= Vector_t(1000.0); // m --> mm } //TODO: This can be completely rewritten with TracerParticle -DW void ParallelCyclotronTracker::singleParticleDump() { IpplTimings::startTimer(DumpTimer_m); if(Ippl::getNodes() > 1 ) { double x; int id; dvector_t tmpr; ivector_t tmpi; int tag = Ippl::Comm->next_tag(IPPL_APP_TAG4, IPPL_APP_CYCLE); // for all nodes, find the location of particle with ID = 0 & 1 in bunch containers int found[2] = {-1, -1}; int counter = 0; for(size_t i = 0; i < itsBunch_m->getLocalNum(); ++i) { if(itsBunch_m->ID[i] == 0) { found[counter] = i; counter++; } if(itsBunch_m->ID[i] == 1) { found[counter] = i; counter++; } } if(myNode_m == 0) { int notReceived = Ippl::getNodes() - 1; int numberOfPart = 0; // receive other nodes while(notReceived > 0) { int node = COMM_ANY_NODE; Message *rmsg = Ippl::Comm->receive_block(node, tag); if(rmsg == nullptr) ERRORMSG("Could not receive from client nodes in main." << endl); --notReceived; rmsg->get(&numberOfPart); for(int i = 0; i < numberOfPart; ++i) { rmsg->get(&id); tmpi.push_back(id); for(int ii = 0; ii < 6; ii++) { rmsg->get(&x); tmpr.push_back(x); } } delete rmsg; } // own node for(int i = 0; i < counter; ++i) { tmpi.push_back(itsBunch_m->ID[found[i]]); for(int j = 0; j < 3; ++j) { tmpr.push_back(itsBunch_m->R[found[i]](j)); tmpr.push_back(itsBunch_m->P[found[i]](j)); } } // store dvector_t::iterator itParameter = tmpr.begin(); for(auto tmpid : tmpi) { outfTrackOrbit_m << "ID" << tmpid; if (tmpid == 0) { // for stat file itsBunch_m->RefPartR_m[0] = *itParameter; itsBunch_m->RefPartR_m[1] = *(itParameter + 2); itsBunch_m->RefPartR_m[2] = *(itParameter + 4); itsBunch_m->RefPartP_m[0] = *(itParameter + 1); itsBunch_m->RefPartP_m[1] = *(itParameter + 3); itsBunch_m->RefPartP_m[2] = *(itParameter + 5); } for(int ii = 0; ii < 6; ii++) { outfTrackOrbit_m << " " << *itParameter; ++itParameter; } outfTrackOrbit_m << std::endl; } } else { // for other nodes Message *smsg = new Message(); smsg->put(counter); for(int i = 0; i < counter; i++) { smsg->put(itsBunch_m->ID[found[i]]); for(int j = 0; j < 3; j++) { smsg->put(itsBunch_m->R[found[i]](j)); smsg->put(itsBunch_m->P[found[i]](j)); } } if(!Ippl::Comm->send(smsg, 0, tag)) ERRORMSG("Ippl::Comm->send(smsg, 0, tag) failed " << endl); } Ippl::Comm->barrier(); } else { for(size_t i = 0; i < itsBunch_m->getLocalNum(); i++) { if(itsBunch_m->ID[i] == 0 || itsBunch_m->ID[i] == 1) { outfTrackOrbit_m << "ID" << itsBunch_m->ID[i] << " "; outfTrackOrbit_m << itsBunch_m->R[i](0) << " " << itsBunch_m->P[i](0) << " "; outfTrackOrbit_m << itsBunch_m->R[i](1) << " " << itsBunch_m->P[i](1) << " "; outfTrackOrbit_m << itsBunch_m->R[i](2) << " " << itsBunch_m->P[i](2) << std::endl; if (itsBunch_m->ID[i] == 0) { // for stat file itsBunch_m->RefPartR_m = itsBunch_m->R[i]; itsBunch_m->RefPartP_m = itsBunch_m->P[i]; } } } } IpplTimings::stopTimer(DumpTimer_m); } void ParallelCyclotronTracker::bunchDumpStatData(){ IpplTimings::startTimer(DumpTimer_m); // dump stat file per bunch in case of multi-bunch mode if (isMultiBunch()) { double phi = 0.0, psi = 0.0; Vector_t meanR = calcMeanR(); // Bunch (global) angle w.r.t. x-axis (cylinder coordinates) double theta = calculateAngle(meanR(0), meanR(1)) * Physics::rad2deg; // fix azimuth_m --> make monotonically increasing dumpAngle(theta, prevAzimuth_m, azimuth_m); updateAzimuthAndRadius(); if(Options::psDumpFrame != Options::GLOBAL) { Vector_t meanP = calcMeanP(); // Bunch (local) azimuth at meanR w.r.t. y-axis phi = calculateAngle(meanP(0), meanP(1)) - 0.5 * pi; // Bunch (local) elevation at meanR w.r.t. xy plane psi = 0.5 * pi - acos(meanP(2) / sqrt(dot(meanP, meanP))); // Rotate so Pmean is in positive y direction. No shift, so that normalized emittance and // unnormalized emittance as well as centroids are calculated correctly in // PartBunch::calcBeamParameters() globalToLocal(itsBunch_m->R, phi, psi, meanR); globalToLocal(itsBunch_m->P, phi, psi); } itsDataSink->writeMultiBunchStatistics(itsBunch_m, mbHandler_m.get()); if(Options::psDumpFrame != Options::GLOBAL) { localToGlobal(itsBunch_m->R, phi, psi, meanR); localToGlobal(itsBunch_m->P, phi, psi); } IpplTimings::stopTimer(DumpTimer_m); return; } /* in case before a bunchDumpPhasespaceData has happen, the calcBeamParameters() is not needed. maybe we have to introduce a dirty-stat-data flag. this flag would be set if the bunch is dirty and reset by calcBeamParameters */ itsBunch_m->R *= Vector_t(0.001); // mm --> m itsBunch_m->calcBeamParameters(); itsBunch_m->R *= Vector_t(1000.0); // m --> mm // --------------------------------- Get some Values ---------------------------------------- // double const temp_t = itsBunch_m->getT() * 1e9; // s -> ns Vector_t meanR; Vector_t meanP; if (Options::psDumpFrame == Options::BUNCH_MEAN) { meanR = calcMeanR(); meanP = calcMeanP(); } else if (itsBunch_m->getLocalNum() > 0) { meanR = itsBunch_m->R[0]; meanP = itsBunch_m->P[0]; } double phi = 0; double psi = 0; // Bunch (global) angle w.r.t. x-axis (cylinder coordinates) double azimuth = calculateAngle(meanR(0), meanR(1)) * Physics::rad2deg; // fix azimuth_m --> make monotonically increasing dumpAngle(azimuth, prevAzimuth_m, azimuth_m); // -------------- Calculate the external fields at the center of the bunch ----------------- // beamline_list::iterator DumpSindex = FieldDimensions.begin(); extE_m = Vector_t(0.0, 0.0, 0.0); extB_m = Vector_t(0.0, 0.0, 0.0); (((*DumpSindex)->second).second)->apply(meanR, meanP, temp_t, extE_m, extB_m); // If we are saving in local frame, bunch and fields at the bunch center have to be rotated // TODO: Make decision if we maybe want to always save statistics data in local frame? -DW if(Options::psDumpFrame != Options::GLOBAL) { // -------------------- ----------- Do Transformations ---------------------------------- // // Bunch (local) azimuth at meanR w.r.t. y-axis phi = calculateAngle(meanP(0), meanP(1)) - 0.5 * pi; // Bunch (local) elevation at meanR w.r.t. xy plane psi = 0.5 * pi - acos(meanP(2) / sqrt(dot(meanP, meanP))); // Rotate so Pmean is in positive y direction. No shift, so that normalized emittance and // unnormalized emittance as well as centroids are calculated correctly in // PartBunch::calcBeamParameters() globalToLocal(extB_m, phi, psi); globalToLocal(extE_m, phi, psi); globalToLocal(itsBunch_m->R, phi, psi); globalToLocal(itsBunch_m->P, phi, psi); } //itsBunch_m->R *= Vector_t(0.001); // mm -> m FDext_m[0] = extB_m / 10.0; // kgauss --> T FDext_m[1] = extE_m; // kV/mm? -DW // Save the stat file itsDataSink->dumpSDDS(itsBunch_m, FDext_m, azimuth_m); //itsBunch_m->R *= Vector_t(1000.0); // m -> mm // If we are in local mode, transform back after saving if(Options::psDumpFrame != Options::GLOBAL) { localToGlobal(itsBunch_m->R, phi, psi); localToGlobal(itsBunch_m->P, phi, psi); } IpplTimings::stopTimer(DumpTimer_m); } void ParallelCyclotronTracker::bunchDumpPhaseSpaceData() { // --------------------------------- Particle dumping --------------------------------------- // // Note: Don't dump when // 1. after one turn // in order to sychronize the dump step for multi-bunch and single bunch for comparison // with each other during post-process phase. // ------------------------------------------------------------------------------------------ // IpplTimings::startTimer(DumpTimer_m); // --------------------------------- Get some Values ---------------------------------------- // double const temp_t = itsBunch_m->getT() * 1.0e9; // s -> ns Vector_t meanR; Vector_t meanP; // in case of multi-bunch mode take always bunch mean (although it takes all bunches) if (Options::psDumpFrame == Options::BUNCH_MEAN || isMultiBunch()) { meanR = calcMeanR(); meanP = calcMeanP(); } else if (itsBunch_m->getLocalNum() > 0) { meanR = itsBunch_m->R[0]; meanP = itsBunch_m->P[0]; } double const betagamma_temp = sqrt(dot(meanP, meanP)); double const E = itsBunch_m->get_meanKineticEnergy(); // Bunch (global) angle w.r.t. x-axis (cylinder coordinates) double const theta = atan2(meanR(1), meanR(0)); // Bunch (local) azimuth at meanR w.r.t. y-axis double const phi = calculateAngle(meanP(0), meanP(1)) - 0.5 * pi; // Bunch (local) elevation at meanR w.r.t. xy plane double const psi = 0.5 * pi - acos(meanP(2) / sqrt(dot(meanP, meanP))); // ---------------- Re-calculate reference values in format of input values ----------------- // // Position: // New OPAL 2.0: Init in m (back to mm just for output) -DW referenceR = computeRadius(meanR); // includes m --> mm conversion referenceTheta = theta / Physics::deg2rad; referenceZ = 1000.0 * meanR(2); // Momentum in Theta-hat, R-hat, Z-hat coordinates at position meanR: referencePtot = betagamma_temp; referencePz = meanP(2); referencePr = meanP(0) * cos(theta) + meanP(1) * sin(theta); referencePt = sqrt(referencePtot * referencePtot - \ referencePz * referencePz - referencePr * referencePr); // ----- Calculate the external fields at the center of the bunch (Cave: Global Frame) ----- // beamline_list::iterator DumpSindex = FieldDimensions.begin(); extE_m = Vector_t(0.0, 0.0, 0.0); extB_m = Vector_t(0.0, 0.0, 0.0); (((*DumpSindex)->second).second)->apply(meanR, meanP, temp_t, extE_m, extB_m); FDext_m[0] = extB_m * 0.1; // kgauss --> T FDext_m[1] = extE_m; // -------------- If flag psDumpFrame is global, dump bunch in global frame ------------- // if (Options::psDumpFreq < std::numeric_limits::max() ){ bool dumpLocalFrame = true; std::string dumpString = "local"; if (Options::psDumpFrame == Options::GLOBAL) { dumpLocalFrame = false; dumpString = "global"; } if (dumpLocalFrame == true) { // ---------------- If flag psDumpFrame is local, dump bunch in local frame ---------------- // // The bunch is transformed into a local coordinate system with meanP in direction of y-axis globalToLocal(itsBunch_m->R, phi, psi, meanR); //globalToLocal(itsBunch_m->R, phi, psi, meanR * Vector_t(0.001)); globalToLocal(itsBunch_m->P, phi, psi); // P should only be rotated globalToLocal(extB_m, phi, psi); globalToLocal(extE_m, phi, psi); } FDext_m[0] = extB_m * 0.1; // kgauss --> T FDext_m[1] = extE_m; lastDumpedStep_m = itsDataSink->dumpH5(itsBunch_m, // Local and in m FDext_m, E, referencePr, referencePt, referencePz, referenceR, referenceTheta, referenceZ, phi / Physics::deg2rad, // P_mean azimuth // at ref. R/Th/Z psi / Physics::deg2rad, // P_mean elevation // at ref. R/Th/Z dumpLocalFrame); // Flag localFrame if (dumpLocalFrame == true) { // Return to global frame localToGlobal(itsBunch_m->R, phi, psi, meanR); //localToGlobal(itsBunch_m->R, phi, psi, meanR * Vector_t(0.001)); localToGlobal(itsBunch_m->P, phi, psi); } // Tell user in which mode we are dumping // New: no longer dumping for num part < 3, omit phase space dump number info if (lastDumpedStep_m == -1){ *gmsg << endl << "* Integration step " << step_m + 1 << " (no phase space dump for <= 2 particles)" << endl; } else { *gmsg << endl << "* Phase space dump " << lastDumpedStep_m << " (" << dumpString << " frame) at integration step " << step_m + 1 << endl; } } // Print dump information on screen *gmsg << "* T = " << temp_t << " ns" << ", Live Particles: " << itsBunch_m->getTotalNum() << endl; *gmsg << "* E = " << E << " MeV" << ", beta * gamma = " << betagamma_temp << endl; *gmsg << "* Bunch position: R = " << referenceR << " mm" << ", Theta = " << referenceTheta << " Deg" << ", Z = " << referenceZ << " mm" << endl; *gmsg << "* Local Azimuth = " << phi / Physics::deg2rad << " Deg" << ", Local Elevation = " << psi / Physics::deg2rad << " Deg" << endl; IpplTimings::stopTimer(DumpTimer_m); } bool ParallelCyclotronTracker::isTurnDone() { return (step_m > 10) && (((step_m + 1) %setup_m.stepsPerTurn) == 0); } void ParallelCyclotronTracker::update_m(double& t, const double& dt, const bool& finishedTurn) { // Reference parameters t += dt; updateTime(dt); itsBunch_m->setLocalTrackStep((step_m + 1)); if (!(step_m + 1 % 1000)) *gmsg << "Step " << step_m + 1 << endl; updatePathLength(dt); // Here is global frame, don't do itsBunch_m->boundp(); if (itsBunch_m->getTotalNum()>0) { // Only dump last step if we have particles left. // Check separately for phase space (ps) and statistics (stat) data dump frequency if ( mode_m != MODE::SEO && ( ((step_m + 1) % Options::psDumpFreq == 0) || (Options::psDumpEachTurn && finishedTurn))) { // Write phase space data to h5 file bunchDumpPhaseSpaceData(); } if ( mode_m != MODE::SEO && ( ((step_m + 1) % Options::statDumpFreq == 0) || (Options::psDumpEachTurn && finishedTurn))) { // Write statistics data to stat file bunchDumpStatData(); } } if (Options::psDumpEachTurn && finishedTurn) for (PluginElement* element : pluginElements_m) element->save(); } std::tuple ParallelCyclotronTracker::initializeTracking_m() { // Read in some control parameters setup_m.scSolveFreq = (spiral_flag) ? 1 : Options::scSolveFreq; setup_m.stepsPerTurn = itsBunch_m->getStepsPerTurn(); // Define 3 special azimuthal angles where dump particle's six parameters // at each turn into 3 ASCII files. only important in single particle tracking azimuth_angle_m.resize(3); azimuth_angle_m[0] = 0.0; azimuth_angle_m[1] = 22.5 * Physics::deg2rad; azimuth_angle_m[2] = 45.0 * Physics::deg2rad; double harm = getHarmonicNumber(); double dt = itsBunch_m->getdT() * 1.0e9 * harm; // time step size (s --> ns) double t = itsBunch_m->getT() * 1.0e9; // current time (s --> ns) double oldReferenceTheta = referenceTheta * Physics::deg2rad; // init here, reset each step setup_m.deltaTheta = pi / (setup_m.stepsPerTurn); // half of the average angle per step //int stepToLastInj = itsBunch_m->getSteptoLastInj(); // TODO: Do we need this? -DW // Record how many bunches have already been injected. ONLY FOR MBM if (isMultiBunch()) mbHandler_m->setNumBunch(itsBunch_m->getNumBunch()); initTrackOrbitFile(); // Get data from h5 file for restart run and reset current step // to last step of previous simulation if(OpalData::getInstance()->inRestartRun()) { restartStep0_m = itsBunch_m->getLocalTrackStep(); step_m = restartStep0_m; turnnumber_m = step_m / setup_m.stepsPerTurn + 1; *gmsg << "* Restart at integration step " << restartStep0_m << " at turn " << turnnumber_m - 1 << endl; initPathLength(); } setup_m.stepsNextCheck = step_m + setup_m.stepsPerTurn; // Steps to next check for transition initDistInGlobalFrame(); if (isMultiBunch()) mbHandler_m->updateParticleBins(itsBunch_m); // --- Output to user --- // *gmsg << "* Beginning of this run is at t = " << t << " [ns]" << endl; *gmsg << "* The time step is set to dt = " << dt << " [ns]" << endl; if ( isMultiBunch() ) { *gmsg << "* MBM: Time interval between neighbour bunches is set to " << setup_m.stepsPerTurn * dt << "[ns]" << endl; *gmsg << "* MBM: The particles energy bin reset frequency is set to " << Options::rebinFreq << endl; } *gmsg << "* Single particle trajectory dump frequency is set to " << Options::sptDumpFreq << endl; *gmsg << "* The frequency to solve space charge fields is set to " << setup_m.scSolveFreq << endl; *gmsg << "* The repartition frequency is set to " << Options::repartFreq << endl; switch ( mode_m ) { case MODE::SEO: *gmsg << endl; *gmsg << "* ------------------------- STATIC EQUILIBRIUM ORBIT MODE ----------------------------- *" << endl << "* Instruction: When the total particle number is equal to 2, SEO mode is triggered *" << endl << "* automatically. This mode does NOT include any RF cavities. The initial distribution *" << endl << "* file must be specified. In the file the first line is for reference particle and the *" << endl << "* second line is for off-center particle. The tune is calculated by FFT routines based *" << endl << "* on these two particles. *" << endl << "* ---------------- NOTE: SEO MODE ONLY WORKS SERIALLY ON SINGLE NODE ------------------ *" << endl; if(Ippl::getNodes() != 1) throw OpalException("Error in ParallelCyclotronTracker::initializeTracking_m", "SEO MODE ONLY WORKS SERIALLY ON SINGLE NODE!"); break; case MODE::SINGLE: *gmsg << endl; *gmsg << "* ------------------------------ SINGLE PARTICLE MODE --------------------------------- *" << endl << "* Instruction: When the total particle number is equal to 1, single particle mode is *" << endl << "* triggered automatically. The initial distribution file must be specified which should *" << endl << "* contain only one line for the single particle *" << endl << "* ---------NOTE: SINGLE PARTICLE MODE ONLY WORKS SERIALLY ON A SINGLE NODE ------------ *" << endl; if(Ippl::getNodes() != 1) throw OpalException("Error in ParallelCyclotronTracker::initializeTracking_m", "SINGLE PARTICLE MODE ONLY WORKS SERIALLY ON A SINGLE NODE!"); // For single particle mode open output files openFiles(azimuth_angle_m.size() + 1, OpalData::getInstance()->getInputBasename()); break; case MODE::BUNCH: break; case MODE::UNDEFINED: default: throw OpalException("ParallelCyclotronTracker::initializeTracking_m()", "No such tracking mode."); } return std::make_tuple(t, dt, oldReferenceTheta); } void ParallelCyclotronTracker::finalizeTracking_m(dvector_t& Ttime, dvector_t& Tdeltr, dvector_t& Tdeltz, ivector_t& TturnNumber) { for(size_t ii = 0; ii < (itsBunch_m->getLocalNum()); ii++) { if(itsBunch_m->ID[ii] == 0) { double FinalMomentum2 = pow(itsBunch_m->P[ii](0), 2.0) + pow(itsBunch_m->P[ii](1), 2.0) + pow(itsBunch_m->P[ii](2), 2.0); double FinalEnergy = (sqrt(1.0 + FinalMomentum2) - 1.0) * itsBunch_m->getM() * 1.0e-6; *gmsg << "* Final energy of reference particle = " << FinalEnergy << " [MeV]" << endl; *gmsg << "* Total phase space dump number(includes the initial distribution) = " << lastDumpedStep_m + 1 << endl; *gmsg << "* One can restart simulation from the last dump step (--restart " << lastDumpedStep_m << ")" << endl; } } Ippl::Comm->barrier(); switch ( mode_m ) { case MODE::SEO: { // Calculate tunes after tracking. *gmsg << endl; *gmsg << "* **************** The result for tune calulation (NO space charge) ******************* *" << endl << "* Number of tracked turns: " << TturnNumber.back() << endl; double nur, nuz; getTunes(Ttime, Tdeltr, Tdeltz, TturnNumber.back(), nur, nuz); break; } case MODE::SINGLE: closeFiles(); // fall through case MODE::BUNCH: // we do nothing case MODE::UNDEFINED: default: { // not for multibunch if (!isMultiBunch()) { *gmsg << "*" << endl; *gmsg << "* Finished during turn " << turnnumber_m << " (" << turnnumber_m - 1 << " turns completed)" << endl; *gmsg << "* Cave: Turn number is not correct for restart mode"<< endl; } } } Ippl::Comm->barrier(); if (myNode_m == 0) outfTrackOrbit_m.close(); *gmsg << endl << "* *********************** Bunch information in global frame: ***********************"; if (itsBunch_m->getTotalNum() > 0){ // Print out the Bunch information at end of the run. itsBunch_m->calcBeamParameters(); *gmsg << *itsBunch_m << endl; } else { *gmsg << endl << "* No Particles left in bunch!" << endl; *gmsg << "* **********************************************************************************" << endl; } } void ParallelCyclotronTracker::seoMode_m(double& t, const double dt, bool& finishedTurn, dvector_t& Ttime, dvector_t& Tdeltr, dvector_t& Tdeltz, ivector_t& TturnNumber) { // 2 particles: Trigger SEO mode // (Switch off cavity and calculate betatron oscillation tuning) double r_tuning[2], z_tuning[2] ; IpplTimings::startTimer(IntegrationTimer_m); for(size_t i = 0; i < (itsBunch_m->getLocalNum()); i++) { if((step_m % Options::sptDumpFreq == 0)) { outfTrackOrbit_m << "ID" << (itsBunch_m->ID[i]); outfTrackOrbit_m << " " << itsBunch_m->R[i](0) << " " << itsBunch_m->P[i](0) << " " << itsBunch_m->R[i](1) << " " << itsBunch_m->P[i](1) << " " << itsBunch_m->R[i](2) << " " << itsBunch_m->P[i](2) << std::endl; } double OldTheta = calculateAngle(itsBunch_m->R[i](0), itsBunch_m->R[i](1)); r_tuning[i] = itsBunch_m->R[i](0) * cos(OldTheta) + itsBunch_m->R[i](1) * sin(OldTheta); z_tuning[i] = itsBunch_m->R[i](2); // Integrate for one step in the lab Cartesian frame (absolute value). itsStepper_mp->advance(itsBunch_m, i, t, dt); if( (i == 0) && isTurnDone() ) turnnumber_m++; } // end for: finish one step tracking for all particles for initialTotalNum_m == 2 mode IpplTimings::stopTimer(IntegrationTimer_m); // store dx and dz for future tune calculation if higher precision needed, reduce freqSample. if(step_m % Options::sptDumpFreq == 0) { Ttime.push_back(t * 1.0e-9); Tdeltz.push_back(z_tuning[1]); Tdeltr.push_back(r_tuning[1] - r_tuning[0]); TturnNumber.push_back(turnnumber_m); } } void ParallelCyclotronTracker::singleMode_m(double& t, const double dt, bool& finishedTurn, double& oldReferenceTheta) { // 1 particle: Trigger single particle mode // ********************************************************************************** // // * This was moved here b/c collision should be tested before the actual * // // * timestep (bgf_main_collision_test() predicts the next step automatically) * // // apply the plugin elements: probe, collimator, stripper, septum bool flagNeedUpdate = applyPluginElements(dt); // check if we lose particles at the boundary bgf_main_collision_test(); // ********************************************************************************** // if (itsBunch_m->getLocalNum() == 0) return; // might happen if particle is in collimator IpplTimings::startTimer(IntegrationTimer_m); unsigned int i = 0; // we only have a single particle if ( step_m % Options::sptDumpFreq == 0 ) { outfTrackOrbit_m << "ID" <ID[i] << " " << itsBunch_m->R[i](0) << " " << itsBunch_m->P[i](0) << " " << itsBunch_m->R[i](1) << " " << itsBunch_m->P[i](1) << " " << itsBunch_m->R[i](2) << " " << itsBunch_m->P[i](2) << std::endl; } double temp_meanTheta = calculateAngle2(itsBunch_m->R[i](0), itsBunch_m->R[i](1)); // [-pi, pi] dumpThetaEachTurn_m(t, itsBunch_m->R[i], itsBunch_m->P[i], temp_meanTheta, finishedTurn); dumpAzimuthAngles_m(t, itsBunch_m->R[i], itsBunch_m->P[i], oldReferenceTheta, temp_meanTheta); oldReferenceTheta = temp_meanTheta; // used for gap crossing checking Vector_t Rold = itsBunch_m->R[i]; // [x,y,z] (mm) Vector_t Pold = itsBunch_m->P[i]; // [px,py,pz] (beta*gamma) // integrate for one step in the lab Cartesian frame (absolute value). itsStepper_mp->advance(itsBunch_m, i, t, dt); // If gap crossing happens, do momenta kicking (not if gap crossing just happened) if (itsBunch_m->cavityGapCrossed[i] == true) itsBunch_m->cavityGapCrossed[i] = false; else gapCrossKick_m(i, t, dt, Rold, Pold); IpplTimings::stopTimer(IntegrationTimer_m); // Destroy particles if they are marked as Bin = -1 in the plugin elements // or out of the global aperture flagNeedUpdate |= (itsBunch_m->Bin[i] < 0); deleteParticle(flagNeedUpdate); } void ParallelCyclotronTracker::bunchMode_m(double& t, const double dt, bool& finishedTurn) { // Flag for transition from single-bunch to multi-bunches mode static bool flagTransition = false; // single particle dumping if(step_m % Options::sptDumpFreq == 0) singleParticleDump(); // Find out if we need to do bunch injection injectBunch(flagTransition); // oldReferenceTheta = calculateAngle2(meanR(0), meanR(1)); // Calculate SC field before each time step and keep constant during integration. // Space Charge effects are included only when total macropaticles number is NOT LESS THAN 1000. if (itsBunch_m->hasFieldSolver()) { if (step_m % setup_m.scSolveFreq == 0) { computeSpaceChargeFields_m(); } else { // If we are not solving for the space charge fields at this time step // we will apply the fields from the previous step and have to rotate them // accordingly. For this we find the quaternion between the previous mean momentum (PreviousMeanP) // and the current mean momentum (meanP) and rotate the fields with this quaternion. Vector_t const meanP = calcMeanP(); Quaternion_t quaternionToNewMeanP; getQuaternionTwoVectors(PreviousMeanP, meanP, quaternionToNewMeanP); // Reset PreviousMeanP. Cave: This HAS to be after the quaternion is calculated! PreviousMeanP = calcMeanP(); // Rotate the fields globalToLocal(itsBunch_m->Ef, quaternionToNewMeanP); globalToLocal(itsBunch_m->Bf, quaternionToNewMeanP); } } // *** This was moved here b/c collision should be tested before the ********************** // *** actual timestep (bgf_main_collision_test() predicts the next step automatically) -DW // Apply the plugin elements: probe, collimator, stripper, septum bool flagNeedUpdate = applyPluginElements(dt); // check if we lose particles at the boundary bgf_main_collision_test(); IpplTimings::startTimer(IntegrationTimer_m); for(size_t i = 0; i < itsBunch_m->getLocalNum(); i++) { // used for gap crossing checking Vector_t Rold = itsBunch_m->R[i]; // [x,y,z] (mm) Vector_t Pold = itsBunch_m->P[i]; // [px,py,pz] (beta*gamma) // Integrate for one step in the lab Cartesian frame (absolute value). itsStepper_mp->advance(itsBunch_m, i, t, dt); // If gap crossing happens, do momenta kicking (not if gap crossing just happened) if (itsBunch_m->cavityGapCrossed[i] == true) itsBunch_m->cavityGapCrossed[i] = false; else gapCrossKick_m(i, t, dt, Rold, Pold); flagNeedUpdate |= (itsBunch_m->Bin[i] < 0); } IpplTimings::stopTimer(IntegrationTimer_m); // Destroy particles if they are marked as Bin = -1 in the plugin elements // or out of global aperture deleteParticle(flagNeedUpdate); // Recalculate bingamma and reset the BinID for each particles according to its current gamma if (isMultiBunch() && step_m % Options::rebinFreq == 0) { mbHandler_m->updateParticleBins(itsBunch_m); } // Some status output for user after each turn if ( isTurnDone() ) { turnnumber_m++; finishedTurn = true; *gmsg << endl; *gmsg << "*** Finished turn " << turnnumber_m - 1 << ", Total number of live particles: " << itsBunch_m->getTotalNum() << endl; } Ippl::Comm->barrier(); } void ParallelCyclotronTracker::gapCrossKick_m(size_t i, double t, double dt, const Vector_t& Rold, const Vector_t& Pold) { for (beamline_list::iterator sindex = ++(FieldDimensions.begin()); sindex != FieldDimensions.end(); ++sindex) { bool tag_crossing = false; double DistOld = 0.0; //mm RFCavity * rfcav; if (((*sindex)->first) == ElementBase::RFCAVITY) { // here check gap cross in the list, if do , set tag_crossing to TRUE rfcav = static_cast(((*sindex)->second).second); tag_crossing = checkGapCross(Rold, itsBunch_m->R[i], rfcav, DistOld); } if ( tag_crossing ) { itsBunch_m->cavityGapCrossed[i] = true; double oldMomentum2 = dot(Pold, Pold); double oldBetgam = sqrt(oldMomentum2); double oldGamma = sqrt(1.0 + oldMomentum2); double oldBeta = oldBetgam / oldGamma; double dt1 = DistOld / (Physics::c * oldBeta * 1.0e-6); // ns double dt2 = dt - dt1; // retrack particle from the old postion to cavity gap point // restore the old coordinates and momenta itsBunch_m->R[i] = Rold; itsBunch_m->P[i] = Pold; if (dt / dt1 < 1.0e9) itsStepper_mp->advance(itsBunch_m, i, t, dt1); // Momentum kick RFkick(rfcav, t, dt1, i); /* Retrack particle from cavity gap point for * the left time to finish the entire timestep */ if (dt / dt2 < 1.0e9) itsStepper_mp->advance(itsBunch_m, i, t, dt2); } } } void ParallelCyclotronTracker::dumpAzimuthAngles_m(const double& t, const Vector_t& R, const Vector_t& P, const double& oldReferenceTheta, const double& temp_meanTheta) { for (unsigned int i=0; i<=2; i++) { if ((oldReferenceTheta < azimuth_angle_m[i] - setup_m.deltaTheta) && ( temp_meanTheta >= azimuth_angle_m[i] - setup_m.deltaTheta)) { *(outfTheta_m[i]) << "#Turn number = " << turnnumber_m << ", Time = " << t << " [ns]" << std::endl << " " << std::hypot(R(0), R(1)) << " " << P(0) * std::cos(temp_meanTheta) + P(1) * std::sin(temp_meanTheta) << " " << temp_meanTheta * Physics::rad2deg << " " << - P(0) * std::sin(temp_meanTheta) + P(1) * std::cos(temp_meanTheta) << " " << R(2) << " " << P(2) << std::endl; } } } void ParallelCyclotronTracker::dumpThetaEachTurn_m(const double& t, const Vector_t& R, const Vector_t& P, const double& temp_meanTheta, bool& finishedTurn) { if ( isTurnDone() ) { ++turnnumber_m; finishedTurn = true; *gmsg << "* SPT: Finished turn " << turnnumber_m - 1 << endl; *(outfTheta_m[3]) << "#Turn number = " << turnnumber_m << ", Time = " << t << " [ns]" << std::endl << " " << std::sqrt(R(0) * R(0) + R(1) * R(1)) << " " << P(0) * std::cos(temp_meanTheta) + P(1) * std::sin(temp_meanTheta) << " " << temp_meanTheta / Physics::deg2rad << " " << - P(0) * std::sin(temp_meanTheta) + P(1) * std::cos(temp_meanTheta) << " " << R(2) << " " << P(2) << std::endl; } } void ParallelCyclotronTracker::computeSpaceChargeFields_m() { // Firstly reset E and B to zero before fill new space charge field data for each track step itsBunch_m->Bf = Vector_t(0.0); itsBunch_m->Ef = Vector_t(0.0); if (spiral_flag && itsBunch_m->getFieldSolverType() == "SAAMG") { // --- Single bunch mode with spiral inflector --- // // If we are doing a spiral inflector simulation and are using the SAAMG solver // we don't rotate or translate the bunch and gamma is 1.0 (non-relativistic). //itsBunch_m->R *= Vector_t(0.001); // mm --> m itsBunch_m->setGlobalMeanR(Vector_t(0.0, 0.0, 0.0)); itsBunch_m->setGlobalToLocalQuaternion(Quaternion_t(1.0, 0.0, 0.0, 0.0)); itsBunch_m->computeSelfFields_cycl(1.0); //itsBunch_m->R *= Vector_t(1000.0); // m --> mm } else { Vector_t const meanR = calcMeanR(); // (m) PreviousMeanP = calcMeanP(); // Since calcMeanP takes into account all particles of all bins (TODO: Check this! -DW) // Using the quaternion method with PreviousMeanP and yaxis should give the correct result Quaternion_t quaternionToYAxis; getQuaternionTwoVectors(PreviousMeanP, yaxis, quaternionToYAxis); globalToLocal(itsBunch_m->R, quaternionToYAxis, meanR); //itsBunch_m->R *= Vector_t(0.001); // mm --> m if ((step_m + 1) % Options::boundpDestroyFreq == 0) itsBunch_m->boundp_destroy(); else itsBunch_m->boundp(); if (hasMultiBunch()) { // --- Multibunch mode --- // // Calculate gamma for each energy bin itsBunch_m->calcGammas_cycl(); repartition(); // Calculate space charge field for each energy bin for(int b = 0; b < itsBunch_m->getLastemittedBin(); b++) { itsBunch_m->setBinCharge(b, itsBunch_m->getChargePerParticle()); //itsBunch_m->setGlobalMeanR(0.001 * meanR); itsBunch_m->setGlobalMeanR(meanR); itsBunch_m->setGlobalToLocalQuaternion(quaternionToYAxis); itsBunch_m->computeSelfFields_cycl(b); } itsBunch_m->Q = itsBunch_m->getChargePerParticle(); } else { // --- Single bunch mode --- // double temp_meangamma = std::sqrt(1.0 + dot(PreviousMeanP, PreviousMeanP)); repartition(); //itsBunch_m->setGlobalMeanR(0.001 * meanR); itsBunch_m->setGlobalMeanR(meanR); itsBunch_m->setGlobalToLocalQuaternion(quaternionToYAxis); itsBunch_m->computeSelfFields_cycl(temp_meangamma); } //scale coordinates back //itsBunch_m->R *= Vector_t(1000.0); // m --> mm // Transform coordinates back to global localToGlobal(itsBunch_m->R, quaternionToYAxis, meanR); // Transform self field back to global frame (rotate only) localToGlobal(itsBunch_m->Ef, quaternionToYAxis); localToGlobal(itsBunch_m->Bf, quaternionToYAxis); } } bool ParallelCyclotronTracker::computeExternalFields_m(const size_t& i, const double& t, Vector_t& Efield, Vector_t& Bfield) { beamline_list::iterator sindex = FieldDimensions.begin(); // Flag whether a particle is out of field bool outOfBound = (((*sindex)->second).second)->apply(i, t, Efield, Bfield); Bfield *= 0.10; // kGauss --> T Efield *= 1.0e6; // kV/mm --> V/m return outOfBound; } void ParallelCyclotronTracker::injectBunch(bool& flagTransition) { if (!isMultiBunch() || step_m != setup_m.stepsNextCheck) { return; } const short result = mbHandler_m->injectBunch(itsBunch_m, itsReference, flagTransition); switch ( result ) { case 0: { // nothing happened break; } case 1: { // bunch got saved saveInjectValues(); setup_m.stepsNextCheck += setup_m.stepsPerTurn; if ( flagTransition ) { *gmsg << "* MBM: Saving beam distribution at turn " << turnnumber_m << endl; *gmsg << "* MBM: After one revolution, Multi-Bunch Mode will be invoked" << endl; } break; } case 2: { // bunch got injected setup_m.stepsNextCheck += setup_m.stepsPerTurn; break; } default: throw OpalException("ParallelCyclotronTracker::injectBunch()", "Unknown return value " + std::to_string(result)); } } void ParallelCyclotronTracker::saveInjectValues() { if ( !isMultiBunch() ) { return; } Vector_t meanR = calcMeanR(); // Bunch (global) angle w.r.t. x-axis (cylinder coordinates) double theta = calculateAngle(meanR(0), meanR(1)) * Physics::rad2deg; // fix azimuth_m --> make monotonically increasing dumpAngle(theta, prevAzimuth_m, azimuth_m); const double radius = computeRadius(meanR); MultiBunchHandler::injection_t& inj = mbHandler_m->getInjectionValues(); inj.time = itsBunch_m->getT() * 1.0e9; /*ns*/ inj.pathlength = itsBunch_m->get_sPos(); inj.azimuth = azimuth_m; inj.radius = radius; } void ParallelCyclotronTracker::updatePathLength(const double& dt) { /* the last element includes all particles, * all other elements are bunch number specific */ std::vector lpaths(1); if ( isMultiBunch() ) { lpaths.resize(mbHandler_m->getNumBunch() + 1); } computePathLengthUpdate(lpaths, dt); pathLength_m += lpaths.back(); itsBunch_m->set_sPos(pathLength_m); if ( isMultiBunch() ) { mbHandler_m->updatePathLength(lpaths); } } void ParallelCyclotronTracker::updateTime(const double& dt) { // t is in seconds double t = itsBunch_m->getT(); // dt: ns --> s itsBunch_m->setT(t + dt * 1.0e-9); if ( isMultiBunch() ) { mbHandler_m->updateTime(dt); } } void ParallelCyclotronTracker::updateAzimuthAndRadius() { if (!isMultiBunch()) { return; } for (short b = 0; b < mbHandler_m->getNumBunch(); ++b) { Vector_t meanR = calcMeanR(b); MultiBunchHandler::beaminfo_t& binfo = mbHandler_m->getBunchInfo(b); binfo.radius = computeRadius(meanR); double azimuth = calculateAngle(meanR(0), meanR(1)) * Physics::rad2deg; dumpAngle(azimuth, binfo.prevAzimuth, binfo.azimuth, b); } } void ParallelCyclotronTracker::initPathLength() { if ( isMultiBunch() ) { // we need to reset the path length of each bunch itsDataSink->setMultiBunchInitialPathLengh(mbHandler_m.get()); } }