Code indexing in gitaly is broken and leads to code not being visible to the user. We work on the issue with highest priority.

Skip to content
Snippets Groups Projects

Opal maps

Merged Opal maps
OPAL-maps into master
Merged snuverink_j requested to merge OPAL-maps into master
Compare
and
Show latest version
2 files
+ 100
− 150
Compare changes
  • Side-by-side
  • Inline
Files
2
src/Algorithms/ThickTracker.cpp
+ 85
− 126
@@ -179,7 +179,7 @@ void ThickTracker::changeDT() {
}
}
/*
//TODO complete and test fringe fields
void ThickTracker::insertFringeField(SBend* pSBend, lstruct_t& mBL,
double& beta0, double& gamma0, double& P0, double& q, std::array<double,2>& paramFringe, std::string e){
@@ -232,7 +232,7 @@ void ThickTracker::insertFringeField(SBend* pSBend, lstruct_t& mBL,
fringeField.elementMap=fieldMap;
mBL.push_back(fringeField);
}
*/
/**
@@ -390,17 +390,18 @@ void ThickTracker::track_m()
map_t transferMap;
int sliceidx=0;
dumpStats(sliceidx, true, true);
fMatrix_t refSigma;
refSigma = (itsBunch_m->getSigmaMatrix());
double position = zstart_m;
double spos = zstart_m;
this->advanceDispersion_m(tFMatrix, dispInitialVal, spos);
this->advanceDispersion_m(tFMatrix, dispInitialVal, position);
std::size_t step = 0;
this->dump_m();
//(1) Loop Beamline
for(beamline_t::const_iterator it = elements_m.cbegin(); it != elements_m.end(); ++it) {
@@ -424,36 +425,54 @@ void ThickTracker::track_m()
this->concatenateMaps_m(map, transferMap);
tFMatrix= transferMap.linearTerms();
position += ds;
this->advanceDispersion_m(tFMatrix, dispInitialVal, spos);
this->advanceDispersion_m(tFMatrix, dispInitialVal, position);
spos += ds;
++step;
this->update_m(spos, step);
bool const psDump = true; //((itsBunch_m->getGlobalTrackStep() % Options::psDumpFreq) + 1 == Options::psDumpFreq);
bool const statDump = true; //((itsBunch_m->getGlobalTrackStep() % Options::statDumpFreq) + 1 == Options::statDumpFreq);
changeDT();
itsBunch_m->set_sPos(position);
this->dump_m();
dumpStats(sliceidx, psDump, statDump);
refSigma = itsBunch_m->getSigmaMatrix();
//printPhaseShift(refSigma ,mapBeamLineit->elementMap.linearTerms(), N);
sliceidx++;
}
}
std::ofstream tmap;
tmap.open ("TransferMap.txt");
tmap << std::setprecision(16);
tmap << transferMap << std::endl;
tmap.close();
this->write_m(transferMap);
IpplTimings::stopTimer(mapTracking_m);
}
void ThickTracker::write_m(const map_t& map) {
if ( Ippl::myNode() == 0 ) {
static bool first = true;
std::string fn = OpalData::getInstance()->getInputBasename() + ".map";
std::ofstream out;
if ( first ) {
first = false;
out.open(fn, std::ios::out);
} else {
out.open(fn, std::ios::app);
}
out << std::setprecision(16)
<< map
<< std::endl;
out.close();
}
}
void ThickTracker::advanceParticles_m(const map_t& map) {
double betagamma = itsBunch_m->getInitialBeta() * itsBunch_m->getInitialGamma();
@@ -557,65 +576,34 @@ void ThickTracker::advanceDispersion_m(fMatrix_t tempMatrix,
}
void ThickTracker::dumpStats(long long step, bool psDump, bool statDump) {
OPALTimer::Timer myt2;
Inform msg("ThickTracker", *gmsg);
std::size_t numParticlesInSimulation_m = itsBunch_m->getTotalNum();
if (itsBunch_m->getGlobalTrackStep() % 10 + 1 == 10) {
msg << level1;
} else if (true){ //itsBunch_m->getGlobalTrackStep() % 1 + 1 == 1) {
msg << level2;
} else {
msg << level3;
}
if (numParticlesInSimulation_m == 0) {
msg << myt2.time() << " "
<< "Step " << std::setw(6) << itsBunch_m->getGlobalTrackStep() << "; "
<< " -- no emission yet -- "
<< "t= " << Util::getTimeString(itsBunch_m->getT())
<< endl;
void ThickTracker::dump_m() {
if ( itsBunch_m->getTotalNum() == 0 )
return;
}
itsBunch_m->calcEMean();
//size_t totalParticles_f = numParticlesInSimulation_m;
if (std::isnan(pathLength_m) || std::isinf(pathLength_m)) {
throw OpalException("ThickTracker::dumpStats()",
"there seems to be something wrong with the position of the bunch!");
} else {
msg << myt2.time() << " "
<< "Step " << std::setw(6) << itsBunch_m->getGlobalTrackStep() << " "
<< "at " << Util::getLengthString(pathLength_m) << ", "
<< "t= " << Util::getTimeString(itsBunch_m->getT()) << ", "
<< "E=" << Util::getEnergyString(itsBunch_m->get_meanKineticEnergy())
<< endl;
writePhaseSpace(step, psDump, statDump);
}
}
void ThickTracker::writePhaseSpace(const long long step, bool psDump, bool statDump) {
extern Inform *gmsg;
Inform msg("OPAL ", *gmsg);
Vector_t externalE, externalB;
Vector_t FDext[2]; // FDext = {BHead, EHead, BRef, ERef, BTail, ETail}.
Inform msg("ThickTracker", *gmsg);
msg << *itsBunch_m << endl;
const std::size_t step = itsBunch_m->getGlobalTrackStep();
bool psDump = (step + 1) % Options::psDumpFreq;
bool statDump = (step + 1) % Options::statDumpFreq;
// Sample fields at (xmin, ymin, zmin), (xmax, ymax, zmax) and the centroid location. We
// are sampling the electric and magnetic fields at the back, front and
// center of the beam.
Vector_t rmin, rmax;
itsBunch_m->get_bounds(rmin, rmax);
Vector_t FDext[2]; // FDext = {BHead, EHead, BRef, ERef, BTail, ETail}.
if (psDump || statDump) {
Vector_t externalE, externalB;
Vector_t rmin, rmax;
itsBunch_m->get_bounds(rmin, rmax);
externalB = Vector_t(0.0);
externalE = Vector_t(0.0);
itsOpalBeamline_m.getFieldAt(referenceToLabCSTrafo_m.transformTo(RefPartR_m),
@@ -626,65 +614,36 @@ void ThickTracker::writePhaseSpace(const long long step, bool psDump, bool statD
FDext[0] = referenceToLabCSTrafo_m.rotateFrom(externalB);
FDext[1] = referenceToLabCSTrafo_m.rotateFrom(externalE * 1e-6);
}
if (statDump) {
if ( psDump ) {
itsDataSink_m->writePhaseSpace(itsBunch_m, FDext);
}
if ( statDump ) {
std::vector<std::pair<std::string, unsigned int> > collimatorLosses;
// Write statistical data.
itsDataSink_m->writeStatData(itsBunch_m, FDext, collimatorLosses);
msg << level3 << "* Wrote beam statistics." << endl;
}
}
if (psDump && (itsBunch_m->getTotalNum() > 0)) {
// Write fields to .h5 file.
const size_t localNum = itsBunch_m->getLocalNum();
double distToLastStop = zStop_m.back() - pathLength_m;
Vector_t driftPerTimeStep = itsBunch_m->getdT() * Physics::c * RefPartP_m / Util::getGamma(RefPartP_m);
bool driftToCorrectPosition = std::abs(distToLastStop) < 0.5 * euclidean_norm(driftPerTimeStep);
Ppos_t stashedR;
if (driftToCorrectPosition) {
const double tau = distToLastStop / euclidean_norm(driftPerTimeStep) * itsBunch_m->getdT();
if (localNum > 0) {
stashedR.create(localNum);
stashedR = itsBunch_m->R;
for (size_t i = 0; i < localNum; ++ i) {
itsBunch_m->R[i] += tau * (Physics::c * itsBunch_m->P[i] / Util::getGamma(itsBunch_m->P[i]) -
driftPerTimeStep / itsBunch_m->getdT());
}
}
itsBunch_m->RefPartR_m = referenceToLabCSTrafo_m.transformTo(RefPartR_m + tau * driftPerTimeStep / itsBunch_m->getdT());
CoordinateSystemTrafo update(tau * driftPerTimeStep / itsBunch_m->getdT(),
Quaternion());
itsBunch_m->toLabTrafo_m = referenceToLabCSTrafo_m * update.inverted();
itsBunch_m->set_sPos(zStop_m.back());
itsBunch_m->calcBeamParameters();
}
if (!statDump && !driftToCorrectPosition) itsBunch_m->calcBeamParameters();
msg << *itsBunch_m << endl;
itsDataSink_m->writePhaseSpace(itsBunch_m, FDext);
if (driftToCorrectPosition) {
if (localNum > 0) {
itsBunch_m->R = stashedR;
stashedR.destroy(localNum, 0);
}
itsBunch_m->RefPartR_m = referenceToLabCSTrafo_m.transformTo(RefPartR_m);
itsBunch_m->set_sPos(pathLength_m);
itsBunch_m->calcBeamParameters();
}
void ThickTracker::update_m(const double& spos,
const std::size_t& step)
{
itsBunch_m->set_sPos(spos);
itsBunch_m->setGlobalTrackStep(step);
itsBunch_m->calcBeamParameters();
itsBunch_m->calcEMean();
// itsBunch_m->setT();
// changeDT();
msg << level2 << "* Wrote beam phase space." << endl;
}
}
void ThickTracker::setTime() {
const unsigned int localNum = itsBunch_m->getLocalNum();
for (unsigned int i = 0; i < localNum; ++i) {