Commit 04ed0139 authored by kraus's avatar kraus
Browse files

- Distribution: always allow to use SIGMAR and CUTOFFR when using Gaussian distribution

- Expressions: adding error function
- LossDataSink: we are using git
- OpalBeamline: let Autophase tracker smoothly (i.e. don't throw a segmentation fault) even if there are no cavities
- ParallelTTracker: don't save phase space after last step if saved at end of last step
parent bbdac8f9
......@@ -51,8 +51,6 @@
using Physics::pi;
using namespace std;
extern Inform *gmsg;
// Class PartBunch
......@@ -486,7 +484,7 @@ double PartBunch::getZ(int i) {
*/
void PartBunch::calcLineDensity() {
// e_dim_tag decomp[3];
list<ListElem> listz;
std::list<ListElem> listz;
// for (int d=0; d < 3; ++d) { // this does not seem to work properly
// decomp[d] = getFieldLayout().getRequestedDistribution(d);
......@@ -575,23 +573,23 @@ void PartBunch::calcLineDensity() {
reduce(&(lineDensity_m[0]), &(lineDensity_m[0]) + nBinsLineDensity_m, &(lineDensity_m[0]), OpAddAssign());
}
void PartBunch::fillArray(double *lineDensity, const list<ListElem> &l) {
void PartBunch::fillArray(double *lineDensity, const std::list<ListElem> &l) {
unsigned int mmax = 0;
unsigned int nmax = 0;
unsigned int count = 0;
for(list<ListElem>::const_iterator it = l.begin(); it != l.end() ; ++it) {
for(std::list<ListElem>::const_iterator it = l.begin(); it != l.end() ; ++it) {
if(it->m > mmax) mmax = it->m;
if(it->n > nmax) nmax = it->n;
}
for(list<ListElem>::const_iterator it = l.begin(); it != l.end(); ++it)
for(std::list<ListElem>::const_iterator it = l.begin(); it != l.end(); ++it)
if((it->m < mmax) && (it->n < nmax)) {
lineDensity[count] = it->den;
count++;
}
}
void PartBunch::getLineDensity(vector<double> &lineDensity) {
void PartBunch::getLineDensity(std::vector<double> &lineDensity) {
if(bool(lineDensity_m)) {
if(lineDensity.size() != nBinsLineDensity_m)
lineDensity.resize(nBinsLineDensity_m, 0.0);
......@@ -622,7 +620,7 @@ void PartBunch::calcGammas() {
reduce(pInBin, pInBin, OpAddAssign());
if(pInBin != 0) {
bingamma_m[i] /= pInBin;
INFOMSG(level2 << "Bin " << std::setw(3) << i << " gamma = " << setw(8) << scientific << setprecision(5) << bingamma_m[i] << "; NpInBin= " << setw(8) << setfill(' ') << pInBin << endl);
INFOMSG(level2 << "Bin " << std::setw(3) << i << " gamma = " << std::setw(8) << std::scientific << std::setprecision(5) << bingamma_m[i] << "; NpInBin= " << std::setw(8) << std::setfill(' ') << pInBin << endl);
} else {
bingamma_m[i] = 1.0;
INFOMSG(level2 << "Bin " << std::setw(3) << i << " has no particles " << endl);
......@@ -1963,7 +1961,7 @@ void PartBunch::calcMomentsInitial() {
for(size_t k = 0; k < pbin_m->getNp(); k++) {
for(int binNumber = 0; binNumber < pbin_m->getNBins(); binNumber++) {
vector<double> p;
std::vector<double> p;
if(pbin_m->getPart(k, binNumber, p)) {
part[0] = p.at(0);
......@@ -2130,7 +2128,7 @@ void PartBunch::calcBeamParametersInitial() {
rrms_m(i) = sqrt(rsqsum(i) / N);
prms_m(i) = sqrt(psqsum(i) / N);
eps_m(i) = sqrt(max(eps2(i), zero));
eps_m(i) = sqrt(std::max(eps2(i), zero));
double tmp = rrms_m(i) * prms_m(i);
fac(i) = (tmp == 0) ? zero : 1.0 / tmp;
}
......@@ -2241,7 +2239,7 @@ void PartBunch::moveBunchToCathode(double &t) {
for(int bin = 0; bin < getNumBins(); ++bin) {
for(size_t i = 0; i < pbin_m->getNp(); ++i) {
vector<double> p;
std::vector<double> p;
if(pbin_m->getPart(i, bin, p)) {
avrg_betagamma += sqrt(1.0 + p[3] * p[3] + p[4] * p[4] + p[5] * p[5]);
if(p[2] > maxspos) maxspos = p[2];
......@@ -2258,7 +2256,7 @@ void PartBunch::moveBunchToCathode(double &t) {
gmsg << "move bunch by " << num_steps *dist_per_step << "; DT = " << num_steps *getdT() << endl;
for(int bin = 0; bin < getNumBins(); ++bin) {
for(size_t i = 0; i < pbin_m->getNp(); ++i) {
vector<double> p;
std::vector<double> p;
if(pbin_m->getPart(i, bin, p)) {
pbin_m->updatePartPos(i, bin, p[2] + num_steps * dist_per_step);
}
......@@ -2283,7 +2281,7 @@ void PartBunch::printBinHist() {
if(Ippl::myNode() == 0) {
for(int bin = 0; bin < getNumBins(); ++bin) {
for(size_t i = 0; i < pbin_m->getNp(); ++i) {
vector<double> p;
std::vector<double> p;
if(pbin_m->getPart(i, bin, p)) {
if(p[2] > maxz) maxz = p[2];
if(p[2] < minz) minz = p[2];
......@@ -2294,7 +2292,7 @@ void PartBunch::printBinHist() {
int minihist[20] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
for(int bin = 0; bin < getNumBins(); ++bin) {
for(size_t i = 0; i < pbin_m->getNp(); ++i) {
vector<double> p;
std::vector<double> p;
if(pbin_m->getPart(i, bin, p)) {
++minihist[(int)floor((p[2] - minz) / dz)];
}
......@@ -2317,27 +2315,27 @@ void PartBunch::printBinHist() {
Inform &PartBunch::print(Inform &os) {
if(this->getTotalNum() != 0) { // to suppress Nan's
Inform::FmtFlags_t ff = os.flags();
os << scientific;
os << std::scientific;
os << level1 << "\n";
os << "* ************** B U N C H ********************************************************* \n";
os << "* NP = " << this->getTotalNum() << "\n";
os << "* Qtot = " << setw(12) << setprecision(5) << abs(sum(Q)) * 1.0e9 << " [nC] "
<< "Qi = " << setw(12) << std::abs(qi_m) * 1e9 << " [nC]" << "\n";
os << "* Ekin = " << setw(12) << setprecision(5) << eKin_m << " [MeV] "
<< "dEkin = " << setw(12) << dE_m << " [MeV]\n";
os << "* rmax = " << setw(12) << setprecision(5) << rmax_m << " [m]\n";
os << "* rmin = " << setw(12) << setprecision(5) << rmin_m << " [m]\n";
os << "* rms beam size = " << setw(12) << setprecision(5) << rrms_m << " [m]\n";
os << "* rms momenta = " << setw(12) << setprecision(5) << prms_m << " [beta gamma]\n";
os << "* mean position = " << setw(12) << setprecision(5) << rmean_m << " [m]\n";
os << "* mean momenta = " << setw(12) << setprecision(5) << pmean_m << " [beta gamma]\n";
os << "* rms emittance = " << setw(12) << setprecision(5) << eps_m << " (not normalized)\n";
os << "* rms correlation = " << setw(12) << setprecision(5) << rprms_m << "\n";
os << "* hr = " << setw(12) << setprecision(5) << hr_m << " [m]\n";
os << "* dh = " << setw(12) << setprecision(5) << dh_m << " [m]\n";
os << "* t = " << setw(12) << setprecision(5) << getT() << " [s] "
<< "dT = " << setw(12) << getdT() << " [s]\n";
os << "* spos = " << setw(12) << setprecision(5) << get_sPos() << " [m]\n";
os << "* Qtot = " << std::setw(12) << std::setprecision(5) << abs(sum(Q)) * 1.0e9 << " [nC] "
<< "Qi = " << std::setw(12) << std::abs(qi_m) * 1e9 << " [nC]" << "\n";
os << "* Ekin = " << std::setw(12) << std::setprecision(5) << eKin_m << " [MeV] "
<< "dEkin = " << std::setw(12) << dE_m << " [MeV]\n";
os << "* rmax = " << std::setw(12) << std::setprecision(5) << rmax_m << " [m]\n";
os << "* rmin = " << std::setw(12) << std::setprecision(5) << rmin_m << " [m]\n";
os << "* rms beam size = " << std::setw(12) << std::setprecision(5) << rrms_m << " [m]\n";
os << "* rms momenta = " << std::setw(12) << std::setprecision(5) << prms_m << " [beta gamma]\n";
os << "* mean position = " << std::setw(12) << std::setprecision(5) << rmean_m << " [m]\n";
os << "* mean momenta = " << std::setw(12) << std::setprecision(5) << pmean_m << " [beta gamma]\n";
os << "* rms emittance = " << std::setw(12) << std::setprecision(5) << eps_m << " (not normalized)\n";
os << "* rms correlation = " << std::setw(12) << std::setprecision(5) << rprms_m << "\n";
os << "* hr = " << std::setw(12) << std::setprecision(5) << hr_m << " [m]\n";
os << "* dh = " << std::setw(12) << std::setprecision(5) << dh_m << " [m]\n";
os << "* t = " << std::setw(12) << std::setprecision(5) << getT() << " [s] "
<< "dT = " << std::setw(12) << getdT() << " [s]\n";
os << "* spos = " << std::setw(12) << std::setprecision(5) << get_sPos() << " [m]\n";
os << "* ********************************************************************************** " << endl;
os.flags(ff);
}
......@@ -2372,8 +2370,8 @@ void PartBunch::calcBeamParameters_cycl() {
for(unsigned int i = 0 ; i < Dim; i++) {
rrms_m(i) = sqrt(rsqsum(i) / TotalNp);
prms_m(i) = sqrt(psqsum(i) / TotalNp);
//eps_m(i) = sqrt( max( eps2(i), zero ) );
eps_norm_m(i) = sqrt(max(eps2(i), zero));
//eps_m(i) = sqrt( std::max( eps2(i), zero ) );
eps_norm_m(i) = sqrt(std::max(eps2(i), zero));
double tmp = rrms_m(i) * prms_m(i);
fac(i) = (tmp == 0) ? zero : 1.0 / tmp;
}
......@@ -2957,4 +2955,4 @@ bool PartBunch::WeHaveEnergyBins() {
Vector_t PartBunch::get_pmean_Distribution() const {
return dist_m->get_pmean();
}
}
\ No newline at end of file
......@@ -74,7 +74,7 @@ void LossDataSink::writeHeaderH5() {
h5_int64_t rc;
// Write file attributes to describe phase space to H5 file.
std::stringstream OPAL_version;
OPAL_version << PACKAGE_NAME << " " << PACKAGE_VERSION << " svn rev. " << GIT_VERSION;
OPAL_version << PACKAGE_NAME << " " << PACKAGE_VERSION << " git rev. " << GIT_VERSION;
rc = H5WriteFileAttribString(H5file_m, "OPAL_version", OPAL_version.str().c_str());
if(rc != H5_SUCCESS)
ERRORMSG("H5 rc= " << rc << " in " << __FILE__ << " @ line " << __LINE__ << endl);
......@@ -434,4 +434,4 @@ void LossDataSink::saveASCII() {
if(! res)
ERRORMSG("LossDataSink Ippl::Comm->send(smsg, 0, tag) failed " << endl;);
}
}
}
\ No newline at end of file
......@@ -539,7 +539,10 @@ void ParallelTTracker::executeDefaultTracker() {
itsBunch->boundp();
numParticlesInSimulation_m = itsBunch->getTotalNum();
}
writePhaseSpace((step + 1), itsBunch->get_sPos(), true, true);
bool const psDump = (itsBunch->getGlobalTrackStep() - 1) % Options::psDumpFreq != 0;
bool const statDump = (itsBunch->getGlobalTrackStep() - 1) % Options::statDumpFreq != 0;
writePhaseSpace((step + 1), itsBunch->get_sPos(), psDump, statDump);
msg << level2 << "Dump phase space of last step" << endl;
OPALTimer::Timer myt3;
itsOpalBeamline_m.switchElementsOff();
......@@ -2952,4 +2955,4 @@ void ParallelTTracker::freeDeviceMemory() {
// mode:c
// c-basic-offset: 4
// indent-tabs-mode:nil
// End:
// End:
\ No newline at end of file
......@@ -4241,13 +4241,13 @@ void Distribution::SetDistParametersGauss(double massIneV) {
sigmaR_m[2] = std::abs(Attributes::getReal(itsAttr[AttributesT::SIGMAZ]));
}
}
if (std::abs(Attributes::getReal(itsAttr[AttributesT::SIGMAR])) > 0.0) {
sigmaR_m[0] = std::abs(Attributes::getReal(itsAttr[AttributesT::SIGMAR]));
sigmaR_m[1] = std::abs(Attributes::getReal(itsAttr[AttributesT::SIGMAR]));
cutoffR_m[0] = Attributes::getReal(itsAttr[AttributesT::CUTOFFR]);
cutoffR_m[1] = Attributes::getReal(itsAttr[AttributesT::CUTOFFR]);
}
if (std::abs(Attributes::getReal(itsAttr[AttributesT::SIGMAR])) > 0.0) {
sigmaR_m[0] = std::abs(Attributes::getReal(itsAttr[AttributesT::SIGMAR]));
sigmaR_m[1] = std::abs(Attributes::getReal(itsAttr[AttributesT::SIGMAR]));
cutoffR_m[0] = Attributes::getReal(itsAttr[AttributesT::CUTOFFR]);
cutoffR_m[1] = Attributes::getReal(itsAttr[AttributesT::CUTOFFR]);
}
}
......@@ -4669,4 +4669,4 @@ void Distribution::WriteOutFileInjection() {
reduce(numberOfParticles, numberOfParticles, OpAddAssign());
}
}
}
}
\ No newline at end of file
......@@ -8,7 +8,9 @@ extern Inform *gmsg;
OpalBeamline::OpalBeamline():
sections_m(),
elements_m(),
prepared_m(false) {
prepared_m(false),
online_secs_m(NULL)
{
online_sections_m = new int[5 * Ippl::getNodes()];
// This is needed to communicate across the nodes which sections are online. We should allocate memory to store
// (# of nodes) * (max # of sections which are online) integers.
......@@ -54,6 +56,9 @@ CompVec &OpalBeamline::getSuccessors(std::shared_ptr<const Component> element) {
size_t index;
bool found = false;
if (sections_m.size() == 0)
return dummy_list_m;
if (element == NULL)
return sections_m[0].getElements();
......@@ -111,6 +116,8 @@ void OpalBeamline::getSectionIndexAt(const Vector_t &pos, long &initial_guess) c
if(initial_guess > -1 && (size_t) initial_guess >= sections_m.size()) initial_guess = static_cast<long>(sections_m.size()) - 1;
if(prepared_m) {
if (sections_m.size() == 0) return;
if(initial_guess == -1 && pos(2) > sections_m[0].getStart(pos(0), pos(1)))
initial_guess = 0;
......@@ -370,6 +377,11 @@ void OpalBeamline::prepareSections() {
list<double>::iterator pos_it, next_it, last_it;
const double tolerance = 1.e-4;
if (elements_m.size() == 0) {
prepared_m = true;
return;
}
/* there might be elements with length zero or extremely short ones.
we delete them such that they don't appear in the simulation
*/
......
......@@ -66,9 +66,6 @@
#endif
#include <vector>
using std::cerr;
using std::endl;
// Namespace Expressions
// ------------------------------------------------------------------------
......@@ -178,23 +175,24 @@ namespace Expressions {
// ----------------------------------------------------------------------
static const TFunction1<double, double> table1[] = {
{ "TRUNC", -1, Truncate },
{ "ROUND", -1, Round },
{ "FLOOR", -1, std::floor },
{ "CEIL", -1, std::ceil },
{ "SIGN", -1, Sign },
{ "SQRT", -1, std::sqrt },
{ "LOG", -1, std::log },
{ "EXP", -1, std::exp },
{ "SIN", -1, std::sin },
{ "COS", -1, std::cos },
{ "ABS", -1, std::abs },
{ "TAN", -1, std::tan },
{ "ASIN", -1, std::asin },
{ "ACOS", -1, std::acos },
{ "ATAN", -1, std::atan },
{ "TGAUSS", -2, Tgauss },
{ 0, -1, 0 }
{ "TRUNC", -1, Truncate },
{ "ROUND", -1, Round },
{ "FLOOR", -1, std::floor },
{ "CEIL", -1, std::ceil },
{ "SIGN", -1, Sign },
{ "SQRT", -1, std::sqrt },
{ "LOG", -1, std::log },
{ "EXP", -1, std::exp },
{ "SIN", -1, std::sin },
{ "COS", -1, std::cos },
{ "ABS", -1, std::abs },
{ "TAN", -1, std::tan },
{ "ASIN", -1, std::asin },
{ "ACOS", -1, std::acos },
{ "ATAN", -1, std::atan },
{ "TGAUSS", -2, Tgauss },
{ "ERF", -1, std::erf },
{ 0, -1, 0 }
};
// Real functions with two arguments.
......@@ -236,8 +234,8 @@ namespace Expressions {
}
return result;
} else if(Options::warn) {
cerr << "\n### Warning ### \"VMIN\" function of empty array.\n"
<< endl;
std::cerr << "\n### Warning ### \"VMIN\" function of empty array.\n"
<< std::endl;
}
return 0.0;
}
......@@ -250,8 +248,8 @@ namespace Expressions {
}
return result;
} else if(Options::warn) {
cerr << "\n### Warning ### \"VMAX\" function of empty array.\n"
<< endl;
std::cerr << "\n### Warning ### \"VMAX\" function of empty array.\n"
<< std::endl;
}
return 0.0;
}
......@@ -264,8 +262,8 @@ namespace Expressions {
}
return sqrt(result / double(array.size()));
} else if(Options::warn) {
cerr << "\n### Warning ### \"VRMS\" function of empty array.\n"
<< endl;
std::cerr << "\n### Warning ### \"VRMS\" function of empty array.\n"
<< std::endl;
}
return 0.0;
}
......@@ -278,8 +276,8 @@ namespace Expressions {
}
return result;
} else if(Options::warn) {
cerr << "\n### Warning ### \"VABSMAX\" function of empty array.\n"
<< endl;
std::cerr << "\n### Warning ### \"VABSMAX\" function of empty array.\n"
<< std::endl;
}
return 0.0;
}
......@@ -1088,4 +1086,4 @@ namespace Expressions {
return currentArray.release();
}
}
}
\ 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