Skip to content

GitLab

Commit 4d4b095e authored by kraus's avatar kraus
Browse files

introducing enum for particle types

parent 3ad092ea
Hide whitespace changes
Inline Side-by-side
Showing with 167 additions and 159 deletions
classic/5.0/src/AbsBeamline/Collimator.cpp
View file @ 4d4b095e
......@@ -111,7 +111,7 @@ Collimator::Collimator(const Collimator &right):
lossDs_m(nullptr),
sphys_m(NULL)
{
setGeom();
setGeom();
}
......@@ -156,8 +156,8 @@ Collimator::Collimator(const std::string &name):
Collimator::~Collimator() {
if(online_m)
goOffline();
if(online_m)
goOffline();
}
......@@ -167,61 +167,61 @@ void Collimator::accept(BeamlineVisitor &visitor) const {
inline bool Collimator::isInColl(Vector_t R, Vector_t P, double recpgamma) {
/**
check if we are in the longitudinal
range of the collimator
*/
const double z = R(2) + P(2) * recpgamma;
if((z > position_m) && (z <= position_m + getElementLength())) {
if(isAPepperPot_m) {
/**
------------
|(0)| |(0)|
---- -----
| a) |
| |
---- -----
|(0)| |(0)|
yL------------
xL
|---| d
|--| pitch
Observation: the area in a) is much larger than the
area(s) (0). In a) particles are lost in (0)
particles they are not lost.
check if we are in the longitudinal
range of the collimator
*/
const double h = pitch_m;
const double xL = - 0.5 * h * (nHolesX_m - 1);
const double yL = - 0.5 * h * (nHolesY_m - 1);
bool alive = false;
for(unsigned int m = 0; (m < nHolesX_m && (!alive)); m++) {
for(unsigned int n = 0; (n < nHolesY_m && (!alive)); n++) {
double x_m = xL + (m * h);
double y_m = yL + (n * h);
/** are we in a) ? */
double rr = std::pow((R(0) - x_m) / rHole_m, 2) + std::pow((R(1) - y_m) / rHole_m, 2);
alive = (rr < 1.0);
}
}
return !alive;
} else if(isASlit_m) {
return (R(0) <= -getXsize() || R(1) <= -getYsize() || R(0) >= getXsize() || R(1) >= getYsize());
} else if (isARColl_m) {
return (R(0) <= -getXsize() || R(1) <= -getYsize() || R(0) >= getXsize() || R(1) >= getYsize());
} else if(isAWire_m) {
ERRORMSG("Not yet implemented");
} else {
// case of an elliptic collimator
const double trm1 = ((R(0)*R(0))/(getXsize()*getXsize()));
const double trm2 = ((R(1)*R(1))/(getYsize()*getYsize()));
return (trm1 + trm2) > 1.0;
const double z = R(2) + P(2) * recpgamma;
if((z > position_m) && (z <= position_m + getElementLength())) {
if(isAPepperPot_m) {
/**
------------
|(0)| |(0)|
---- -----
| a) |
| |
---- -----
|(0)| |(0)|
yL------------
xL
|---| d
|--| pitch
Observation: the area in a) is much larger than the
area(s) (0). In a) particles are lost in (0)
particles they are not lost.
*/
const double h = pitch_m;
const double xL = - 0.5 * h * (nHolesX_m - 1);
const double yL = - 0.5 * h * (nHolesY_m - 1);
bool alive = false;
for(unsigned int m = 0; (m < nHolesX_m && (!alive)); m++) {
for(unsigned int n = 0; (n < nHolesY_m && (!alive)); n++) {
double x_m = xL + (m * h);
double y_m = yL + (n * h);
/** are we in a) ? */
double rr = std::pow((R(0) - x_m) / rHole_m, 2) + std::pow((R(1) - y_m) / rHole_m, 2);
alive = (rr < 1.0);
}
}
return !alive;
} else if(isASlit_m) {
return (R(0) <= -getXsize() || R(1) <= -getYsize() || R(0) >= getXsize() || R(1) >= getYsize());
} else if (isARColl_m) {
return (R(0) <= -getXsize() || R(1) <= -getYsize() || R(0) >= getXsize() || R(1) >= getYsize());
} else if(isAWire_m) {
ERRORMSG("Not yet implemented");
} else {
// case of an elliptic collimator
const double trm1 = ((R(0)*R(0))/(getXsize()*getXsize()));
const double trm2 = ((R(1)*R(1))/(getYsize()*getYsize()));
return (trm1 + trm2) > 1.0;
}
}
}
return false;
return false;
}
bool Collimator::apply(const size_t &i, const double &t, double E[], double B[]) {
......@@ -237,9 +237,9 @@ bool Collimator::apply(const size_t &i, const double &t, Vector_t &E, Vector_t &
pdead = isInColl(R,P,recpgamma);
if(pdead) {
double frac = (R(2) - position_m) / P(2) * recpgamma;
lossDs_m->addParticle(R,P, RefPartBunch_m->ID[i], t + frac * RefPartBunch_m->getdT(), 0);
++losses_m;
double frac = (R(2) - position_m) / P(2) * recpgamma;
lossDs_m->addParticle(R,P, RefPartBunch_m->ID[i], t + frac * RefPartBunch_m->getdT(), 0);
++losses_m;
}
return pdead;
}
......@@ -250,19 +250,19 @@ bool Collimator::apply(const Vector_t &R, const Vector_t &centroid, const double
bool Collimator::checkCollimator(Vector_t r, Vector_t rmin, Vector_t rmax) {
double r_start = sqrt(xstart_m * xstart_m + ystart_m * ystart_m);
double r_end = sqrt(xend_m * xend_m + yend_m * yend_m);
double r1 = sqrt(rmax(0) * rmax(0) + rmax(1) * rmax(1));
bool isDead = false;
if(rmax(2) >= zstart_m && rmin(2) <= zend_m) {
if( r1 > r_start - 10.0 && r1 < r_end + 10.0 ){
if(r(2) < zend_m && r(2) > zstart_m ) {
int pflag = checkPoint(r(0), r(1));
isDead = (pflag != 0);
}
double r_start = sqrt(xstart_m * xstart_m + ystart_m * ystart_m);
double r_end = sqrt(xend_m * xend_m + yend_m * yend_m);
double r1 = sqrt(rmax(0) * rmax(0) + rmax(1) * rmax(1));
bool isDead = false;
if(rmax(2) >= zstart_m && rmin(2) <= zend_m) {
if( r1 > r_start - 10.0 && r1 < r_end + 10.0 ){
if(r(2) < zend_m && r(2) > zstart_m ) {
int pflag = checkPoint(r(0), r(1));
isDead = (pflag != 0);
}
}
}
}
return isDead;
return isDead;
}
......@@ -279,25 +279,25 @@ bool Collimator::checkCollimator(PartBunch &bunch, const int turnnumber, const d
double r1 = sqrt(rmax(0) * rmax(0) + rmax(1) * rmax(1));
if(rmax(2) >= zstart_m && rmin(2) <= zend_m) {
if( r1 > r_start - 10.0 && r1 < r_end + 10.0 ){
size_t tempnum = bunch.getLocalNum();
int pflag = 0;
for(unsigned int i = 0; i < tempnum; ++i) {
if(bunch.PType[i] == 0 && bunch.R[i](2) < zend_m && bunch.R[i](2) > zstart_m ) {
pflag = checkPoint(bunch.R[i](0), bunch.R[i](1));
if(pflag != 0) {
if (!sphys_m)
lossDs_m->addParticle(bunch.R[i], bunch.P[i], bunch.ID[i]);
bunch.Bin[i] = -1;
flagNeedUpdate = true;
}
}
}
}
if( r1 > r_start - 10.0 && r1 < r_end + 10.0 ){
size_t tempnum = bunch.getLocalNum();
int pflag = 0;
for(unsigned int i = 0; i < tempnum; ++i) {
if(bunch.PType[i] == ParticleType::REGULAR && bunch.R[i](2) < zend_m && bunch.R[i](2) > zstart_m ) {
pflag = checkPoint(bunch.R[i](0), bunch.R[i](1));
if(pflag != 0) {
if (!sphys_m)
lossDs_m->addParticle(bunch.R[i], bunch.P[i], bunch.ID[i]);
bunch.Bin[i] = -1;
flagNeedUpdate = true;
}
}
}
}
}
reduce(&flagNeedUpdate, &flagNeedUpdate + 1, &flagNeedUpdate, OpBitwiseOrAssign());
if (flagNeedUpdate && sphys_m) {
sphys_m->apply(bunch);
sphys_m->apply(bunch);
}
return flagNeedUpdate;
}
......@@ -310,10 +310,10 @@ void Collimator::initialise(PartBunch *bunch, double &startField, double &endFie
sphys_m = getSurfacePhysics();
if (!sphys_m) {
if (filename_m == std::string(""))
lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(getName(), !Options::asciidump));
else
lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(filename_m.substr(0, filename_m.rfind(".")), !Options::asciidump));
if (filename_m == std::string(""))
lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(getName(), !Options::asciidump));
else
lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(filename_m.substr(0, filename_m.rfind(".")), !Options::asciidump));
}
goOnline(-1e6);
......@@ -325,10 +325,10 @@ void Collimator::initialise(PartBunch *bunch, const double &scaleFactor) {
sphys_m = getSurfacePhysics();
if (!sphys_m) {
if (filename_m == std::string(""))
lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(getName(), !Options::asciidump));
else
lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(filename_m.substr(0, filename_m.rfind(".")), !Options::asciidump));
if (filename_m == std::string(""))
lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(getName(), !Options::asciidump));
else
lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(filename_m.substr(0, filename_m.rfind(".")), !Options::asciidump));
}
goOnline(-1e6);
......@@ -336,9 +336,9 @@ void Collimator::initialise(PartBunch *bunch, const double &scaleFactor) {
void Collimator::finalise()
{
if(online_m)
goOffline();
*gmsg << "* Finalize probe" << endl;
if(online_m)
goOffline();
*gmsg << "* Finalize probe" << endl;
}
void Collimator::goOnline(const double &) {
......@@ -392,7 +392,7 @@ void Collimator::print() {
if(!informed_m) {
std::string errormsg = Fieldmap::typeset_msg("BUNCH SIZE NOT SET", "warning");
*gmsg << errormsg << "\n"
<< endl;
<< endl;
if(Ippl::myNode() == 0) {
ofstream omsg("errormsg.txt", ios_base::app);
omsg << errormsg << endl;
......@@ -419,9 +419,9 @@ void Collimator::print() {
}
void Collimator::goOffline() {
if (online_m && lossDs_m)
lossDs_m->save();
online_m = false;
if (online_m && lossDs_m)
lossDs_m->save();
online_m = false;
}
bool Collimator::bends() const {
......@@ -597,9 +597,9 @@ void Collimator::setGeom() {
double slope;
if (xend_m == xstart_m)
slope = 1.0e12;
slope = 1.0e12;
else
slope = (yend_m - ystart_m) / (xend_m - xstart_m);
slope = (yend_m - ystart_m) / (xend_m - xstart_m);
double coeff2 = sqrt(1 + slope * slope);
double coeff1 = slope / coeff2;
......@@ -620,10 +620,10 @@ void Collimator::setGeom() {
geom_m[4].y = geom_m[0].y;
if (zstart_m > zend_m){
double tempz = 0.0;
tempz = zstart_m;
zstart_m = zend_m;
zend_m = tempz;
double tempz = 0.0;
tempz = zstart_m;
zstart_m = zend_m;
zend_m = tempz;
}
}
......@@ -641,4 +641,4 @@ int Collimator::checkPoint(const double &x, const double &y) {
}
}
return (cn & 1); // 0 if even (out), and 1 if odd (in)
}
\ No newline at end of file
}
classic/5.0/src/AbsBeamline/Stripper.cpp
View file @ 4d4b095e
......@@ -298,7 +298,7 @@ bool Stripper::checkStripper(PartBunch &bunch, const int turnnumber, const doub
setGeom(Swidth);
for(unsigned int i = 0; i < tempnum; ++i) {
if(bunch.PType[i] == 0) {
if(bunch.PType[i] == ParticleType::REGULAR) {
pflag = checkPoint(bunch.R[i](0), bunch.R[i](1));
if(pflag != 0) {
// dist1 > 0, right hand, dt > 0; dist1 < 0, left hand, dt < 0
......@@ -341,7 +341,7 @@ bool Stripper::checkStripper(PartBunch &bunch, const int turnnumber, const doub
// change the mass and charge
bunch.M[i] = opmass_m;
bunch.Q[i] = opcharge_m * q_e;
bunch.PType[i] = 1;
bunch.PType[i] = ParticleType::STRIPPED;
int j = 1;
//create new particles
......@@ -352,7 +352,7 @@ bool Stripper::checkStripper(PartBunch &bunch, const int turnnumber, const doub
bunch.Q[tempnum+count] = bunch.Q[i];
bunch.M[tempnum+count] = bunch.M[i];
// once the particle is stripped, change PType from 0 to 1 as a flag so as to avoid repetitive stripping.
bunch.PType[tempnum+count] = 1;
bunch.PType[tempnum+count] = ParticleType::STRIPPED;
count++;
j++;
}
......
classic/5.0/src/Algorithms/PBunchDefs.h
View file @ 4d4b095e
......@@ -44,4 +44,13 @@ typedef Field<dcomplex, 3, Mesh_t, Center_t> CxField_t;
typedef FFT<RCTransform, 3, double> FFT_t;
typedef FFT<SineTransform, 3, double> SINE_t;
namespace ParticleType {
enum { REGULAR,
FIELDEMISSION,
SECONDARY,
NEWSECONDARY,
STRIPPED,
PROBE};
}
#endif
classic/5.0/src/Algorithms/PartBunch.cpp
View file @ 4d4b095e
......@@ -2852,12 +2852,12 @@ void PartBunch::pop() {
double PartBunch::getZPos() {
if(sum(PType != 0)) {
if(sum(PType != ParticleType::REGULAR)) {
size_t numberOfPrimaryParticles = 0;
double zAverage = 0.0;
if(getLocalNum() != 0) {
for(size_t partIndex = 0; partIndex < getLocalNum(); partIndex++) {
if(PType[partIndex] == 0) {
if(PType[partIndex] == ParticleType::REGULAR) {
zAverage += X[partIndex](2);
numberOfPrimaryParticles++;
}
......@@ -2964,4 +2964,4 @@ bool PartBunch::WeHaveEnergyBins() {
Vector_t PartBunch::get_pmean_Distribution() const {
return dist_m->get_pmean();
}
\ No newline at end of file
}
classic/5.0/src/Algorithms/PartBunch.h
View file @ 4d4b095e
......@@ -87,7 +87,7 @@ public:
void makHistograms();
/** \brief After each Schottky scan we delete
all the particles.
all the particles.
*/
void cleanUpParticles();
......@@ -124,7 +124,7 @@ public:
/*
Energy bins related functions
Energy bins related functions
*/
......@@ -187,7 +187,7 @@ public:
/*
Mesh and Field Layout related functions
Mesh and Field Layout related functions
*/
......@@ -523,10 +523,10 @@ private:
double calculateAngle(double x, double y);
double calculateAngle2(double x, double y);
public:
public:
void calcEMean(); // update eKin_m;
void correctEnergy(double avrgp);
private:
private:
/*
Member variables starts here
......@@ -689,7 +689,7 @@ private:
Distribution *dist_m;
// flag to tell if we are a DC-beam
// flag to tell if we are a DC-beam
bool dcBeam_m;
......@@ -1012,22 +1012,22 @@ double PartBunch::getGaBeLa() const { // used in Distribution
inline
double PartBunch::get_sPos() {
if(sum(PType != 0)) {
if(sum(PType != ParticleType::REGULAR)) {
const size_t n = getLocalNum();
size_t numPrimBeamParts = 0;
double z = 0.0;
if(n != 0) {
for(size_t i = 0; i < n; i++) {
if(PType[i] == 0) {
z += R[i](2);
numPrimBeamParts++;
}
}
for(size_t i = 0; i < n; i++) {
if(PType[i] == ParticleType::REGULAR) {
z += R[i](2);
numPrimBeamParts++;
}
}
}
reduce(z, z, OpAddAssign());
reduce(numPrimBeamParts, numPrimBeamParts, OpAddAssign());
if(numPrimBeamParts != 0)
z = z / numPrimBeamParts;
z = z / numPrimBeamParts;
return z;
} else {
const size_t n = getTotalNum();
......@@ -1151,10 +1151,10 @@ size_t PartBunch::getLoadBalance(int p) const {
inline
size_t PartBunch::getMinLocalNum() {
/// Get the minimal number of particles per node
if (minLocNum_m < 0)
gatherLoadBalanceStatistics();
return minLocNum_m;
/// Get the minimal number of particles per node
if (minLocNum_m < 0)
gatherLoadBalanceStatistics();
return minLocNum_m;
}
inline
......
src/Algorithms/AutophaseTracker.cpp
View file @ 4d4b095e
......@@ -31,7 +31,7 @@ AutophaseTracker::AutophaseTracker(const Beamline &beamline,
itsBunch_m.P[0] = Vector_t(0.0, 0.0, initialP);
itsBunch_m.Bin[0] = 0;
itsBunch_m.Q[0] = itsBunch_m.getChargePerParticle();
itsBunch_m.PType[0] = 0;
itsBunch_m.PType[0] = ParticleType::REGULAR;
itsBunch_m.LastSection[0] = 0;
}
......@@ -541,4 +541,4 @@ void AutophaseTracker::save(const std::string &fname) {
}
out.close();
}
\ No newline at end of file
}
src/Algorithms/ParallelCyclotronTracker.cpp
View file @ 4d4b095e
......@@ -4171,7 +4171,7 @@ bool ParallelCyclotronTracker::readOneBunchFromFile(const size_t BinID) {
itsBunch->P[localNum] = tmpBunch.P[ii];
itsBunch->M[localNum] = tmpBunch.M[ii];
itsBunch->Q[localNum] = tmpBunch.Q[ii];
itsBunch->PType[localNum] = 0;
itsBunch->PType[localNum] = ParticleType::REGULAR;
itsBunch->Bin[localNum] = BinID;
r_mb[ii] = itsBunch->R[localNum];
......
src/Algorithms/ParallelTTracker.cpp
View file @ 4d4b095e
......@@ -459,10 +459,10 @@ void ParallelTTracker::executeDefaultTracker() {
//allocate memory on device for particle
allocateDeviceMemory();
//page lock itsBunch->X, itsBunch->R, itsBunch-P
registerHostMemory();
//write R, P and X data to device
dksbase.writeDataAsync<Vector_t>(r_ptr, &itsBunch->R[0], itsBunch->getLocalNum());
dksbase.writeDataAsync<Vector_t>(p_ptr, &itsBunch->P[0], itsBunch->getLocalNum());
......@@ -552,7 +552,7 @@ void ParallelTTracker::executeDefaultTracker() {
//unregister page lock itsBunch->X, itsBunch->R, itsBunch-P
unregisterHostMemory();
#endif
}
......@@ -1439,9 +1439,8 @@ void ParallelTTracker::bgf_main_collision_test() {
double seyNum = 0;
if(secondaryFlg_m != 1) {
// mark secondaries generated in last step as old
if(itsBunch->PType[i] == 3)
itsBunch->PType[i] = 2;
if(itsBunch->PType[i] == ParticleType::NEWSECONDARY)
itsBunch->PType[i] = ParticleType::SECONDARY;
}
/*
......@@ -1621,7 +1620,7 @@ void ParallelTTracker::timeIntegration1_bgf(BorisPusher & pusher) {
Vector_t intecoords = outr;
int triId = 0;
if(itsBunch->PType[i] == 3) { // only test newly generated secondaries
if(itsBunch->PType[i] == ParticleType::NEWSECONDARY) {
particleHitBoundary = bgf_m->PartInside(itsBunch->R[i], itsBunch->P[i], 0.5 * itsBunch->dt[i], itsBunch->PType[i], itsBunch->Q[i], intecoords, triId) == 0;
}
......@@ -1686,7 +1685,7 @@ void ParallelTTracker::timeIntegration2(BorisPusher & pusher) {
IpplTimings::startTimer(timeIntegrationTimer2Loop1_m);
#ifdef OPAL_DKS
//check if number of particles in bunch has changed then write
if (surfaceStatus_m || itsBunch->getLocalNum() != numDeviceElements) {
......@@ -1697,7 +1696,7 @@ void ParallelTTracker::timeIntegration2(BorisPusher & pusher) {
allocateDeviceMemory();
registerHostMemory();
}
//write R to device
dksbase.writeDataAsync<Vector_t>(r_ptr, &itsBunch->R[0], itsBunch->getLocalNum(), stream1);
//write P to device
......@@ -1724,7 +1723,7 @@ void ParallelTTracker::timeIntegration2(BorisPusher & pusher) {
dt_ptr, itsBunch->getdT(), Physics::c, true, stream2);
//read R from device
dksbase.readDataAsync<Vector_t>(x_ptr, &itsBunch->X[0], itsBunch->getLocalNum(), stream2);
//sync and wait till all tasks and reads are complete
dksbase.syncDevice();
......@@ -2819,7 +2818,7 @@ void ParallelTTracker::unregisterHostMemory() {
void ParallelTTracker::allocateDeviceMemory() {
//allocate memory on device
//allocate memory on device
r_ptr = dksbase.allocateMemory<Vector_t>(itsBunch->getLocalNum(), ierr);
p_ptr = dksbase.allocateMemory<Vector_t>(itsBunch->getLocalNum(), ierr);
x_ptr = dksbase.allocateMemory<Vector_t>(itsBunch->getLocalNum(), ierr);
......
src/Distribution/Distribution.cpp
View file @ 4d4b095e
......@@ -527,7 +527,7 @@ void Distribution::CreatePriPart(PartBunch *beam, BoundaryGeometry &bg) {
beam->P[lowMark + count] = betagamma * bg.getMomenta(count);
}
beam->Bin[lowMark + count] = 0;
beam->PType[lowMark + count] = 0; // create primary particle bunch;
beam->PType[lowMark + count] = ParticleType::REGULAR;
beam->TriID[lowMark + count] = 0;
beam->Q[lowMark + count] = beam->getChargePerParticle();
beam->LastSection[lowMark + count] = 0;
......@@ -569,7 +569,7 @@ void Distribution::CreatePriPart(PartBunch *beam, BoundaryGeometry &bg) {
beam->R[lowMark + count] = bg.getCooridinate(count); // for node0 the particle number N_mean = N_mean + N_extra
beam->P[lowMark + count] = Vector_t(0.0);
beam->Bin[lowMark + count] = 0;
beam->PType[lowMark + count] = 0; // create primary particle bunch.
beam->PType[lowMark + count] = ParticleType::REGULAR;
beam->TriID[lowMark + count] = 0;
beam->Q[lowMark + count] = beam->getChargePerParticle();
beam->LastSection[lowMark + count] = 0;
......@@ -1428,7 +1428,7 @@ void Distribution::CreateBoundaryGeometry(PartBunch &beam, BoundaryGeometry &bg
beam.R[lowMark + count] = bg.getCooridinate(count);
beam.P[lowMark + count] = Vector_t(0.0);
beam.Bin[lowMark + count] = 0;
beam.PType[lowMark + count] = 1; // create darkcurrent particles.
beam.PType[lowMark + count] = ParticleType::FIELDEMISSION;