Commit c0da8d5d authored by frey_m's avatar frey_m
Browse files

add multi-bunch functions

parent 617459ce
......@@ -223,12 +223,12 @@ bool MultiBunchHandler::readBunch(PartBunchBase<double, 3> *beam,
}
int MultiBunchHandler::injectBunch(PartBunchBase<double, 3> *beam,
const PartData& ref,
bool& flagTransition,
const double& azimuth)
short MultiBunchHandler::injectBunch(PartBunchBase<double, 3> *beam,
const PartData& ref,
bool& flagTransition,
const double& azimuth)
{
int result = 0;
short result = 0;
if ((bunchCount_m == 1) && (mode_m == MB_MODE::AUTO) && (!flagTransition)) {
// If all of the following conditions are met, this code will be executed
......
......@@ -48,10 +48,10 @@ public:
* 1 - if bunch got saved
* 2 - if bunch got injected
*/
int injectBunch(PartBunchBase<double, 3> *beam,
const PartData& ref,
bool& flagTransition,
const double& azimuth);
short injectBunch(PartBunchBase<double, 3> *beam,
const PartData& ref,
bool& flagTransition,
const double& azimuth);
void updateParticleBins(PartBunchBase<double, 3> *beam);
......@@ -64,24 +64,36 @@ public:
void setRadiusTurns(const double& radius);
/// set total number of tracked bunches
void setTotalNumBunch(int n);
void setTotalNumBunch(short n);
/// get total number of tracked bunches
int getTotalNumBunch() const;
short getTotalNumBunch() const;
void setNumBunch(int n);
void setNumBunch(short n);
int getNumBunch() const;
short getNumBunch() const;
bool isForceMode() const;
void setPathLength(const double& pathlength, short bunchNr);
const double& getPathLength(short bunchNr) const;
void setAzimuth(const double& azimuth, short bunchNr);
const double& getAzimuth(short bunchNr) const;
void setTime(const double& time, short bunchNr);
const double& getTime(short bunchNr) const;
private:
// store the data of the beam which are required for injecting a
// new bunch for multibunch filename
std::string onebunch_m;
/// The number of bunches specified in TURNS of RUN commond
int numBunch_m;
short numBunch_m;
// parameter for reset bin in multi-bunch run
double eta_m;
......@@ -105,30 +117,30 @@ private:
DataSink& ds_m;
// record how many bunches has already been injected.
int bunchCount_m;
short bunchCount_m;
};
inline
void MultiBunchHandler::setTotalNumBunch(int n) {
void MultiBunchHandler::setTotalNumBunch(short n) {
numBunch_m = n;
}
inline
int MultiBunchHandler::getTotalNumBunch() const {
short MultiBunchHandler::getTotalNumBunch() const {
return numBunch_m;
}
inline
void MultiBunchHandler::setNumBunch(int n) {
void MultiBunchHandler::setNumBunch(short n) {
bunchCount_m = n;
}
inline
int MultiBunchHandler::getNumBunch() const {
short MultiBunchHandler::getNumBunch() const {
return bunchCount_m;
}
......@@ -138,4 +150,46 @@ bool MultiBunchHandler::isForceMode() const {
return (mode_m == MB_MODE::FORCE);
}
inline
void MultiBunchHandler::setPathLength(const double& pathlength, short bunchNr) {
MultiBunchDump *mbd = ds_m.getMultiBunchWriter(bunchNr);
mbd->setPathLength(pathlength);
}
inline
const double& MultiBunchHandler::getPathLength(short bunchNr) const {
MultiBunchDump *mbd = ds_m.getMultiBunchWriter(bunchNr);
return mbd->getPathLength();
}
inline
void MultiBunchHandler::setAzimuth(const double& azimuth, short bunchNr) {
MultiBunchDump *mbd = ds_m.getMultiBunchWriter(bunchNr);
mbd->setAzimuth(azimuth);
}
inline
const double& MultiBunchHandler::getAzimuth(short bunchNr) const {
MultiBunchDump *mbd = ds_m.getMultiBunchWriter(bunchNr);
return mbd->getAzimuth();
}
inline
void MultiBunchHandler::setTime(const double& time, short bunchNr) {
MultiBunchDump *mbd = ds_m.getMultiBunchWriter(bunchNr);
mbd->setTime(time);
}
inline
const double& MultiBunchHandler::getTime(short bunchNr) const {
MultiBunchDump *mbd = ds_m.getMultiBunchWriter(bunchNr);
return mbd->getTime();
}
#endif
......@@ -207,20 +207,23 @@ void ParallelCyclotronTracker::bgf_main_collision_test() {
}
// only used for dumping into stat file
void ParallelCyclotronTracker::dumpAngle_m(const double& theta) {
if ( prevAzimuth_m < 0.0 ) { // only at first occurrence
azimuth_m = theta;
void ParallelCyclotronTracker::dumpAngle(const double& theta,
double& prevAzimuth,
double& azimuth)
{
if ( prevAzimuth < 0.0 ) { // only at first occurrence
azimuth = theta;
} else {
double dtheta = theta - prevAzimuth_m;
double dtheta = theta - prevAzimuth;
if ( dtheta < 0 ) {
dtheta += 360.0;
}
if ( dtheta > 180 ) { // rotating clockwise, reduce angle
dtheta -= 360;
}
azimuth_m += dtheta;
azimuth += dtheta;
}
prevAzimuth_m = theta;
prevAzimuth = theta;
}
......@@ -1596,17 +1599,27 @@ double ParallelCyclotronTracker::getHarmonicNumber() const {
}
Vector_t ParallelCyclotronTracker::calcMeanR() const {
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());
return meanR / Vector_t(itsBunch_m->getTotalNum());
size_t n = itsBunch_m->getTotalNum();
if ( bunchNr > -1 )
n = itsBunch_m->getTotalNumPerBunch(bunchNr);
return meanR / Vector_t(n);
}
Vector_t ParallelCyclotronTracker::calcMeanP() const {
......@@ -1958,11 +1971,10 @@ bool ParallelCyclotronTracker::push(double h) {
}
flagNeedUpdate |= (itsBunch_m->Bin[i] < 0);
}
itsBunch_m->setT(itsBunch_m->getT() + h);
// Path length update
pathLength_m += computePathLengthUpdate(h * 1.0e9 /*s --> ns*/);
itsBunch_m->setLPath(pathLength_m);
updateTime(h * 1.0e9 /*s --> ns*/);
updatePathLength(h * 1.0e9 /*s --> ns*/);
IpplTimings::stopTimer(IntegrationTimer_m);
return flagNeedUpdate;
......@@ -2461,7 +2473,7 @@ void ParallelCyclotronTracker::bunchDumpStatData(){
double theta = calculateAngle(meanR(0), meanR(1)) * Physics::rad2deg;
// fix azimuth_m --> make monotonically increasing
dumpAngle_m(theta);
dumpAngle(theta, prevAzimuth_m, azimuth_m);
if(Options::psDumpFrame != Options::GLOBAL) {
Vector_t meanP = calcMeanP();
......@@ -2519,7 +2531,7 @@ void ParallelCyclotronTracker::bunchDumpStatData(){
double azimuth = calculateAngle(meanR(0), meanR(1)) * Physics::rad2deg;
// fix azimuth_m --> make monotonically increasing
dumpAngle_m(azimuth);
dumpAngle(azimuth, prevAzimuth_m, azimuth_m);
// -------------- Calculate the external fields at the center of the bunch ----------------- //
beamline_list::iterator DumpSindex = FieldDimensions.begin();
......@@ -2705,14 +2717,14 @@ void ParallelCyclotronTracker::update_m(double& t, const double& dt,
{
// Reference parameters
t += dt;
itsBunch_m->setT(t * 1.0e-9);
updateTime(dt);
itsBunch_m->setLocalTrackStep((step_m + 1));
if (!(step_m + 1 % 1000))
*gmsg << "Step " << step_m + 1 << endl;
pathLength_m += computePathLengthUpdate(dt); // unit: m
itsBunch_m->setLPath(pathLength_m);
updatePathLength(dt);
// Here is global frame, don't do itsBunch_m->boundp();
......@@ -3324,12 +3336,12 @@ bool ParallelCyclotronTracker::computeExternalFields_m(const size_t& i, const do
void ParallelCyclotronTracker::injectBunch(bool& flagTransition) {
if ( !isMultiBunch() && step_m != setup_m.stepsNextCheck ) ) {
if (!isMultiBunch() && step_m != setup_m.stepsNextCheck) {
return;
}
const int result = mbHandler_m->injectBunch(itsBunch_m, itsReference,
flagTransition, azimuth_m);
const short result = mbHandler_m->injectBunch(itsBunch_m, itsReference,
flagTransition, azimuth_m);
switch ( result ) {
case 0: {
......@@ -3355,3 +3367,48 @@ void ParallelCyclotronTracker::injectBunch(bool& flagTransition) {
"Unknown return value " + std::to_string(result));
}
}
void ParallelCyclotronTracker::updatePathLength(const double& dt) {
pathLength_m += computePathLengthUpdate(dt);
itsBunch_m->setLPath(pathLength_m);
if ( isMultiBunch() ) {
for (short b = 0; b < mbHandler_m->getNumBunch(); ++b) {
double lpath = mbHandler_m->getPathLength(b);
lpath += computePathLengthUpdate(dt, b);
mbHandler_m->setPathLength(lpath, b);
}
}
}
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() ) {
for (short b = 0; b < mbHandler_m->getNumBunch(); ++b) {
double time = mbHandler_m->getTime(b);
mbHandler_m->setTime(time + dt, b);
}
}
}
void ParallelCyclotronTracker::updateAzimuth() {
for (short b = 0; b < mbHandler_m->getNumBunch(); ++b) {
Vector_t meanR = calcMeanR(b);
double azimuth = calculateAngle(meanR(0), meanR(1)) * Physics::rad2deg;
const double& prevAzimuth = mbHandler_m->getAzimuth();
double newAzimuth = prevAzimuth;
dumpAngle(azimuth, prevAzimuth, newAzimuth);
mbHandler_m->setAzimuth(newAzimuth);
}
}
......@@ -275,8 +275,14 @@ private:
double azimuth_m;
double prevAzimuth_m;
// only for dumping to stat-file --> make azimuth monotonically increasing
void dumpAngle_m(const double& theta);
/* only for dumping to stat-file --> make azimuth monotonically increasing
* @param theta computed azimuth [deg]
* @param prevAzimuth previous angle [deg]
* @param azimuth to be updated [deg]
*/
void dumpAngle(const double& theta,
double& prevAzimuth,
double& azimuth);
// take all particles if bunchNr <= -1
double computePathLengthUpdate(const double& dt, short bunchNr = -1);
......@@ -318,7 +324,11 @@ private:
IpplTimings::TimerRef PluginElemTimer_m;
IpplTimings::TimerRef DelParticleTimer_m;
Vector_t calcMeanR() const;
/*
* @param bunchNr if <= -1 --> take all particles in simulation, if bunchNr > -1,
* take only particles of *this* bunch
*/
Vector_t calcMeanR(short bunchNr = -1) const;
Vector_t calcMeanP() const;
......@@ -492,6 +502,18 @@ private:
bool isMultiBunch() const;
bool hasMultiBunch() const;
/*
* @param dt time step in ns
*/
void updatePathLength(const double& dt);
/*
* @param dt time step in ns
*/
void updateTime(const double& dt);
void updateAzimuth();
};
/**
......
......@@ -9,6 +9,7 @@ public:
struct beaminfo_t {
double time;
double azimuth;
double pathlength;
long unsigned int nParticles;
double ekin;
double dEkin;
......@@ -29,6 +30,36 @@ public:
bool calcBunchBeamParameters(PartBunchBase<double, 3>* beam);
/*!
* @param azimuth in deg
*/
void setAzimuth(const double& azimuth);
/*!
* @returns azimuth in deg
*/
const double& getAzimuth() const;
/*!
* @param time in ns
*/
void setTime(const double& time);
/*!
* @preturns time in ns
*/
const double& getTime() const;
/*!
* @param pathlength in m
*/
void setPathLength(const double& pathlength);
/*!
* @returns pathlength in m
*/
const double& getPathLength() const;
private:
bool isNotFirst_m;
......@@ -38,4 +69,40 @@ private:
beaminfo_t binfo_m;
};
inline
void MultiBunchDump::setAzimuth(const double& azimuth) {
binfo_m.azimuth = azimuth;
}
inline
const double& MultiBunchDump::getAzimuth() const {
return binfo_m.azimuth;
}
inline
void MultiBunchDump::setTime(const double& time) {
binfo_m.time = time;
}
inline
const double& MultiBunchDump::getTime() const {
return binfo_m.time;
}
inline
void MultiBunchDump::setPathLength(const double& pathlength) {
binfo_m.pathlength = pathlength;
}
inline
const double& MultiBunchDump::getPathLength() const {
return binfo_m.pathlength;
}
#endif
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment