Commit 848b22ba authored by kraus's avatar kraus
Browse files

- remove attribute DISTRIBUTIONS and instead make attribute DISTRIBUTION an array of strings

- fix output of number of particles (suppress warning about number of particles on Beam command where there shouldn't be one)
parent ce7458c4
......@@ -396,7 +396,9 @@ void PartBunch::setDistribution(Distribution *d,
bool scan) {
Inform m("setDistribution " );
dist_m = d;
dist_m->CreateOpalT(*this, addedDistributions, np, scan);
// if (Options::cZero)
// dist_m->Create(*this, addedDistributions, np / 2, scan);
// else
......@@ -769,18 +771,18 @@ void PartBunch::computeSelfFields(int binNumber) {
int mz2 = (int)nr_m[2] / 2;
for (int i=0; i<mx; i++ )
*gmsg << "Bin " << binNumber
<< ", Self Field along x axis E = " << eg_m[i][my2][mz2]
*gmsg << "Bin " << binNumber
<< ", Self Field along x axis E = " << eg_m[i][my2][mz2]
<< ", Pot = " << rho_m[i][my2][mz2] << endl;
for (int i=0; i<my; i++ )
*gmsg << "Bin " << binNumber
<< ", Self Field along y axis E = " << eg_m[mx2][i][mz2]
*gmsg << "Bin " << binNumber
<< ", Self Field along y axis E = " << eg_m[mx2][i][mz2]
<< ", Pot = " << rho_m[mx2][i][mz2] << endl;
for (int i=0; i<mz; i++ )
*gmsg << "Bin " << binNumber
<< ", Self Field along z axis E = " << eg_m[mx2][my2][i]
*gmsg << "Bin " << binNumber
<< ", Self Field along z axis E = " << eg_m[mx2][my2][i]
<< ", Pot = " << rho_m[mx2][my2][i] << endl;
#endif
......@@ -847,18 +849,18 @@ void PartBunch::computeSelfFields(int binNumber) {
//int mz2 = (int)nr_m[2] / 2;
for (int i=0; i<mx; i++ )
*gmsg << "Bin " << binNumber
<< ", Image Field along x axis E = " << eg_m[i][my2][mz2]
*gmsg << "Bin " << binNumber
<< ", Image Field along x axis E = " << eg_m[i][my2][mz2]
<< ", Pot = " << rho_m[i][my2][mz2] << endl;
for (int i=0; i<my; i++ )
*gmsg << "Bin " << binNumber
<< ", Image Field along y axis E = " << eg_m[mx2][i][mz2]
*gmsg << "Bin " << binNumber
<< ", Image Field along y axis E = " << eg_m[mx2][i][mz2]
<< ", Pot = " << rho_m[mx2][i][mz2] << endl;
for (int i=0; i<mz; i++ )
*gmsg << "Bin " << binNumber
<< ", Image Field along z axis E = " << eg_m[mx2][my2][i]
*gmsg << "Bin " << binNumber
<< ", Image Field along z axis E = " << eg_m[mx2][my2][i]
<< ", Pot = " << rho_m[mx2][my2][i] << endl;
#endif
......@@ -1505,18 +1507,18 @@ void PartBunch::computeSelfFields_cycl(int bin) {
int mz2 = (int)nr_m[2] / 2;
for (int i=0; i<mx; i++ )
*gmsg << "Bin " << bin
<< ", Field along x axis Ex = " << eg_m[i][my2][mz2]
*gmsg << "Bin " << bin
<< ", Field along x axis Ex = " << eg_m[i][my2][mz2]
<< ", Pot = " << rho_m[i][my2][mz2] << endl;
for (int i=0; i<my; i++ )
*gmsg << "Bin " << bin
<< ", Field along y axis Ey = " << eg_m[mx2][i][mz2]
*gmsg << "Bin " << bin
<< ", Field along y axis Ey = " << eg_m[mx2][i][mz2]
<< ", Pot = " << rho_m[mx2][i][mz2] << endl;
for (int i=0; i<mz; i++ )
*gmsg << "Bin " << bin
<< ", Field along z axis Ez = " << eg_m[mx2][my2][i]
*gmsg << "Bin " << bin
<< ", Field along z axis Ez = " << eg_m[mx2][my2][i]
<< ", Pot = " << rho_m[mx2][my2][i] << endl;
#endif
......@@ -1870,7 +1872,7 @@ void PartBunch::beamEllipsoid(FVector<double, 6> &centroid,
void PartBunch::gatherLoadBalanceStatistics() {
minLocNum_m = std::numeric_limits<size_t>::max();
for(int i = 0; i < Ippl::getNodes(); i++)
partPerNode_m[i] = globalPartPerNode_m[i] = 0.0;
......@@ -2467,7 +2469,7 @@ size_t PartBunch::boundp_destroyT() {
if((Bin[i] < 0) && ((this->getLocalNum()-i)>1)) { // need in minimum 1 particle per node
ne++;
this->destroy(1, i);
}
}
}
lowParticleCount_m = ((this->getLocalNum()-ne) <= 1);
reduce(lowParticleCount_m, lowParticleCount_m, OpOr());
......@@ -2940,4 +2942,4 @@ bool PartBunch::WeHaveEnergyBins() {
return dist_m->GetNumberOfEnergyBins() > 0;
else
return false;
}
}
\ No newline at end of file
......@@ -196,6 +196,7 @@ Distribution::Distribution():
Definition(AttributesT::SIZE + LegacyAttributesT::SIZE, "DISTRIBUTION",
"The DISTRIBUTION statement defines data for the 6D particle distribution."),
distrTypeT_m(DistrTypeT::NODIST),
numberOfDistributions_m(1),
emitting_m(false),
scan_m(false),
emissionModel_m(EmissionModelT::NONE),
......@@ -269,6 +270,7 @@ Distribution::Distribution(const std::string &name, Distribution *parent):
Definition(name, parent),
distT_m(parent->distT_m),
distrTypeT_m(DistrTypeT::NODIST),
numberOfDistributions_m(parent->numberOfDistributions_m),
emitting_m(parent->emitting_m),
scan_m(parent->scan_m),
particleRefData_m(parent->particleRefData_m),
......@@ -1471,27 +1473,22 @@ void Distribution::ApplyEmissModelNonEquil(double eZ,
void Distribution::CalcPartPerDist(size_t numberOfParticles) {
double totalWeight = GetWeight();
size_t numberOfDistributions = 1;
std::vector<Distribution *>::iterator distributionIt;
for (distributionIt = addedDistributions_m.begin();
distributionIt != addedDistributions_m.end(); distributionIt++) {
totalWeight += (*distributionIt)->GetWeight();
numberOfDistributions++;
}
typedef std::vector<Distribution *>::iterator iterator;
if (numberOfDistributions == 1)
if (numberOfDistributions_m == 1)
particlesPerDist_m.push_back(numberOfParticles);
else {
double totalWeight = GetWeight();
for (iterator it = addedDistributions_m.begin(); it != addedDistributions_m.end(); it++) {
totalWeight += (*it)->GetWeight();
}
particlesPerDist_m.push_back(0);
size_t numberOfCommittedPart = 0;
for (distributionIt = addedDistributions_m.begin();
distributionIt != addedDistributions_m.end(); distributionIt++) {
particlesPerDist_m.push_back(numberOfParticles
* (*distributionIt)->GetWeight()
/ totalWeight);
numberOfCommittedPart += particlesPerDist_m.back();
for (iterator it = addedDistributions_m.begin(); it != addedDistributions_m.end(); it++) {
size_t particlesCurrentDist = numberOfParticles * (*it)->GetWeight() / totalWeight;
particlesPerDist_m.push_back(particlesCurrentDist);
numberOfCommittedPart += particlesCurrentDist;
}
// Remaining particles go into main distribution.
......@@ -1545,24 +1542,25 @@ void Distribution::CheckIfEmitted() {
void Distribution::CheckParticleNumber(size_t &numberOfParticles) {
size_t numberOfDistParticles = tOrZDist_m.size();
reduce(numberOfDistParticles, numberOfDistParticles, OpAddAssign());
if (numberOfDistParticles != numberOfParticles) {
size_t numberOfDistParticles[] = {tOrZDist_m.size(), numberOfParticles};
reduce(numberOfDistParticles, numberOfDistParticles + 2, numberOfDistParticles, OpAddAssign());
if (numberOfDistParticles[0] != numberOfDistParticles[1]) {
*gmsg << "\n--------------------------------------------------" << endl
<< "Warning!! The number of particles in the initial" << endl
<< "distribution is " << numberOfDistParticles << "." << endl << endl
<< "distribution is " << numberOfDistParticles[0] << "." << endl << endl
<< "This is different from the number of particles" << endl
<< "defined by the BEAM command: " << numberOfParticles << endl << endl
<< "defined by the BEAM command: " << numberOfDistParticles[1] << endl << endl
<< "This is often happens when using a FROMFILE type" << endl
<< "distribution and not matching the number of" << endl
<< "particles in the particles file(s) with the number" << endl
<< "given in the BEAM command." << endl << endl
<< "The number of particles in the initial distribution" << endl
<< "(" << numberOfDistParticles << ") "
<< "(" << numberOfDistParticles[0] << ") "
<< "will take precedence." << endl
<< "---------------------------------------------------\n" << endl;
}
numberOfParticles = numberOfDistParticles;
numberOfParticles = numberOfDistParticles[0];
}
void Distribution::ChooseInputMomentumUnits(InputMomentumUnitsT::InputMomentumUnitsT inputMoUnits) {
......@@ -1711,9 +1709,9 @@ void Distribution::CreateDistributionFromFile(size_t numberOfParticles, double m
void Distribution::CreateMatchedGaussDistribution(size_t numberOfParticles, double massIneV) {
/*
/*
ToDo:
- add flag in order to calculate tunes from FMLOWE to FMHIGHE
- add flag in order to calculate tunes from FMLOWE to FMHIGHE
- eliminate physics and error
*/
......@@ -1724,7 +1722,7 @@ void Distribution::CreateMatchedGaussDistribution(size_t numberOfParticles, doub
SpecificElementVisitor<Cyclotron> CyclotronVisitor(*LineSequence->fetchLine());
CyclotronVisitor.execute();
size_t NumberOfCyclotrons = CyclotronVisitor.size();
if (NumberOfCyclotrons > 1) {
throw OpalException("Distribution::CreateMatchedGaussDistribution",
"I am confused, found more than one Cyclotron element in line");
......@@ -1734,7 +1732,7 @@ void Distribution::CreateMatchedGaussDistribution(size_t numberOfParticles, doub
"didn't find any Cyclotron element in line");
}
const Cyclotron* CyclotronElement = dynamic_cast<const Cyclotron*>(CyclotronVisitor.front());
*gmsg << "----------------------------------------------------" << endl;
*gmsg << "About to find closed orbit and matched distribution " << endl;
*gmsg << "I= " << I_m*1E3 << " (mA) E= " << E_m*1E-6 << " (MeV)" << endl;
......@@ -1755,7 +1753,7 @@ void Distribution::CreateMatchedGaussDistribution(size_t numberOfParticles, doub
ERRORMSG("FMHIGHE of FMLOWE not set propperly" << endl);
exit(1);
}
int nr, nth, nsc;
double rmin, dr, dth;
......@@ -1787,7 +1785,7 @@ void Distribution::CreateMatchedGaussDistribution(size_t numberOfParticles, doub
}
INFOMSG(endl);
}
/*
Now setup the distribution generator
......@@ -1801,11 +1799,11 @@ void Distribution::CreateMatchedGaussDistribution(size_t numberOfParticles, doub
sigmaR_m[0] = std::sqrt(1E-3*siggen.getSigma()(0,0)); // rms x
sigmaR_m[1] = std::sqrt(1E-3*siggen.getSigma()(2,2)); // rms y
sigmaR_m[2] = std::sqrt(1E-3*siggen.getSigma()(4,4)); // rms z
sigmaP_m[0] = std::sqrt(1E-3*siggen.getSigma()(1,1)); // UNITS ...
sigmaP_m[1] = std::sqrt(1E-3*siggen.getSigma()(3,3));
sigmaP_m[2] = std::sqrt(1E-3*siggen.getSigma()(5,5));
distCorr_m.push_back(1E-3*siggen.getSigma()(0,1)); // R12
distCorr_m.push_back(1E-3*siggen.getSigma()(2,3)); // R34
distCorr_m.push_back(1E-3*siggen.getSigma()(4,5)); // R56
......@@ -1813,11 +1811,11 @@ void Distribution::CreateMatchedGaussDistribution(size_t numberOfParticles, doub
distCorr_m.push_back(1E-3*siggen.getSigma()(1,5)); // R26
distCorr_m.push_back(1E-3*siggen.getSigma()(0,4)); // R15
distCorr_m.push_back(1E-3*siggen.getSigma()(1,4)); // R25
} else {
*gmsg << "Not converged." << endl;
}
CreateDistributionGauss(numberOfParticles, massIneV);
}
}
......@@ -1835,7 +1833,7 @@ void Distribution::CreateDistributionGauss(size_t numberOfParticles, double mass
GenerateTransverseGauss(numberOfParticles);
GenerateLongFlattopT(numberOfParticles);
} else {
GenerateGaussZ(numberOfParticles);
GenerateGaussZ(numberOfParticles);
}
}
......@@ -3075,6 +3073,7 @@ void Distribution::GenerateLongFlattopT(size_t numberOfParticles) {
gsl_rng_env_setup();
gsl_rng *randGenGSL = gsl_rng_alloc(gsl_rng_default);
int saveProcessor = -1;
for (size_t partIndex = 0; partIndex < numFall; partIndex++) {
double t = 0.0;
......@@ -3477,9 +3476,9 @@ void Distribution::PrintDist(Inform &os, size_t numberOfParticles) const {
if (numberOfParticles > 0) {
os << "Number of particles: "
<< numberOfParticles
<< endl
<< endl;
<< numberOfParticles * (Options::cZero ? 2: 1)
<< endl
<< endl;
}
switch (distrTypeT_m) {
......@@ -5053,4 +5052,4 @@ void Distribution::WriteOutFileInjection() {
reduce(numberOfParticles, numberOfParticles, OpAddAssign());
}
}
}
}
\ No newline at end of file
......@@ -193,6 +193,7 @@ public:
double GetPsi() {return psi_m;}
bool GetPreviousH5Local() {return previousH5Local_m;}
void setNumberOfDistributions(unsigned int n) { numberOfDistributions_m = n; }
private:
Distribution(const std::string &name, Distribution *parent);
......@@ -275,6 +276,8 @@ private:
std::ofstream os_m; /// Output file to write distribution.
void writeToFileCycl(PartBunch &beam, size_t Np);
unsigned int numberOfDistributions_m;
bool emitting_m; /// Distribution is an emitted, and is currently
/// emitting, rather than an injected, beam.
......@@ -338,7 +341,7 @@ private:
std::vector<size_t> binWrite_m;
// for compatibility reasons
double avrgpz_m;
double avrgpz_m;
......@@ -427,7 +430,7 @@ private:
double referenceTheta_m;
double referenceZ_m;
double bega_m;
// Cyclotron for restart in local mode
double phi_m;
double psi_m;
......@@ -439,4 +442,4 @@ inline Inform &operator<<(Inform &os, const Distribution &d) {
return d.printInfo(os);
}
#endif // OPAL_Distribution_HH
#endif // OPAL_Distribution_HH
\ No newline at end of file
......@@ -74,7 +74,6 @@ namespace {
FIELDSOLVER, // The field solver attached
BOUNDARYGEOMETRY, // The boundary geometry
DISTRIBUTION, // The particle distribution
DISTRIBUTIONS, // A list of particle distributions
MULTIPACTING, // MULTIPACTING flag
// THE INTEGRATION TIMESTEP IN SEC
SIZE
......@@ -106,10 +105,8 @@ TrackRun::TrackRun():
("FIELDSOLVER", "Field solver to be used ", "FIELDSOLVER");
itsAttr[BOUNDARYGEOMETRY] = Attributes::makeString
("BOUNDARYGEOMETRY", "Boundary geometry to be used NONE (default)", "NONE");
itsAttr[DISTRIBUTION] = Attributes::makeString
("DISTRIBUTION", "Particle distribution to be used ", "DISTRIBUTION");
itsAttr[DISTRIBUTIONS] = Attributes::makeStringArray
("DISTRIBUTIONS", "List of particle distributions to be used ");
itsAttr[DISTRIBUTION] = Attributes::makeStringArray
("DISTRIBUTION", "List of particle distributions to be used ");
itsAttr[MULTIPACTING] = Attributes::makeBool
("MULTIPACTING", "Multipacting flag, default: false. Set true to initialize primary particles according to BoundaryGeometry", false);
OPAL = OpalData::getInstance();
......@@ -163,25 +160,31 @@ void TrackRun::execute() {
}
Beam *beam = Beam::find(Attributes::getString(itsAttr[BEAM]));
dist = Distribution::find(Attributes::getString(itsAttr[DISTRIBUTION]));
std::vector<std::string> distr_str = Attributes::getStringArray(itsAttr[DISTRIBUTION]);
const size_t numberOfDistributions = distr_str.size();
dist = Distribution::find(distr_str.at(0));
dist->setNumberOfDistributions(numberOfDistributions);
fs = FieldSolver::find(Attributes::getString(itsAttr[FIELDSOLVER]));
fs->initCartesianFields();
double charge = 0.0;
std::vector<std::string> distr_str = Attributes::getStringArray(itsAttr[DISTRIBUTIONS]);
if(distr_str.size() > 0) {
if(numberOfDistributions > 1) {
*gmsg << "Found more than one distribution: ";
for(std::vector<std::string>::const_iterator dit = distr_str.begin(); dit != distr_str.end(); ++ dit) {
Distribution *d = Distribution::find(*dit);
*gmsg << " " << *dit;
for(size_t i = 1; i < numberOfDistributions; ++ i) {
Distribution *d = Distribution::find(distr_str.at(i));
d->setNumberOfDistributions(numberOfDistributions);
distrs_m.push_back(d);
*gmsg << " " << distr_str.at(i);
}
*gmsg << endl;
}
fs = FieldSolver::find(Attributes::getString(itsAttr[FIELDSOLVER]));
fs->initCartesianFields();
double charge = 0.0;
if(!OPAL->hasSLBunchAllocated()) {
if(!OPAL->inRestartRun()) {
......@@ -457,7 +460,6 @@ void TrackRun::execute() {
// }
// }
#endif
if(!OPAL->inRestartRun()) {
if(!OPAL->hasDataSinkAllocated() && !Options::scan) {
OPAL->setDataSink(new DataSink());
......@@ -532,7 +534,8 @@ void TrackRun::execute() {
Track::block->bunch->PType = 0;
dist = Distribution::find(Attributes::getString(itsAttr[DISTRIBUTION]));
std::vector<std::string> distr_str = Attributes::getStringArray(itsAttr[DISTRIBUTION]);
dist = Distribution::find(distr_str.at(0));
// set macromass and charge for simulation particles
double macromass = 0.0;
......@@ -543,8 +546,8 @@ void TrackRun::execute() {
if(beam->getNumberOfParticles() < 3 || beam->getCurrent() == 0.0) {
macrocharge = beam->getCharge() * q_e;
macromass = beam->getMass();
dist->CreateOpalCycl(*Track::block->bunch,
beam->getNumberOfParticles(),
dist->CreateOpalCycl(*Track::block->bunch,
beam->getNumberOfParticles(),
beam->getCurrent(),*Track::block->use->fetchLine(),
Options::scan);
......@@ -799,34 +802,32 @@ double TrackRun::SetDistributionParallelT(Beam *beam) {
* it to create the full distribution.
*/
std::vector<std::string> distributionArray
= Attributes::getStringArray(itsAttr[DISTRIBUTIONS]);
= Attributes::getStringArray(itsAttr[DISTRIBUTION]);
const size_t numberOfDistributions = distributionArray.size();
if (distributionArray.size() > 0) {
dist = Distribution::find(distributionArray.at(0));
dist->setNumberOfDistributions(numberOfDistributions);
if (numberOfDistributions > 1) {
*gmsg << endl
<< "---------------------------------" << endl
<< "Found more than one distribution:" << endl << endl;
*gmsg << "Main Distribution" << endl
<< "---------------------------------" << endl
<< distributionArray.at(0) << endl << endl
<< "Secondary Distribution(s)" << endl
<< "---------------------------------" << endl;
for (std::vector<std::string>::const_iterator distIterator
= distributionArray.begin();
distIterator != distributionArray.end(); ++distIterator) {
if (distIterator == distributionArray.begin()) {
dist = Distribution::find(*distIterator);
*gmsg << "Main Distribution" << endl
<< "-----------------" << endl
<< *distIterator << endl << endl
<< "Secondary Distribution(s)" << endl
<< "-------------------------" << endl;
} else {
Distribution *distribution = Distribution::find(*distIterator);
*gmsg << *distIterator << endl;
distrs_m.push_back(distribution);
}
for (size_t i = 1; i < numberOfDistributions; ++ i) {
Distribution *distribution = Distribution::find(distributionArray.at(i));
distribution->setNumberOfDistributions(numberOfDistributions);
distrs_m.push_back(distribution);
*gmsg << distributionArray.at(i) << endl;
}
*gmsg << endl
<< "---------------------------------" << endl << endl;
} else
dist = Distribution::find(Attributes::getString(itsAttr[DISTRIBUTION]));
}
}
......@@ -897,7 +898,8 @@ ParallelTTracker *TrackRun::setupForAutophase() {
Track::block->bunch->setSolver(fs);
dist = Distribution::find(Attributes::getString(itsAttr[DISTRIBUTION]));
std::vector<std::string> distr_str = Attributes::getStringArray(itsAttr[DISTRIBUTION]);
dist = Distribution::find(distr_str.at(0));
double charge = beam->getCharge() * beam->getCurrent() / beam->getFrequency();
......@@ -986,4 +988,4 @@ std::pair<Box,unsigned int> TrackRun::getBlGrids(std::string str){
return std::pair<Box,unsigned int>(bx,theGrid);
}
#endif
#endif
\ No newline at end of file
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